pkg/mom_vecinv/mom_vecinv.F

cec2469d72 Alis*0001 #include "MOM_VECINV_OPTIONS.h"
88b07ffa37 Jean*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
f0501c53d1 Jean*0005 #ifdef ALLOW_MOM_COMMON
                0006 # include "MOM_COMMON_OPTIONS.h"
                0007 #endif
9293d3c672 Hajo*0008 #ifdef ALLOW_GGL90
                0009 # include "GGL90_OPTIONS.h"
                0010 #endif
aea29c8517 Alis*0011
039fe61d35 Jean*0012       SUBROUTINE MOM_VECINV(
07e17fa184 Jean*0013      I        bi,bj,k,iMin,iMax,jMin,jMax,
1833b564cb Jean*0014      I        kappaRU, kappaRV,
07e17fa184 Jean*0015      I        fVerUkm, fVerVkm,
                0016      O        fVerUkp, fVerVkp,
4667e4066d Jean*0017      O        guDiss, gvDiss,
07e17fa184 Jean*0018      I        myTime, myIter, myThid )
                0019 C     *==========================================================*
66d5e1666c Alis*0020 C     | S/R MOM_VECINV                                           |
aea29c8517 Alis*0021 C     | o Form the right hand-side of the momentum equation.     |
07e17fa184 Jean*0022 C     *==========================================================*
aea29c8517 Alis*0023 C     | Terms are evaluated one layer at a time working from     |
                0024 C     | the bottom to the top. The vertically integrated         |
                0025 C     | barotropic flow tendency term is evluated by summing the |
                0026 C     | tendencies.                                              |
                0027 C     | Notes:                                                   |
                0028 C     | We have not sorted out an entirely satisfactory formula  |
                0029 C     | for the diffusion equation bc with lopping. The present  |
                0030 C     | form produces a diffusive flux that does not scale with  |
                0031 C     | open-area. Need to do something to solidfy this and to   |
                0032 C     | deal "properly" with thin walls.                         |
07e17fa184 Jean*0033 C     *==========================================================*
aea29c8517 Alis*0034       IMPLICIT NONE
                0035
                0036 C     == Global variables ==
                0037 #include "SIZE.h"
                0038 #include "EEPARAMS.h"
                0039 #include "PARAMS.h"
                0040 #include "GRID.h"
01f860d49e Jean*0041 #include "SURFACE.h"
f0501c53d1 Jean*0042 #include "DYNVARS.h"
79074ef66b Jean*0043 #include "FFIELDS.h"
f0501c53d1 Jean*0044 #ifdef ALLOW_MOM_COMMON
                0045 # include "MOM_VISC.h"
                0046 #endif
9293d3c672 Hajo*0047 #if ( defined ALLOW_GGL90 && defined ALLOW_GGL90_LANGMUIR )
                0048 # include "GGL90.h"
                0049 #endif
cab1a69b8a Jean*0050 #ifdef ALLOW_TIMEAVE
f0501c53d1 Jean*0051 # include "TIMEAVE_STATV.h"
                0052 #endif
                0053 #ifdef ALLOW_MNC
                0054 # include "MNC_PARAMS.h"
cab1a69b8a Jean*0055 #endif
cd3c16aeda Patr*0056 #ifdef ALLOW_AUTODIFF_TAMC
                0057 # include "tamc.h"
                0058 #endif
aea29c8517 Alis*0059
                0060 C     == Routine arguments ==
07e17fa184 Jean*0061 C     bi,bj   :: current tile indices
                0062 C     k       :: current vertical level
                0063 C     iMin,iMax,jMin,jMax :: loop ranges
                0064 C     fVerU   :: Flux of momentum in the vertical direction, out of the upper
f1a4cec01a Jean*0065 C     fVerV   :: face of a cell k ( flux into the cell above ).
07e17fa184 Jean*0066 C     fVerUkm :: vertical viscous flux of U, interface above (k-1/2)
                0067 C     fVerVkm :: vertical viscous flux of V, interface above (k-1/2)
                0068 C     fVerUkp :: vertical viscous flux of U, interface below (k+1/2)
                0069 C     fVerVkp :: vertical viscous flux of V, interface below (k+1/2)
                0070
                0071 C     guDiss  :: dissipation tendency (all explicit terms), u component
                0072 C     gvDiss  :: dissipation tendency (all explicit terms), v component
                0073 C     myTime  :: current time
                0074 C     myIter  :: current time-step number
                0075 C     myThid  :: my Thread Id number
                0076       INTEGER bi,bj,k
                0077       INTEGER iMin,iMax,jMin,jMax
1833b564cb Jean*0078       _RL kappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
                0079       _RL kappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
07e17fa184 Jean*0080       _RL fVerUkm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0081       _RL fVerVkm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0082       _RL fVerUkp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0083       _RL fVerVkp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
4667e4066d Jean*0084       _RL guDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0085       _RL gvDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
f7d48db11c Jean*0086       _RL     myTime
3a279374db Alis*0087       INTEGER myIter
                0088       INTEGER myThid
aea29c8517 Alis*0089
99bc607d7a Ed H*0090 #ifdef ALLOW_MOM_VECINV
cab1a69b8a Jean*0091
3a279374db Alis*0092 C     == Functions ==
94a46dfe0d Jean*0093       LOGICAL  DIFFERENT_MULTIPLE
                0094       EXTERNAL DIFFERENT_MULTIPLE
3a279374db Alis*0095
aea29c8517 Alis*0096 C     == Local variables ==
ed2dbbe83d Jean*0097 C     strainBC :: same as strain but account for no-slip BC
                0098 C     vort3BC  :: same as vort3  but account for no-slip BC
aea29c8517 Alis*0099       _RL      vF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
df1ec59c8a Jean*0100       _RL      vrF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0101       _RL      uCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0102       _RL      vCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0103       _RS hFacZ   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
01f860d49e Jean*0104       _RS h0FacZ  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
df1ec59c8a Jean*0105       _RS r_hFacZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0106       _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0107       _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
9293d3c672 Hajo*0108 #if ( defined ALLOW_GGL90 && defined ALLOW_GGL90_LANGMUIR )
                0109       _RL uRes    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0110       _RL vRes    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0111 #endif
df1ec59c8a Jean*0112       _RL del2u   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0113       _RL del2v   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0114       _RL dStar   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0115       _RL zStar   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0116       _RL tension (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0117       _RL strain  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
ed2dbbe83d Jean*0118       _RL strainBC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
f59d76b0dd Ed D*0119       _RL stretching(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
06244a5e4f Jean*0120 #ifdef ALLOW_LEITH_QG
f59d76b0dd Ed D*0121       _RL Nsquare (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
06244a5e4f Jean*0122 #endif
79074ef66b Jean*0123       _RL cDrag   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
df1ec59c8a Jean*0124       _RL KE      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0125       _RL omega3  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0126       _RL vort3   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
ed2dbbe83d Jean*0127       _RL vort3BC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
df1ec59c8a Jean*0128       _RL hDiv    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0129       _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0130       _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0131       _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0132       _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
07e17fa184 Jean*0133 C     i,j    :: Loop counters
                0134       INTEGER i,j
                0135 C     xxxFac :: On-off tracer parameters used for switching terms off.
aea29c8517 Alis*0136       _RL  ArDudrFac
                0137       _RL  ArDvdrFac
df1ec59c8a Jean*0138       _RL  sideMaskFac
aea29c8517 Alis*0139       LOGICAL bottomDragTerms
f7d48db11c Jean*0140       LOGICAL writeDiag
cd3c16aeda Patr*0141 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0142 C     kkey :: tape key (depends on level and tile indices)
                0143       INTEGER kkey
cd3c16aeda Patr*0144 #endif
aea29c8517 Alis*0145
cb356b4c5f Ed H*0146 #ifdef ALLOW_MNC
                0147       INTEGER offsets(9)
b22b541fe9 Ed H*0148       CHARACTER*(1) pf
cb356b4c5f Ed H*0149 #endif
                0150
88b07ffa37 Jean*0151 #ifdef ALLOW_AUTODIFF
7d3b758ca2 Patr*0152 C--   only the kDown part of fverU/V is set in this subroutine
                0153 C--   the kUp is still required
f1a4cec01a Jean*0154 C--   In the case of mom_fluxform kUp is set as well
7d3b758ca2 Patr*0155 C--   (at least in part)
07e17fa184 Jean*0156       fVerUkm(1,1) = fVerUkm(1,1)
                0157       fVerVkm(1,1) = fVerVkm(1,1)
7d3b758ca2 Patr*0158 #endif
                0159
cd3c16aeda Patr*0160 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0161       kkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0162       kkey = k  + (kkey-1)*Nr
cd3c16aeda Patr*0163 #endif /* ALLOW_AUTODIFF_TAMC */
                0164
94a46dfe0d Jean*0165       writeDiag = DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock)
aea29c8517 Alis*0166
ef1e764710 Ed H*0167 #ifdef ALLOW_MNC
                0168       IF (useMNC .AND. snapshot_mnc .AND. writeDiag) THEN
b22b541fe9 Ed H*0169         IF ( writeBinaryPrec .EQ. precFloat64 ) THEN
                0170           pf(1:1) = 'D'
                0171         ELSE
                0172           pf(1:1) = 'R'
                0173         ENDIF
cb356b4c5f Ed H*0174         IF ((bi .EQ. 1).AND.(bj .EQ. 1).AND.(k .EQ. 1)) THEN
                0175           CALL MNC_CW_SET_UDIM('mom_vi', -1, myThid)
987ff12cb6 Ed H*0176           CALL MNC_CW_RL_W_S('D','mom_vi',0,0,'T',myTime,myThid)
cb356b4c5f Ed H*0177           CALL MNC_CW_SET_UDIM('mom_vi', 0, myThid)
987ff12cb6 Ed H*0178           CALL MNC_CW_I_W_S('I','mom_vi',0,0,'iter',myIter,myThid)
cb356b4c5f Ed H*0179         ENDIF
                0180         DO i = 1,9
                0181           offsets(i) = 0
                0182         ENDDO
                0183         offsets(3) = k
9c98e82bbb Jean*0184 c       write(*,*) 'offsets = ',(offsets(i),i=1,9)
ef1e764710 Ed H*0185       ENDIF
                0186 #endif /*  ALLOW_MNC  */
                0187
9c98e82bbb Jean*0188 C--   Initialise intermediate terms
                0189       DO j=1-OLy,sNy+OLy
                0190        DO i=1-OLx,sNx+OLx
4667e4066d Jean*0191         vF(i,j)    = 0.
                0192         vrF(i,j)   = 0.
aea29c8517 Alis*0193         uCf(i,j)   = 0.
                0194         vCf(i,j)   = 0.
                0195         del2u(i,j) = 0.
                0196         del2v(i,j) = 0.
                0197         dStar(i,j) = 0.
                0198         zStar(i,j) = 0.
4667e4066d Jean*0199         guDiss(i,j)= 0.
                0200         gvDiss(i,j)= 0.
9293d3c672 Hajo*0201 #if ( defined ALLOW_GGL90 && defined ALLOW_GGL90_LANGMUIR )
                0202 c       uRes(i,j)  = 0.
                0203 c       vRes(i,j)  = 0.
                0204 #endif
aea29c8517 Alis*0205         vort3(i,j) = 0.
4667e4066d Jean*0206         omega3(i,j)= 0.
df1ec59c8a Jean*0207         KE(i,j)    = 0.
9c98e82bbb Jean*0208 C-    need to initialise hDiv for MOM_VI_DEL2UV(call FILL_CS_CORNER_TR_RL)
                0209         hDiv(i,j)  = 0.
34e137f064 Jean*0210 c       viscAh_Z(i,j) = 0.
                0211 c       viscAh_D(i,j) = 0.
                0212 c       viscA4_Z(i,j) = 0.
                0213 c       viscA4_D(i,j) = 0.
629d440dd4 Patr*0214         strain(i,j)  = 0. _d 0
ed2dbbe83d Jean*0215         strainBC(i,j)= 0. _d 0
f59d76b0dd Ed D*0216         stretching(i,j) = 0. _d 0
06244a5e4f Jean*0217 #ifdef ALLOW_LEITH_QG
f59d76b0dd Ed D*0218         Nsquare(i,j) = 0. _d 0
06244a5e4f Jean*0219 #endif
629d440dd4 Patr*0220         tension(i,j) = 0. _d 0
88b07ffa37 Jean*0221 #ifdef ALLOW_AUTODIFF
cdc9f269ae Patr*0222         hFacZ(i,j)   = 0. _d 0
629d440dd4 Patr*0223 #endif
aea29c8517 Alis*0224        ENDDO
                0225       ENDDO
                0226
                0227 C--   Term by term tracer parmeters
                0228 C     o U momentum equation
                0229       ArDudrFac    = vfFacMom*1.
                0230 C     o V momentum equation
                0231       ArDvdrFac    = vfFacMom*1.
                0232
df1ec59c8a Jean*0233 C note: using standard stencil (no mask) results in under-estimating
                0234 C       vorticity at a no-slip boundary by a factor of 2 = sideDragFactor
                0235       IF ( no_slip_sides ) THEN
                0236         sideMaskFac = sideDragFactor
                0237       ELSE
                0238         sideMaskFac = 0. _d 0
                0239       ENDIF
                0240
99731c717f Jean*0241       IF ( selectImplicitDrag.EQ.0 .AND.
                0242      &      (  no_slip_bottom
932b38363b Jean*0243      &    .OR. selectBotDragQuadr.GE.0
99731c717f Jean*0244      &    .OR. bottomDragLinear.NE.0. ) ) THEN
aea29c8517 Alis*0245        bottomDragTerms=.TRUE.
                0246       ELSE
                0247        bottomDragTerms=.FALSE.
                0248       ENDIF
                0249
                0250 C--   Calculate open water fraction at vorticity points
                0251       CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
                0252
                0253 C     Make local copies of horizontal flow field
                0254       DO j=1-OLy,sNy+OLy
                0255        DO i=1-OLx,sNx+OLx
                0256         uFld(i,j) = uVel(i,j,k,bi,bj)
                0257         vFld(i,j) = vVel(i,j,k,bi,bj)
                0258        ENDDO
                0259       ENDDO
                0260
d09cd10c60 Gael*0261 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0262 CADJ STORE hFacZ(:,:)   = comlev1_bibj_k, key = kkey, byte = isbyte
                0263 CADJ STORE r_hFacZ(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0264 CADJ STORE fverukm(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0265 CADJ STORE fvervkm(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
d09cd10c60 Gael*0266 #endif
                0267
cab1a69b8a Jean*0268 C note (jmc) : Dissipation and Vort3 advection do not necesary
                0269 C              use the same maskZ (and hFacZ)  => needs 2 call(s)
                0270 c     CALL MOM_VI_HFACZ_DISS(bi,bj,k,hFacZ,r_hFacZ,myThid)
                0271
b8082fc644 Jean*0272       CALL MOM_CALC_KE(bi,bj,k,selectKEscheme,uFld,vFld,KE,myThid)
aea29c8517 Alis*0273
7c7b0b4a46 Alis*0274       CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
aea29c8517 Alis*0275
ed2dbbe83d Jean*0276 C-    mask vort3 and account for no-slip / free-slip BC in vort3BC:
                0277       DO j=1-OLy,sNy+OLy
                0278        DO i=1-OLx,sNx+OLx
                0279          vort3BC(i,j) = vort3(i,j)
                0280          IF ( hFacZ(i,j).EQ.zeroRS ) THEN
                0281            vort3BC(i,j) = sideMaskFac*vort3BC(i,j)
                0282            vort3(i,j)   = 0.
                0283          ENDIF
                0284        ENDDO
                0285       ENDDO
                0286
d09cd10c60 Gael*0287 #ifdef ALLOW_AUTODIFF_TAMC
7448700841 Mart*0288 CADJ STORE KE(:,:)      = comlev1_bibj_k, key = kkey, byte = isbyte
edb6656069 Mart*0289 CADJ STORE vort3(:,:)   = comlev1_bibj_k, key = kkey, byte = isbyte
                0290 CADJ STORE vort3bc(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
d09cd10c60 Gael*0291 #endif
                0292
aea29c8517 Alis*0293       IF (momViscosity) THEN
039fe61d35 Jean*0294 C--    For viscous term, compute horizontal divergence, tension & strain
df1ec59c8a Jean*0295 C      and mask relative vorticity (free-slip case):
                0296
01f860d49e Jean*0297        DO j=1-OLy,sNy+OLy
                0298         DO i=1-OLx,sNx+OLx
                0299           h0FacZ(i,j) = hFacZ(i,j)
                0300         ENDDO
                0301        ENDDO
                0302 #ifdef NONLIN_FRSURF
                0303        IF ( no_slip_sides .AND. nonlinFreeSurf.GT.0 ) THEN
                0304         DO j=2-OLy,sNy+OLy
                0305          DO i=2-OLx,sNx+OLx
                0306           h0FacZ(i,j) = MIN(
                0307      &       MIN( h0FacW(i,j,k,bi,bj), h0FacW(i,j-1,k,bi,bj) ),
                0308      &       MIN( h0FacS(i,j,k,bi,bj), h0FacS(i-1,j,k,bi,bj) ) )
                0309          ENDDO
                0310         ENDDO
                0311        ENDIF
7448700841 Mart*0312 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0313 CADJ STORE h0FacZ(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
7448700841 Mart*0314 # endif
                0315 #endif /* NONLIN_FRSURF */
d09cd10c60 Gael*0316
df1ec59c8a Jean*0317        CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
                0318
ed2dbbe83d Jean*0319        IF ( useVariableVisc .OR. useStrainTensionVisc ) THEN
                0320         CALL MOM_CALC_TENSION( bi,bj,k,uFld,vFld,tension,myThid )
                0321         CALL MOM_CALC_STRAIN( bi,bj,k,uFld,vFld,hFacZ,strain,myThid )
                0322 C-    mask strain and account for no-slip / free-slip BC in strainBC:
                0323         DO j=1-OLy,sNy+OLy
                0324          DO i=1-OLx,sNx+OLx
                0325            strainBC(i,j) = strain(i,j)
                0326            IF ( hFacZ(i,j).EQ.zeroRS ) THEN
                0327              strainBC(i,j) = sideMaskFac*strainBC(i,j)
                0328              strain(i,j)   = 0.
                0329            ENDIF
                0330          ENDDO
df1ec59c8a Jean*0331         ENDDO
f59d76b0dd Ed D*0332 #ifdef ALLOW_LEITH_QG
a125ab7be7 Jean*0333         IF ( viscC2LeithQG.NE.zeroRL ) THEN
f59d76b0dd Ed D*0334           CALL MOM_VISC_QGL_STRETCH(bi,bj,k,
                0335      O                            stretching, Nsquare,
                0336      I                            myTime, myIter, myThid )
                0337           CALL MOM_VISC_QGL_LIMIT(bi,bj,k,
                0338      O                          stretching,
                0339      I                          Nsquare, uFld,vFld,vort3,
                0340      I                          myTime, myIter, myThid )
                0341         ENDIF
                0342 #endif /* ALLOW_LEITH_QG */
ed2dbbe83d Jean*0343        ENDIF
df1ec59c8a Jean*0344
d09cd10c60 Gael*0345 #ifdef ALLOW_AUTODIFF_TAMC
7448700841 Mart*0346 CADJ STORE hDiv(:,:)     = comlev1_bibj_k, key = kkey, byte = isbyte
edb6656069 Mart*0347 CADJ STORE tension(:,:)  = comlev1_bibj_k, key = kkey, byte = isbyte
                0348 CADJ STORE strain(:,:)   = comlev1_bibj_k, key = kkey, byte = isbyte
7448700841 Mart*0349 CADJ STORE strainBC(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
d09cd10c60 Gael*0350 #endif
                0351
34e137f064 Jean*0352 C--    Calculate Lateral Viscosities
                0353        DO j=1-OLy,sNy+OLy
                0354         DO i=1-OLx,sNx+OLx
                0355          viscAh_D(i,j) = viscAhD
                0356          viscAh_Z(i,j) = viscAhZ
                0357          viscA4_D(i,j) = viscA4D
                0358          viscA4_Z(i,j) = viscA4Z
                0359         ENDDO
                0360        ENDDO
                0361        IF ( useVariableVisc ) THEN
ed2dbbe83d Jean*0362 C-     uses vort3BC & strainBC which account for no-slip / free-slip BC
79074ef66b Jean*0363 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0364 CADJ STORE KE(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
79074ef66b Jean*0365 #endif
34e137f064 Jean*0366          CALL MOM_CALC_VISC( bi, bj, k,
                0367      O            viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
f59d76b0dd Ed D*0368      I            hDiv, vort3BC, tension, strainBC, stretching,
                0369      I            KE, hfacZ, myThid )
d09cd10c60 Gael*0370 #ifdef ALLOW_AUTODIFF_TAMC
4240547d2d Mart*0371 # ifndef AUTODIFF_ALLOW_VISCFACADJ
                0372 C     These store directive must not be here, if you want to recompute
                0373 C     these viscosity coefficients with a modified viscFacAdj = viscFacInAd
                0374 C     because the store directives intentionally prevent the recomputation.
edb6656069 Mart*0375 CADJ STORE viscAh_Z(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0376 CADJ STORE viscAh_D(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0377 CADJ STORE viscA4_Z(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0378 CADJ STORE viscA4_D(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
4240547d2d Mart*0379 # endif /* AUTODIFF_ALLOW_VISCFACADJ */
                0380 #endif /* ALLOW_AUTODIFF_TAMC */
                0381        ENDIF
d09cd10c60 Gael*0382
                0383 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0384 CADJ STORE hDiv(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
d09cd10c60 Gael*0385 #endif
                0386
aea29c8517 Alis*0387 C      Calculate del^2 u and del^2 v for bi-harmonic term
f0501c53d1 Jean*0388        IF (useBiharmonicVisc) THEN
3a279374db Alis*0389          CALL MOM_VI_DEL2UV(bi,bj,k,hDiv,vort3,hFacZ,
                0390      O                      del2u,del2v,
ed2dbbe83d Jean*0391      I                      myThid)
88e5e58941 Jean*0392          CALL MOM_CALC_HDIV(bi,bj,k,2,del2u,del2v,dStar,myThid)
                0393          CALL MOM_CALC_RELVORT3(bi,bj,k,
                0394      &                          del2u,del2v,hFacZ,zStar,myThid)
df1ec59c8a Jean*0395        ENDIF
                0396
d09cd10c60 Gael*0397 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0398 CADJ STORE del2u(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0399 CADJ STORE del2v(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0400 CADJ STORE dStar(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0401 CADJ STORE zStar(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
d09cd10c60 Gael*0402 #endif
                0403
df1ec59c8a Jean*0404 C---   Calculate dissipation terms for U and V equations
b0c3bd7ab0 Bayl*0405
ed2dbbe83d Jean*0406 C-     in terms of tension and strain
b0c3bd7ab0 Bayl*0407        IF (useStrainTensionVisc) THEN
ed2dbbe83d Jean*0408 C      use masked strain as if free-slip since side-drag is computed separately
f0501c53d1 Jean*0409          CALL MOM_HDISSIP( bi, bj, k,
ed2dbbe83d Jean*0410      I            tension, strain, hFacZ,
f0501c53d1 Jean*0411      I            viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
                0412      I            useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
                0413      O            guDiss, gvDiss,
                0414      I            myThid )
b0c3bd7ab0 Bayl*0415        ELSE
ed2dbbe83d Jean*0416 C-     in terms of vorticity and divergence
f0501c53d1 Jean*0417          CALL MOM_VI_HDISSIP( bi, bj, k,
ed2dbbe83d Jean*0418      I            hDiv, vort3, dStar, zStar, hFacZ,
f0501c53d1 Jean*0419      I            viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
                0420      I            useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
                0421      O            guDiss, gvDiss,
ed2dbbe83d Jean*0422      I            myThid )
07cc642809 Alis*0423        ENDIF
cab1a69b8a Jean*0424
df1ec59c8a Jean*0425 C---  Other dissipation terms in Zonal momentum equation
aea29c8517 Alis*0426
                0427 C--   Vertical flux (fVer is at upper face of "u" cell)
                0428 C     Eddy component of vertical flux (interior component only) -> vrF
79074ef66b Jean*0429        IF ( .NOT.implicitViscosity ) THEN
                0430         CALL MOM_U_RVISCFLUX(bi,bj,k+1,uVel,kappaRU,vrF,myThid)
aea29c8517 Alis*0431 C     Combine fluxes
79074ef66b Jean*0432         DO j=jMin,jMax
                0433          DO i=iMin,iMax
                0434           fVerUkp(i,j) = ArDudrFac*vrF(i,j)
                0435          ENDDO
4667e4066d Jean*0436         ENDDO
d09cd10c60 Gael*0437
                0438 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0439 CADJ STORE fVerUkp(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
d09cd10c60 Gael*0440 #endif
                0441
4667e4066d Jean*0442 C--   Tendency is minus divergence of the fluxes
f1a4cec01a Jean*0443 C     vert.visc.flx is scaled by deepFac2F (deep-atmos) and rhoFac (anelastic)
79074ef66b Jean*0444         DO j=jMin,jMax
                0445          DO i=iMin,iMax
                0446           guDiss(i,j) = guDiss(i,j)
                0447      &    -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
                0448      &    *recip_rAw(i,j,bi,bj)
                0449      &    *( fVerUkp(i,j) - fVerUkm(i,j) )*rkSign
                0450      &    *recip_deepFac2C(k)*recip_rhoFacC(k)
                0451          ENDDO
4667e4066d Jean*0452         ENDDO
79074ef66b Jean*0453        ENDIF
aea29c8517 Alis*0454
039fe61d35 Jean*0455 C-- No-slip and drag BCs appear as body forces in cell abutting topography
79074ef66b Jean*0456        IF ( no_slip_sides ) THEN
aea29c8517 Alis*0457 C-     No-slip BCs impose a drag at walls...
79074ef66b Jean*0458         CALL MOM_U_SIDEDRAG( bi, bj, k,
                0459      I           uFld, del2u, h0FacZ,
                0460      I           viscAh_Z, viscA4_Z,
                0461      I           useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
                0462      O           vF,
                0463      I           myThid )
                0464         DO j=jMin,jMax
                0465          DO i=iMin,iMax
                0466           guDiss(i,j) = guDiss(i,j)+vF(i,j)
                0467          ENDDO
aea29c8517 Alis*0468         ENDDO
79074ef66b Jean*0469        ENDIF
34e137f064 Jean*0470
aea29c8517 Alis*0471 C-    No-slip BCs impose a drag at bottom
79074ef66b Jean*0472        IF ( bottomDragTerms ) THEN
                0473 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0474 CADJ STORE KE(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
79074ef66b Jean*0475 #endif
a125ab7be7 Jean*0476         CALL MOM_U_BOTDRAG_COEFF( bi, bj, k, .TRUE.,
79074ef66b Jean*0477      I                  uFld, vFld, kappaRU, KE,
                0478      O                  cDrag,
                0479      I                  myIter, myThid )
7448700841 Mart*0480 #ifdef ALLOW_AUTODIFF_TAMC
                0481 CADJ STORE cDrag(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0482 #endif
79074ef66b Jean*0483         DO j=jMin,jMax
                0484          DO i=iMin,iMax
                0485            gUdiss(i,j) = gUdiss(i,j)
                0486      &                 - cDrag(i,j)*uFld(i,j)
                0487      &                 *_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
                0488          ENDDO
aea29c8517 Alis*0489         ENDDO
79074ef66b Jean*0490         IF ( useDiagnostics ) THEN
                0491          DO j=jMin,jMax
                0492           DO i=iMin,iMax
                0493            botDragU(i,j,bi,bj) = botDragU(i,j,bi,bj)
                0494      &                         - cDrag(i,j)*uFld(i,j)*rUnit2mass
                0495           ENDDO
                0496          ENDDO
                0497         ENDIF
                0498        ENDIF
dd49642a1d Mart*0499 #ifdef ALLOW_SHELFICE
e2d988bd46 Jean*0500        IF ( useShelfIce .AND. selectImplicitDrag.EQ.0 ) THEN
79074ef66b Jean*0501 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0502 CADJ STORE KE(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
79074ef66b Jean*0503 #endif
e2d988bd46 Jean*0504         CALL SHELFICE_U_DRAG_COEFF( bi, bj, k, .TRUE.,
                0505      I                  uFld, vFld, kappaRU, KE,
                0506      O                  cDrag,
                0507      I                  myIter, myThid )
7448700841 Mart*0508 #ifdef ALLOW_AUTODIFF_TAMC
                0509 CADJ STORE cDrag(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0510 #endif
79074ef66b Jean*0511         DO j=jMin,jMax
                0512          DO i=iMin,iMax
e2d988bd46 Jean*0513            gUdiss(i,j) = gUdiss(i,j)
                0514      &                 - cDrag(i,j)*uFld(i,j)
                0515      &                 *_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
79074ef66b Jean*0516          ENDDO
dd49642a1d Mart*0517         ENDDO
79074ef66b Jean*0518        ENDIF
dd49642a1d Mart*0519 #endif /* ALLOW_SHELFICE */
                0520
df1ec59c8a Jean*0521 C---  Other dissipation terms in Meridional momentum equation
aea29c8517 Alis*0522
                0523 C--   Vertical flux (fVer is at upper face of "v" cell)
                0524 C     Eddy component of vertical flux (interior component only) -> vrF
79074ef66b Jean*0525        IF ( .NOT.implicitViscosity ) THEN
                0526         CALL MOM_V_RVISCFLUX(bi,bj,k+1,vVel,kappaRV,vrF,myThid)
aea29c8517 Alis*0527 C     Combine fluxes -> fVerV
79074ef66b Jean*0528         DO j=jMin,jMax
                0529          DO i=iMin,iMax
                0530           fVerVkp(i,j) = ArDvdrFac*vrF(i,j)
                0531          ENDDO
4667e4066d Jean*0532         ENDDO
d09cd10c60 Gael*0533 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0534 CADJ STORE fVerVkp(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
d09cd10c60 Gael*0535 #endif
4667e4066d Jean*0536 C--   Tendency is minus divergence of the fluxes
f1a4cec01a Jean*0537 C     vert.visc.flx is scaled by deepFac2F (deep-atmos) and rhoFac (anelastic)
79074ef66b Jean*0538         DO j=jMin,jMax
                0539          DO i=iMin,iMax
                0540           gvDiss(i,j) = gvDiss(i,j)
                0541      &    -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
                0542      &    *recip_rAs(i,j,bi,bj)
                0543      &    *( fVerVkp(i,j) - fVerVkm(i,j) )*rkSign
                0544      &    *recip_deepFac2C(k)*recip_rhoFacC(k)
                0545          ENDDO
4667e4066d Jean*0546         ENDDO
79074ef66b Jean*0547        ENDIF
aea29c8517 Alis*0548
039fe61d35 Jean*0549 C-- No-slip and drag BCs appear as body forces in cell abutting topography
79074ef66b Jean*0550        IF ( no_slip_sides ) THEN
aea29c8517 Alis*0551 C-     No-slip BCs impose a drag at walls...
79074ef66b Jean*0552         CALL MOM_V_SIDEDRAG( bi, bj, k,
                0553      I           vFld, del2v, h0FacZ,
                0554      I           viscAh_Z, viscA4_Z,
                0555      I           useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
                0556      O           vF,
                0557      I           myThid )
                0558         DO j=jMin,jMax
                0559          DO i=iMin,iMax
                0560           gvDiss(i,j) = gvDiss(i,j)+vF(i,j)
                0561          ENDDO
aea29c8517 Alis*0562         ENDDO
79074ef66b Jean*0563        ENDIF
34e137f064 Jean*0564
aea29c8517 Alis*0565 C-    No-slip BCs impose a drag at bottom
79074ef66b Jean*0566        IF ( bottomDragTerms ) THEN
                0567 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0568 CADJ STORE KE(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
79074ef66b Jean*0569 #endif
a125ab7be7 Jean*0570         CALL MOM_V_BOTDRAG_COEFF( bi, bj, k, .TRUE.,
79074ef66b Jean*0571      I                  uFld, vFld, kappaRV, KE,
                0572      O                  cDrag,
                0573      I                  myIter, myThid )
7448700841 Mart*0574 #ifdef ALLOW_AUTODIFF_TAMC
                0575 CADJ STORE cDrag(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0576 #endif
79074ef66b Jean*0577         DO j=jMin,jMax
                0578          DO i=iMin,iMax
                0579            gvDiss(i,j) = gvDiss(i,j)
                0580      &                 - cDrag(i,j)*vFld(i,j)
                0581      &                 *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
                0582          ENDDO
aea29c8517 Alis*0583         ENDDO
79074ef66b Jean*0584         IF ( useDiagnostics ) THEN
                0585          DO j=jMin,jMax
                0586           DO i=iMin,iMax
                0587            botDragV(i,j,bi,bj) = botDragV(i,j,bi,bj)
                0588      &                         - cDrag(i,j)*vFld(i,j)*rUnit2mass
                0589           ENDDO
                0590          ENDDO
                0591         ENDIF
                0592        ENDIF
dd49642a1d Mart*0593 #ifdef ALLOW_SHELFICE
e2d988bd46 Jean*0594        IF ( useShelfIce .AND. selectImplicitDrag.EQ.0 ) THEN
79074ef66b Jean*0595 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0596 CADJ STORE KE(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
79074ef66b Jean*0597 #endif
e2d988bd46 Jean*0598         CALL SHELFICE_V_DRAG_COEFF( bi, bj, k, .TRUE.,
                0599      I                  uFld, vFld, kappaRV, KE,
                0600      O                  cDrag,
                0601      I                  myIter, myThid )
7448700841 Mart*0602 #ifdef ALLOW_AUTODIFF_TAMC
                0603 CADJ STORE cDrag(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0604 #endif
79074ef66b Jean*0605         DO j=jMin,jMax
                0606          DO i=iMin,iMax
e2d988bd46 Jean*0607            gvDiss(i,j) = gvDiss(i,j)
                0608      &                 - cDrag(i,j)*vFld(i,j)
                0609      &                 *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
79074ef66b Jean*0610          ENDDO
932b38363b Jean*0611         ENDDO
79074ef66b Jean*0612        ENDIF
dd49642a1d Mart*0613 #endif /* ALLOW_SHELFICE */
                0614
34e137f064 Jean*0615 C--   if (momViscosity) end of block.
                0616       ENDIF
                0617
                0618 C-    Return to standard hfacZ (min-4) and mask vort3 accordingly:
                0619 c     CALL MOM_VI_MASK_VORT3(bi,bj,k,hFacZ,r_hFacZ,vort3,myThid)
aea29c8517 Alis*0620
df1ec59c8a Jean*0621 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0622
                0623 C---  Prepare for Advection & Coriolis terms:
ed2dbbe83d Jean*0624 C-    calculate absolute vorticity
df1ec59c8a Jean*0625       IF (useAbsVorticity)
                0626      &  CALL MOM_CALC_ABSVORT3(bi,bj,k,vort3,omega3,myThid)
aea29c8517 Alis*0627
d09cd10c60 Gael*0628 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0629 CADJ STORE omega3(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
d09cd10c60 Gael*0630 #endif
                0631
aea29c8517 Alis*0632 C--   Horizontal Coriolis terms
a75214ff78 Jean*0633 c     IF (useCoriolis .AND. .NOT.useCDscheme
                0634 c    &    .AND. .NOT. useAbsVorticity) THEN
                0635 C- jmc: change it to keep the Coriolis terms when useAbsVorticity=T & momAdvection=F
711329b07b Jean*0636       IF ( useCoriolis .AND.
a75214ff78 Jean*0637      &     .NOT.( useCDscheme .OR. useAbsVorticity.AND.momAdvection )
                0638      &   ) THEN
                0639        IF (useAbsVorticity) THEN
7c3c2cec96 Jean*0640         CALL MOM_VI_U_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
                0641      &                         vFld,omega3,hFacZ,r_hFacZ,
a75214ff78 Jean*0642      &                         uCf,myThid)
7c3c2cec96 Jean*0643         CALL MOM_VI_V_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
                0644      &                         uFld,omega3,hFacZ,r_hFacZ,
a75214ff78 Jean*0645      &                         vCf,myThid)
9293d3c672 Hajo*0646 #if ( defined ALLOW_GGL90 && defined ALLOW_GGL90_LANGMUIR )
                0647        ELSEIF ( useLANGMUIR ) THEN
                0648         CALL GGL90_ADD_STOKESDRIFT(
                0649      O                 uRes, vRes,
                0650      I                 uFld, vFld, k, bi, bj, myThid )
                0651         CALL MOM_VI_CORIOLIS(bi,bj,k,uRes,vRes,hFacZ,r_hFacZ,
                0652      &                       uCf,vCf,myThid)
                0653 #endif
a75214ff78 Jean*0654        ELSE
                0655         CALL MOM_VI_CORIOLIS(bi,bj,k,uFld,vFld,hFacZ,r_hFacZ,
                0656      &                       uCf,vCf,myThid)
                0657        ENDIF
481f592dad Jean*0658        DO j=jMin,jMax
                0659         DO i=iMin,iMax
1aff67ca82 Jean*0660          gU(i,j,k,bi,bj) = uCf(i,j)
                0661          gV(i,j,k,bi,bj) = vCf(i,j)
481f592dad Jean*0662         ENDDO
aea29c8517 Alis*0663        ENDDO
f7d48db11c Jean*0664        IF ( writeDiag ) THEN
ef1e764710 Ed H*0665          IF (snapshot_mdsio) THEN
                0666            CALL WRITE_LOCAL_RL('fV','I10',1,uCf,bi,bj,k,myIter,myThid)
                0667            CALL WRITE_LOCAL_RL('fU','I10',1,vCf,bi,bj,k,myIter,myThid)
                0668          ENDIF
                0669 #ifdef ALLOW_MNC
                0670          IF (useMNC .AND. snapshot_mnc) THEN
b22b541fe9 Ed H*0671            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'fV', uCf,
cb356b4c5f Ed H*0672      &          offsets, myThid)
b22b541fe9 Ed H*0673            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'fU', vCf,
cb356b4c5f Ed H*0674      &          offsets, myThid)
ef1e764710 Ed H*0675          ENDIF
                0676 #endif /*  ALLOW_MNC  */
f7d48db11c Jean*0677        ENDIF
711329b07b Jean*0678 #ifdef ALLOW_DIAGNOSTICS
                0679        IF ( useDiagnostics ) THEN
                0680          CALL DIAGNOSTICS_FILL(uCf,'Um_Cori ',k,1,2,bi,bj,myThid)
                0681          CALL DIAGNOSTICS_FILL(vCf,'Vm_Cori ',k,1,2,bi,bj,myThid)
                0682        ENDIF
                0683 #endif /* ALLOW_DIAGNOSTICS */
4667e4066d Jean*0684       ELSE
                0685        DO j=jMin,jMax
                0686         DO i=iMin,iMax
1aff67ca82 Jean*0687          gU(i,j,k,bi,bj) = 0. _d 0
                0688          gV(i,j,k,bi,bj) = 0. _d 0
4667e4066d Jean*0689         ENDDO
                0690        ENDDO
481f592dad Jean*0691       ENDIF
                0692
d09cd10c60 Gael*0693 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0694 CADJ STORE ucf(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0695 CADJ STORE vcf(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
d09cd10c60 Gael*0696 #endif
                0697
481f592dad Jean*0698       IF (momAdvection) THEN
3add9696e1 Jean*0699 C--   Horizontal advection of relative (or absolute) vorticity
7fe6343684 Jean*0700        IF ( (highOrderVorticity.OR.upwindVorticity)
                0701      &     .AND.useAbsVorticity ) THEN
79074ef66b Jean*0702         CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,selectVortScheme,
3370e71df9 Mart*0703      &                         highOrderVorticity,upwindVorticity,
                0704      &                         vFld,omega3,r_hFacZ,
5d5fff3a81 Alis*0705      &                         uCf,myThid)
7fe6343684 Jean*0706        ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
79074ef66b Jean*0707         CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,selectVortScheme,
3370e71df9 Mart*0708      &                         highOrderVorticity, upwindVorticity,
                0709      &                         vFld,vort3, r_hFacZ,
3add9696e1 Jean*0710      &                         uCf,myThid)
7fe6343684 Jean*0711        ELSEIF ( useAbsVorticity ) THEN
3370e71df9 Mart*0712         CALL MOM_VI_U_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
                0713      &                         vFld,omega3,hFacZ,r_hFacZ,
d4c99033f5 Jean*0714      &                         uCf,myThid)
5d5fff3a81 Alis*0715        ELSE
3370e71df9 Mart*0716         CALL MOM_VI_U_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
                0717      &                         vFld,vort3, hFacZ,r_hFacZ,
5d5fff3a81 Alis*0718      &                         uCf,myThid)
                0719        ENDIF
481f592dad Jean*0720        DO j=jMin,jMax
                0721         DO i=iMin,iMax
                0722          gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
                0723         ENDDO
aea29c8517 Alis*0724        ENDDO
7fe6343684 Jean*0725        IF ( (highOrderVorticity.OR.upwindVorticity)
                0726      &     .AND.useAbsVorticity ) THEN
79074ef66b Jean*0727         CALL MOM_VI_V_CORIOLIS_C4(bi,bj,k,selectVortScheme,
3370e71df9 Mart*0728      &                         highOrderVorticity, upwindVorticity,
                0729      &                         uFld,omega3,r_hFacZ,
5d5fff3a81 Alis*0730      &                         vCf,myThid)
7fe6343684 Jean*0731        ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
79074ef66b Jean*0732         CALL MOM_VI_V_CORIOLIS_C4(bi,bj,k,selectVortScheme,
3370e71df9 Mart*0733      &                         highOrderVorticity, upwindVorticity,
                0734      &                         uFld,vort3, r_hFacZ,
3add9696e1 Jean*0735      &                         vCf,myThid)
7fe6343684 Jean*0736        ELSEIF ( useAbsVorticity ) THEN
3370e71df9 Mart*0737         CALL MOM_VI_V_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
                0738      &                         uFld,omega3,hFacZ,r_hFacZ,
d4c99033f5 Jean*0739      &                         vCf,myThid)
5d5fff3a81 Alis*0740        ELSE
3370e71df9 Mart*0741         CALL MOM_VI_V_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
                0742      &                         uFld,vort3, hFacZ,r_hFacZ,
5d5fff3a81 Alis*0743      &                         vCf,myThid)
                0744        ENDIF
481f592dad Jean*0745        DO j=jMin,jMax
                0746         DO i=iMin,iMax
                0747          gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
                0748         ENDDO
aea29c8517 Alis*0749        ENDDO
                0750
d09cd10c60 Gael*0751 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0752 CADJ STORE ucf(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
                0753 CADJ STORE vcf(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
d09cd10c60 Gael*0754 #endif
                0755
f7d48db11c Jean*0756        IF ( writeDiag ) THEN
ef1e764710 Ed H*0757          IF (snapshot_mdsio) THEN
                0758            CALL WRITE_LOCAL_RL('zV','I10',1,uCf,bi,bj,k,myIter,myThid)
                0759            CALL WRITE_LOCAL_RL('zU','I10',1,vCf,bi,bj,k,myIter,myThid)
                0760          ENDIF
                0761 #ifdef ALLOW_MNC
                0762          IF (useMNC .AND. snapshot_mnc) THEN
b22b541fe9 Ed H*0763            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'zV', uCf,
cb356b4c5f Ed H*0764      &          offsets, myThid)
b22b541fe9 Ed H*0765            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'zU', vCf,
cb356b4c5f Ed H*0766      &          offsets, myThid)
ef1e764710 Ed H*0767          ENDIF
                0768 #endif /*  ALLOW_MNC  */
f7d48db11c Jean*0769        ENDIF
ef1e764710 Ed H*0770
cab1a69b8a Jean*0771 #ifdef ALLOW_TIMEAVE
                0772        IF (taveFreq.GT.0.) THEN
                0773          CALL TIMEAVE_CUMUL_1K1T(uZetatave,vCf,deltaTClock,
                0774      &                           Nr, k, bi, bj, myThid)
                0775          CALL TIMEAVE_CUMUL_1K1T(vZetatave,uCf,deltaTClock,
                0776      &                           Nr, k, bi, bj, myThid)
                0777        ENDIF
5751fc37e0 Jean*0778 #endif /* ALLOW_TIMEAVE */
711329b07b Jean*0779 #ifdef ALLOW_DIAGNOSTICS
                0780        IF ( useDiagnostics ) THEN
                0781          CALL DIAGNOSTICS_FILL(uCf,'Um_AdvZ3',k,1,2,bi,bj,myThid)
                0782          CALL DIAGNOSTICS_FILL(vCf,'Vm_AdvZ3',k,1,2,bi,bj,myThid)
                0783        ENDIF
                0784 #endif /* ALLOW_DIAGNOSTICS */
cab1a69b8a Jean*0785
481f592dad Jean*0786 C--   Vertical shear terms (-w*du/dr & -w*dv/dr)
6d4b75df68 Jean*0787        IF ( .NOT. momImplVertAdv ) THEN
f1a4cec01a Jean*0788         CALL MOM_VI_U_VERTSHEAR(bi,bj,k,uVel,wVel,uCf,myThid)
6d4b75df68 Jean*0789         DO j=jMin,jMax
                0790          DO i=iMin,iMax
                0791           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
                0792          ENDDO
481f592dad Jean*0793         ENDDO
f1a4cec01a Jean*0794         CALL MOM_VI_V_VERTSHEAR(bi,bj,k,vVel,wVel,vCf,myThid)
6d4b75df68 Jean*0795         DO j=jMin,jMax
                0796          DO i=iMin,iMax
                0797           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
                0798          ENDDO
481f592dad Jean*0799         ENDDO
711329b07b Jean*0800 #ifdef ALLOW_DIAGNOSTICS
                0801         IF ( useDiagnostics ) THEN
                0802          CALL DIAGNOSTICS_FILL(uCf,'Um_AdvRe',k,1,2,bi,bj,myThid)
                0803          CALL DIAGNOSTICS_FILL(vCf,'Vm_AdvRe',k,1,2,bi,bj,myThid)
                0804         ENDIF
                0805 #endif /* ALLOW_DIAGNOSTICS */
6d4b75df68 Jean*0806        ENDIF
aea29c8517 Alis*0807
                0808 C--   Bernoulli term
f1a4cec01a Jean*0809        CALL MOM_VI_U_GRAD_KE(bi,bj,k,KE,uCf,myThid)
481f592dad Jean*0810        DO j=jMin,jMax
                0811         DO i=iMin,iMax
                0812          gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
                0813         ENDDO
aea29c8517 Alis*0814        ENDDO
f1a4cec01a Jean*0815        CALL MOM_VI_V_GRAD_KE(bi,bj,k,KE,vCf,myThid)
481f592dad Jean*0816        DO j=jMin,jMax
                0817         DO i=iMin,iMax
                0818          gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
                0819         ENDDO
                0820        ENDDO
f7d48db11c Jean*0821        IF ( writeDiag ) THEN
ef1e764710 Ed H*0822          IF (snapshot_mdsio) THEN
                0823            CALL WRITE_LOCAL_RL('KEx','I10',1,uCf,bi,bj,k,myIter,myThid)
                0824            CALL WRITE_LOCAL_RL('KEy','I10',1,vCf,bi,bj,k,myIter,myThid)
                0825          ENDIF
                0826 #ifdef ALLOW_MNC
                0827          IF (useMNC .AND. snapshot_mnc) THEN
b22b541fe9 Ed H*0828            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'KEx', uCf,
cb356b4c5f Ed H*0829      &          offsets, myThid)
b22b541fe9 Ed H*0830            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'KEy', vCf,
cb356b4c5f Ed H*0831      &          offsets, myThid)
df1ec59c8a Jean*0832          ENDIF
ef1e764710 Ed H*0833 #endif /*  ALLOW_MNC  */
f7d48db11c Jean*0834        ENDIF
                0835
481f592dad Jean*0836 C--   end if momAdvection
                0837       ENDIF
                0838
3daafce20b Jean*0839 C--   3.D Coriolis term (horizontal momentum, Eastward component: -fprime*w)
427e24e121 Jean*0840       IF ( select3dCoriScheme.GE.1 ) THEN
039fe61d35 Jean*0841         CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,uCf,myThid)
                0842         DO j=jMin,jMax
                0843          DO i=iMin,iMax
                0844           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
                0845          ENDDO
                0846         ENDDO
427e24e121 Jean*0847        IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
039fe61d35 Jean*0848 C-     presently, non zero angleSinC array only supported with Curvilinear-Grid
                0849         CALL MOM_V_CORIOLIS_NH(bi,bj,k,wVel,vCf,myThid)
                0850         DO j=jMin,jMax
                0851          DO i=iMin,iMax
                0852           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
                0853          ENDDO
                0854         ENDDO
                0855        ENDIF
                0856       ENDIF
                0857
                0858 C--   Non-Hydrostatic (spherical) metric terms
                0859       IF ( useNHMTerms ) THEN
                0860        CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,uCf,myThid)
                0861        DO j=jMin,jMax
                0862         DO i=iMin,iMax
                0863          gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
                0864         ENDDO
                0865        ENDDO
                0866        CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,vCf,myThid)
                0867        DO j=jMin,jMax
                0868         DO i=iMin,iMax
                0869          gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
                0870         ENDDO
                0871        ENDDO
                0872       ENDIF
df1ec59c8a Jean*0873
481f592dad Jean*0874 C--   Set du/dt & dv/dt on boundaries to zero
aea29c8517 Alis*0875       DO j=jMin,jMax
                0876        DO i=iMin,iMax
481f592dad Jean*0877         gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
                0878         gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
aea29c8517 Alis*0879        ENDDO
                0880       ENDDO
481f592dad Jean*0881
5751fc37e0 Jean*0882 #ifdef ALLOW_DEBUG
8830b8f970 Jean*0883       IF ( debugLevel .GE. debLevC
5751fc37e0 Jean*0884      &   .AND. k.EQ.4 .AND. myIter.EQ.nIter0
                0885      &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
                0886      &   .AND. useCubedSphereExchange ) THEN
e85db1faec Jean*0887         CALL DEBUG_CS_CORNER_UV( ' uDiss,vDiss from MOM_VECINV',
4667e4066d Jean*0888      &             guDiss,gvDiss, k, standardMessageUnit,bi,bj,myThid )
5751fc37e0 Jean*0889       ENDIF
                0890 #endif /* ALLOW_DEBUG */
aea29c8517 Alis*0891
f7d48db11c Jean*0892       IF ( writeDiag ) THEN
ed2dbbe83d Jean*0893         IF (useBiharmonicVisc) THEN
                0894          CALL WRITE_LOCAL_RL( 'del2u', 'I10', 1, del2u,
                0895      &                        bi,bj,k, myIter, myThid )
                0896          CALL WRITE_LOCAL_RL( 'del2v', 'I10', 1, del2v,
                0897      &                        bi,bj,k, myIter, myThid )
                0898          CALL WRITE_LOCAL_RL( 'dStar', 'I10', 1, dStar,
                0899      &                        bi,bj,k, myIter, myThid )
                0900          CALL WRITE_LOCAL_RL( 'zStar', 'I10', 1, zStar,
                0901      &                        bi,bj,k, myIter, myThid )
                0902         ENDIF
ef1e764710 Ed H*0903         IF (snapshot_mdsio) THEN
df1ec59c8a Jean*0904          CALL WRITE_LOCAL_RL('W3','I10',1,omega3, bi,bj,k,myIter,myThid)
ed2dbbe83d Jean*0905          CALL WRITE_LOCAL_RL('Z3','I10',1,vort3BC,bi,bj,k,myIter,myThid)
df1ec59c8a Jean*0906          CALL WRITE_LOCAL_RL('KE','I10',1,KE,     bi,bj,k,myIter,myThid)
                0907          CALL WRITE_LOCAL_RL('D', 'I10',1,hDiv,   bi,bj,k,myIter,myThid)
                0908          CALL WRITE_LOCAL_RL('Dt','I10',1,tension,bi,bj,k,myIter,myThid)
ed2dbbe83d Jean*0909          CALL WRITE_LOCAL_RL( 'Ds', 'I10', 1, strainBC,
                0910      &                        bi,bj,k, myIter, myThid )
df1ec59c8a Jean*0911          CALL WRITE_LOCAL_RL('Du','I10',1,guDiss, bi,bj,k,myIter,myThid)
                0912          CALL WRITE_LOCAL_RL('Dv','I10',1,gvDiss, bi,bj,k,myIter,myThid)
ef1e764710 Ed H*0913         ENDIF
                0914 #ifdef ALLOW_MNC
                0915         IF (useMNC .AND. snapshot_mnc) THEN
df1ec59c8a Jean*0916           CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'W3',omega3,
                0917      &          offsets, myThid)
ed2dbbe83d Jean*0918           CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Z3',vort3BC,
                0919      &          offsets, myThid)
df1ec59c8a Jean*0920           CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'KE',KE,
                0921      &          offsets, myThid)
                0922           CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'D', hDiv,
cb356b4c5f Ed H*0923      &          offsets, myThid)
b22b541fe9 Ed H*0924           CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dt',tension,
cb356b4c5f Ed H*0925      &          offsets, myThid)
ed2dbbe83d Jean*0926           CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Ds',strainBC,
                0927      &          offsets, myThid)
b22b541fe9 Ed H*0928           CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Du',guDiss,
cb356b4c5f Ed H*0929      &          offsets, myThid)
b22b541fe9 Ed H*0930           CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dv',gvDiss,
cb356b4c5f Ed H*0931      &          offsets, myThid)
ef1e764710 Ed H*0932         ENDIF
                0933 #endif /*  ALLOW_MNC  */
3a279374db Alis*0934       ENDIF
3add9696e1 Jean*0935
711329b07b Jean*0936 #ifdef ALLOW_DIAGNOSTICS
                0937       IF ( useDiagnostics ) THEN
ed2dbbe83d Jean*0938         CALL DIAGNOSTICS_FILL(vort3BC,'momVort3',k,1,2,bi,bj,myThid)
df1ec59c8a Jean*0939         CALL DIAGNOSTICS_FILL(KE,     'momKE   ',k,1,2,bi,bj,myThid)
711329b07b Jean*0940        IF (momViscosity) THEN
df1ec59c8a Jean*0941         CALL DIAGNOSTICS_FILL(hDiv,   'momHDiv ',k,1,2,bi,bj,myThid)
711329b07b Jean*0942        ENDIF
ed2dbbe83d Jean*0943        IF ( useVariableVisc .OR. useStrainTensionVisc ) THEN
                0944         CALL DIAGNOSTICS_FILL(tension, 'Tension ',k,1,2,bi,bj,myThid)
                0945         CALL DIAGNOSTICS_FILL(strainBC,'Strain  ',k,1,2,bi,bj,myThid)
f59d76b0dd Ed D*0946 C         stretching will be zero, unless using QG Leith
a125ab7be7 Jean*0947         IF ( viscC2LeithQG.NE.zeroRL ) THEN
f59d76b0dd Ed D*0948           CALL DIAGNOSTICS_FILL(stretching,
                0949      I                          'Stretch ',k,1,2,bi,bj,myThid)
                0950         ENDIF
ed2dbbe83d Jean*0951        ENDIF
07e17fa184 Jean*0952         CALL DIAGNOSTICS_FILL(gU(1-OLx,1-OLy,k,bi,bj),
df1ec59c8a Jean*0953      &                                'Um_Advec',k,1,2,bi,bj,myThid)
07e17fa184 Jean*0954         CALL DIAGNOSTICS_FILL(gV(1-OLx,1-OLy,k,bi,bj),
df1ec59c8a Jean*0955      &                                'Vm_Advec',k,1,2,bi,bj,myThid)
711329b07b Jean*0956       ENDIF
                0957 #endif /* ALLOW_DIAGNOSTICS */
                0958
99bc607d7a Ed H*0959 #endif /* ALLOW_MOM_VECINV */
cab1a69b8a Jean*0960
aea29c8517 Alis*0961       RETURN
                0962       END