[MITgcm-devel] potential problem in nonhydrostatic code

[Martin Losch]

I have this suggestion for the lateral boundary conditions for the 
viscosity terms in calc_gw.F.
right now the boundary conditions seem OK to me for vertical wall, but 
with steps and
in particular with partial (looped) cells they will give the wrong 
fluxes. I suggest this:

for example the x-direction (y should be analogous):
instead of
>      &    -viscAh*_recip_dxC(I,J,bi,bj)
>      &      *(1. _d 0 + slipSideFac*
>      &         (maskW(I,J,K-1,bi,bj)+maskW(I,J,K,bi,bj)-2. _d 0))
>      &                   *(wVel(I,J,K,bi,bj)-wVel(I-1,J,K,bi,bj))
with slipSideFac = +/-0.5
put
      &    -viscAh*_recip_dxC(I,J,bi,bj)
      &      *(1. _d 0 + slipSideFac*
      &         (maskW(I,J,K-1,bi,bj)*maskW(I,J,K,bi,bj)-1. _d 0))
      &         
*(max(hFacW(I,J,K-1,bi,bj)-HALF,0)+min(hFacW(I,J,K,bi,bj),HALF)
      &                   *(wVel(I,J,K,bi,bj)-wVel(I-1,J,K,bi,bj))
with slipSideFac = +/-1.0

The advective terms look OK to me.
What do you think?
Martin

On Friday, May 14, 2004, at 04:12 PM, Martin Losch wrote:

> 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
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://dev.mitgcm.org/mailman/listinfo/mitgcm-devel