[MITgcm-devel] potential problem in nonhydrostatic code

Tue May 18 08:08:56 EDT 2004

Hi,

are there any comments on the boundary value handling? Otherwise I am 
going to check it in like this.
The bug fix in calc_gw.F for w^2 near the bottom changes the 
verification experiment plume_on_slope. Because of this, it makes sense 
to change the (probably ad-hoc) half slip condition for vertical 
velocities to the default lateral slip conditions (which is 
no_slip_sides=.true.). I'll do that, too. Any objections?

Martin

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

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