[MITgcm-devel] potential problem in nonhydrostatic code

[Martin Losch]

I thought that the vertical viscous term is correct because there 
w(Nr+1)=wOverRide*w(kp1)=0, which is correct (no w through the 
boundary).

Currently I don't need any lateral boundary conditions for 
non-hydrostatic w, so I may not notice in (the likely) case I implement 
a zillion bugs.
I have a hard time understanding the lateral boundary conditions right 
now, probably because it's Friday afternoon, but I'll give it some more 
thought and then check in any changes next week (including the vertical 
fix, because all of these changes will break plume_on_slope).

Martin

On Friday, May 14, 2004, at 03:12 PM, Alistair Adcroft wrote:

> I think you're right - the term w^2 is over-estimated by 2 near the
> boundary. I think the viscous term may be under-estimated by a factor 
> of 2
> which you didn't fix. While you're at it, can you implement lateral 
> boundary
> conditions too? I never put any boundary conditions in any of the 
> terms in
> calc_gw() .
>
> A.
> --
> Dr Alistair Adcroft            http://www.mit.edu/~adcroft
> MIT Climate Modeling Initiative        tel: (617) 253-5938
> EAPS 54-1624,  77 Massachusetts Ave,  Cambridge,  MA,  USA
>
> -----Original Message-----
> From: mitgcm-devel-bounces at mitgcm.org
> [mailto:mitgcm-devel-bounces at mitgcm.org] On Behalf Of Martin Losch
> Sent: Friday, May 14, 2004 7:47 AM
> To: MITgcm-devel at mitgcm.org
> Subject: [MITgcm-devel] potential problem in nonhydrostatic code
>
>
> Hi,
>
> I may have stumbled over a bug in the routine calc_gw.F
> The following snipped is the original:
> C Flux on Lower face
>           DO J=J0,Jn
>            DO I=I0,In
>             Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)
>             tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)+wVel(I,J,Kp1,bi,bj))
>             Flx_Dn(I,J,bi,bj)=
>       &     tmp_WbarZ*tmp_WbarZ
>       &    -viscAr*recip_drF(K)
>       &       *( wVel(I,J,K,bi,bj)-wOverRide*wVel(I,J,Kp1,bi,bj) )
>            ENDDO
>           ENDDO
>
> kp1=k+1, except when k=Nr, then kp1=k, in which case wOverRide=0
> instead of the default value of 1
> In this formulation tmp_WbarZ = wVel(I,J,Nr,bi,bj) at the bottom for
> k=Nr. I have the impression that this is not correct, at least not for
> a flat bottom, where w=0 (no flux through the bottom) at bottom, so
> that tmp_WbarZ (at the lower interface of the w-control volume) should
> be .5*wVel(I,J,Nr,bi,bj).
>
> In fact I get numerical instabilities with the original formulation
> after 15000 timesteps, when convective plumes start hitting the bottom
> in my experiment with 10m resolution (vertical and horizontal). These
> instabilities don't seem to occur, when I replace the above
> snippet by this:
> C Flux on Lower face
>           DO J=J0,Jn
>            DO I=I0,In
>             Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)
> CML(
>             tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)+
>       &          wOverRide*wVel(I,J,Kp1,bi,bj))
> CML
>             Flx_Dn(I,J,bi,bj)=
>       &     tmp_WbarZ*tmp_WbarZ
>       &    -viscAr*recip_drF(K)
>       &       *( wVel(I,J,K,bi,bj)-wOverRide*wVel(I,J,Kp1,bi,bj) )
>            ENDDO
>           ENDDO
>
>
> What do you think?
>
> Martin
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org 
> http://dev.mitgcm.org/mailman/listinfo/mitgcm-devel
>
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://dev.mitgcm.org/mailman/listinfo/mitgcm-devel