[MITgcm-devel] potential problem in nonhydrostatic code

Wed May 19 09:23:10 EDT 2004

Alistair,
exp5 is also affected, so I modified it's output.txt, too.

Martin

On Wednesday, May 19, 2004, at 03:07 PM, Alistair Adcroft wrote:

> OK by me. Is the gravity plume expt the only one that changes? Why not
> internal wave test?
>
> Cc'ing Legg who will probably be affected by this but only for the 
> better.
>
> 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: Tuesday, May 18, 2004 8:09 AM
> To: MITgcm-devel at mitgcm.org
> Subject: Re: [MITgcm-devel] potential problem in nonhydrostatic code
>
>
> 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
>
> _______________________________________________
> 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