[MITgcm-devel] potential problem in nonhydrostatic code
Alistair Adcroft
adcroft at MIT.EDU
Wed May 19 09:07:04 EDT 2004
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
More information about the MITgcm-devel
mailing list