pkg/mom_fluxform/mom_fluxform.F

71207ba845 Alis*0001 CBOI
                0002 C !TITLE: pkg/mom\_advdiff
                0003 C !AUTHORS: adcroft@mit.edu
cfdc763dc5 Alis*0004 C !INTRODUCTION: Flux-form Momentum Equations Package
71207ba845 Alis*0005 C
                0006 C Package "mom\_fluxform" provides methods for calculating explicit terms
                0007 C in the momentum equation cast in flux-form:
                0008 C \begin{eqnarray*}
                0009 C G^u & = & -\frac{1}{\rho} \partial_x \phi_h
                0010 C           -\nabla \cdot {\bf v} u
                0011 C           -fv
                0012 C           +\frac{1}{\rho} \nabla \cdot {\bf \tau}^x
                0013 C           + \mbox{metrics}
                0014 C \\
                0015 C G^v & = & -\frac{1}{\rho} \partial_y \phi_h
                0016 C           -\nabla \cdot {\bf v} v
                0017 C           +fu
                0018 C           +\frac{1}{\rho} \nabla \cdot {\bf \tau}^y
                0019 C           + \mbox{metrics}
                0020 C \end{eqnarray*}
                0021 C where ${\bf v}=(u,v,w)$ and $\tau$, the stress tensor, includes surface
                0022 C stresses as well as internal viscous stresses.
                0023 CEOI
                0024
6d54cf9ca1 Ed H*0025 #include "MOM_FLUXFORM_OPTIONS.h"
88b07ffa37 Jean*0026 #ifdef ALLOW_AUTODIFF
                0027 # include "AUTODIFF_OPTIONS.h"
                0028 #endif
f0501c53d1 Jean*0029 #ifdef ALLOW_MOM_COMMON
                0030 # include "MOM_COMMON_OPTIONS.h"
                0031 #endif
9293d3c672 Hajo*0032 #ifdef ALLOW_GGL90
                0033 # include "GGL90_OPTIONS.h"
                0034 #endif
cda1c18f72 Jean*0035 C- to switch back to (very) old accounting of geometry/Area factor
                0036 #undef OLD_UV_GEOM
aea29c8517 Alis*0037
71207ba845 Alis*0038 CBOP
                0039 C !ROUTINE: MOM_FLUXFORM
                0040
                0041 C !INTERFACE: ==========================================================
9496c6c9ef Jean*0042       SUBROUTINE MOM_FLUXFORM(
07e17fa184 Jean*0043      I        bi,bj,k,iMin,iMax,jMin,jMax,
1833b564cb Jean*0044      I        kappaRU, kappaRV,
07e17fa184 Jean*0045      U        fVerUkm, fVerVkm,
                0046      O        fVerUkp, fVerVkp,
722a76e285 Jean*0047      O        guDiss, gvDiss,
07e17fa184 Jean*0048      I        myTime, myIter, myThid )
71207ba845 Alis*0049
                0050 C !DESCRIPTION:
                0051 C Calculates all the horizontal accelerations except for the implicit surface
eaba2fd266 Jean*0052 C pressure gradient and implicit vertical viscosity.
aea29c8517 Alis*0053
71207ba845 Alis*0054 C !USES: ===============================================================
aea29c8517 Alis*0055 C     == Global variables ==
71207ba845 Alis*0056       IMPLICIT NONE
aea29c8517 Alis*0057 #include "SIZE.h"
                0058 #include "EEPARAMS.h"
                0059 #include "PARAMS.h"
                0060 #include "GRID.h"
f0501c53d1 Jean*0061 #include "DYNVARS.h"
                0062 #include "FFIELDS.h"
aea29c8517 Alis*0063 #include "SURFACE.h"
f0501c53d1 Jean*0064 #ifdef ALLOW_MOM_COMMON
                0065 # include "MOM_VISC.h"
                0066 #endif
9293d3c672 Hajo*0067 #if ( defined ALLOW_GGL90 && defined ALLOW_GGL90_LANGMUIR )
                0068 # include "GGL90.h"
                0069 #endif
88b07ffa37 Jean*0070 #ifdef ALLOW_AUTODIFF
7c50f07931 Mart*0071 # ifdef ALLOW_AUTODIFF_TAMC
                0072 #  include "tamc.h"
                0073 # endif
aa2d1573fa Patr*0074 # include "MOM_FLUXFORM.h"
                0075 #endif
aea29c8517 Alis*0076
71207ba845 Alis*0077 C !INPUT PARAMETERS: ===================================================
07e17fa184 Jean*0078 C  bi,bj                :: current tile indices
                0079 C  k                    :: current vertical level
                0080 C  iMin,iMax,jMin,jMax  :: loop ranges
1833b564cb Jean*0081 C  kappaRU              :: vertical viscosity
                0082 C  kappaRV              :: vertical viscosity
07e17fa184 Jean*0083 C  fVerUkm              :: vertical advective flux of U, interface above (k-1/2)
                0084 C  fVerVkm              :: vertical advective flux of V, interface above (k-1/2)
                0085 C  fVerUkp              :: vertical advective flux of U, interface below (k+1/2)
                0086 C  fVerVkp              :: vertical advective flux of V, interface below (k+1/2)
722a76e285 Jean*0087 C  guDiss               :: dissipation tendency (all explicit terms), u component
                0088 C  gvDiss               :: dissipation tendency (all explicit terms), v component
bd2e80b12f Jean*0089 C  myTime               :: current time
71207ba845 Alis*0090 C  myIter               :: current time-step number
07e17fa184 Jean*0091 C  myThid               :: my Thread Id number
                0092       INTEGER bi,bj,k
                0093       INTEGER iMin,iMax,jMin,jMax
1833b564cb Jean*0094       _RL kappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
                0095       _RL kappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
07e17fa184 Jean*0096       _RL fVerUkm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0097       _RL fVerVkm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0098       _RL fVerUkp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0099       _RL fVerVkp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
722a76e285 Jean*0100       _RL guDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0101       _RL gvDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
bd2e80b12f Jean*0102       _RL     myTime
3a279374db Alis*0103       INTEGER myIter
                0104       INTEGER myThid
aea29c8517 Alis*0105
71207ba845 Alis*0106 C !OUTPUT PARAMETERS: ==================================================
                0107 C None - updates gU() and gV() in common blocks
                0108
                0109 C !LOCAL VARIABLES: ====================================================
                0110 C  i,j                  :: loop indices
31fb0e0e6d Jean*0111 C  gAdd                 :: individual momentum tendency term to add
                0112 C  del2                 :: horizontal second derivative (Laplacian)
9293d3c672 Hajo*0113 C  uCf,vCf              :: Coriolis acceleration
31fb0e0e6d Jean*0114 C  vMT                  :: vertical/Non-Hyd. metric terms
71207ba845 Alis*0115 C  fZon                 :: zonal fluxes
                0116 C  fMer                 :: meridional fluxes
07e17fa184 Jean*0117 C  fVrUp,fVrDw          :: vertical viscous fluxes at interface k & k+1
71207ba845 Alis*0118       INTEGER i,j
9496c6c9ef Jean*0119 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0120 C     kkey :: tape key (depends on level and tile indices)
                0121       INTEGER kkey
aa2d1573fa Patr*0122 #endif
31fb0e0e6d Jean*0123       _RL gAdd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0124       _RL del2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0125       _RL  uCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0126       _RL  vCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0127       _RL  vMT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71207ba845 Alis*0128       _RL fZon(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0129       _RL fMer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
722a76e285 Jean*0130       _RL fVrUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0131       _RL fVrDw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
4d2713b160 Jean*0132 C     afFacMom     :: Tracer parameters for turning terms on and off.
9496c6c9ef Jean*0133 C     vfFacMom
                0134 C     pfFacMom        afFacMom - Advective terms
aea29c8517 Alis*0135 C     cfFacMom        vfFacMom - Eddy viscosity terms
4d2713b160 Jean*0136 C     mtFacMom        pfFacMom - Pressure terms
aea29c8517 Alis*0137 C                     cfFacMom - Coriolis terms
                0138 C                     foFacMom - Forcing
4d2713b160 Jean*0139 C                     mtFacMom - Metric term
722a76e285 Jean*0140 C     uDudxFac, AhDudxFac, etc ... individual term parameters for switching terms off
aea29c8517 Alis*0141       _RS    hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
fd362e9f7c Jean*0142       _RS   h0FacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
aea29c8517 Alis*0143       _RS  r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
9293d3c672 Hajo*0144       _RS       xA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0145       _RS       yA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0146       _RL   uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0147       _RL   vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0148       _RL     uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0149       _RL     vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
bd2e80b12f Jean*0150       _RL  rTransU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0151       _RL  rTransV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
9293d3c672 Hajo*0152 #if ( defined ALLOW_GGL90 && defined ALLOW_GGL90_LANGMUIR )
                0153       _RL     uRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0154       _RL     vRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0155 #endif
79074ef66b Jean*0156       _RL       KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0157       _RL    cDrag(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
05b9f17ae6 Bayl*0158       _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0159       _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0160       _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0161       _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79074ef66b Jean*0162       _RL    vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0163       _RL     hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0164       _RL   strain(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0165       _RL  tension(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
f59d76b0dd Ed D*0166       _RL stretching(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
06244a5e4f Jean*0167 #ifdef ALLOW_LEITH_QG
2ff762fc08 Jean*0168       _RL  Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
06244a5e4f Jean*0169 #endif
aea29c8517 Alis*0170       _RL  uDudxFac
                0171       _RL  AhDudxFac
                0172       _RL  vDudyFac
                0173       _RL  AhDudyFac
31fb0e0e6d Jean*0174       _RL  wDudrFac
aea29c8517 Alis*0175       _RL  ArDudrFac
                0176       _RL  fuFac
                0177       _RL  mtFacU
4d2713b160 Jean*0178       _RL  mtNHFacU
aea29c8517 Alis*0179       _RL  uDvdxFac
                0180       _RL  AhDvdxFac
                0181       _RL  vDvdyFac
                0182       _RL  AhDvdyFac
31fb0e0e6d Jean*0183       _RL  wDvdrFac
aea29c8517 Alis*0184       _RL  ArDvdrFac
                0185       _RL  fvFac
                0186       _RL  mtFacV
4d2713b160 Jean*0187       _RL  mtNHFacV
ea2dd4993a Jean*0188       _RL  sideMaskFac
31fb0e0e6d Jean*0189       _RL  rAdvDeepFac
                0190       LOGICAL bottomDragTerms, diagFlag
71207ba845 Alis*0191 CEOP
229b6d70e7 Jean*0192 #ifdef MOM_BOUNDARY_CONSERVE
                0193       COMMON / MOM_FLUXFORM_LOCAL / uBnd, vBnd
                0194       _RL  uBnd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0195       _RL  vBnd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0196 #endif /* MOM_BOUNDARY_CONSERVE */
aea29c8517 Alis*0197
aa2d1573fa Patr*0198 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0199       kkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0200       kkey = k  + (kkey-1)*Nr
aa2d1573fa Patr*0201 #endif /* ALLOW_AUTODIFF_TAMC */
                0202
aea29c8517 Alis*0203 C     Initialise intermediate terms
722a76e285 Jean*0204       DO j=1-OLy,sNy+OLy
                0205        DO i=1-OLx,sNx+OLx
31fb0e0e6d Jean*0206         gAdd(i,j) = 0.
                0207         del2(i,j) = 0.
9293d3c672 Hajo*0208         uCf(i,j)  = 0.
                0209         vCf(i,j)  = 0.
31fb0e0e6d Jean*0210         vMT(i,j)  = 0.
aea29c8517 Alis*0211         fZon(i,j) = 0.
                0212         fMer(i,j) = 0.
722a76e285 Jean*0213         fVrUp(i,j)= 0.
                0214         fVrDw(i,j)= 0.
                0215         rTransU(i,j)= 0.
                0216         rTransV(i,j)= 0.
54ecf2c8fe Jean*0217 c       KE(i,j)     = 0.
cc254a4876 Patr*0218         hDiv(i,j)   = 0.
54ecf2c8fe Jean*0219         vort3(i,j)  = 0.
7c7b0b4a46 Alis*0220         strain(i,j) = 0.
722a76e285 Jean*0221         tension(i,j)= 0.
f59d76b0dd Ed D*0222         stretching(i,j) = 0.
06244a5e4f Jean*0223 #ifdef ALLOW_LEITH_QG
                0224         Nsquare(i,j) = 0.
                0225 #endif
722a76e285 Jean*0226         guDiss(i,j) = 0.
                0227         gvDiss(i,j) = 0.
9293d3c672 Hajo*0228 #if ( defined ALLOW_GGL90 && defined ALLOW_GGL90_LANGMUIR )
                0229 c       uRes(i,j)   = 0.
                0230 c       vRes(i,j)   = 0.
                0231 #endif
aea29c8517 Alis*0232        ENDDO
                0233       ENDDO
                0234
                0235 C--   Term by term tracer parmeters
                0236 C     o U momentum equation
                0237       uDudxFac     = afFacMom*1.
                0238       AhDudxFac    = vfFacMom*1.
                0239       vDudyFac     = afFacMom*1.
                0240       AhDudyFac    = vfFacMom*1.
31fb0e0e6d Jean*0241       wDudrFac     = afFacMom*1.
aea29c8517 Alis*0242       ArDudrFac    = vfFacMom*1.
4d2713b160 Jean*0243       mtFacU       = mtFacMom*1.
                0244       mtNHFacU     = 1.
aea29c8517 Alis*0245       fuFac        = cfFacMom*1.
                0246 C     o V momentum equation
                0247       uDvdxFac     = afFacMom*1.
                0248       AhDvdxFac    = vfFacMom*1.
                0249       vDvdyFac     = afFacMom*1.
                0250       AhDvdyFac    = vfFacMom*1.
31fb0e0e6d Jean*0251       wDvdrFac     = afFacMom*1.
aea29c8517 Alis*0252       ArDvdrFac    = vfFacMom*1.
4d2713b160 Jean*0253       mtFacV       = mtFacMom*1.
                0254       mtNHFacV     = 1.
aea29c8517 Alis*0255       fvFac        = cfFacMom*1.
722a76e285 Jean*0256
                0257       IF (implicitViscosity) THEN
                0258         ArDudrFac  = 0.
                0259         ArDvdrFac  = 0.
                0260       ENDIF
aea29c8517 Alis*0261
ea2dd4993a Jean*0262 C note: using standard stencil (no mask) results in under-estimating
                0263 C       vorticity at a no-slip boundary by a factor of 2 = sideDragFactor
                0264       IF ( no_slip_sides ) THEN
                0265         sideMaskFac = sideDragFactor
                0266       ELSE
                0267         sideMaskFac = 0. _d 0
                0268       ENDIF
                0269
99731c717f Jean*0270       IF ( selectImplicitDrag.EQ.0 .AND.
                0271      &      (  no_slip_bottom
ca594b8231 Jean*0272      &    .OR. selectBotDragQuadr.GE.0
99731c717f Jean*0273      &    .OR. bottomDragLinear.NE.0. ) ) THEN
aea29c8517 Alis*0274        bottomDragTerms=.TRUE.
                0275       ELSE
                0276        bottomDragTerms=.FALSE.
                0277       ENDIF
                0278
                0279 C--   Calculate open water fraction at vorticity points
fd362e9f7c Jean*0280       CALL MOM_CALC_HFACZ( bi,bj,k,hFacZ,r_hFacZ,myThid )
7448700841 Mart*0281 #ifdef ALLOW_AUTODIFF_TAMC
                0282 CADJ STORE hFacZ(:,:)  = comlev1_bibj_k, key = kkey, byte = isbyte
                0283 #endif
aea29c8517 Alis*0284
                0285 C---- Calculate common quantities used in both U and V equations
                0286 C     Calculate tracer cell face open areas
                0287       DO j=1-OLy,sNy+OLy
                0288        DO i=1-OLx,sNx+OLx
eaba2fd266 Jean*0289         xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
                0290      &          *drF(k)*_hFacW(i,j,k,bi,bj)
                0291         yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
                0292      &          *drF(k)*_hFacS(i,j,k,bi,bj)
fd362e9f7c Jean*0293         h0FacZ(i,j) = hFacZ(i,j)
aea29c8517 Alis*0294        ENDDO
                0295       ENDDO
fd362e9f7c Jean*0296 #ifdef NONLIN_FRSURF
                0297       IF ( momViscosity .AND. no_slip_sides
                0298      &                  .AND. nonlinFreeSurf.GT.0 ) THEN
                0299         DO j=2-OLy,sNy+OLy
                0300          DO i=2-OLx,sNx+OLx
                0301           h0FacZ(i,j) = MIN(
                0302      &       MIN( h0FacW(i,j,k,bi,bj), h0FacW(i,j-1,k,bi,bj) ),
                0303      &       MIN( h0FacS(i,j,k,bi,bj), h0FacS(i-1,j,k,bi,bj) ) )
                0304          ENDDO
                0305         ENDDO
7448700841 Mart*0306       ENDIF
                0307 # ifdef ALLOW_AUTODIFF_TAMC
                0308 CADJ STORE h0FacZ(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0309 # endif
fd362e9f7c Jean*0310 #endif /* NONLIN_FRSURF */
aea29c8517 Alis*0311
                0312 C     Make local copies of horizontal flow field
                0313       DO j=1-OLy,sNy+OLy
                0314        DO i=1-OLx,sNx+OLx
                0315         uFld(i,j) = uVel(i,j,k,bi,bj)
                0316         vFld(i,j) = vVel(i,j,k,bi,bj)
                0317        ENDDO
                0318       ENDDO
                0319
                0320 C     Calculate velocity field "volume transports" through tracer cell faces.
eaba2fd266 Jean*0321 C     anelastic: transports are scaled by rhoFacC (~ mass transport)
aea29c8517 Alis*0322       DO j=1-OLy,sNy+OLy
                0323        DO i=1-OLx,sNx+OLx
eaba2fd266 Jean*0324         uTrans(i,j) = uFld(i,j)*xA(i,j)*rhoFacC(k)
                0325         vTrans(i,j) = vFld(i,j)*yA(i,j)*rhoFacC(k)
aea29c8517 Alis*0326        ENDDO
                0327       ENDDO
                0328
fd362e9f7c Jean*0329       CALL MOM_CALC_KE( bi,bj,k,2,uFld,vFld,KE,myThid )
3c031bdc7b Jean*0330       IF ( useVariableVisc ) THEN
fd362e9f7c Jean*0331         CALL MOM_CALC_HDIV( bi,bj,k,2,uFld,vFld,hDiv,myThid )
                0332         CALL MOM_CALC_RELVORT3( bi,bj,k,uFld,vFld,hFacZ,vort3,myThid )
                0333         CALL MOM_CALC_TENSION( bi,bj,k,uFld,vFld,tension,myThid )
                0334         CALL MOM_CALC_STRAIN( bi,bj,k,uFld,vFld,hFacZ,strain,myThid )
f59d76b0dd Ed D*0335 #ifdef ALLOW_LEITH_QG
                0336         IF ( viscC2LeithQG.NE.zeroRL ) THEN
                0337           CALL MOM_VISC_QGL_STRETCH(bi,bj,k,
                0338      O                            stretching, Nsquare,
                0339      I                            myTime, myIter, myThid)
                0340           CALL MOM_VISC_QGL_LIMIT(bi,bj,k,
                0341      O                          stretching,
                0342      I                          Nsquare, uFld,vFld,vort3,
                0343      I                          myTime, myIter, myThid )
                0344         ENDIF
06244a5e4f Jean*0345 #endif /* ALLOW_LEITH_QG */
9e4f1da9cf Jean*0346         DO j=1-OLy,sNy+OLy
                0347          DO i=1-OLx,sNx+OLx
ea2dd4993a Jean*0348            IF ( hFacZ(i,j).EQ.0. ) THEN
                0349              vort3(i,j)  = sideMaskFac*vort3(i,j)
                0350              strain(i,j) = sideMaskFac*strain(i,j)
                0351            ENDIF
                0352          ENDDO
                0353         ENDDO
                0354 #ifdef ALLOW_DIAGNOSTICS
                0355         IF ( useDiagnostics ) THEN
                0356           CALL DIAGNOSTICS_FILL(hDiv,   'momHDiv ',k,1,2,bi,bj,myThid)
                0357           CALL DIAGNOSTICS_FILL(vort3,  'momVort3',k,1,2,bi,bj,myThid)
                0358           CALL DIAGNOSTICS_FILL(tension,'Tension ',k,1,2,bi,bj,myThid)
                0359           CALL DIAGNOSTICS_FILL(strain, 'Strain  ',k,1,2,bi,bj,myThid)
f59d76b0dd Ed D*0360 C     stretching will be zero, unless using QG Leith
                0361           IF ( viscC2LeithQG.NE.zeroRL ) THEN
                0362             CALL DIAGNOSTICS_FILL(stretching,
                0363      I                            'Stretch ',k,1,2,bi,bj,myThid)
                0364           ENDIF
ea2dd4993a Jean*0365         ENDIF
                0366 #endif
                0367       ENDIF
7448700841 Mart*0368 #ifdef ALLOW_AUTODIFF_TAMC
                0369 C     KE is used a lot, so we save here in order to avoid multiple
                0370 C     recomputations of this variable.
                0371 CADJ STORE KE(:,:)      = comlev1_bibj_k, key = kkey, byte = isbyte
                0372 C     Variables tension and strain are relatively cheap to recompute, so
                0373 C     we do not save them here (but saving them would clearly avoid
                0374 C     the recomputations).
                0375 cCADJ STORE tension(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0376 cCADJ STORE strain(:,:)  = comlev1_bibj_k, key = kkey, byte = isbyte
                0377 C     hDiv and vort3 are relatively expensive to recompute, so we save
                0378 C     them here.
                0379 CADJ STORE hDiv(:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
                0380 CADJ STORE vort3(:,:)   = comlev1_bibj_k, key = kkey, byte = isbyte
                0381 #endif
                0382
07e17fa184 Jean*0383 C---  First call (k=1): compute vertical adv. flux fVerUkm & fVerVkm
bd2e80b12f Jean*0384       IF (momAdvection.AND.k.EQ.1) THEN
                0385
229b6d70e7 Jean*0386 #ifdef MOM_BOUNDARY_CONSERVE
                0387         CALL MOM_UV_BOUNDARY( bi, bj, k,
                0388      I                        uVel, vVel,
                0389      O                        uBnd(1-OLx,1-OLy,k,bi,bj),
                0390      O                        vBnd(1-OLx,1-OLy,k,bi,bj),
                0391      I                        myTime, myIter, myThid )
                0392 #endif /* MOM_BOUNDARY_CONSERVE */
                0393
bd2e80b12f Jean*0394 C-    Calculate vertical transports above U & V points (West & South face):
aa2d1573fa Patr*0395
                0396 #ifdef ALLOW_AUTODIFF_TAMC
a377feeea9 Patr*0397 # ifdef NONLIN_FRSURF
                0398 #  ifndef DISABLE_RSTAR_CODE
edb6656069 Mart*0399 CADJ STORE dwtransc(:,:,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
                0400 CADJ STORE dwtransu(:,:,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
                0401 CADJ STORE dwtransv(:,:,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
a377feeea9 Patr*0402 #  endif
                0403 # endif /* NONLIN_FRSURF */
aa2d1573fa Patr*0404 #endif /* ALLOW_AUTODIFF_TAMC */
722a76e285 Jean*0405         CALL MOM_CALC_RTRANS( k, bi, bj,
                0406      O                        rTransU, rTransV,
fd362e9f7c Jean*0407      I                        myTime, myIter, myThid )
bd2e80b12f Jean*0408
                0409 C-    Free surface correction term (flux at k=1)
31fb0e0e6d Jean*0410         CALL MOM_U_ADV_WU( bi, bj, k, deepFacAdv,
                0411      I                     uVel, wVel, rTransU,
07e17fa184 Jean*0412      O                     fVerUkm, myThid )
bd2e80b12f Jean*0413
31fb0e0e6d Jean*0414         CALL MOM_V_ADV_WV( bi, bj, k, deepFacAdv,
                0415      I                     vVel, wVel, rTransV,
07e17fa184 Jean*0416      O                     fVerVkm, myThid )
bd2e80b12f Jean*0417
                0418 C---  endif momAdvection & k=1
                0419       ENDIF
                0420
                0421 C---  Calculate vertical transports (at k+1) below U & V points :
                0422       IF (momAdvection) THEN
722a76e285 Jean*0423         CALL MOM_CALC_RTRANS( k+1, bi, bj,
                0424      O                        rTransU, rTransV,
fd362e9f7c Jean*0425      I                        myTime, myIter, myThid )
bd2e80b12f Jean*0426       ENDIF
7448700841 Mart*0427 #ifdef ALLOW_AUTODIFF_TAMC
                0428 C     This saves us from calling the relatively expensive
                0429 C     mom_calc_rtrans in AD mode.
                0430 CADJ STORE rTransU(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0431 CADJ STORE rTransV(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0432 #endif
bd2e80b12f Jean*0433
229b6d70e7 Jean*0434 #ifdef MOM_BOUNDARY_CONSERVE
                0435       IF ( momAdvection .AND. k.LT.Nr ) THEN
                0436         CALL MOM_UV_BOUNDARY( bi, bj, k+1,
                0437      I                        uVel, vVel,
                0438      O                        uBnd(1-OLx,1-OLy,k+1,bi,bj),
                0439      O                        vBnd(1-OLx,1-OLy,k+1,bi,bj),
                0440      I                        myTime, myIter, myThid )
                0441       ENDIF
                0442 #endif /* MOM_BOUNDARY_CONSERVE */
                0443
05b9f17ae6 Bayl*0444       IF (momViscosity) THEN
b639c0f848 Jean*0445        DO j=1-OLy,sNy+OLy
                0446         DO i=1-OLx,sNx+OLx
                0447          viscAh_D(i,j) = viscAhD
                0448          viscAh_Z(i,j) = viscAhZ
                0449          viscA4_D(i,j) = viscA4D
                0450          viscA4_Z(i,j) = viscA4Z
                0451         ENDDO
                0452        ENDDO
                0453        IF ( useVariableVisc ) THEN
                0454         CALL MOM_CALC_VISC( bi, bj, k,
                0455      O           viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
f59d76b0dd Ed D*0456      I           hDiv, vort3, tension, strain, stretching, KE, hFacZ,
b639c0f848 Jean*0457      I           myThid )
                0458        ENDIF
7448700841 Mart*0459 #ifdef ALLOW_AUTODIFF_TAMC
                0460 # ifndef AUTODIFF_ALLOW_VISCFACADJ
                0461 C     These store directive must not be here, if you want to recompute
                0462 C     these viscosity coefficients with a modified viscFacAdj = viscFacInAd
                0463 C     because the store directives intentionally prevent the recomputation.
                0464 CADJ STORE viscAh_Z(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0465 CADJ STORE viscAh_D(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0466 CADJ STORE viscA4_Z(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0467 CADJ STORE viscA4_D(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0468 # endif /* AUTODIFF_ALLOW_VISCFACADJ */
                0469 #endif /* ALLOW_AUTODIFF_TAMC */
05b9f17ae6 Bayl*0470       ENDIF
bd2e80b12f Jean*0471
722a76e285 Jean*0472 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
aea29c8517 Alis*0473
722a76e285 Jean*0474 C---- Zonal momentum equation starts here
aea29c8517 Alis*0475
722a76e285 Jean*0476       IF (momAdvection) THEN
                0477 C---  Calculate mean fluxes (advection)   between cells for zonal flow.
aea29c8517 Alis*0478
31fb0e0e6d Jean*0479 C--   For Deep-Model (deepAtmosphere=T), no need for NHMTerm w*u/r since
                0480 C     vertical advective flux (in MOM_U_ADV_WU) is changed to W*(u*deepFacC)
                0481 C     with rAdvDeepFac factor in flux divergence
                0482
229b6d70e7 Jean*0483 #ifdef MOM_BOUNDARY_CONSERVE
                0484         CALL MOM_U_ADV_UU( bi,bj,k,uTrans,uBnd(1-OLx,1-OLy,k,bi,bj),
                0485      O                     fZon,myThid )
                0486         CALL MOM_U_ADV_VU( bi,bj,k,vTrans,uBnd(1-OLx,1-OLy,k,bi,bj),
                0487      O                     fMer,myThid )
31fb0e0e6d Jean*0488         CALL MOM_U_ADV_WU( bi, bj, k+1, deepFacAdv,
                0489      I                     uBnd, wVel, rTransU,
07e17fa184 Jean*0490      O                     fVerUkp, myThid )
229b6d70e7 Jean*0491 #else /* MOM_BOUNDARY_CONSERVE */
aea29c8517 Alis*0492 C--   Zonal flux (fZon is at east face of "u" cell)
722a76e285 Jean*0493 C     Mean flow component of zonal flux -> fZon
fd362e9f7c Jean*0494         CALL MOM_U_ADV_UU( bi,bj,k,uTrans,uFld,fZon,myThid )
aea29c8517 Alis*0495
                0496 C--   Meridional flux (fMer is at south face of "u" cell)
722a76e285 Jean*0497 C     Mean flow component of meridional flux -> fMer
fd362e9f7c Jean*0498         CALL MOM_U_ADV_VU( bi,bj,k,vTrans,uFld,fMer,myThid )
aea29c8517 Alis*0499
                0500 C--   Vertical flux (fVer is at upper face of "u" cell)
722a76e285 Jean*0501 C     Mean flow component of vertical flux (at k+1) -> fVer
31fb0e0e6d Jean*0502         CALL MOM_U_ADV_WU( bi, bj, k+1, deepFacAdv,
                0503      I                     uVel, wVel, rTransU,
07e17fa184 Jean*0504      O                     fVerUkp, myThid )
229b6d70e7 Jean*0505 #endif /* MOM_BOUNDARY_CONSERVE */
aea29c8517 Alis*0506
                0507 C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
31fb0e0e6d Jean*0508         rAdvDeepFac = wDudrFac*rkSign
                0509 #ifndef MOM_USE_OLD_DEEP_VERT_ADV
                0510         IF ( useNHMTerms ) rAdvDeepFac = rAdvDeepFac*recip_deepFacC(k)
                0511 #endif
                0512         IF ( selectMetricTerms.EQ.3 ) THEN
                0513 C-      In this case fMer contains the advective flux vTrans*(u*dxC)
                0514          DO j=jMin,jMax
                0515           DO i=iMin,iMax
                0516            gU(i,j,k,bi,bj) =
                0517      &      -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
                0518      &      *recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
                0519      &      *( ( fZon( i, j ) - fZon(i-1,j ) )*uDudxFac
                0520      &        +( fMer( i,j+1) - fMer( i, j ) )*vDudyFac
                0521      &                                        *recip_dxC(i,j,bi,bj)
                0522      &        +( fVerUkp(i,j) - fVerUkm(i,j) )*rAdvDeepFac
                0523      &       )
                0524           ENDDO
                0525          ENDDO
                0526         ELSE
                0527          DO j=jMin,jMax
                0528           DO i=iMin,iMax
                0529            gU(i,j,k,bi,bj) =
aea29c8517 Alis*0530 #ifdef OLD_UV_GEOM
31fb0e0e6d Jean*0531      &      -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
722a76e285 Jean*0532      &      ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
aea29c8517 Alis*0533 #else
31fb0e0e6d Jean*0534      &      -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
                0535      &      *recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
aea29c8517 Alis*0536 #endif
31fb0e0e6d Jean*0537      &      *( ( fZon( i, j ) - fZon(i-1,j ) )*uDudxFac
                0538      &        +( fMer( i,j+1) - fMer( i, j ) )*vDudyFac
                0539      &        +( fVerUkp(i,j) - fVerUkm(i,j) )*rAdvDeepFac
                0540      &       )
                0541           ENDDO
722a76e285 Jean*0542          ENDDO
31fb0e0e6d Jean*0543         ENDIF
aea29c8517 Alis*0544
711329b07b Jean*0545 #ifdef ALLOW_DIAGNOSTICS
                0546         IF ( useDiagnostics ) THEN
07e17fa184 Jean*0547           CALL DIAGNOSTICS_FILL( fZon,  'ADVx_Um ',k,1,2,bi,bj,myThid)
                0548           CALL DIAGNOSTICS_FILL( fMer,  'ADVy_Um ',k,1,2,bi,bj,myThid)
                0549           CALL DIAGNOSTICS_FILL(fVerUkm,'ADVrE_Um',k,1,2,bi,bj,myThid)
711329b07b Jean*0550         ENDIF
                0551 #endif
                0552
bd2e80b12f Jean*0553 #ifdef NONLIN_FRSURF
                0554 C-- account for 3.D divergence of the flow in rStar coordinate:
cdc9f269ae Patr*0555 # ifndef DISABLE_RSTAR_CODE
722a76e285 Jean*0556         IF ( select_rStar.GT.0 ) THEN
                0557          DO j=jMin,jMax
                0558           DO i=iMin,iMax
                0559            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
88b07ffa37 Jean*0560      &     - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTFreeSurf
bd2e80b12f Jean*0561      &       *uVel(i,j,k,bi,bj)
722a76e285 Jean*0562           ENDDO
                0563          ENDDO
                0564         ENDIF
                0565         IF ( select_rStar.LT.0 ) THEN
                0566          DO j=jMin,jMax
                0567           DO i=iMin,iMax
                0568            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
bd2e80b12f Jean*0569      &     - rStarDhWDt(i,j,bi,bj)*uVel(i,j,k,bi,bj)
722a76e285 Jean*0570           ENDDO
                0571          ENDDO
                0572         ENDIF
cdc9f269ae Patr*0573 # endif /* DISABLE_RSTAR_CODE */
722a76e285 Jean*0574 #endif /* NONLIN_FRSURF */
                0575
9e4f1da9cf Jean*0576 #ifdef ALLOW_ADDFLUID
                0577         IF ( selectAddFluid.GE.1 ) THEN
                0578          DO j=jMin,jMax
                0579           DO i=iMin,iMax
                0580            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
                0581      &     + uVel(i,j,k,bi,bj)*mass2rUnit*0.5 _d 0
                0582      &      *( addMass(i-1,j,k,bi,bj) + addMass(i,j,k,bi,bj) )
                0583      &      *_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)*recip_rhoFacC(k)
                0584      &      * recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)
                0585           ENDDO
                0586          ENDDO
                0587         ENDIF
                0588 #endif /* ALLOW_ADDFLUID */
                0589
722a76e285 Jean*0590       ELSE
                0591 C-    if momAdvection / else
                0592         DO j=1-OLy,sNy+OLy
                0593          DO i=1-OLx,sNx+OLx
                0594            gU(i,j,k,bi,bj) = 0. _d 0
                0595          ENDDO
bd2e80b12f Jean*0596         ENDDO
722a76e285 Jean*0597
                0598 C-    endif momAdvection.
bd2e80b12f Jean*0599       ENDIF
722a76e285 Jean*0600
                0601       IF (momViscosity) THEN
                0602 C---  Calculate eddy fluxes (dissipation) between cells for zonal flow.
                0603
31fb0e0e6d Jean*0604 C     Bi-harmonic term del^2 U -> del2
f0501c53d1 Jean*0605         IF ( useBiharmonicVisc )
fd362e9f7c Jean*0606      &  CALL MOM_U_DEL2U( bi, bj, k, uFld, hFacZ, h0FacZ,
31fb0e0e6d Jean*0607      O                    del2, myThid )
7448700841 Mart*0608 #ifdef ALLOW_AUTODIFF_TAMC
31fb0e0e6d Jean*0609 CADJ STORE del2(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
7448700841 Mart*0610 #endif
722a76e285 Jean*0611
                0612 C     Laplacian and bi-harmonic terms, Zonal  Fluxes -> fZon
31fb0e0e6d Jean*0613         CALL MOM_U_XVISCFLUX( bi,bj,k,uFld,del2,fZon,
fd362e9f7c Jean*0614      I                        viscAh_D,viscA4_D,myThid )
722a76e285 Jean*0615
                0616 C     Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
31fb0e0e6d Jean*0617         CALL MOM_U_YVISCFLUX( bi,bj,k,uFld,del2,hFacZ,fMer,
fd362e9f7c Jean*0618      I                        viscAh_Z,viscA4_Z,myThid )
722a76e285 Jean*0619
                0620 C     Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
                0621        IF (.NOT.implicitViscosity) THEN
1833b564cb Jean*0622         CALL MOM_U_RVISCFLUX( bi,bj, k, uVel,kappaRU,fVrUp,myThid )
                0623         CALL MOM_U_RVISCFLUX( bi,bj,k+1,uVel,kappaRU,fVrDw,myThid )
722a76e285 Jean*0624        ENDIF
                0625
                0626 C--   Tendency is minus divergence of the fluxes
eaba2fd266 Jean*0627 C     anelastic: hor.visc.fluxes are not scaled by rhoFac (by vert.visc.flx is)
722a76e285 Jean*0628         DO j=jMin,jMax
                0629          DO i=iMin,iMax
                0630           guDiss(i,j) =
                0631 #ifdef OLD_UV_GEOM
                0632      &     -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
                0633      &      ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
                0634 #else
                0635      &     -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
eaba2fd266 Jean*0636      &     *recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)
722a76e285 Jean*0637 #endif
eaba2fd266 Jean*0638      &     *( ( fZon(i,j  ) - fZon(i-1,j) )*AhDudxFac
                0639      &       +( fMer(i,j+1) - fMer(i,  j) )*AhDudyFac
                0640      &       +( fVrDw(i,j)  - fVrUp(i,j)  )*rkSign*ArDudrFac
                0641      &                                     *recip_rhoFacC(k)
722a76e285 Jean*0642      &     )
                0643          ENDDO
                0644         ENDDO
bd2e80b12f Jean*0645
711329b07b Jean*0646 #ifdef ALLOW_DIAGNOSTICS
                0647         IF ( useDiagnostics ) THEN
                0648           CALL DIAGNOSTICS_FILL(fZon, 'VISCx_Um',k,1,2,bi,bj,myThid)
                0649           CALL DIAGNOSTICS_FILL(fMer, 'VISCy_Um',k,1,2,bi,bj,myThid)
                0650           IF (.NOT.implicitViscosity)
                0651      &    CALL DIAGNOSTICS_FILL(fVrUp,'VISrE_Um',k,1,2,bi,bj,myThid)
                0652         ENDIF
                0653 #endif
                0654
9496c6c9ef Jean*0655 C-- No-slip and drag BCs appear as body forces in cell abutting topography
722a76e285 Jean*0656         IF (no_slip_sides) THEN
aea29c8517 Alis*0657 C-     No-slip BCs impose a drag at walls...
f0501c53d1 Jean*0658          CALL MOM_U_SIDEDRAG( bi, bj, k,
31fb0e0e6d Jean*0659      I        uFld, del2, h0FacZ,
f0501c53d1 Jean*0660      I        viscAh_Z, viscA4_Z,
                0661      I        useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
31fb0e0e6d Jean*0662      O        gAdd,
f0501c53d1 Jean*0663      I        myThid )
722a76e285 Jean*0664          DO j=jMin,jMax
                0665           DO i=iMin,iMax
31fb0e0e6d Jean*0666            guDiss(i,j) = guDiss(i,j) + gAdd(i,j)
722a76e285 Jean*0667           ENDDO
                0668          ENDDO
                0669         ENDIF
aea29c8517 Alis*0670 C-    No-slip BCs impose a drag at bottom
e2d988bd46 Jean*0671         IF ( bottomDragTerms ) THEN
79074ef66b Jean*0672 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0673 CADJ STORE KE(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
79074ef66b Jean*0674 #endif
a125ab7be7 Jean*0675          CALL MOM_U_BOTDRAG_COEFF( bi, bj, k, .TRUE.,
79074ef66b Jean*0676      I                   uFld, vFld, kappaRU, KE,
                0677      O                   cDrag,
                0678      I                   myIter, myThid )
7448700841 Mart*0679 #ifdef ALLOW_AUTODIFF_TAMC
                0680 CADJ STORE cDrag(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0681 #endif
722a76e285 Jean*0682          DO j=jMin,jMax
                0683           DO i=iMin,iMax
79074ef66b Jean*0684             guDiss(i,j) = guDiss(i,j)
                0685      &                  - cDrag(i,j)*uFld(i,j)
                0686      &                  *_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
722a76e285 Jean*0687           ENDDO
                0688          ENDDO
79074ef66b Jean*0689          IF ( useDiagnostics ) THEN
                0690           DO j=jMin,jMax
                0691            DO i=iMin,iMax
                0692             botDragU(i,j,bi,bj) = botDragU(i,j,bi,bj)
                0693      &                          - cDrag(i,j)*uFld(i,j)*rUnit2mass
                0694            ENDDO
                0695           ENDDO
                0696          ENDIF
722a76e285 Jean*0697         ENDIF
                0698
dd49642a1d Mart*0699 #ifdef ALLOW_SHELFICE
e2d988bd46 Jean*0700         IF ( useShelfIce .AND. selectImplicitDrag.EQ.0 ) THEN
                0701 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0702 CADJ STORE KE(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
e2d988bd46 Jean*0703 #endif
                0704          CALL SHELFICE_U_DRAG_COEFF( bi, bj, k, .TRUE.,
                0705      I                   uFld, vFld, kappaRU, KE,
                0706      O                   cDrag,
                0707      I                   myIter, myThid )
7448700841 Mart*0708 #ifdef ALLOW_AUTODIFF_TAMC
                0709 CADJ STORE cDrag(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0710 #endif
dd49642a1d Mart*0711          DO j=jMin,jMax
                0712           DO i=iMin,iMax
e2d988bd46 Jean*0713             guDiss(i,j) = guDiss(i,j)
                0714      &                  - cDrag(i,j)*uFld(i,j)
                0715      &                  *_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
dd49642a1d Mart*0716           ENDDO
                0717          ENDDO
                0718         ENDIF
                0719 #endif /* ALLOW_SHELFICE */
                0720
722a76e285 Jean*0721 C-    endif momViscosity
aea29c8517 Alis*0722       ENDIF
                0723
7fc445a0ef Jean*0724 C--   Forcing term (moved to timestep.F)
                0725 c     IF (momForcing)
                0726 c    &  CALL EXTERNAL_FORCING_U(
                0727 c    I     iMin,iMax,jMin,jMax,bi,bj,k,
                0728 c    I     myTime,myThid)
aea29c8517 Alis*0729
                0730 C--   Metric terms for curvilinear grid systems
31fb0e0e6d Jean*0731       diagFlag = .FALSE.
                0732 #ifdef MOM_USE_OLD_DEEP_VERT_ADV
                0733       IF ( useNHMTerms ) THEN
                0734 #else
                0735       IF ( useNHMTerms .AND. .NOT.deepAtmosphere ) THEN
                0736 #endif
                0737 C      o Non-Hydrostatic/vertical (spherical) metric terms
                0738        CALL MOM_U_METRIC_NH( bi,bj,k,uFld,wVel,vMT,myThid )
                0739        diagFlag = .TRUE.
aea29c8517 Alis*0740        DO j=jMin,jMax
                0741         DO i=iMin,iMax
31fb0e0e6d Jean*0742          gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtNHFacU*vMT(i,j)
aea29c8517 Alis*0743         ENDDO
                0744        ENDDO
31fb0e0e6d Jean*0745 #ifdef ALLOW_DIAGNOSTICS
                0746        IF ( useDiagnostics .AND.
                0747      &     ( selectMetricTerms.EQ.0 .OR. selectMetricTerms.EQ.3 ) ) THEN
                0748         CALL DIAGNOSTICS_FILL( vMT, 'Um_Metr ', k,1,2, bi,bj, myThid )
                0749        ENDIF
                0750 #endif /* ALLOW_DIAGNOSTICS */
3121bb922d Alis*0751       ENDIF
31fb0e0e6d Jean*0752       IF ( selectMetricTerms.GE.1 .AND. selectMetricTerms.NE.3 ) THEN
                0753        IF ( usingSphericalPolarGrid ) THEN
4d2713b160 Jean*0754 C      o Spherical polar grid metric terms
31fb0e0e6d Jean*0755         CALL MOM_U_METRIC_SPHERE( bi,bj,k,uFld,vFld,gAdd,myThid )
                0756         DO j=jMin,jMax
                0757          DO i=iMin,iMax
                0758           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*gAdd(i,j)
                0759          ENDDO
aea29c8517 Alis*0760         ENDDO
31fb0e0e6d Jean*0761        ENDIF
                0762        IF ( usingCylindricalGrid ) THEN
4d2713b160 Jean*0763 C      o Cylindrical grid metric terms
31fb0e0e6d Jean*0764         CALL MOM_U_METRIC_CYLINDER( bi,bj,k,uFld,vFld,gAdd,myThid )
                0765         DO j=jMin,jMax
                0766          DO i=iMin,iMax
                0767           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*gAdd(i,j)
                0768          ENDDO
4d2713b160 Jean*0769         ENDDO
31fb0e0e6d Jean*0770        ENDIF
                0771 #ifdef ALLOW_DIAGNOSTICS
                0772        IF ( useDiagnostics ) THEN
                0773         IF ( diagFlag ) THEN
                0774          DO j=jMin,jMax
                0775           DO i=iMin,iMax
                0776            gAdd(i,j) = gAdd(i,j) + vMT(i,j)
                0777           ENDDO
                0778          ENDDO
                0779         ENDIF
                0780         CALL DIAGNOSTICS_FILL( gAdd, 'Um_Metr ', k,1,2, bi,bj, myThid )
                0781        ENDIF
                0782 #endif /* ALLOW_DIAGNOSTICS */
aea29c8517 Alis*0783       ENDIF
                0784
722a76e285 Jean*0785 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
aea29c8517 Alis*0786
                0787 C---- Meridional momentum equation starts here
                0788
722a76e285 Jean*0789       IF (momAdvection) THEN
229b6d70e7 Jean*0790
31fb0e0e6d Jean*0791 C--   For Deep-Model (deepAtmosphere=T), no need for NHMTerm w*v/r since
                0792 C     vertical advective flux (in MOM_V_ADV_WV) is changed to W*(v*deepFacC)
                0793 C     with rAdvDeepFac factor in flux divergence
                0794
229b6d70e7 Jean*0795 #ifdef MOM_BOUNDARY_CONSERVE
                0796         CALL MOM_V_ADV_UV( bi,bj,k,uTrans,vBnd(1-OLx,1-OLy,k,bi,bj),
                0797      O                     fZon,myThid )
                0798         CALL MOM_V_ADV_VV( bi,bj,k,vTrans,vBnd(1-OLx,1-OLy,k,bi,bj),
                0799      O                     fMer,myThid )
31fb0e0e6d Jean*0800         CALL MOM_V_ADV_WV( bi, bj, k+1, deepFacAdv,
                0801      I                     vBnd, wVel, rTransV,
07e17fa184 Jean*0802      O                     fVerVkp, myThid )
229b6d70e7 Jean*0803 #else /* MOM_BOUNDARY_CONSERVE */
722a76e285 Jean*0804 C---  Calculate mean fluxes (advection)   between cells for meridional flow.
                0805 C     Mean flow component of zonal flux -> fZon
07e17fa184 Jean*0806         CALL MOM_V_ADV_UV( bi,bj,k,uTrans,vFld,fZon,myThid )
aea29c8517 Alis*0807
                0808 C--   Meridional flux (fMer is at north face of "v" cell)
722a76e285 Jean*0809 C     Mean flow component of meridional flux -> fMer
07e17fa184 Jean*0810         CALL MOM_V_ADV_VV( bi,bj,k,vTrans,vFld,fMer,myThid )
aea29c8517 Alis*0811
                0812 C--   Vertical flux (fVer is at upper face of "v" cell)
722a76e285 Jean*0813 C     Mean flow component of vertical flux (at k+1) -> fVerV
31fb0e0e6d Jean*0814         CALL MOM_V_ADV_WV( bi, bj, k+1, deepFacAdv,
                0815      I                     vVel, wVel, rTransV,
07e17fa184 Jean*0816      O                     fVerVkp, myThid )
229b6d70e7 Jean*0817 #endif /* MOM_BOUNDARY_CONSERVE */
aea29c8517 Alis*0818
                0819 C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
31fb0e0e6d Jean*0820         rAdvDeepFac = wDvdrFac*rkSign
                0821 #ifndef MOM_USE_OLD_DEEP_VERT_ADV
                0822         IF ( useNHMTerms ) rAdvDeepFac = rAdvDeepFac*recip_deepFacC(k)
                0823 #endif
722a76e285 Jean*0824         DO j=jMin,jMax
                0825          DO i=iMin,iMax
                0826           gV(i,j,k,bi,bj) =
aea29c8517 Alis*0827 #ifdef OLD_UV_GEOM
722a76e285 Jean*0828      &     -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
                0829      &      ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
aea29c8517 Alis*0830 #else
722a76e285 Jean*0831      &     -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
eaba2fd266 Jean*0832      &     *recip_rAs(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
aea29c8517 Alis*0833 #endif
07e17fa184 Jean*0834      &     *( ( fZon(i+1,j)  - fZon(i,j  )  )*uDvdxFac
                0835      &       +( fMer(i,  j)  - fMer(i,j-1)  )*vDvdyFac
31fb0e0e6d Jean*0836      &       +( fVerVkp(i,j) - fVerVkm(i,j) )*rAdvDeepFac
722a76e285 Jean*0837      &     )
711329b07b Jean*0838          ENDDO
                0839         ENDDO
                0840
                0841 #ifdef ALLOW_DIAGNOSTICS
                0842         IF ( useDiagnostics ) THEN
07e17fa184 Jean*0843           CALL DIAGNOSTICS_FILL( fZon,  'ADVx_Vm ',k,1,2,bi,bj,myThid)
                0844           CALL DIAGNOSTICS_FILL( fMer,  'ADVy_Vm ',k,1,2,bi,bj,myThid)
                0845           CALL DIAGNOSTICS_FILL(fVerVkm,'ADVrE_Vm',k,1,2,bi,bj,myThid)
711329b07b Jean*0846         ENDIF
                0847 #endif
aea29c8517 Alis*0848
bd2e80b12f Jean*0849 #ifdef NONLIN_FRSURF
                0850 C-- account for 3.D divergence of the flow in rStar coordinate:
cdc9f269ae Patr*0851 # ifndef DISABLE_RSTAR_CODE
722a76e285 Jean*0852         IF ( select_rStar.GT.0 ) THEN
                0853          DO j=jMin,jMax
                0854           DO i=iMin,iMax
                0855            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
88b07ffa37 Jean*0856      &     - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTFreeSurf
bd2e80b12f Jean*0857      &       *vVel(i,j,k,bi,bj)
722a76e285 Jean*0858           ENDDO
                0859          ENDDO
                0860         ENDIF
                0861         IF ( select_rStar.LT.0 ) THEN
                0862          DO j=jMin,jMax
                0863           DO i=iMin,iMax
                0864            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
bd2e80b12f Jean*0865      &     - rStarDhSDt(i,j,bi,bj)*vVel(i,j,k,bi,bj)
722a76e285 Jean*0866           ENDDO
                0867          ENDDO
                0868         ENDIF
cdc9f269ae Patr*0869 # endif /* DISABLE_RSTAR_CODE */
722a76e285 Jean*0870 #endif /* NONLIN_FRSURF */
                0871
9e4f1da9cf Jean*0872 #ifdef ALLOW_ADDFLUID
                0873         IF ( selectAddFluid.GE.1 ) THEN
                0874          DO j=jMin,jMax
                0875           DO i=iMin,iMax
                0876            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
                0877      &     + vVel(i,j,k,bi,bj)*mass2rUnit*0.5 _d 0
                0878      &      *( addMass(i,j-1,k,bi,bj) + addMass(i,j,k,bi,bj) )
                0879      &      *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)*recip_rhoFacC(k)
                0880      &      * recip_rAs(i,j,bi,bj)*recip_deepFac2C(k)
                0881           ENDDO
                0882          ENDDO
                0883         ENDIF
                0884 #endif /* ALLOW_ADDFLUID */
                0885
722a76e285 Jean*0886       ELSE
                0887 C-    if momAdvection / else
                0888         DO j=1-OLy,sNy+OLy
                0889          DO i=1-OLx,sNx+OLx
                0890            gV(i,j,k,bi,bj) = 0. _d 0
                0891          ENDDO
bd2e80b12f Jean*0892         ENDDO
722a76e285 Jean*0893
                0894 C-    endif momAdvection.
bd2e80b12f Jean*0895       ENDIF
722a76e285 Jean*0896
                0897       IF (momViscosity) THEN
                0898 C---  Calculate eddy fluxes (dissipation) between cells for meridional flow.
31fb0e0e6d Jean*0899 C     Bi-harmonic term del^2 V -> del2
f0501c53d1 Jean*0900         IF ( useBiharmonicVisc )
fd362e9f7c Jean*0901      &  CALL MOM_V_DEL2V( bi, bj, k, vFld, hFacZ, h0FacZ,
31fb0e0e6d Jean*0902      O                    del2, myThid )
7448700841 Mart*0903 #ifdef ALLOW_AUTODIFF_TAMC
31fb0e0e6d Jean*0904 CADJ STORE del2(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
7448700841 Mart*0905 #endif
722a76e285 Jean*0906
                0907 C     Laplacian and bi-harmonic terms, Zonal  Fluxes -> fZon
31fb0e0e6d Jean*0908         CALL MOM_V_XVISCFLUX( bi,bj,k,vFld,del2,hFacZ,fZon,
fd362e9f7c Jean*0909      I                        viscAh_Z,viscA4_Z,myThid )
722a76e285 Jean*0910
                0911 C     Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
31fb0e0e6d Jean*0912         CALL MOM_V_YVISCFLUX( bi,bj,k,vFld,del2,fMer,
fd362e9f7c Jean*0913      I                        viscAh_D,viscA4_D,myThid )
722a76e285 Jean*0914
                0915 C     Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
                0916        IF (.NOT.implicitViscosity) THEN
fd362e9f7c Jean*0917         CALL MOM_V_RVISCFLUX( bi,bj, k, vVel,KappaRV,fVrUp,myThid )
                0918         CALL MOM_V_RVISCFLUX( bi,bj,k+1,vVel,KappaRV,fVrDw,myThid )
722a76e285 Jean*0919        ENDIF
                0920
                0921 C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
eaba2fd266 Jean*0922 C     anelastic: hor.visc.fluxes are not scaled by rhoFac (by vert.visc.flx is)
722a76e285 Jean*0923         DO j=jMin,jMax
                0924          DO i=iMin,iMax
                0925           gvDiss(i,j) =
                0926 #ifdef OLD_UV_GEOM
                0927      &     -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
                0928      &      ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
                0929 #else
                0930      &     -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
eaba2fd266 Jean*0931      &      *recip_rAs(i,j,bi,bj)*recip_deepFac2C(k)
722a76e285 Jean*0932 #endif
eaba2fd266 Jean*0933      &     *( ( fZon(i+1,j)  - fZon(i,j  ) )*AhDvdxFac
                0934      &       +( fMer(i,  j)  - fMer(i,j-1) )*AhDvdyFac
                0935      &       +( fVrDw(i,j)   - fVrUp(i,j) )*rkSign*ArDvdrFac
                0936      &                                     *recip_rhoFacC(k)
722a76e285 Jean*0937      &     )
                0938          ENDDO
                0939         ENDDO
bd2e80b12f Jean*0940
711329b07b Jean*0941 #ifdef ALLOW_DIAGNOSTICS
                0942         IF ( useDiagnostics ) THEN
                0943           CALL DIAGNOSTICS_FILL(fZon, 'VISCx_Vm',k,1,2,bi,bj,myThid)
                0944           CALL DIAGNOSTICS_FILL(fMer, 'VISCy_Vm',k,1,2,bi,bj,myThid)
                0945           IF (.NOT.implicitViscosity)
                0946      &    CALL DIAGNOSTICS_FILL(fVrUp,'VISrE_Vm',k,1,2,bi,bj,myThid)
                0947         ENDIF
                0948 #endif
                0949
9496c6c9ef Jean*0950 C-- No-slip and drag BCs appear as body forces in cell abutting topography
dd49642a1d Mart*0951         IF (no_slip_sides) THEN
aea29c8517 Alis*0952 C-     No-slip BCs impose a drag at walls...
f0501c53d1 Jean*0953          CALL MOM_V_SIDEDRAG( bi, bj, k,
31fb0e0e6d Jean*0954      I        vFld, del2, h0FacZ,
fd362e9f7c Jean*0955      I        viscAh_Z, viscA4_Z,
f0501c53d1 Jean*0956      I        useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
31fb0e0e6d Jean*0957      O        gAdd,
f0501c53d1 Jean*0958      I        myThid )
722a76e285 Jean*0959          DO j=jMin,jMax
                0960           DO i=iMin,iMax
31fb0e0e6d Jean*0961            gvDiss(i,j) = gvDiss(i,j) + gAdd(i,j)
722a76e285 Jean*0962           ENDDO
                0963          ENDDO
                0964         ENDIF
aea29c8517 Alis*0965 C-    No-slip BCs impose a drag at bottom
e2d988bd46 Jean*0966         IF ( bottomDragTerms ) THEN
79074ef66b Jean*0967 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0968 CADJ STORE KE(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
79074ef66b Jean*0969 #endif
a125ab7be7 Jean*0970          CALL MOM_V_BOTDRAG_COEFF( bi, bj, k, .TRUE.,
79074ef66b Jean*0971      I                   uFld, vFld, kappaRV, KE,
                0972      O                   cDrag,
                0973      I                   myIter, myThid )
7448700841 Mart*0974 #ifdef ALLOW_AUTODIFF_TAMC
                0975 CADJ STORE cDrag(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0976 #endif
722a76e285 Jean*0977          DO j=jMin,jMax
                0978           DO i=iMin,iMax
79074ef66b Jean*0979             gvDiss(i,j) = gvDiss(i,j)
                0980      &                  - cDrag(i,j)*vFld(i,j)
                0981      &                  *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
722a76e285 Jean*0982           ENDDO
                0983          ENDDO
79074ef66b Jean*0984          IF ( useDiagnostics ) THEN
                0985           DO j=jMin,jMax
                0986            DO i=iMin,iMax
                0987             botDragV(i,j,bi,bj) = botDragV(i,j,bi,bj)
                0988      &                          - cDrag(i,j)*vFld(i,j)*rUnit2mass
                0989            ENDDO
                0990           ENDDO
                0991          ENDIF
722a76e285 Jean*0992         ENDIF
                0993
dd49642a1d Mart*0994 #ifdef ALLOW_SHELFICE
e2d988bd46 Jean*0995         IF ( useShelfIce .AND. selectImplicitDrag.EQ.0 ) THEN
                0996 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0997 CADJ STORE KE(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
e2d988bd46 Jean*0998 #endif
                0999          CALL SHELFICE_V_DRAG_COEFF( bi, bj, k, .TRUE.,
                1000      I                   uFld, vFld, kappaRV, KE,
                1001      O                   cDrag,
                1002      I                   myIter, myThid )
7448700841 Mart*1003 #ifdef ALLOW_AUTODIFF_TAMC
                1004 CADJ STORE cDrag(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                1005 #endif
dd49642a1d Mart*1006          DO j=jMin,jMax
                1007           DO i=iMin,iMax
e2d988bd46 Jean*1008             gvDiss(i,j) = gvDiss(i,j)
                1009      &                  - cDrag(i,j)*vFld(i,j)
                1010      &                  *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
dd49642a1d Mart*1011           ENDDO
                1012          ENDDO
                1013         ENDIF
                1014 #endif /* ALLOW_SHELFICE */
                1015
722a76e285 Jean*1016 C-    endif momViscosity
aea29c8517 Alis*1017       ENDIF
                1018
7fc445a0ef Jean*1019 C--   Forcing term (moved to timestep.F)
                1020 c     IF (momForcing)
                1021 c    & CALL EXTERNAL_FORCING_V(
                1022 c    I     iMin,iMax,jMin,jMax,bi,bj,k,
                1023 c    I     myTime,myThid)
aea29c8517 Alis*1024
                1025 C--   Metric terms for curvilinear grid systems
31fb0e0e6d Jean*1026       diagFlag = .FALSE.
                1027 #ifdef MOM_USE_OLD_DEEP_VERT_ADV
                1028       IF ( useNHMTerms ) THEN
                1029 #else
                1030       IF ( useNHMTerms .AND. .NOT.deepAtmosphere ) THEN
                1031 #endif
                1032 C      o Non-Hydrostatic/vertical (spherical) metric terms
                1033        CALL MOM_V_METRIC_NH( bi,bj,k,vFld,wVel,vMT,myThid )
                1034        diagFlag = .TRUE.
aea29c8517 Alis*1035        DO j=jMin,jMax
                1036         DO i=iMin,iMax
31fb0e0e6d Jean*1037          gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtNHFacV*vMT(i,j)
aea29c8517 Alis*1038         ENDDO
                1039        ENDDO
31fb0e0e6d Jean*1040 #ifdef ALLOW_DIAGNOSTICS
                1041        IF ( useDiagnostics .AND. selectMetricTerms.EQ.0 ) THEN
                1042         CALL DIAGNOSTICS_FILL( vMT, 'Vm_Metr ', k,1,2, bi,bj, myThid )
                1043        ENDIF
                1044 #endif /* ALLOW_DIAGNOSTICS */
3121bb922d Alis*1045       ENDIF
31fb0e0e6d Jean*1046       IF ( selectMetricTerms.GE.1 ) THEN
                1047        IF ( usingSphericalPolarGrid ) THEN
4d2713b160 Jean*1048 C      o Spherical polar grid metric terms
31fb0e0e6d Jean*1049         CALL MOM_V_METRIC_SPHERE( bi,bj,k,uFld,gAdd,myThid )
                1050         DO j=jMin,jMax
                1051          DO i=iMin,iMax
                1052           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*gAdd(i,j)
                1053          ENDDO
aea29c8517 Alis*1054         ENDDO
31fb0e0e6d Jean*1055        ENDIF
                1056        IF ( usingCylindricalGrid ) THEN
4d2713b160 Jean*1057 C      o Cylindrical grid metric terms
31fb0e0e6d Jean*1058         CALL MOM_V_METRIC_CYLINDER( bi,bj,k,uFld,vFld,gAdd,myThid )
                1059         DO j=jMin,jMax
                1060          DO i=iMin,iMax
                1061           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*gAdd(i,j)
                1062          ENDDO
4d2713b160 Jean*1063         ENDDO
31fb0e0e6d Jean*1064        ENDIF
                1065 #ifdef ALLOW_DIAGNOSTICS
                1066        IF ( useDiagnostics ) THEN
                1067         IF ( diagFlag ) THEN
                1068          DO j=jMin,jMax
                1069           DO i=iMin,iMax
                1070            gAdd(i,j) = gAdd(i,j) + vMT(i,j)
                1071           ENDDO
                1072          ENDDO
                1073         ENDIF
                1074         CALL DIAGNOSTICS_FILL( gAdd, 'Vm_Metr ', k,1,2, bi,bj, myThid )
                1075        ENDIF
                1076 #endif /* ALLOW_DIAGNOSTICS */
0ac260a803 Andr*1077       ENDIF
aea29c8517 Alis*1078
722a76e285 Jean*1079 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
aea29c8517 Alis*1080
c831f9444b Jean*1081 C--   Coriolis term (call to CD_CODE_SCHEME has been moved to timestep.F)
9293d3c672 Hajo*1082       IF ( .NOT.useCDscheme ) THEN
                1083 #if ( defined ALLOW_GGL90 && defined ALLOW_GGL90_LANGMUIR )
                1084        IF ( useLANGMUIR ) THEN
                1085         CALL GGL90_ADD_STOKESDRIFT(
                1086      O                 uRes, vRes,
                1087      I                 uFld, vFld, k, bi, bj, myThid )
                1088         CALL MOM_U_CORIOLIS( bi,bj,k,vRes,uCf,myThid )
                1089         CALL MOM_V_CORIOLIS( bi,bj,k,uRes,vCf,myThid )
                1090        ELSE
                1091 #endif
                1092         CALL MOM_U_CORIOLIS( bi,bj,k,vFld,uCf,myThid )
                1093         CALL MOM_V_CORIOLIS( bi,bj,k,uFld,vCf,myThid )
                1094 #if ( defined ALLOW_GGL90 && defined ALLOW_GGL90_LANGMUIR )
                1095        ENDIF
711329b07b Jean*1096 #endif
7fc445a0ef Jean*1097         DO j=jMin,jMax
                1098          DO i=iMin,iMax
9293d3c672 Hajo*1099           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + fuFac*uCf(i,j)
                1100           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + fvFac*vCf(i,j)
7fc445a0ef Jean*1101          ENDDO
                1102         ENDDO
711329b07b Jean*1103 #ifdef ALLOW_DIAGNOSTICS
9293d3c672 Hajo*1104         IF ( useDiagnostics ) THEN
                1105           CALL DIAGNOSTICS_FILL( uCf,'Um_Cori ',k,1,2,bi,bj,myThid )
                1106           CALL DIAGNOSTICS_FILL( vCf,'Vm_Cori ',k,1,2,bi,bj,myThid )
                1107         ENDIF
711329b07b Jean*1108 #endif
7fc445a0ef Jean*1109       ENDIF
                1110
3daafce20b Jean*1111 C--   3.D Coriolis term (horizontal momentum, Eastward component: -fprime*w)
427e24e121 Jean*1112       IF ( select3dCoriScheme.GE.1 ) THEN
9293d3c672 Hajo*1113         CALL MOM_U_CORIOLIS_NH( bi,bj,k,wVel,uCf,myThid )
ac8c612649 Jean*1114         DO j=jMin,jMax
                1115          DO i=iMin,iMax
9293d3c672 Hajo*1116           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + fuFac*uCf(i,j)
ac8c612649 Jean*1117          ENDDO
c5990018f4 Alis*1118         ENDDO
427e24e121 Jean*1119        IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
ac8c612649 Jean*1120 C-     presently, non zero angleSinC array only supported with Curvilinear-Grid
9293d3c672 Hajo*1121         CALL MOM_V_CORIOLIS_NH( bi,bj,k,wVel,vCf,myThid )
ac8c612649 Jean*1122         DO j=jMin,jMax
                1123          DO i=iMin,iMax
9293d3c672 Hajo*1124           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + fvFac*vCf(i,j)
ac8c612649 Jean*1125          ENDDO
                1126         ENDDO
                1127        ENDIF
c5990018f4 Alis*1128       ENDIF
aea29c8517 Alis*1129
722a76e285 Jean*1130 C--   Set du/dt & dv/dt on boundaries to zero
                1131       DO j=jMin,jMax
                1132        DO i=iMin,iMax
                1133         gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
                1134         guDiss(i,j)     = guDiss(i,j)    *_maskW(i,j,k,bi,bj)
                1135         gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
                1136         gvDiss(i,j)     = gvDiss(i,j)    *_maskS(i,j,k,bi,bj)
                1137        ENDDO
                1138       ENDDO
                1139
711329b07b Jean*1140 #ifdef ALLOW_DIAGNOSTICS
                1141       IF ( useDiagnostics ) THEN
b114b87ee5 Bayl*1142         CALL DIAGNOSTICS_FILL(KE,    'momKE   ',k,1,2,bi,bj,myThid)
9e4f1da9cf Jean*1143         CALL DIAGNOSTICS_FILL(gU(1-OLx,1-OLy,k,bi,bj),
711329b07b Jean*1144      &                               'Um_Advec',k,1,2,bi,bj,myThid)
9e4f1da9cf Jean*1145         CALL DIAGNOSTICS_FILL(gV(1-OLx,1-OLy,k,bi,bj),
711329b07b Jean*1146      &                               'Vm_Advec',k,1,2,bi,bj,myThid)
                1147       ENDIF
                1148 #endif /* ALLOW_DIAGNOSTICS */
                1149
aea29c8517 Alis*1150       RETURN
                1151       END