[MITgcm-devel] potential problem in nonhydrostatic code

[Martin Losch]

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