[MITgcm-devel] potential problem in nonhydrostatic code

[Alistair Adcroft]

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