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
                0111 C  vF                   :: viscous flux
                0112 C  v4F                  :: bi-harmonic viscous flux
9293d3c672 Hajo*0113 C  uCf,vCf              :: Coriolis acceleration
71207ba845 Alis*0114 C  mT                   :: Metric terms
                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
9293d3c672 Hajo*0123       _RL  vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71207ba845 Alis*0124       _RL v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
9293d3c672 Hajo*0125       _RL uCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0126       _RL vCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0127       _RL  mT(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
                0174       _RL  rVelDudrFac
                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
                0183       _RL  rVelDvdrFac
                0184       _RL  ArDvdrFac
                0185       _RL  fvFac
                0186       _RL  mtFacV
4d2713b160 Jean*0187       _RL  mtNHFacV
ea2dd4993a Jean*0188       _RL  sideMaskFac
427e24e121 Jean*0189       LOGICAL bottomDragTerms, metricTerms
71207ba845 Alis*0190 CEOP
229b6d70e7 Jean*0191 #ifdef MOM_BOUNDARY_CONSERVE
                0192       COMMON / MOM_FLUXFORM_LOCAL / uBnd, vBnd
                0193       _RL  uBnd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0194       _RL  vBnd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0195 #endif /* MOM_BOUNDARY_CONSERVE */
aea29c8517 Alis*0196
aa2d1573fa Patr*0197 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0198       kkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0199       kkey = k  + (kkey-1)*Nr
aa2d1573fa Patr*0200 #endif /* ALLOW_AUTODIFF_TAMC */
                0201
aea29c8517 Alis*0202 C     Initialise intermediate terms
722a76e285 Jean*0203       DO j=1-OLy,sNy+OLy
                0204        DO i=1-OLx,sNx+OLx
aea29c8517 Alis*0205         vF(i,j)   = 0.
                0206         v4F(i,j)  = 0.
9293d3c672 Hajo*0207         uCf(i,j)  = 0.
                0208         vCf(i,j)  = 0.
aea29c8517 Alis*0209         mT(i,j)   = 0.
                0210         fZon(i,j) = 0.
                0211         fMer(i,j) = 0.
722a76e285 Jean*0212         fVrUp(i,j)= 0.
                0213         fVrDw(i,j)= 0.
                0214         rTransU(i,j)= 0.
                0215         rTransV(i,j)= 0.
54ecf2c8fe Jean*0216 c       KE(i,j)     = 0.
cc254a4876 Patr*0217         hDiv(i,j)   = 0.
54ecf2c8fe Jean*0218         vort3(i,j)  = 0.
7c7b0b4a46 Alis*0219         strain(i,j) = 0.
722a76e285 Jean*0220         tension(i,j)= 0.
f59d76b0dd Ed D*0221         stretching(i,j) = 0.
06244a5e4f Jean*0222 #ifdef ALLOW_LEITH_QG
                0223         Nsquare(i,j) = 0.
                0224 #endif
722a76e285 Jean*0225         guDiss(i,j) = 0.
                0226         gvDiss(i,j) = 0.
9293d3c672 Hajo*0227 #if ( defined ALLOW_GGL90 && defined ALLOW_GGL90_LANGMUIR )
                0228 c       uRes(i,j)   = 0.
                0229 c       vRes(i,j)   = 0.
                0230 #endif
aea29c8517 Alis*0231        ENDDO
                0232       ENDDO
                0233
                0234 C--   Term by term tracer parmeters
                0235 C     o U momentum equation
                0236       uDudxFac     = afFacMom*1.
                0237       AhDudxFac    = vfFacMom*1.
                0238       vDudyFac     = afFacMom*1.
                0239       AhDudyFac    = vfFacMom*1.
                0240       rVelDudrFac  = afFacMom*1.
                0241       ArDudrFac    = vfFacMom*1.
4d2713b160 Jean*0242       mtFacU       = mtFacMom*1.
                0243       mtNHFacU     = 1.
aea29c8517 Alis*0244       fuFac        = cfFacMom*1.
                0245 C     o V momentum equation
                0246       uDvdxFac     = afFacMom*1.
                0247       AhDvdxFac    = vfFacMom*1.
                0248       vDvdyFac     = afFacMom*1.
                0249       AhDvdyFac    = vfFacMom*1.
                0250       rVelDvdrFac  = afFacMom*1.
                0251       ArDvdrFac    = vfFacMom*1.
4d2713b160 Jean*0252       mtFacV       = mtFacMom*1.
                0253       mtNHFacV     = 1.
aea29c8517 Alis*0254       fvFac        = cfFacMom*1.
722a76e285 Jean*0255
427e24e121 Jean*0256       metricTerms  = selectMetricTerms.GE.1
722a76e285 Jean*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)
722a76e285 Jean*0410         CALL MOM_U_ADV_WU( bi,bj,k,uVel,wVel,rTransU,
07e17fa184 Jean*0411      O                     fVerUkm, myThid )
bd2e80b12f Jean*0412
722a76e285 Jean*0413         CALL MOM_V_ADV_WV( bi,bj,k,vVel,wVel,rTransV,
07e17fa184 Jean*0414      O                     fVerVkm, myThid )
bd2e80b12f Jean*0415
                0416 C---  endif momAdvection & k=1
                0417       ENDIF
                0418
                0419 C---  Calculate vertical transports (at k+1) below U & V points :
                0420       IF (momAdvection) THEN
722a76e285 Jean*0421         CALL MOM_CALC_RTRANS( k+1, bi, bj,
                0422      O                        rTransU, rTransV,
fd362e9f7c Jean*0423      I                        myTime, myIter, myThid )
bd2e80b12f Jean*0424       ENDIF
7448700841 Mart*0425 #ifdef ALLOW_AUTODIFF_TAMC
                0426 C     This saves us from calling the relatively expensive
                0427 C     mom_calc_rtrans in AD mode.
                0428 CADJ STORE rTransU(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0429 CADJ STORE rTransV(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0430 #endif
bd2e80b12f Jean*0431
229b6d70e7 Jean*0432 #ifdef MOM_BOUNDARY_CONSERVE
                0433       IF ( momAdvection .AND. k.LT.Nr ) THEN
                0434         CALL MOM_UV_BOUNDARY( bi, bj, k+1,
                0435      I                        uVel, vVel,
                0436      O                        uBnd(1-OLx,1-OLy,k+1,bi,bj),
                0437      O                        vBnd(1-OLx,1-OLy,k+1,bi,bj),
                0438      I                        myTime, myIter, myThid )
                0439       ENDIF
                0440 #endif /* MOM_BOUNDARY_CONSERVE */
                0441
05b9f17ae6 Bayl*0442       IF (momViscosity) THEN
b639c0f848 Jean*0443        DO j=1-OLy,sNy+OLy
                0444         DO i=1-OLx,sNx+OLx
                0445          viscAh_D(i,j) = viscAhD
                0446          viscAh_Z(i,j) = viscAhZ
                0447          viscA4_D(i,j) = viscA4D
                0448          viscA4_Z(i,j) = viscA4Z
                0449         ENDDO
                0450        ENDDO
                0451        IF ( useVariableVisc ) THEN
                0452         CALL MOM_CALC_VISC( bi, bj, k,
                0453      O           viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
f59d76b0dd Ed D*0454      I           hDiv, vort3, tension, strain, stretching, KE, hFacZ,
b639c0f848 Jean*0455      I           myThid )
                0456        ENDIF
7448700841 Mart*0457 #ifdef ALLOW_AUTODIFF_TAMC
                0458 # ifndef AUTODIFF_ALLOW_VISCFACADJ
                0459 C     These store directive must not be here, if you want to recompute
                0460 C     these viscosity coefficients with a modified viscFacAdj = viscFacInAd
                0461 C     because the store directives intentionally prevent the recomputation.
                0462 CADJ STORE viscAh_Z(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0463 CADJ STORE viscAh_D(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0464 CADJ STORE viscA4_Z(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0465 CADJ STORE viscA4_D(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0466 # endif /* AUTODIFF_ALLOW_VISCFACADJ */
                0467 #endif /* ALLOW_AUTODIFF_TAMC */
05b9f17ae6 Bayl*0468       ENDIF
bd2e80b12f Jean*0469
722a76e285 Jean*0470 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
aea29c8517 Alis*0471
722a76e285 Jean*0472 C---- Zonal momentum equation starts here
aea29c8517 Alis*0473
722a76e285 Jean*0474       IF (momAdvection) THEN
                0475 C---  Calculate mean fluxes (advection)   between cells for zonal flow.
aea29c8517 Alis*0476
229b6d70e7 Jean*0477 #ifdef MOM_BOUNDARY_CONSERVE
                0478         CALL MOM_U_ADV_UU( bi,bj,k,uTrans,uBnd(1-OLx,1-OLy,k,bi,bj),
                0479      O                     fZon,myThid )
                0480         CALL MOM_U_ADV_VU( bi,bj,k,vTrans,uBnd(1-OLx,1-OLy,k,bi,bj),
                0481      O                     fMer,myThid )
                0482         CALL MOM_U_ADV_WU(
                0483      I                     bi,bj,k+1,uBnd,wVel,rTransU,
07e17fa184 Jean*0484      O                     fVerUkp, myThid )
229b6d70e7 Jean*0485 #else /* MOM_BOUNDARY_CONSERVE */
aea29c8517 Alis*0486 C--   Zonal flux (fZon is at east face of "u" cell)
722a76e285 Jean*0487 C     Mean flow component of zonal flux -> fZon
fd362e9f7c Jean*0488         CALL MOM_U_ADV_UU( bi,bj,k,uTrans,uFld,fZon,myThid )
aea29c8517 Alis*0489
                0490 C--   Meridional flux (fMer is at south face of "u" cell)
722a76e285 Jean*0491 C     Mean flow component of meridional flux -> fMer
fd362e9f7c Jean*0492         CALL MOM_U_ADV_VU( bi,bj,k,vTrans,uFld,fMer,myThid )
aea29c8517 Alis*0493
                0494 C--   Vertical flux (fVer is at upper face of "u" cell)
722a76e285 Jean*0495 C     Mean flow component of vertical flux (at k+1) -> fVer
                0496         CALL MOM_U_ADV_WU(
                0497      I                     bi,bj,k+1,uVel,wVel,rTransU,
07e17fa184 Jean*0498      O                     fVerUkp, myThid )
229b6d70e7 Jean*0499 #endif /* MOM_BOUNDARY_CONSERVE */
aea29c8517 Alis*0500
                0501 C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
722a76e285 Jean*0502         DO j=jMin,jMax
                0503          DO i=iMin,iMax
                0504           gU(i,j,k,bi,bj) =
aea29c8517 Alis*0505 #ifdef OLD_UV_GEOM
722a76e285 Jean*0506      &     -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
                0507      &      ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
aea29c8517 Alis*0508 #else
722a76e285 Jean*0509      &     -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
eaba2fd266 Jean*0510      &     *recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
aea29c8517 Alis*0511 #endif
07e17fa184 Jean*0512      &     *( ( fZon(i,j  )  - fZon(i-1,j)  )*uDudxFac
                0513      &       +( fMer(i,j+1)  - fMer(i,  j)  )*vDudyFac
                0514      &       +( fVerUkp(i,j) - fVerUkm(i,j) )*rkSign*rVelDudrFac
722a76e285 Jean*0515      &     )
                0516          ENDDO
                0517         ENDDO
aea29c8517 Alis*0518
711329b07b Jean*0519 #ifdef ALLOW_DIAGNOSTICS
                0520         IF ( useDiagnostics ) THEN
07e17fa184 Jean*0521           CALL DIAGNOSTICS_FILL( fZon,  'ADVx_Um ',k,1,2,bi,bj,myThid)
                0522           CALL DIAGNOSTICS_FILL( fMer,  'ADVy_Um ',k,1,2,bi,bj,myThid)
                0523           CALL DIAGNOSTICS_FILL(fVerUkm,'ADVrE_Um',k,1,2,bi,bj,myThid)
711329b07b Jean*0524         ENDIF
                0525 #endif
                0526
bd2e80b12f Jean*0527 #ifdef NONLIN_FRSURF
                0528 C-- account for 3.D divergence of the flow in rStar coordinate:
cdc9f269ae Patr*0529 # ifndef DISABLE_RSTAR_CODE
722a76e285 Jean*0530         IF ( select_rStar.GT.0 ) THEN
                0531          DO j=jMin,jMax
                0532           DO i=iMin,iMax
                0533            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
88b07ffa37 Jean*0534      &     - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTFreeSurf
bd2e80b12f Jean*0535      &       *uVel(i,j,k,bi,bj)
722a76e285 Jean*0536           ENDDO
                0537          ENDDO
                0538         ENDIF
                0539         IF ( select_rStar.LT.0 ) THEN
                0540          DO j=jMin,jMax
                0541           DO i=iMin,iMax
                0542            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
bd2e80b12f Jean*0543      &     - rStarDhWDt(i,j,bi,bj)*uVel(i,j,k,bi,bj)
722a76e285 Jean*0544           ENDDO
                0545          ENDDO
                0546         ENDIF
cdc9f269ae Patr*0547 # endif /* DISABLE_RSTAR_CODE */
722a76e285 Jean*0548 #endif /* NONLIN_FRSURF */
                0549
9e4f1da9cf Jean*0550 #ifdef ALLOW_ADDFLUID
                0551         IF ( selectAddFluid.GE.1 ) THEN
                0552          DO j=jMin,jMax
                0553           DO i=iMin,iMax
                0554            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
                0555      &     + uVel(i,j,k,bi,bj)*mass2rUnit*0.5 _d 0
                0556      &      *( addMass(i-1,j,k,bi,bj) + addMass(i,j,k,bi,bj) )
                0557      &      *_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)*recip_rhoFacC(k)
                0558      &      * recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)
                0559           ENDDO
                0560          ENDDO
                0561         ENDIF
                0562 #endif /* ALLOW_ADDFLUID */
                0563
722a76e285 Jean*0564       ELSE
                0565 C-    if momAdvection / else
                0566         DO j=1-OLy,sNy+OLy
                0567          DO i=1-OLx,sNx+OLx
                0568            gU(i,j,k,bi,bj) = 0. _d 0
                0569          ENDDO
bd2e80b12f Jean*0570         ENDDO
722a76e285 Jean*0571
                0572 C-    endif momAdvection.
bd2e80b12f Jean*0573       ENDIF
722a76e285 Jean*0574
                0575       IF (momViscosity) THEN
                0576 C---  Calculate eddy fluxes (dissipation) between cells for zonal flow.
                0577
                0578 C     Bi-harmonic term del^2 U -> v4F
f0501c53d1 Jean*0579         IF ( useBiharmonicVisc )
fd362e9f7c Jean*0580      &  CALL MOM_U_DEL2U( bi, bj, k, uFld, hFacZ, h0FacZ,
                0581      O                    v4f, myThid )
7448700841 Mart*0582 #ifdef ALLOW_AUTODIFF_TAMC
                0583 CADJ STORE v4f(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0584 #endif
722a76e285 Jean*0585
                0586 C     Laplacian and bi-harmonic terms, Zonal  Fluxes -> fZon
fd362e9f7c Jean*0587         CALL MOM_U_XVISCFLUX( bi,bj,k,uFld,v4F,fZon,
                0588      I                        viscAh_D,viscA4_D,myThid )
722a76e285 Jean*0589
                0590 C     Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
fd362e9f7c Jean*0591         CALL MOM_U_YVISCFLUX( bi,bj,k,uFld,v4F,hFacZ,fMer,
                0592      I                        viscAh_Z,viscA4_Z,myThid )
722a76e285 Jean*0593
                0594 C     Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
                0595        IF (.NOT.implicitViscosity) THEN
1833b564cb Jean*0596         CALL MOM_U_RVISCFLUX( bi,bj, k, uVel,kappaRU,fVrUp,myThid )
                0597         CALL MOM_U_RVISCFLUX( bi,bj,k+1,uVel,kappaRU,fVrDw,myThid )
722a76e285 Jean*0598        ENDIF
                0599
                0600 C--   Tendency is minus divergence of the fluxes
eaba2fd266 Jean*0601 C     anelastic: hor.visc.fluxes are not scaled by rhoFac (by vert.visc.flx is)
722a76e285 Jean*0602         DO j=jMin,jMax
                0603          DO i=iMin,iMax
                0604           guDiss(i,j) =
                0605 #ifdef OLD_UV_GEOM
                0606      &     -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
                0607      &      ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
                0608 #else
                0609      &     -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
eaba2fd266 Jean*0610      &     *recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)
722a76e285 Jean*0611 #endif
eaba2fd266 Jean*0612      &     *( ( fZon(i,j  ) - fZon(i-1,j) )*AhDudxFac
                0613      &       +( fMer(i,j+1) - fMer(i,  j) )*AhDudyFac
                0614      &       +( fVrDw(i,j)  - fVrUp(i,j)  )*rkSign*ArDudrFac
                0615      &                                     *recip_rhoFacC(k)
722a76e285 Jean*0616      &     )
                0617          ENDDO
                0618         ENDDO
bd2e80b12f Jean*0619
711329b07b Jean*0620 #ifdef ALLOW_DIAGNOSTICS
                0621         IF ( useDiagnostics ) THEN
                0622           CALL DIAGNOSTICS_FILL(fZon, 'VISCx_Um',k,1,2,bi,bj,myThid)
                0623           CALL DIAGNOSTICS_FILL(fMer, 'VISCy_Um',k,1,2,bi,bj,myThid)
                0624           IF (.NOT.implicitViscosity)
                0625      &    CALL DIAGNOSTICS_FILL(fVrUp,'VISrE_Um',k,1,2,bi,bj,myThid)
                0626         ENDIF
                0627 #endif
                0628
9496c6c9ef Jean*0629 C-- No-slip and drag BCs appear as body forces in cell abutting topography
722a76e285 Jean*0630         IF (no_slip_sides) THEN
aea29c8517 Alis*0631 C-     No-slip BCs impose a drag at walls...
f0501c53d1 Jean*0632          CALL MOM_U_SIDEDRAG( bi, bj, k,
fd362e9f7c Jean*0633      I        uFld, v4f, h0FacZ,
f0501c53d1 Jean*0634      I        viscAh_Z, viscA4_Z,
                0635      I        useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
998681995e Bayl*0636      O        vF,
f0501c53d1 Jean*0637      I        myThid )
722a76e285 Jean*0638          DO j=jMin,jMax
                0639           DO i=iMin,iMax
79074ef66b Jean*0640            guDiss(i,j) = guDiss(i,j) + vF(i,j)
722a76e285 Jean*0641           ENDDO
                0642          ENDDO
                0643         ENDIF
aea29c8517 Alis*0644 C-    No-slip BCs impose a drag at bottom
e2d988bd46 Jean*0645         IF ( bottomDragTerms ) THEN
79074ef66b Jean*0646 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0647 CADJ STORE KE(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
79074ef66b Jean*0648 #endif
a125ab7be7 Jean*0649          CALL MOM_U_BOTDRAG_COEFF( bi, bj, k, .TRUE.,
79074ef66b Jean*0650      I                   uFld, vFld, kappaRU, KE,
                0651      O                   cDrag,
                0652      I                   myIter, myThid )
7448700841 Mart*0653 #ifdef ALLOW_AUTODIFF_TAMC
                0654 CADJ STORE cDrag(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0655 #endif
722a76e285 Jean*0656          DO j=jMin,jMax
                0657           DO i=iMin,iMax
79074ef66b Jean*0658             guDiss(i,j) = guDiss(i,j)
                0659      &                  - cDrag(i,j)*uFld(i,j)
                0660      &                  *_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
722a76e285 Jean*0661           ENDDO
                0662          ENDDO
79074ef66b Jean*0663          IF ( useDiagnostics ) THEN
                0664           DO j=jMin,jMax
                0665            DO i=iMin,iMax
                0666             botDragU(i,j,bi,bj) = botDragU(i,j,bi,bj)
                0667      &                          - cDrag(i,j)*uFld(i,j)*rUnit2mass
                0668            ENDDO
                0669           ENDDO
                0670          ENDIF
722a76e285 Jean*0671         ENDIF
                0672
dd49642a1d Mart*0673 #ifdef ALLOW_SHELFICE
e2d988bd46 Jean*0674         IF ( useShelfIce .AND. selectImplicitDrag.EQ.0 ) THEN
                0675 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0676 CADJ STORE KE(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
e2d988bd46 Jean*0677 #endif
                0678          CALL SHELFICE_U_DRAG_COEFF( bi, bj, k, .TRUE.,
                0679      I                   uFld, vFld, kappaRU, KE,
                0680      O                   cDrag,
                0681      I                   myIter, myThid )
7448700841 Mart*0682 #ifdef ALLOW_AUTODIFF_TAMC
                0683 CADJ STORE cDrag(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0684 #endif
dd49642a1d Mart*0685          DO j=jMin,jMax
                0686           DO i=iMin,iMax
e2d988bd46 Jean*0687             guDiss(i,j) = guDiss(i,j)
                0688      &                  - cDrag(i,j)*uFld(i,j)
                0689      &                  *_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
dd49642a1d Mart*0690           ENDDO
                0691          ENDDO
                0692         ENDIF
                0693 #endif /* ALLOW_SHELFICE */
                0694
722a76e285 Jean*0695 C-    endif momViscosity
aea29c8517 Alis*0696       ENDIF
                0697
7fc445a0ef Jean*0698 C--   Forcing term (moved to timestep.F)
                0699 c     IF (momForcing)
                0700 c    &  CALL EXTERNAL_FORCING_U(
                0701 c    I     iMin,iMax,jMin,jMax,bi,bj,k,
                0702 c    I     myTime,myThid)
aea29c8517 Alis*0703
                0704 C--   Metric terms for curvilinear grid systems
3121bb922d Alis*0705       IF (useNHMTerms) THEN
4d2713b160 Jean*0706 C      o Non-Hydrostatic (spherical) metric terms
fd362e9f7c Jean*0707        CALL MOM_U_METRIC_NH( bi,bj,k,uFld,wVel,mT,myThid )
aea29c8517 Alis*0708        DO j=jMin,jMax
                0709         DO i=iMin,iMax
4d2713b160 Jean*0710          gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtNHFacU*mT(i,j)
aea29c8517 Alis*0711         ENDDO
                0712        ENDDO
3121bb922d Alis*0713       ENDIF
4d2713b160 Jean*0714       IF ( usingSphericalPolarGrid .AND. metricTerms ) THEN
                0715 C      o Spherical polar grid metric terms
fd362e9f7c Jean*0716        CALL MOM_U_METRIC_SPHERE( bi,bj,k,uFld,vFld,mT,myThid )
aea29c8517 Alis*0717        DO j=jMin,jMax
                0718         DO i=iMin,iMax
4d2713b160 Jean*0719          gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*mT(i,j)
aea29c8517 Alis*0720         ENDDO
                0721        ENDDO
0de9c264f0 Andr*0722       ENDIF
4d2713b160 Jean*0723       IF ( usingCylindricalGrid .AND. metricTerms ) THEN
                0724 C      o Cylindrical grid metric terms
fd362e9f7c Jean*0725        CALL MOM_U_METRIC_CYLINDER( bi,bj,k,uFld,vFld,mT,myThid )
4d2713b160 Jean*0726        DO j=jMin,jMax
                0727         DO i=iMin,iMax
                0728          gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*mT(i,j)
                0729         ENDDO
0ac260a803 Andr*0730        ENDDO
aea29c8517 Alis*0731       ENDIF
                0732
722a76e285 Jean*0733 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
aea29c8517 Alis*0734
                0735 C---- Meridional momentum equation starts here
                0736
722a76e285 Jean*0737       IF (momAdvection) THEN
229b6d70e7 Jean*0738
                0739 #ifdef MOM_BOUNDARY_CONSERVE
                0740         CALL MOM_V_ADV_UV( bi,bj,k,uTrans,vBnd(1-OLx,1-OLy,k,bi,bj),
                0741      O                     fZon,myThid )
                0742         CALL MOM_V_ADV_VV( bi,bj,k,vTrans,vBnd(1-OLx,1-OLy,k,bi,bj),
                0743      O                     fMer,myThid )
07e17fa184 Jean*0744         CALL MOM_V_ADV_WV( bi,bj,k+1,vBnd,wVel,rTransV,
                0745      O                     fVerVkp, myThid )
229b6d70e7 Jean*0746 #else /* MOM_BOUNDARY_CONSERVE */
722a76e285 Jean*0747 C---  Calculate mean fluxes (advection)   between cells for meridional flow.
                0748 C     Mean flow component of zonal flux -> fZon
07e17fa184 Jean*0749         CALL MOM_V_ADV_UV( bi,bj,k,uTrans,vFld,fZon,myThid )
aea29c8517 Alis*0750
                0751 C--   Meridional flux (fMer is at north face of "v" cell)
722a76e285 Jean*0752 C     Mean flow component of meridional flux -> fMer
07e17fa184 Jean*0753         CALL MOM_V_ADV_VV( bi,bj,k,vTrans,vFld,fMer,myThid )
aea29c8517 Alis*0754
                0755 C--   Vertical flux (fVer is at upper face of "v" cell)
722a76e285 Jean*0756 C     Mean flow component of vertical flux (at k+1) -> fVerV
07e17fa184 Jean*0757         CALL MOM_V_ADV_WV( bi,bj,k+1,vVel,wVel,rTransV,
                0758      O                     fVerVkp, myThid )
229b6d70e7 Jean*0759 #endif /* MOM_BOUNDARY_CONSERVE */
aea29c8517 Alis*0760
                0761 C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
722a76e285 Jean*0762         DO j=jMin,jMax
                0763          DO i=iMin,iMax
                0764           gV(i,j,k,bi,bj) =
aea29c8517 Alis*0765 #ifdef OLD_UV_GEOM
722a76e285 Jean*0766      &     -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
                0767      &      ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
aea29c8517 Alis*0768 #else
722a76e285 Jean*0769      &     -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
eaba2fd266 Jean*0770      &     *recip_rAs(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
aea29c8517 Alis*0771 #endif
07e17fa184 Jean*0772      &     *( ( fZon(i+1,j)  - fZon(i,j  )  )*uDvdxFac
                0773      &       +( fMer(i,  j)  - fMer(i,j-1)  )*vDvdyFac
                0774      &       +( fVerVkp(i,j) - fVerVkm(i,j) )*rkSign*rVelDvdrFac
722a76e285 Jean*0775      &     )
711329b07b Jean*0776          ENDDO
                0777         ENDDO
                0778
                0779 #ifdef ALLOW_DIAGNOSTICS
                0780         IF ( useDiagnostics ) THEN
07e17fa184 Jean*0781           CALL DIAGNOSTICS_FILL( fZon,  'ADVx_Vm ',k,1,2,bi,bj,myThid)
                0782           CALL DIAGNOSTICS_FILL( fMer,  'ADVy_Vm ',k,1,2,bi,bj,myThid)
                0783           CALL DIAGNOSTICS_FILL(fVerVkm,'ADVrE_Vm',k,1,2,bi,bj,myThid)
711329b07b Jean*0784         ENDIF
                0785 #endif
aea29c8517 Alis*0786
bd2e80b12f Jean*0787 #ifdef NONLIN_FRSURF
                0788 C-- account for 3.D divergence of the flow in rStar coordinate:
cdc9f269ae Patr*0789 # ifndef DISABLE_RSTAR_CODE
722a76e285 Jean*0790         IF ( select_rStar.GT.0 ) THEN
                0791          DO j=jMin,jMax
                0792           DO i=iMin,iMax
                0793            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
88b07ffa37 Jean*0794      &     - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTFreeSurf
bd2e80b12f Jean*0795      &       *vVel(i,j,k,bi,bj)
722a76e285 Jean*0796           ENDDO
                0797          ENDDO
                0798         ENDIF
                0799         IF ( select_rStar.LT.0 ) THEN
                0800          DO j=jMin,jMax
                0801           DO i=iMin,iMax
                0802            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
bd2e80b12f Jean*0803      &     - rStarDhSDt(i,j,bi,bj)*vVel(i,j,k,bi,bj)
722a76e285 Jean*0804           ENDDO
                0805          ENDDO
                0806         ENDIF
cdc9f269ae Patr*0807 # endif /* DISABLE_RSTAR_CODE */
722a76e285 Jean*0808 #endif /* NONLIN_FRSURF */
                0809
9e4f1da9cf Jean*0810 #ifdef ALLOW_ADDFLUID
                0811         IF ( selectAddFluid.GE.1 ) THEN
                0812          DO j=jMin,jMax
                0813           DO i=iMin,iMax
                0814            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
                0815      &     + vVel(i,j,k,bi,bj)*mass2rUnit*0.5 _d 0
                0816      &      *( addMass(i,j-1,k,bi,bj) + addMass(i,j,k,bi,bj) )
                0817      &      *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)*recip_rhoFacC(k)
                0818      &      * recip_rAs(i,j,bi,bj)*recip_deepFac2C(k)
                0819           ENDDO
                0820          ENDDO
                0821         ENDIF
                0822 #endif /* ALLOW_ADDFLUID */
                0823
722a76e285 Jean*0824       ELSE
                0825 C-    if momAdvection / else
                0826         DO j=1-OLy,sNy+OLy
                0827          DO i=1-OLx,sNx+OLx
                0828            gV(i,j,k,bi,bj) = 0. _d 0
                0829          ENDDO
bd2e80b12f Jean*0830         ENDDO
722a76e285 Jean*0831
                0832 C-    endif momAdvection.
bd2e80b12f Jean*0833       ENDIF
722a76e285 Jean*0834
                0835       IF (momViscosity) THEN
                0836 C---  Calculate eddy fluxes (dissipation) between cells for meridional flow.
                0837 C     Bi-harmonic term del^2 V -> v4F
f0501c53d1 Jean*0838         IF ( useBiharmonicVisc )
fd362e9f7c Jean*0839      &  CALL MOM_V_DEL2V( bi, bj, k, vFld, hFacZ, h0FacZ,
                0840      O                    v4f, myThid )
7448700841 Mart*0841 #ifdef ALLOW_AUTODIFF_TAMC
                0842 CADJ STORE v4f(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0843 #endif
722a76e285 Jean*0844
                0845 C     Laplacian and bi-harmonic terms, Zonal  Fluxes -> fZon
fd362e9f7c Jean*0846         CALL MOM_V_XVISCFLUX( bi,bj,k,vFld,v4f,hFacZ,fZon,
                0847      I                        viscAh_Z,viscA4_Z,myThid )
722a76e285 Jean*0848
                0849 C     Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
fd362e9f7c Jean*0850         CALL MOM_V_YVISCFLUX( bi,bj,k,vFld,v4f,fMer,
                0851      I                        viscAh_D,viscA4_D,myThid )
722a76e285 Jean*0852
                0853 C     Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
                0854        IF (.NOT.implicitViscosity) THEN
fd362e9f7c Jean*0855         CALL MOM_V_RVISCFLUX( bi,bj, k, vVel,KappaRV,fVrUp,myThid )
                0856         CALL MOM_V_RVISCFLUX( bi,bj,k+1,vVel,KappaRV,fVrDw,myThid )
722a76e285 Jean*0857        ENDIF
                0858
                0859 C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
eaba2fd266 Jean*0860 C     anelastic: hor.visc.fluxes are not scaled by rhoFac (by vert.visc.flx is)
722a76e285 Jean*0861         DO j=jMin,jMax
                0862          DO i=iMin,iMax
                0863           gvDiss(i,j) =
                0864 #ifdef OLD_UV_GEOM
                0865      &     -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
                0866      &      ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
                0867 #else
                0868      &     -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
eaba2fd266 Jean*0869      &      *recip_rAs(i,j,bi,bj)*recip_deepFac2C(k)
722a76e285 Jean*0870 #endif
eaba2fd266 Jean*0871      &     *( ( fZon(i+1,j)  - fZon(i,j  ) )*AhDvdxFac
                0872      &       +( fMer(i,  j)  - fMer(i,j-1) )*AhDvdyFac
                0873      &       +( fVrDw(i,j)   - fVrUp(i,j) )*rkSign*ArDvdrFac
                0874      &                                     *recip_rhoFacC(k)
722a76e285 Jean*0875      &     )
                0876          ENDDO
                0877         ENDDO
bd2e80b12f Jean*0878
711329b07b Jean*0879 #ifdef ALLOW_DIAGNOSTICS
                0880         IF ( useDiagnostics ) THEN
                0881           CALL DIAGNOSTICS_FILL(fZon, 'VISCx_Vm',k,1,2,bi,bj,myThid)
                0882           CALL DIAGNOSTICS_FILL(fMer, 'VISCy_Vm',k,1,2,bi,bj,myThid)
                0883           IF (.NOT.implicitViscosity)
                0884      &    CALL DIAGNOSTICS_FILL(fVrUp,'VISrE_Vm',k,1,2,bi,bj,myThid)
                0885         ENDIF
                0886 #endif
                0887
9496c6c9ef Jean*0888 C-- No-slip and drag BCs appear as body forces in cell abutting topography
dd49642a1d Mart*0889         IF (no_slip_sides) THEN
aea29c8517 Alis*0890 C-     No-slip BCs impose a drag at walls...
f0501c53d1 Jean*0891          CALL MOM_V_SIDEDRAG( bi, bj, k,
fd362e9f7c Jean*0892      I        vFld, v4f, h0FacZ,
                0893      I        viscAh_Z, viscA4_Z,
f0501c53d1 Jean*0894      I        useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
998681995e Bayl*0895      O        vF,
f0501c53d1 Jean*0896      I        myThid )
722a76e285 Jean*0897          DO j=jMin,jMax
                0898           DO i=iMin,iMax
                0899            gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
                0900           ENDDO
                0901          ENDDO
                0902         ENDIF
aea29c8517 Alis*0903 C-    No-slip BCs impose a drag at bottom
e2d988bd46 Jean*0904         IF ( bottomDragTerms ) THEN
79074ef66b Jean*0905 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0906 CADJ STORE KE(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
79074ef66b Jean*0907 #endif
a125ab7be7 Jean*0908          CALL MOM_V_BOTDRAG_COEFF( bi, bj, k, .TRUE.,
79074ef66b Jean*0909      I                   uFld, vFld, kappaRV, KE,
                0910      O                   cDrag,
                0911      I                   myIter, myThid )
7448700841 Mart*0912 #ifdef ALLOW_AUTODIFF_TAMC
                0913 CADJ STORE cDrag(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0914 #endif
722a76e285 Jean*0915          DO j=jMin,jMax
                0916           DO i=iMin,iMax
79074ef66b Jean*0917             gvDiss(i,j) = gvDiss(i,j)
                0918      &                  - cDrag(i,j)*vFld(i,j)
                0919      &                  *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
722a76e285 Jean*0920           ENDDO
                0921          ENDDO
79074ef66b Jean*0922          IF ( useDiagnostics ) THEN
                0923           DO j=jMin,jMax
                0924            DO i=iMin,iMax
                0925             botDragV(i,j,bi,bj) = botDragV(i,j,bi,bj)
                0926      &                          - cDrag(i,j)*vFld(i,j)*rUnit2mass
                0927            ENDDO
                0928           ENDDO
                0929          ENDIF
722a76e285 Jean*0930         ENDIF
                0931
dd49642a1d Mart*0932 #ifdef ALLOW_SHELFICE
e2d988bd46 Jean*0933         IF ( useShelfIce .AND. selectImplicitDrag.EQ.0 ) THEN
                0934 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0935 CADJ STORE KE(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
e2d988bd46 Jean*0936 #endif
                0937          CALL SHELFICE_V_DRAG_COEFF( bi, bj, k, .TRUE.,
                0938      I                   uFld, vFld, kappaRV, KE,
                0939      O                   cDrag,
                0940      I                   myIter, myThid )
7448700841 Mart*0941 #ifdef ALLOW_AUTODIFF_TAMC
                0942 CADJ STORE cDrag(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0943 #endif
dd49642a1d Mart*0944          DO j=jMin,jMax
                0945           DO i=iMin,iMax
e2d988bd46 Jean*0946             gvDiss(i,j) = gvDiss(i,j)
                0947      &                  - cDrag(i,j)*vFld(i,j)
                0948      &                  *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
dd49642a1d Mart*0949           ENDDO
                0950          ENDDO
                0951         ENDIF
                0952 #endif /* ALLOW_SHELFICE */
                0953
722a76e285 Jean*0954 C-    endif momViscosity
aea29c8517 Alis*0955       ENDIF
                0956
7fc445a0ef Jean*0957 C--   Forcing term (moved to timestep.F)
                0958 c     IF (momForcing)
                0959 c    & CALL EXTERNAL_FORCING_V(
                0960 c    I     iMin,iMax,jMin,jMax,bi,bj,k,
                0961 c    I     myTime,myThid)
aea29c8517 Alis*0962
                0963 C--   Metric terms for curvilinear grid systems
3121bb922d Alis*0964       IF (useNHMTerms) THEN
4d2713b160 Jean*0965 C      o Non-Hydrostatic (spherical) metric terms
fd362e9f7c Jean*0966        CALL MOM_V_METRIC_NH( bi,bj,k,vFld,wVel,mT,myThid )
aea29c8517 Alis*0967        DO j=jMin,jMax
                0968         DO i=iMin,iMax
4d2713b160 Jean*0969          gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtNHFacV*mT(i,j)
aea29c8517 Alis*0970         ENDDO
                0971        ENDDO
3121bb922d Alis*0972       ENDIF
4d2713b160 Jean*0973       IF ( usingSphericalPolarGrid .AND. metricTerms ) THEN
                0974 C      o Spherical polar grid metric terms
fd362e9f7c Jean*0975        CALL MOM_V_METRIC_SPHERE( bi,bj,k,uFld,mT,myThid )
aea29c8517 Alis*0976        DO j=jMin,jMax
                0977         DO i=iMin,iMax
4d2713b160 Jean*0978          gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*mT(i,j)
aea29c8517 Alis*0979         ENDDO
                0980        ENDDO
                0981       ENDIF
4d2713b160 Jean*0982       IF ( usingCylindricalGrid .AND. metricTerms ) THEN
                0983 C      o Cylindrical grid metric terms
fd362e9f7c Jean*0984        CALL MOM_V_METRIC_CYLINDER( bi,bj,k,uFld,vFld,mT,myThid )
4d2713b160 Jean*0985        DO j=jMin,jMax
                0986         DO i=iMin,iMax
                0987          gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*mT(i,j)
                0988         ENDDO
                0989        ENDDO
0ac260a803 Andr*0990       ENDIF
aea29c8517 Alis*0991
722a76e285 Jean*0992 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
aea29c8517 Alis*0993
c831f9444b Jean*0994 C--   Coriolis term (call to CD_CODE_SCHEME has been moved to timestep.F)
9293d3c672 Hajo*0995       IF ( .NOT.useCDscheme ) THEN
                0996 #if ( defined ALLOW_GGL90 && defined ALLOW_GGL90_LANGMUIR )
                0997        IF ( useLANGMUIR ) THEN
                0998         CALL GGL90_ADD_STOKESDRIFT(
                0999      O                 uRes, vRes,
                1000      I                 uFld, vFld, k, bi, bj, myThid )
                1001         CALL MOM_U_CORIOLIS( bi,bj,k,vRes,uCf,myThid )
                1002         CALL MOM_V_CORIOLIS( bi,bj,k,uRes,vCf,myThid )
                1003        ELSE
                1004 #endif
                1005         CALL MOM_U_CORIOLIS( bi,bj,k,vFld,uCf,myThid )
                1006         CALL MOM_V_CORIOLIS( bi,bj,k,uFld,vCf,myThid )
                1007 #if ( defined ALLOW_GGL90 && defined ALLOW_GGL90_LANGMUIR )
                1008        ENDIF
711329b07b Jean*1009 #endif
7fc445a0ef Jean*1010         DO j=jMin,jMax
                1011          DO i=iMin,iMax
9293d3c672 Hajo*1012           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + fuFac*uCf(i,j)
                1013           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + fvFac*vCf(i,j)
7fc445a0ef Jean*1014          ENDDO
                1015         ENDDO
711329b07b Jean*1016 #ifdef ALLOW_DIAGNOSTICS
9293d3c672 Hajo*1017         IF ( useDiagnostics ) THEN
                1018           CALL DIAGNOSTICS_FILL( uCf,'Um_Cori ',k,1,2,bi,bj,myThid )
                1019           CALL DIAGNOSTICS_FILL( vCf,'Vm_Cori ',k,1,2,bi,bj,myThid )
                1020         ENDIF
711329b07b Jean*1021 #endif
7fc445a0ef Jean*1022       ENDIF
                1023
3daafce20b Jean*1024 C--   3.D Coriolis term (horizontal momentum, Eastward component: -fprime*w)
427e24e121 Jean*1025       IF ( select3dCoriScheme.GE.1 ) THEN
9293d3c672 Hajo*1026         CALL MOM_U_CORIOLIS_NH( bi,bj,k,wVel,uCf,myThid )
ac8c612649 Jean*1027         DO j=jMin,jMax
                1028          DO i=iMin,iMax
9293d3c672 Hajo*1029           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + fuFac*uCf(i,j)
ac8c612649 Jean*1030          ENDDO
c5990018f4 Alis*1031         ENDDO
427e24e121 Jean*1032        IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
ac8c612649 Jean*1033 C-     presently, non zero angleSinC array only supported with Curvilinear-Grid
9293d3c672 Hajo*1034         CALL MOM_V_CORIOLIS_NH( bi,bj,k,wVel,vCf,myThid )
ac8c612649 Jean*1035         DO j=jMin,jMax
                1036          DO i=iMin,iMax
9293d3c672 Hajo*1037           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + fvFac*vCf(i,j)
ac8c612649 Jean*1038          ENDDO
                1039         ENDDO
                1040        ENDIF
c5990018f4 Alis*1041       ENDIF
aea29c8517 Alis*1042
722a76e285 Jean*1043 C--   Set du/dt & dv/dt on boundaries to zero
                1044       DO j=jMin,jMax
                1045        DO i=iMin,iMax
                1046         gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
                1047         guDiss(i,j)     = guDiss(i,j)    *_maskW(i,j,k,bi,bj)
                1048         gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
                1049         gvDiss(i,j)     = gvDiss(i,j)    *_maskS(i,j,k,bi,bj)
                1050        ENDDO
                1051       ENDDO
                1052
711329b07b Jean*1053 #ifdef ALLOW_DIAGNOSTICS
                1054       IF ( useDiagnostics ) THEN
b114b87ee5 Bayl*1055         CALL DIAGNOSTICS_FILL(KE,    'momKE   ',k,1,2,bi,bj,myThid)
9e4f1da9cf Jean*1056         CALL DIAGNOSTICS_FILL(gU(1-OLx,1-OLy,k,bi,bj),
711329b07b Jean*1057      &                               'Um_Advec',k,1,2,bi,bj,myThid)
9e4f1da9cf Jean*1058         CALL DIAGNOSTICS_FILL(gV(1-OLx,1-OLy,k,bi,bj),
711329b07b Jean*1059      &                               'Vm_Advec',k,1,2,bi,bj,myThid)
                1060       ENDIF
                1061 #endif /* ALLOW_DIAGNOSTICS */
                1062
aea29c8517 Alis*1063       RETURN
                1064       END