C $Header: /u/gcmpack/MITgcm/model/src/calc_gw.F,v 1.44 2010/05/11 20:54:11 jmc Exp $
C $Name: $
#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, km1, kp1
_RL mskM1, mskP1
_RL tmp_WbarZ
_RL uTrans, vTrans, rTrans
_RL viscLoc
_RL halfRL
_RS halfRS, zeroRS
PARAMETER( halfRL = 0.5 _d 0 )
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
#endif /* ALLOW_DIAGNOSTICS */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#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 (vertical advection of vertical momentum):
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
flxAdvUp(i,j) = 0.
c flxDisUp(i,j) = 0.
ENDDO
ENDDO
C--- Sweep down column
DO k=1,Nr
km1 = MAX( k-1, 1 )
kp1 = MIN( k+1,Nr )
mskM1 = 1.
mskP1 = 1.
IF ( k.EQ. 1 ) mskM1 = 0.
IF ( k.EQ.Nr ) mskP1 = 0.
IF ( k.GT.1 ) THEN
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 */
ELSEIF ( selectNHfreeSurf.GE.1 ) THEN
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
recip_rThickC(i,j) = recip_drC(k)
c rThickC_C(i,j) = drC(k)
c rThickC_W(i,j) = drC(k)
c rThickC_S(i,j) = drC(k)
c xA(i,j) = _dyG(i,j,bi,bj)*drC(k)
c yA(i,j) = _dxG(i,j,bi,bj)*drC(k)
ENDDO
ENDDO
ENDIF
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-- horizontal bi-harmonic dissipation
IF ( momViscosity .AND. k.GT.1 .AND. viscA4W.NE.0. ) THEN
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 .AND. k.GT.1 ) 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)*mskP1
& -wVel(i,j, k ,bi,bj) )*rkSign
& *recip_drF(k)*rA(i,j,bi,bj)
& *deepFac2C(k)*rhoFacC(k)
ENDDO
ENDDO
IF ( k.EQ.2 ) THEN
C Viscous Flux on Upper face of W-Cell (= at tracer-cell center, level k-1)
DO j=jMin,jMax
DO i=iMin,iMax
C Interpolate horizontally (but not vertically) vert viscosity to center:
C Although background visc. might be defined at k=1, this is not
C generally true when using variable visc. (from vertical mixing scheme).
C Therefore, no vert. interp. and only horizontal interpolation.
viscLoc = ( KappaRU(i,j,k) + KappaRU(i+1,j,k)
& +KappaRV(i,j,k) + KappaRV(i,j+1,k)
& )*0.25 _d 0
flxDisUp(i,j) =
& - viscLoc*( wVel(i,j, k ,bi,bj)
& -wVel(i,j,k-1,bi,bj) )*rkSign
& *recip_drF(k-1)*rA(i,j,bi,bj)
& *deepFac2C(k-1)*rhoFacC(k-1)
C to recover old (before 2009/11/30) results (since flxDisUp(k=2) was zero)
c flxDisUp(i,j) = 0.
ENDDO
ENDDO
ENDIF
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. k.GT.1 .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
IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
C Advective Flux on Western face
DO j=jMin,jMax
DO i=iMin,iMax+1
C transport through Western face area:
uTrans = (
& drF(km1)*_hFacW(i,j,km1,bi,bj)*uVel(i,j,km1,bi,bj)
& *rhoFacC(km1)*mskM1
& + 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*(wFld(i,j)+wFld(i-1,j))*halfRL
c flx_EW(i,j)=
c & 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(km1)*_hFacS(i,j,km1,bi,bj)*vVel(i,j,km1,bi,bj)
& *rhoFacC(km1)*mskM1
& +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*(wFld(i,j)+wFld(i,j-1))*halfRL
c flx_NS(i,j)=
c & vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL
ENDDO
ENDDO
ENDIF
C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)
c IF (.TRUE.) THEN
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)*mskP1 )
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)
& *mskP1
& )*rA(i,j,bi,bj)
flx_Dn(i,j) = rTrans*tmp_WbarZ
ENDDO
ENDDO
c ENDIF
IF ( k.EQ.1 .AND. selectNHfreeSurf.GE.1 ) THEN
C Advective Flux on Upper face of W-Cell (= at surface)
DO j=jMin,jMax
DO i=iMin,iMax
tmp_WbarZ = wVel(i,j,k,bi,bj)*rVel2wUnit(k)
rTrans = wVel(i,j,k,bi,bj)*deepFac2F(k)*rhoFacF(k)
& *rA(i,j,bi,bj)
flxAdvUp(i,j) = rTrans*tmp_WbarZ
c flxAdvUp(i,j) = 0.
ENDDO
ENDDO
ENDIF
IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
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
C to recover old (before 2009/11/30) results (since flxAdvUp(k=2) was zero)
c IF (k.EQ.2) flxAdvUp(i,j) = 0.
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)
ENDDO
ENDDO
ENDIF
DO j=jMin,jMax
DO i=iMin,iMax
C-- prepare for next level (k+1)
flxAdvUp(i,j)=flx_Dn(i,j)
ENDDO
ENDDO
c ELSE
C- if momAdvection / else
c DO j=1-OLy,sNy+OLy
c DO i=1-OLx,sNx+OLx
c gW(i,j,k,bi,bj) = 0. _d 0
c ENDDO
c ENDDO
C- endif momAdvection.
ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
IF ( useNHMTerms .AND. k.GT.1 ) 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 .AND. k.GT.1 ) 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)
c 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 )
c 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.
#ifdef ALLOW_ADAMSBASHFORTH_3
CALL ADAMS_BASHFORTH3(
I bi, bj, k,
U gW, gwNm,
I nHydStartAB, myIter, myThid )
#else /* ALLOW_ADAMSBASHFORTH_3 */
CALL ADAMS_BASHFORTH2(
I bi, bj, k,
U gW, gwNm1,
I nHydStartAB, myIter, myThid )
#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