C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/calc_gw.F,v 1.35 2006/12/05 05:25:08 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: CALC_GW
C     !INTERFACE:
      SUBROUTINE CALC_GW(
     I               bi, bj, KappaRU, KappaRV,
     I               myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R CALC_GW
C     | o Calculate vertical velocity tendency terms
C     |   ( Non-Hydrostatic only )
C     *==========================================================*
C     | In NH, the vertical momentum tendency must be
C     | calculated explicitly and included as a source term
C     | for a 3d pressure eqn. Calculate that term here.
C     | This routine is not used in HYD calculations.
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "DYNVARS.h"
#include "NH_VARS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     bi,bj   :: current tile indices
C     KappaRU :: vertical viscosity at U points
C     KappaRV :: vertical viscosity at V points
C     myTime  :: Current time in simulation
C     myIter  :: Current iteration number in simulation
C     myThid  :: Thread number for this instance of the routine.
      INTEGER bi,bj
      _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL     myTime
      INTEGER myIter
      INTEGER myThid

#ifdef ALLOW_NONHYDROSTATIC

C     !LOCAL VARIABLES:
C     == Local variables ==
C     iMin,iMax
C     jMin,jMax
C     xA            :: W-Cell face area normal to X
C     yA            :: W-Cell face area normal to Y
C     rThickC_W     :: thickness (in r-units) of W-Cell at Western Edge
C     rThickC_S     :: thickness (in r-units) of W-Cell at Southern Edge
C     rThickC_C     :: thickness (in r-units) of W-Cell (centered on W pt)
C     recip_rThickC :: reciprol thickness of W-Cell (centered on W-point)
C     flx_NS        :: vertical momentum flux, meridional direction
C     flx_EW        :: vertical momentum flux, zonal direction
C     flxAdvUp      :: vertical mom. advective   flux, vertical direction (@ level k-1)
C     flxDisUp      :: vertical mom. dissipation flux, vertical direction (@ level k-1)
C     flx_Dn        :: vertical momentum flux, vertical direction (@ level k)
C     gwDiss        :: vertical momentum dissipation tendency
C     i,j,k         :: Loop counters
      INTEGER iMin,iMax,jMin,jMax
      _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL    rThickC_C    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL    recip_rThickC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL    flxAdvUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL    flxDisUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL    gwDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL    gwAdd (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL    del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER i,j,k, kp1
      _RL  wOverride
      _RL  tmp_WbarZ
      _RL  uTrans, vTrans, rTrans
      _RL  viscLoc
      _RL  halfRL
      _RS  halfRS, zeroRS
      PARAMETER( halfRL = 0.5D0 )
      PARAMETER( halfRS = 0.5 , zeroRS = 0. )
CEOP

C Catch barotropic mode
      IF ( Nr .LT. 2 ) RETURN

      iMin = 1
      iMax = sNx
      jMin = 1
      jMax = sNy

C--   Initialise gW to zero
      DO k=1,Nr
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
           gW(i,j,k,bi,bj) = 0.
         ENDDO
        ENDDO
      ENDDO
C-    Initialise gwDiss to zero
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
         gwDiss(i,j) = 0.
       ENDDO
      ENDDO

C--   Boundaries condition at top
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
         flxAdvUp(i,j) = 0.
         flxDisUp(i,j) = 0.
       ENDDO
      ENDDO

C---  Sweep down column
      DO k=2,Nr
        kp1=k+1
        wOverRide=1.
        IF (k.EQ.Nr) THEN
          kp1=Nr
          wOverRide=0.
        ENDIF
C--   Compute grid factor arround a W-point:
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
C-    note: assume fluid @ smaller k than bottom: does not work in p-coordinate !
           IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN
             recip_rThickC(i,j) = 0.
           ELSE
             recip_rThickC(i,j) = 1. _d 0 /
     &        ( drF(k-1)*halfRS
     &        + drF( k )*MIN( _hFacC(i,j, k ,bi,bj), halfRS )
     &        )
           ENDIF
c          IF (momViscosity) THEN
#ifdef NONLIN_FRSURF
           rThickC_C(i,j) =
     &          drF(k-1)*MAX( h0FacC(i,j,k-1,bi,bj)-halfRS, zeroRS )
     &        + drF( k )*MIN( h0FacC(i,j,k  ,bi,bj), halfRS )
#else
           rThickC_C(i,j) =
     &          drF(k-1)*MAX( _hFacC(i,j,k-1,bi,bj)-halfRS, zeroRS )
     &        + drF( k )*MIN( _hFacC(i,j,k  ,bi,bj), halfRS )
#endif
           rThickC_W(i,j) =
     &          drF(k-1)*MAX( _hFacW(i,j,k-1,bi,bj)-halfRS, zeroRS )
     &        + drF( k )*MIN( _hFacW(i,j,k  ,bi,bj), halfRS )
           rThickC_S(i,j) =
     &          drF(k-1)*MAX( _hFacS(i,j,k-1,bi,bj)-halfRS, zeroRS )
     &        + drF( k )*MIN( _hFacS(i,j, k ,bi,bj), halfRS )
C     W-Cell Western face area:
           xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
c    &                              *deepFacF(k)
C     W-Cell Southern face area:
           yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
c    &                              *deepFacF(k)
C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
c          ENDIF
         ENDDO
        ENDDO

C--   horizontal bi-harmonic dissipation
        IF (momViscosity .AND. viscA4W.NE.0. ) THEN
C-    calculate the horizontal Laplacian of vertical flow
C     Zonal flux d/dx W
          DO j=1-Oly,sNy+Oly
           flx_EW(1-Olx,j)=0.
           DO i=1-Olx+1,sNx+Olx
            flx_EW(i,j) =
     &              (wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
     &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
#ifdef COSINEMETH_III
     &              *sqCosFacU(j,bi,bj)
#endif
           ENDDO
          ENDDO
C     Meridional flux d/dy W
          DO i=1-Olx,sNx+Olx
           flx_NS(i,1-Oly)=0.
          ENDDO
          DO j=1-Oly+1,sNy+Oly
           DO i=1-Olx,sNx+Olx
            flx_NS(i,j) =
     &              (wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
     &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
#ifdef ISOTROPIC_COS_SCALING
#ifdef COSINEMETH_III
     &              *sqCosFacV(j,bi,bj)
#endif
#endif
           ENDDO
          ENDDO

C     del^2 W
C     Difference of zonal fluxes ...
          DO j=1-Oly,sNy+Oly
           DO i=1-Olx,sNx+Olx-1
            del2w(i,j)=flx_EW(i+1,j)-flx_EW(i,j)
           ENDDO
           del2w(sNx+Olx,j)=0.
          ENDDO

C     ... add difference of meridional fluxes and divide by volume
          DO j=1-Oly,sNy+Oly-1
           DO i=1-Olx,sNx+Olx
            del2w(i,j) = ( del2w(i,j)
     &                    +(flx_NS(i,j+1)-flx_NS(i,j))
     &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
     &                    *recip_deepFac2F(k)
           ENDDO
          ENDDO
C-- No-slip BCs impose a drag at walls...
CML ************************************************************
CML   No-slip Boundary conditions for bi-harmonic dissipation
CML   need to be implemented here!
CML ************************************************************
        ELSE
C-    Initialize del2w to zero:
          DO j=1-Oly,sNy+Oly
           DO i=1-Olx,sNx+Olx
            del2w(i,j) = 0. _d 0
           ENDDO
          ENDDO
        ENDIF

        IF (momViscosity) THEN
C Viscous Flux on Western face
          DO j=jMin,jMax
           DO i=iMin,iMax+1
             flx_EW(i,j)=
     &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*halfRL
     &              *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
     &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
cOld &              *_recip_dxC(i,j,bi,bj)*rThickC_W(i,j)
     &              *cosFacU(j,bi,bj)
     &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL
     &              *(del2w(i,j)-del2w(i-1,j))
     &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
cOld &              *_recip_dxC(i,j,bi,bj)*drC(k)
#ifdef COSINEMETH_III
     &              *sqCosFacU(j,bi,bj)
#else
     &              *cosFacU(j,bi,bj)
#endif
           ENDDO
          ENDDO
C Viscous Flux on Southern face
          DO j=jMin,jMax+1
           DO i=iMin,iMax
             flx_NS(i,j)=
     &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*halfRL
     &              *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
     &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
cOld &              *_recip_dyC(i,j,bi,bj)*rThickC_S(i,j)
#ifdef ISOTROPIC_COS_SCALING
     &              *cosFacV(j,bi,bj)
#endif
     &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL
     &              *(del2w(i,j)-del2w(i,j-1))
     &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
cOld &              *_recip_dyC(i,j,bi,bj)*drC(k)
#ifdef ISOTROPIC_COS_SCALING
#ifdef COSINEMETH_III
     &              *sqCosFacV(j,bi,bj)
#else
     &              *cosFacV(j,bi,bj)
#endif
#endif
           ENDDO
          ENDDO
C Viscous Flux on Lower face of W-Cell (= at tracer-cell center, level k)
          DO j=jMin,jMax
           DO i=iMin,iMax
C     Interpolate vert viscosity to center of tracer-cell (level k):
             viscLoc = ( KappaRU(i,j,k)  +KappaRU(i+1,j,k)
     &                  +KappaRU(i,j,kp1)+KappaRU(i+1,j,kp1)
     &                  +KappaRV(i,j,k)  +KappaRV(i,j+1,k)
     &                  +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)
     &                 )*0.125 _d 0
             flx_Dn(i,j) =
     &          - viscLoc*( wVel(i,j,kp1,bi,bj)*wOverRide
     &                     -wVel(i,j, k ,bi,bj) )*rkSign
     &                   *recip_drF(k)*rA(i,j,bi,bj)
     &                   *deepFac2C(k)*rhoFacC(k)
cOld &                   *recip_drF(k)
           ENDDO
          ENDDO
C     Tendency is minus divergence of viscous fluxes:
C     anelastic: vert.visc.flx is scaled by rhoFac but hor.visc.fluxes are not
          DO j=jMin,jMax
           DO i=iMin,iMax
             gwDiss(i,j) =
     &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
     &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
     &           + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign
     &                                          *recip_rhoFacF(k)
     &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
     &          *recip_deepFac2F(k)
cOld         gwDiss(i,j) =
cOld &        -(
cOld &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )
cOld &          +_recip_dyF(i,j,bi,bj)*( flx_NS(i,j+1)-flx_NS(i,j) )
cOld &          +                      ( flxDisUp(i,j)-flx_Dn(i,j) )
c    &         )*recip_rThickC(i,j)
cOld &         )*recip_drC(k)
C--        prepare for next level (k+1)
             flxDisUp(i,j)=flx_Dn(i,j)
           ENDDO
          ENDDO
        ENDIF

        IF (no_slip_sides) THEN
C-     No-slip BCs impose a drag at walls...
          CALL MOM_W_SIDEDRAG(
     I               bi,bj,k,
     I               wVel, del2w,
     I               rThickC_C, recip_rThickC,
     I               viscAh_W, viscA4_W,
     O               gwAdd,
     I               myThid )
          DO j=jMin,jMax
           DO i=iMin,iMax
            gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
           ENDDO
          ENDDO
        ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

        IF ( momAdvection ) THEN
C Advective Flux on Western face
          DO j=jMin,jMax
           DO i=iMin,iMax+1
C     transport through Western face area:
             uTrans = (
     &          drF(k-1)*_hFacW(i,j,k-1,bi,bj)*uVel(i,j,k-1,bi,bj)
     &                  *rhoFacC(k-1)
     &        + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj)
     &                  *rhoFacC(k)
     &                )*halfRL*_dyG(i,j,bi,bj)*deepFacF(k)
cOld &                )*halfRL
             flx_EW(i,j)=
     &         uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL
           ENDDO
          ENDDO
C Advective Flux on Southern face
          DO j=jMin,jMax+1
           DO i=iMin,iMax
C     transport through Southern face area:
             vTrans = (
     &          drF(k-1)*_hFacS(i,j,k-1,bi,bj)*vVel(i,j,k-1,bi,bj)
     &                  *rhoFacC(k-1)
     &         +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj)
     &                  *rhoFacC(k)
     &                )*halfRL*_dxG(i,j,bi,bj)*deepFacF(k)
cOld &                )*halfRL
             flx_NS(i,j)=
     &         vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL
           ENDDO
          ENDDO
C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)
          DO j=jMin,jMax
           DO i=iMin,iMax
             tmp_WbarZ = halfRL*( wVel(i,j, k ,bi,bj)
     &                           +wVel(i,j,kp1,bi,bj)*wOverRide )
C     transport through Lower face area:
             rTrans = halfRL*
     &              ( wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k )
     &               +wVel(i,j,kp1,bi,bj)*deepFac2F(kp1)*rhoFacF(kp1)
     &                                   *wOverRide
     &              )*rA(i,j,bi,bj)
             flx_Dn(i,j) = rTrans*tmp_WbarZ
cOld         flx_Dn(i,j) = tmp_WbarZ*tmp_WbarZ
           ENDDO
          ENDDO
C     Tendency is minus divergence of advective fluxes:
C     anelastic: all transports & advect. fluxes are scaled by rhoFac
          DO j=jMin,jMax
           DO i=iMin,iMax
             gW(i,j,k,bi,bj) =
     &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
     &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
     &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign
     &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
     &          *recip_deepFac2F(k)*recip_rhoFacF(k)
cOld         gW(i,j,k,bi,bj) =
cOld &        -(
cOld &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )
cOld &          +_recip_dyF(i,j,bi,bj)*( flx_NS(i,j+1)-flx_NS(i,j) )
cOld &          +                      ( flxAdvUp(i,j)-flx_Dn(i,j) )
c    &         )*recip_rThickC(i,j)
cOld &         )*recip_drC(k)
C--          prepare for next level (k+1)
             flxAdvUp(i,j)=flx_Dn(i,j)
           ENDDO
          ENDDO
        ENDIF

        IF ( useNHMTerms ) THEN
          CALL MOM_W_METRIC_NH(
     I               bi,bj,k,
     I               uVel, vVel,
     O               gwAdd,
     I               myThid )
          DO j=jMin,jMax
           DO i=iMin,iMax
             gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
           ENDDO
          ENDDO
        ENDIF
        IF ( use3dCoriolis ) THEN
          CALL MOM_W_CORIOLIS_NH(
     I               bi,bj,k,
     I               uVel, vVel,
     O               gwAdd,
     I               myThid )
          DO j=jMin,jMax
           DO i=iMin,iMax
             gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
           ENDDO
          ENDDO
        ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--   Dissipation term inside the Adams-Bashforth:
        IF ( momViscosity .AND. momDissip_In_AB) THEN
          DO j=jMin,jMax
           DO i=iMin,iMax
             gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
           ENDDO
          ENDDO
        ENDIF

C-    Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)
C     and save gW_[n] into gwNm1 for the next time step.
c#ifdef ALLOW_ADAMSBASHFORTH_3
c       CALL ADAMS_BASHFORTH3(
c    I                         bi, bj, k,
c    U                         gW, gwNm,
c    I                         momStartAB, myIter, myThid )
c#else /* ALLOW_ADAMSBASHFORTH_3 */
        CALL ADAMS_BASHFORTH2(
     I                         bi, bj, k,
     U                         gW, gwNm1,
     I                         myIter, myThid )
c#endif /* ALLOW_ADAMSBASHFORTH_3 */

C--   Dissipation term outside the Adams-Bashforth:
        IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
          DO j=jMin,jMax
           DO i=iMin,iMax
             gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
           ENDDO
          ENDDO
        ENDIF

C-    end of the k loop
      ENDDO

#endif /* ALLOW_NONHYDROSTATIC */

      RETURN
      END