C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/calc_gw.F,v 1.40 2009/02/27 01:47:41 jmc Exp $
C $Name: checkpoint61o $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
#define CALC_GW_NEW_THICK

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 "RESTART.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     gwAdd         :: other tendencies (Coriolis, Metric-terms)
C     del2w         :: laplacian of wVel
C     wFld          :: local copy of wVel
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)
      _RL    wFld  (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. )
      PARAMETER( iMin = 1 , iMax = sNx )
      PARAMETER( jMin = 1 , jMax = sNy )
CEOP
#ifdef ALLOW_DIAGNOSTICS
      LOGICAL diagDiss, diagAdvec
      LOGICAL  DIAGNOSTICS_IS_ON
      EXTERNAL DIAGNOSTICS_IS_ON
#endif /* ALLOW_DIAGNOSTICS */

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

#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics ) THEN
        diagDiss  = DIAGNOSTICS_IS_ON( 'Wm_Diss ', myThid )
        diagAdvec = DIAGNOSTICS_IS_ON( 'Wm_Advec', myThid )
      ELSE
        diagDiss  = .FALSE.
        diagAdvec = .FALSE.
      ENDIF
#endif /* ALLOW_DIAGNOSTICS */

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
      IF (momViscosity) THEN
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

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:
#ifdef CALC_GW_NEW_THICK
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
           IF ( maskC(i,j,k-1,bi,bj).EQ.0. .OR.
     &          maskC(i,j, k ,bi,bj).EQ.0. ) THEN
             recip_rThickC(i,j) = 0.
           ELSE
C-    valid in z & p coord.; also accurate if Interface @ middle between 2 centers
             recip_rThickC(i,j) = 1. _d 0 /
     &        (  MIN( Ro_surf(i,j,bi,bj),rC(k-1) )
     &         - MAX( R_low(i,j,bi,bj),  rC(k)   )
     &        )
           ENDIF
         ENDDO
        ENDDO
        IF (momViscosity) THEN
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           rThickC_C(i,j) = MAX( zeroRS,
     &                           MIN( Ro_surf(i,j,bi,bj), rC(k-1) )
     &                          -MAX(   R_low(i,j,bi,bj),  rC(k)  )
     &                         )
          ENDDO
         ENDDO
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx+1,sNx+Olx
           rThickC_W(i,j) = MAX( zeroRS,
     &                           MIN( rSurfW(i,j,bi,bj), rC(k-1) )
     &                          -MAX(  rLowW(i,j,bi,bj), rC(k)   )
     &                         )
C     W-Cell Western face area:
           xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
c    &                              *deepFacF(k)
          ENDDO
         ENDDO
         DO j=1-Oly+1,sNy+Oly
          DO i=1-Olx,sNx+Olx
           rThickC_S(i,j) = MAX( zeroRS,
     &                           MIN( rSurfS(i,j,bi,bj), rC(k-1) )
     &                          -MAX(  rLowS(i,j,bi,bj), rC(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)
          ENDDO
         ENDDO
        ENDIF
#else /* CALC_GW_NEW_THICK */
        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
#endif /* CALC_GW_NEW_THICK */

C--   horizontal bi-harmonic dissipation
        IF (momViscosity .AND. viscA4W.NE.0. ) THEN

C-    local copy of wVel:
          DO j=1-Oly,sNy+Oly
           DO i=1-Olx,sNx+Olx
             wFld(i,j) = wVel(i,j,k,bi,bj)
           ENDDO
          ENDDO
C-    calculate the horizontal Laplacian of vertical flow
C     Zonal flux d/dx W
          IF ( useCubedSphereExchange ) THEN
C     to compute d/dx(W), fill corners with appropriate values:
            CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
     &                                 wFld, bi,bj, myThid )
          ENDIF
          DO j=1-Oly,sNy+Oly
           flx_EW(1-Olx,j)=0.
           DO i=1-Olx+1,sNx+Olx
            flx_EW(i,j) =
     &               ( wFld(i,j) - wFld(i-1,j) )
     &              *_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
          IF ( useCubedSphereExchange ) THEN
C     to compute d/dy(W), fill corners with appropriate values:
            CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
     &                                 wFld, bi,bj, myThid )
          ENDIF
          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) =
     &               ( wFld(i,j) - wFld(i,j-1) )
     &              *_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     Divergence of horizontal fluxes
          DO j=1-Oly,sNy+Oly-1
           DO i=1-Olx,sNx+Olx-1
            del2w(i,j) = ( ( flx_EW(i+1,j)-flx_EW(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 end if biharmonic viscosity
        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)
     &              *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)
#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)
#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)
#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)
           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)
C--        prepare for next level (k+1)
             flxDisUp(i,j)=flx_Dn(i,j)
           ENDDO
          ENDDO
        ENDIF

        IF ( momViscosity .AND. 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)
             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)
             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
C     NH in p-coord.: advect wSpeed [m/s] with rTrans
             tmp_WbarZ = halfRL*
     &              ( wVel(i,j, k ,bi,bj)*rVel2wUnit(k)
     &               +wVel(i,j,kp1,bi,bj)*rVel2wUnit(kp1)*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
           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*wUnit2rVel(k)
     &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
     &          *recip_deepFac2F(k)*recip_rhoFacF(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-|--+----|

#ifdef ALLOW_DIAGNOSTICS
        IF ( diagDiss )  THEN
          CALL DIAGNOSTICS_FILL( gwDiss, 'Wm_Diss ',
     &                           k, 1, 2, bi,bj, myThid )
C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
C        does it only if k=1 (never the case here)
          IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Diss ',bi,bj,myThid)
        ENDIF
        IF ( diagAdvec ) THEN
          CALL DIAGNOSTICS_FILL( gW, 'Wm_Advec',
     &                           k,Nr, 1, bi,bj, myThid )
          IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Advec',bi,bj,myThid)
        ENDIF
#endif /* ALLOW_DIAGNOSTICS */

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                         nHydStartAB, myIter, myThid )
c#else /* ALLOW_ADAMSBASHFORTH_3 */
        CALL ADAMS_BASHFORTH2(
     I                         bi, bj, k,
     U                         gW, gwNm1,
     I                         nHydStartAB, 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

#ifdef ALLOW_DIAGNOSTICS
      IF (useDiagnostics) THEN
        CALL DIAGNOSTICS_FILL(viscAh_W,'VISCAHW ',0,Nr,1,bi,bj,myThid)
        CALL DIAGNOSTICS_FILL(viscA4_W,'VISCA4W ',0,Nr,1,bi,bj,myThid)
      ENDIF
#endif /* ALLOW_DIAGNOSTICS */

#endif /* ALLOW_NONHYDROSTATIC */

      RETURN
      END