/[MITgcm]/MITgcm/pkg/mom_vecinv/mom_vecinv.F
ViewVC logotype

Diff of /MITgcm/pkg/mom_vecinv/mom_vecinv.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.4 by jmc, Sat Feb 8 02:10:57 2003 UTC revision 1.80 by mlosch, Fri Apr 28 17:17:14 2017 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "MOM_VECINV_OPTIONS.h"
5    #ifdef ALLOW_AUTODIFF
6        SUBROUTINE MOM_VECINV(  # include "AUTODIFF_OPTIONS.h"
7       I        bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,  #endif
8       I        dPhiHydX,dPhiHydY,KappaRU,KappaRV,  #ifdef ALLOW_MOM_COMMON
9       U        fVerU, fVerV,  # include "MOM_COMMON_OPTIONS.h"
10       I        myCurrentTime, myIter, myThid)  #endif
11  C     /==========================================================\  
12          SUBROUTINE MOM_VECINV(
13         I        bi,bj,k,iMin,iMax,jMin,jMax,
14         I        kappaRU, kappaRV,
15         I        fVerUkm, fVerVkm,
16         O        fVerUkp, fVerVkp,
17         O        guDiss, gvDiss,
18         I        myTime, myIter, myThid )
19    C     *==========================================================*
20  C     | S/R MOM_VECINV                                           |  C     | S/R MOM_VECINV                                           |
21  C     | o Form the right hand-side of the momentum equation.     |  C     | o Form the right hand-side of the momentum equation.     |
22  C     |==========================================================|  C     *==========================================================*
23  C     | Terms are evaluated one layer at a time working from     |  C     | Terms are evaluated one layer at a time working from     |
24  C     | the bottom to the top. The vertically integrated         |  C     | the bottom to the top. The vertically integrated         |
25  C     | barotropic flow tendency term is evluated by summing the |  C     | barotropic flow tendency term is evluated by summing the |
# Line 22  C     | for the diffusion equation bc wi Line 30  C     | for the diffusion equation bc wi
30  C     | form produces a diffusive flux that does not scale with  |  C     | form produces a diffusive flux that does not scale with  |
31  C     | open-area. Need to do something to solidfy this and to   |  C     | open-area. Need to do something to solidfy this and to   |
32  C     | deal "properly" with thin walls.                         |  C     | deal "properly" with thin walls.                         |
33  C     \==========================================================/  C     *==========================================================*
34        IMPLICIT NONE        IMPLICIT NONE
35    
36  C     == Global variables ==  C     == Global variables ==
37  #include "SIZE.h"  #include "SIZE.h"
 #include "DYNVARS.h"  
38  #include "EEPARAMS.h"  #include "EEPARAMS.h"
39  #include "PARAMS.h"  #include "PARAMS.h"
40  #include "GRID.h"  #include "GRID.h"
41    #include "SURFACE.h"
42    #include "DYNVARS.h"
43    #ifdef ALLOW_MOM_COMMON
44    # include "MOM_VISC.h"
45    #endif
46    #ifdef ALLOW_TIMEAVE
47    # include "TIMEAVE_STATV.h"
48    #endif
49    #ifdef ALLOW_MNC
50    # include "MNC_PARAMS.h"
51    #endif
52    #ifdef ALLOW_AUTODIFF_TAMC
53    # include "tamc.h"
54    # include "tamc_keys.h"
55    #endif
56    
57  C     == Routine arguments ==  C     == Routine arguments ==
58  C     fVerU   - Flux of momentum in the vertical  C     bi,bj   :: current tile indices
59  C     fVerV     direction out of the upper face of a cell K  C     k       :: current vertical level
60  C               ( flux into the cell above ).  C     iMin,iMax,jMin,jMax :: loop ranges
61  C     dPhiHydX,Y :: Gradient (X & Y dir.) of Hydrostatic Potential  C     fVerU   :: Flux of momentum in the vertical direction, out of the upper
62  C     bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation  C     fVerV   :: face of a cell k ( flux into the cell above ).
63  C                                      results will be set.  C     fVerUkm :: vertical viscous flux of U, interface above (k-1/2)
64  C     kUp, kDown                     - Index for upper and lower layers.  C     fVerVkm :: vertical viscous flux of V, interface above (k-1/2)
65  C     myThid - Instance number for this innvocation of CALC_MOM_RHS  C     fVerUkp :: vertical viscous flux of U, interface below (k+1/2)
66        _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)  C     fVerVkp :: vertical viscous flux of V, interface below (k+1/2)
67        _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)  
68        _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     guDiss  :: dissipation tendency (all explicit terms), u component
69        _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     gvDiss  :: dissipation tendency (all explicit terms), v component
70        _RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  C     myTime  :: current time
71        _RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  C     myIter  :: current time-step number
72        INTEGER kUp,kDown  C     myThid  :: my Thread Id number
73        _RL     myCurrentTime        INTEGER bi,bj,k
74          INTEGER iMin,iMax,jMin,jMax
75          _RL kappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
76          _RL kappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
77          _RL fVerUkm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78          _RL fVerVkm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79          _RL fVerUkp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80          _RL fVerVkp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81          _RL guDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82          _RL gvDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83          _RL     myTime
84        INTEGER myIter        INTEGER myIter
85        INTEGER myThid        INTEGER myThid
86        INTEGER bi,bj,iMin,iMax,jMin,jMax  
87    #ifdef ALLOW_MOM_VECINV
88    
89  C     == Functions ==  C     == Functions ==
90        LOGICAL  DIFFERENT_MULTIPLE        LOGICAL  DIFFERENT_MULTIPLE
91        EXTERNAL DIFFERENT_MULTIPLE        EXTERNAL DIFFERENT_MULTIPLE
92    
93  C     == Local variables ==  C     == Local variables ==
94        _RL      aF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     strainBC :: same as strain but account for no-slip BC
95    C     vort3BC  :: same as vort3  but account for no-slip BC
96        _RL      vF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL      vF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL      vrF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL      vrF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RL      uCf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL      uCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99        _RL      vCf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL      vCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RL      mT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS hFacZ   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL      pF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS h0FacZ  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL del2u(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS r_hFacZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103        _RL del2v(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104        _RL tension(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105        _RL strain(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL del2u   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106        _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL del2v   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107        _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dStar   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108        _RS xA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL zStar   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109        _RS yA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL tension (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL strain  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL strainBC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112        _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KE      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113        _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL omega3  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114        _RL dStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vort3   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115        _RL zStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vort3BC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116        _RL uDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL hDiv    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117        _RL vDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118  C     I,J,K - Loop counters        _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119        INTEGER i,j,k        _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120  C     rVelMaskOverride - Factor for imposing special surface boundary conditions        _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121  C                        ( set according to free-surface condition ).  C     i,j    :: Loop counters
122  C     hFacROpen        - Lopped cell factos used tohold fraction of open        INTEGER i,j
123  C     hFacRClosed        and closed cell wall.  C     xxxFac :: On-off tracer parameters used for switching terms off.
       _RL  rVelMaskOverride  
 C     xxxFac - On-off tracer parameters used for switching terms off.  
       _RL  uDudxFac  
       _RL  AhDudxFac  
       _RL  A4DuxxdxFac  
       _RL  vDudyFac  
       _RL  AhDudyFac  
       _RL  A4DuyydyFac  
       _RL  rVelDudrFac  
124        _RL  ArDudrFac        _RL  ArDudrFac
       _RL  fuFac  
       _RL  phxFac  
       _RL  mtFacU  
       _RL  uDvdxFac  
       _RL  AhDvdxFac  
       _RL  A4DvxxdxFac  
       _RL  vDvdyFac  
       _RL  AhDvdyFac  
       _RL  A4DvyydyFac  
       _RL  rVelDvdrFac  
125        _RL  ArDvdrFac        _RL  ArDvdrFac
126        _RL  fvFac        _RL  sideMaskFac
       _RL  phyFac  
       _RL  vForcFac  
       _RL  mtFacV  
       INTEGER km1,kp1  
       _RL wVelBottomOverride  
127        LOGICAL bottomDragTerms        LOGICAL bottomDragTerms
128        _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        LOGICAL writeDiag
129        _RL omega3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #ifdef ALLOW_AUTODIFF_TAMC
130        _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        INTEGER imomkey
131        _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #endif
132    
133        km1=MAX(1,k-1)  #ifdef ALLOW_MNC
134        kp1=MIN(Nr,k+1)        INTEGER offsets(9)
135        rVelMaskOverride=1.        CHARACTER*(1) pf
136        IF ( k .EQ. 1 ) rVelMaskOverride=freeSurfFac  #endif
137        wVelBottomOverride=1.  
138        IF (k.EQ.Nr) wVelBottomOverride=0.  #ifdef ALLOW_AUTODIFF
139    C--   only the kDown part of fverU/V is set in this subroutine
140  C     Initialise intermediate terms  C--   the kUp is still required
141        DO J=1-OLy,sNy+OLy  C--   In the case of mom_fluxform kUp is set as well
142         DO I=1-OLx,sNx+OLx  C--   (at least in part)
143          aF(i,j)   = 0.        fVerUkm(1,1) = fVerUkm(1,1)
144          vF(i,j)   = 0.        fVerVkm(1,1) = fVerVkm(1,1)
145          vrF(i,j)  = 0.  #endif
146    
147    #ifdef ALLOW_AUTODIFF_TAMC
148              act0 = k - 1
149              max0 = Nr
150              act1 = bi - myBxLo(myThid)
151              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
152              act2 = bj - myByLo(myThid)
153              max2 = myByHi(myThid) - myByLo(myThid) + 1
154              act3 = myThid - 1
155              max3 = nTx*nTy
156              act4 = ikey_dynamics - 1
157              imomkey = (act0 + 1)
158         &                    + act1*max0
159         &                    + act2*max0*max1
160         &                    + act3*max0*max1*max2
161         &                    + act4*max0*max1*max2*max3
162    #endif /* ALLOW_AUTODIFF_TAMC */
163    
164          writeDiag = DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock)
165    
166    #ifdef ALLOW_MNC
167          IF (useMNC .AND. snapshot_mnc .AND. writeDiag) THEN
168            IF ( writeBinaryPrec .EQ. precFloat64 ) THEN
169              pf(1:1) = 'D'
170            ELSE
171              pf(1:1) = 'R'
172            ENDIF
173            IF ((bi .EQ. 1).AND.(bj .EQ. 1).AND.(k .EQ. 1)) THEN
174              CALL MNC_CW_SET_UDIM('mom_vi', -1, myThid)
175              CALL MNC_CW_RL_W_S('D','mom_vi',0,0,'T',myTime,myThid)
176              CALL MNC_CW_SET_UDIM('mom_vi', 0, myThid)
177              CALL MNC_CW_I_W_S('I','mom_vi',0,0,'iter',myIter,myThid)
178            ENDIF
179            DO i = 1,9
180              offsets(i) = 0
181            ENDDO
182            offsets(3) = k
183    c       write(*,*) 'offsets = ',(offsets(i),i=1,9)
184          ENDIF
185    #endif /*  ALLOW_MNC  */
186    
187    C--   Initialise intermediate terms
188          DO j=1-OLy,sNy+OLy
189           DO i=1-OLx,sNx+OLx
190            vF(i,j)    = 0.
191            vrF(i,j)   = 0.
192          uCf(i,j)   = 0.          uCf(i,j)   = 0.
193          vCf(i,j)   = 0.          vCf(i,j)   = 0.
         mT(i,j)   = 0.  
         pF(i,j)   = 0.  
194          del2u(i,j) = 0.          del2u(i,j) = 0.
195          del2v(i,j) = 0.          del2v(i,j) = 0.
196          dStar(i,j) = 0.          dStar(i,j) = 0.
197          zStar(i,j) = 0.          zStar(i,j) = 0.
198          uDiss(i,j) = 0.          guDiss(i,j)= 0.
199          vDiss(i,j) = 0.          gvDiss(i,j)= 0.
200          vort3(i,j) = 0.          vort3(i,j) = 0.
201          omega3(i,j) = 0.          omega3(i,j)= 0.
202          ke(i,j) = 0.          KE(i,j)    = 0.
203    C-    need to initialise hDiv for MOM_VI_DEL2UV(call FILL_CS_CORNER_TR_RL)
204            hDiv(i,j)  = 0.
205    c       viscAh_Z(i,j) = 0.
206    c       viscAh_D(i,j) = 0.
207    c       viscA4_Z(i,j) = 0.
208    c       viscA4_D(i,j) = 0.
209            strain(i,j)  = 0. _d 0
210            strainBC(i,j)= 0. _d 0
211            tension(i,j) = 0. _d 0
212    #ifdef ALLOW_AUTODIFF
213            hFacZ(i,j)   = 0. _d 0
214    #endif
215         ENDDO         ENDDO
216        ENDDO        ENDDO
217    
218  C--   Term by term tracer parmeters  C--   Term by term tracer parmeters
219  C     o U momentum equation  C     o U momentum equation
       uDudxFac     = afFacMom*1.  
       AhDudxFac    = vfFacMom*1.  
       A4DuxxdxFac  = vfFacMom*1.  
       vDudyFac     = afFacMom*1.  
       AhDudyFac    = vfFacMom*1.  
       A4DuyydyFac  = vfFacMom*1.  
       rVelDudrFac  = afFacMom*1.  
220        ArDudrFac    = vfFacMom*1.        ArDudrFac    = vfFacMom*1.
       mTFacU       = mtFacMom*1.  
       fuFac        = cfFacMom*1.  
       phxFac       = pfFacMom*1.  
221  C     o V momentum equation  C     o V momentum equation
       uDvdxFac     = afFacMom*1.  
       AhDvdxFac    = vfFacMom*1.  
       A4DvxxdxFac  = vfFacMom*1.  
       vDvdyFac     = afFacMom*1.  
       AhDvdyFac    = vfFacMom*1.  
       A4DvyydyFac  = vfFacMom*1.  
       rVelDvdrFac  = afFacMom*1.  
222        ArDvdrFac    = vfFacMom*1.        ArDvdrFac    = vfFacMom*1.
223        mTFacV       = mtFacMom*1.  
224        fvFac        = cfFacMom*1.  C note: using standard stencil (no mask) results in under-estimating
225        phyFac       = pfFacMom*1.  C       vorticity at a no-slip boundary by a factor of 2 = sideDragFactor
226        vForcFac     = foFacMom*1.        IF ( no_slip_sides ) THEN
227            sideMaskFac = sideDragFactor
       IF (     no_slip_bottom  
      &    .OR. bottomDragQuadratic.NE.0.  
      &    .OR. bottomDragLinear.NE.0.) THEN  
        bottomDragTerms=.TRUE.  
228        ELSE        ELSE
229         bottomDragTerms=.FALSE.          sideMaskFac = 0. _d 0
230        ENDIF        ENDIF
231    
232  C-- with stagger time stepping, grad Phi_Hyp is directly incoporated in TIMESTEP        IF ( selectImplicitDrag.EQ.0 .AND.
233        IF (staggerTimeStep) THEN       &      (  no_slip_bottom
234          phxFac = 0.       &    .OR. selectBotDragQuadr.GE.0
235          phyFac = 0.       &    .OR. bottomDragLinear.NE.0. ) ) THEN
236           bottomDragTerms=.TRUE.
237          ELSE
238           bottomDragTerms=.FALSE.
239        ENDIF        ENDIF
240    
241  C--   Calculate open water fraction at vorticity points  C--   Calculate open water fraction at vorticity points
242        CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)        CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
243    
 C---- Calculate common quantities used in both U and V equations  
 C     Calculate tracer cell face open areas  
       DO j=1-OLy,sNy+OLy  
        DO i=1-OLx,sNx+OLx  
         xA(i,j) = _dyG(i,j,bi,bj)  
      &   *drF(k)*_hFacW(i,j,k,bi,bj)  
         yA(i,j) = _dxG(i,j,bi,bj)  
      &   *drF(k)*_hFacS(i,j,k,bi,bj)  
        ENDDO  
       ENDDO  
   
244  C     Make local copies of horizontal flow field  C     Make local copies of horizontal flow field
245        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
246         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
# Line 212  C     Make local copies of horizontal fl Line 249  C     Make local copies of horizontal fl
249         ENDDO         ENDDO
250        ENDDO        ENDDO
251    
252  C     Calculate velocity field "volume transports" through tracer cell faces.  #ifdef ALLOW_AUTODIFF_TAMC
253    CADJ STORE ufld(:,:) =
254    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
255    CADJ STORE vfld(:,:) =
256    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
257    CADJ STORE hFacZ(:,:) =
258    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
259    CADJ STORE r_hFacZ(:,:) =
260    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
261    CADJ STORE fverukm(:,:) =
262    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
263    CADJ STORE fvervkm(:,:) =
264    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
265    #endif
266    
267    C note (jmc) : Dissipation and Vort3 advection do not necesary
268    C              use the same maskZ (and hFacZ)  => needs 2 call(s)
269    c     CALL MOM_VI_HFACZ_DISS(bi,bj,k,hFacZ,r_hFacZ,myThid)
270    
271          CALL MOM_CALC_KE(bi,bj,k,selectKEscheme,uFld,vFld,KE,myThid)
272    
273          CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
274    
275    #ifdef ALLOW_AUTODIFF_TAMC
276    CADJ STORE ke(:,:) =
277    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
278    CADJ STORE vort3(:,:) =
279    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
280    CADJ STORE vort3bc(:,:) =
281    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
282    #endif
283    
284    C-    mask vort3 and account for no-slip / free-slip BC in vort3BC:
285        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
286         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
287          uTrans(i,j) = uFld(i,j)*xA(i,j)           vort3BC(i,j) = vort3(i,j)
288          vTrans(i,j) = vFld(i,j)*yA(i,j)           IF ( hFacZ(i,j).EQ.zeroRS ) THEN
289               vort3BC(i,j) = sideMaskFac*vort3BC(i,j)
290               vort3(i,j)   = 0.
291             ENDIF
292         ENDDO         ENDDO
293        ENDDO        ENDDO
294    
295        CALL MOM_VI_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid)  #ifdef ALLOW_AUTODIFF_TAMC
296    CADJ STORE vort3(:,:) =
297    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
298    CADJ STORE vort3bc(:,:) =
299    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
300    #endif
301    
302          IF (momViscosity) THEN
303    C--    For viscous term, compute horizontal divergence, tension & strain
304    C      and mask relative vorticity (free-slip case):
305    
306        CALL MOM_VI_CALC_HDIV(bi,bj,k,uFld,vFld,hDiv,myThid)         DO j=1-OLy,sNy+OLy
307            DO i=1-OLx,sNx+OLx
308              h0FacZ(i,j) = hFacZ(i,j)
309            ENDDO
310           ENDDO
311    #ifdef NONLIN_FRSURF
312           IF ( no_slip_sides .AND. nonlinFreeSurf.GT.0 ) THEN
313            DO j=2-OLy,sNy+OLy
314             DO i=2-OLx,sNx+OLx
315              h0FacZ(i,j) = MIN(
316         &       MIN( h0FacW(i,j,k,bi,bj), h0FacW(i,j-1,k,bi,bj) ),
317         &       MIN( h0FacS(i,j,k,bi,bj), h0FacS(i-1,j,k,bi,bj) ) )
318             ENDDO
319            ENDDO
320           ENDIF
321    #endif /* NONLIN_FRSURF */
322    
323        CALL MOM_VI_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)  #ifdef ALLOW_AUTODIFF_TAMC
324    CADJ STORE h0FacZ(:,:) =
325    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
326    CADJ STORE hFacZ(:,:) =
327    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
328    #endif
329    
330           CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
331    
332           IF ( useVariableVisc .OR. useStrainTensionVisc ) THEN
333            CALL MOM_CALC_TENSION( bi,bj,k,uFld,vFld,tension,myThid )
334            CALL MOM_CALC_STRAIN( bi,bj,k,uFld,vFld,hFacZ,strain,myThid )
335    C-    mask strain and account for no-slip / free-slip BC in strainBC:
336            DO j=1-OLy,sNy+OLy
337             DO i=1-OLx,sNx+OLx
338               strainBC(i,j) = strain(i,j)
339               IF ( hFacZ(i,j).EQ.zeroRS ) THEN
340                 strainBC(i,j) = sideMaskFac*strainBC(i,j)
341                 strain(i,j)   = 0.
342               ENDIF
343             ENDDO
344            ENDDO
345           ENDIF
346    
347        CALL MOM_VI_CALC_ABSVORT3(bi,bj,k,vort3,omega3,myThid)  #ifdef ALLOW_AUTODIFF_TAMC
348    CADJ STORE hdiv(:,:) =
349    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
350    CADJ STORE tension(:,:) =
351    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
352    CADJ STORE strain(:,:) =
353    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
354    CADJ STORE strainbc(:,:) =
355    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
356    #endif
357    
358    C--    Calculate Lateral Viscosities
359           DO j=1-OLy,sNy+OLy
360            DO i=1-OLx,sNx+OLx
361             viscAh_D(i,j) = viscAhD
362             viscAh_Z(i,j) = viscAhZ
363             viscA4_D(i,j) = viscA4D
364             viscA4_Z(i,j) = viscA4Z
365            ENDDO
366           ENDDO
367           IF ( useVariableVisc ) THEN
368    C-     uses vort3BC & strainBC which account for no-slip / free-slip BC
369             CALL MOM_CALC_VISC( bi, bj, k,
370         O            viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
371         I            hDiv, vort3BC, tension, strainBC, KE, hfacZ,
372         I            myThid )
373           ENDIF
374    
375    #ifdef ALLOW_AUTODIFF_TAMC
376    CADJ STORE viscAh_Z(:,:) =
377    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
378    CADJ STORE viscAh_D(:,:) =
379    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
380    CADJ STORE viscA4_Z(:,:) =
381    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
382    CADJ STORE viscA4_D(:,:) =
383    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
384    #endif
385    
386    #ifdef ALLOW_AUTODIFF_TAMC
387    CADJ STORE hDiv(:,:) =
388    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
389    CADJ STORE vort3(:,:) =
390    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
391    CADJ STORE hFacZ(:,:) =
392    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
393    #endif
394    
       IF (momViscosity) THEN  
395  C      Calculate del^2 u and del^2 v for bi-harmonic term  C      Calculate del^2 u and del^2 v for bi-harmonic term
396         IF (viscA4.NE.0.) THEN         IF (useBiharmonicVisc) THEN
397           CALL MOM_VI_DEL2UV(bi,bj,k,hDiv,vort3,hFacZ,           CALL MOM_VI_DEL2UV(bi,bj,k,hDiv,vort3,hFacZ,
398       O                      del2u,del2v,       O                      del2u,del2v,
399       &                      myThid)       I                      myThid)
400           CALL MOM_VI_CALC_HDIV(bi,bj,k,del2u,del2v,dStar,myThid)  #ifdef ALLOW_AUTODIFF_TAMC
401           CALL MOM_VI_CALC_RELVORT3(  CADJ STORE del2u(:,:) =
402       &                         bi,bj,k,del2u,del2v,hFacZ,zStar,myThid)  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
403         ENDIF  CADJ STORE del2v(:,:) =
404  C      Calculate dissipation terms for U and V equations  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
405  C      in terms of vorticity and divergence  #endif
406         IF (viscAh.NE.0. .OR. viscA4.NE.0.) THEN           CALL MOM_CALC_HDIV(bi,bj,k,2,del2u,del2v,dStar,myThid)
407           CALL MOM_VI_HDISSIP(bi,bj,k,hDiv,vort3,hFacZ,dStar,zStar,           CALL MOM_CALC_RELVORT3(bi,bj,k,
408       O                       uDiss,vDiss,       &                          del2u,del2v,hFacZ,zStar,myThid)
409       &                       myThid)         ENDIF
410         ENDIF  
411  C      or in terms of tension and strain  #ifdef ALLOW_AUTODIFF_TAMC
412         IF (viscAstrain.NE.0. .OR. viscAtension.NE.0.) THEN  CADJ STORE del2u(:,:) =
413           CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
414       O                         tension,  CADJ STORE del2v(:,:) =
415       I                         myThid)  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
416           CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,  CADJ STORE dStar(:,:) =
417       O                        strain,  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
418       I                        myThid)  CADJ STORE zStar(:,:) =
419           CALL MOM_HDISSIP(bi,bj,k,  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
420       I                    tension,strain,hFacZ,viscAtension,viscAstrain,  #endif
421       O                    uDiss,vDiss,  
422       I                    myThid)  C---   Calculate dissipation terms for U and V equations
423    
424    C-     in terms of tension and strain
425           IF (useStrainTensionVisc) THEN
426    C      use masked strain as if free-slip since side-drag is computed separately
427             CALL MOM_HDISSIP( bi, bj, k,
428         I            tension, strain, hFacZ,
429         I            viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
430         I            useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
431         O            guDiss, gvDiss,
432         I            myThid )
433           ELSE
434    C-     in terms of vorticity and divergence
435             CALL MOM_VI_HDISSIP( bi, bj, k,
436         I            hDiv, vort3, dStar, zStar, hFacZ,
437         I            viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
438         I            useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
439         O            guDiss, gvDiss,
440         I            myThid )
441         ENDIF         ENDIF
       ENDIF  
442    
443  C---- Zonal momentum equation starts here  C---  Other dissipation terms in Zonal momentum equation
444    
445  C--   Vertical flux (fVer is at upper face of "u" cell)  C--   Vertical flux (fVer is at upper face of "u" cell)
   
446  C     Eddy component of vertical flux (interior component only) -> vrF  C     Eddy component of vertical flux (interior component only) -> vrF
447        IF (momViscosity.AND..NOT.implicitViscosity)        IF ( .NOT.implicitViscosity ) THEN
448       & CALL MOM_U_RVISCFLUX(bi,bj,k,uVel,KappaRU,vrF,myThid)         CALL MOM_U_RVISCFLUX(bi,bj,k+1,uVel,kappaRU,vrF,myThid)
   
449  C     Combine fluxes  C     Combine fluxes
450        DO j=jMin,jMax         DO j=jMin,jMax
451         DO i=iMin,iMax          DO i=iMin,iMax
452          fVerU(i,j,kDown) = ArDudrFac*vrF(i,j)           fVerUkp(i,j) = ArDudrFac*vrF(i,j)
453            ENDDO
454         ENDDO         ENDDO
       ENDDO  
455    
456  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term  #ifdef ALLOW_AUTODIFF_TAMC
457        DO j=2-Oly,sNy+Oly-1  CADJ STORE fVerUkp(:,:) =
458         DO i=2-Olx,sNx+Olx-1  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
459          gU(i,j,k,bi,bj) = uDiss(i,j)  #endif
460    
461    C--   Tendency is minus divergence of the fluxes
462    C     vert.visc.flx is scaled by deepFac2F (deep-atmos) and rhoFac (anelastic)
463           DO j=jMin,jMax
464            DO i=iMin,iMax
465             guDiss(i,j) = guDiss(i,j)
466       &   -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
467       &   *recip_rAw(i,j,bi,bj)       &   *recip_rAw(i,j,bi,bj)
468       &  *(       &   *( fVerUkp(i,j) - fVerUkm(i,j) )*rkSign
469       &   +fVerU(i,j,kUp)*rkFac - fVerU(i,j,kDown)*rkFac       &   *recip_deepFac2C(k)*recip_rhoFacC(k)
470       &   )          ENDDO
      &  - phxFac*dPhiHydX(i,j)  
471         ENDDO         ENDDO
472        ENDDO        ENDIF
473    
474  C-- No-slip and drag BCs appear as body forces in cell abutting topography  C-- No-slip and drag BCs appear as body forces in cell abutting topography
475        IF (momViscosity.AND.no_slip_sides) THEN        IF ( no_slip_sides ) THEN
476  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
477         CALL MOM_U_SIDEDRAG(bi,bj,k,uFld,del2u,hFacZ,vF,myThid)         CALL MOM_U_SIDEDRAG( bi, bj, k,
478         I          uFld, del2u, h0FacZ,
479         I          viscAh_Z, viscA4_Z,
480         I          useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
481         O          vF,
482         I          myThid )
483         DO j=jMin,jMax         DO j=jMin,jMax
484          DO i=iMin,iMax          DO i=iMin,iMax
485           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+vF(i,j)           guDiss(i,j) = guDiss(i,j)+vF(i,j)
486          ENDDO          ENDDO
487         ENDDO         ENDDO
488        ENDIF        ENDIF
489    
490  C-    No-slip BCs impose a drag at bottom  C-    No-slip BCs impose a drag at bottom
491        IF (momViscosity.AND.bottomDragTerms) THEN        IF ( bottomDragTerms ) THEN
492         CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)         CALL MOM_U_BOTTOMDRAG( bi, bj, k,
493         I            uFld, vFld, KE, kappaRU,
494         O            vF,
495         I            myThid )
496         DO j=jMin,jMax         DO j=jMin,jMax
497          DO i=iMin,iMax          DO i=iMin,iMax
498           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+vF(i,j)           guDiss(i,j) = guDiss(i,j)+vF(i,j)
499          ENDDO          ENDDO
500         ENDDO         ENDDO
501        ENDIF        ENDIF
502    #ifdef ALLOW_SHELFICE
503  C--   Forcing term        IF ( useShelfIce ) THEN
504        IF (momForcing)         CALL SHELFICE_U_DRAG( bi, bj, k,
505       &  CALL EXTERNAL_FORCING_U(       I            uFld, vFld, KE, kappaRU,
506       I     iMin,iMax,jMin,jMax,bi,bj,k,       O            vF,
507       I     myCurrentTime,myThid)       I            myThid )
508           DO j=jMin,jMax
509  C--   Metric terms for curvilinear grid systems          DO i=iMin,iMax
510  c     IF (usingSphericalPolarMTerms) THEN           guDiss(i,j) = guDiss(i,j) + vF(i,j)
511  C      o Spherical polar grid metric terms          ENDDO
 c      CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)  
 c      DO j=jMin,jMax  
 c       DO i=iMin,iMax  
 c        gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)  
 c       ENDDO  
 c      ENDDO  
 c     ENDIF  
   
 C--   Set du/dt on boundaries to zero  
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)  
512         ENDDO         ENDDO
513        ENDDO        ENDIF
514    #endif /* ALLOW_SHELFICE */
515    
516  C---- Meridional momentum equation starts here  C---  Other dissipation terms in Meridional momentum equation
517    
518  C--   Vertical flux (fVer is at upper face of "v" cell)  C--   Vertical flux (fVer is at upper face of "v" cell)
   
519  C     Eddy component of vertical flux (interior component only) -> vrF  C     Eddy component of vertical flux (interior component only) -> vrF
520        IF (momViscosity.AND..NOT.implicitViscosity)        IF ( .NOT.implicitViscosity ) THEN
521       & CALL MOM_V_RVISCFLUX(bi,bj,k,vVel,KappaRV,vrf,myThid)         CALL MOM_V_RVISCFLUX(bi,bj,k+1,vVel,kappaRV,vrF,myThid)
   
522  C     Combine fluxes -> fVerV  C     Combine fluxes -> fVerV
523        DO j=jMin,jMax         DO j=jMin,jMax
524         DO i=iMin,iMax          DO i=iMin,iMax
525          fVerV(i,j,kDown) = ArDvdrFac*vrF(i,j)           fVerVkp(i,j) = ArDvdrFac*vrF(i,j)
526            ENDDO
527         ENDDO         ENDDO
528        ENDDO  #ifdef ALLOW_AUTODIFF_TAMC
529    CADJ STORE fVerVkp(:,:) =
530  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
531        DO j=jMin,jMax  #endif
532         DO i=iMin,iMax  C--   Tendency is minus divergence of the fluxes
533          gV(i,j,k,bi,bj) = vDiss(i,j)  C     vert.visc.flx is scaled by deepFac2F (deep-atmos) and rhoFac (anelastic)
534           DO j=jMin,jMax
535            DO i=iMin,iMax
536             gvDiss(i,j) = gvDiss(i,j)
537       &   -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
538       &    *recip_rAs(i,j,bi,bj)       &   *recip_rAs(i,j,bi,bj)
539       &  *(       &   *( fVerVkp(i,j) - fVerVkm(i,j) )*rkSign
540       &   +fVerV(i,j,kUp)*rkFac - fVerV(i,j,kDown)*rkFac       &   *recip_deepFac2C(k)*recip_rhoFacC(k)
541       &   )          ENDDO
      &  - phyFac*dPhiHydY(i,j)  
542         ENDDO         ENDDO
543        ENDDO        ENDIF
544    
545  C-- No-slip and drag BCs appear as body forces in cell abutting topography  C-- No-slip and drag BCs appear as body forces in cell abutting topography
546        IF (momViscosity.AND.no_slip_sides) THEN        IF ( no_slip_sides ) THEN
547  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
548         CALL MOM_V_SIDEDRAG(bi,bj,k,vFld,del2v,hFacZ,vF,myThid)         CALL MOM_V_SIDEDRAG( bi, bj, k,
549         I          vFld, del2v, h0FacZ,
550         I          viscAh_Z, viscA4_Z,
551         I          useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
552         O          vF,
553         I          myThid )
554         DO j=jMin,jMax         DO j=jMin,jMax
555          DO i=iMin,iMax          DO i=iMin,iMax
556           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vF(i,j)           gvDiss(i,j) = gvDiss(i,j)+vF(i,j)
557          ENDDO          ENDDO
558         ENDDO         ENDDO
559        ENDIF        ENDIF
560    
561  C-    No-slip BCs impose a drag at bottom  C-    No-slip BCs impose a drag at bottom
562        IF (momViscosity.AND.bottomDragTerms) THEN        IF ( bottomDragTerms ) THEN
563         CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)         CALL MOM_V_BOTTOMDRAG( bi, bj, k,
564         I            uFld, vFld, KE, kappaRV,
565         O            vF,
566         I            myThid )
567         DO j=jMin,jMax         DO j=jMin,jMax
568          DO i=iMin,iMax          DO i=iMin,iMax
569           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vF(i,j)           gvDiss(i,j) = gvDiss(i,j)+vF(i,j)
570          ENDDO          ENDDO
571         ENDDO         ENDDO
572        ENDIF        ENDIF
573    #ifdef ALLOW_SHELFICE
574          IF ( useShelfIce ) THEN
575           CALL SHELFICE_V_DRAG( bi, bj, k,
576         I            uFld, vFld, KE, kappaRV,
577         O            vF,
578         I            myThid )
579           DO j=jMin,jMax
580            DO i=iMin,iMax
581             gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
582            ENDDO
583           ENDDO
584          ENDIF
585    #endif /* ALLOW_SHELFICE */
586    
587  C--   Forcing term  C--   if (momViscosity) end of block.
588        IF (momForcing)        ENDIF
      & CALL EXTERNAL_FORCING_V(  
      I     iMin,iMax,jMin,jMax,bi,bj,k,  
      I     myCurrentTime,myThid)  
589    
590  C--   Metric terms for curvilinear grid systems  C-    Return to standard hfacZ (min-4) and mask vort3 accordingly:
591  c     IF (usingSphericalPolarMTerms) THEN  c     CALL MOM_VI_MASK_VORT3(bi,bj,k,hFacZ,r_hFacZ,vort3,myThid)
 C      o Spherical polar grid metric terms  
 c      CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)  
 c      DO j=jMin,jMax  
 c       DO i=iMin,iMax  
 c        gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)  
 c       ENDDO  
 c      ENDDO  
 c     ENDIF  
592    
593  C--   Set dv/dt on boundaries to zero  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
594        DO j=jMin,jMax  
595         DO i=iMin,iMax  C---  Prepare for Advection & Coriolis terms:
596          gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)  C-    calculate absolute vorticity
597         ENDDO        IF (useAbsVorticity)
598        ENDDO       &  CALL MOM_CALC_ABSVORT3(bi,bj,k,vort3,omega3,myThid)
599    
600    #ifdef ALLOW_AUTODIFF_TAMC
601    CADJ STORE omega3(:,:) =
602    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
603    #endif
604    
605  C--   Horizontal Coriolis terms  C--   Horizontal Coriolis terms
606        CALL MOM_VI_CORIOLIS(bi,bj,K,uFld,vFld,omega3,r_hFacZ,  c     IF (useCoriolis .AND. .NOT.useCDscheme
607       &                     uCf,vCf,myThid)  c    &    .AND. .NOT. useAbsVorticity) THEN
608        DO j=jMin,jMax  C- jmc: change it to keep the Coriolis terms when useAbsVorticity=T & momAdvection=F
609         DO i=iMin,iMax        IF ( useCoriolis .AND.
610          gU(i,j,k,bi,bj) = (gU(i,j,k,bi,bj)+uCf(i,j))       &     .NOT.( useCDscheme .OR. useAbsVorticity.AND.momAdvection )
611       &                    *_maskW(i,j,k,bi,bj)       &   ) THEN
612          gV(i,j,k,bi,bj) = (gV(i,j,k,bi,bj)+vCf(i,j))         IF (useAbsVorticity) THEN
613       &                    *_maskS(i,j,k,bi,bj)          CALL MOM_VI_U_CORIOLIS(bi,bj,k,vFld,omega3,hFacZ,r_hFacZ,
614         ENDDO       &                         uCf,myThid)
615        ENDDO          CALL MOM_VI_V_CORIOLIS(bi,bj,k,uFld,omega3,hFacZ,r_hFacZ,
616  c     CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,r_hFacZ,uCf,myThid)       &                         vCf,myThid)
617        CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,vort3,r_hFacZ,uCf,myThid)         ELSE
618  c     CALL MOM_VI_U_CORIOLIS_C4(bi,bj,K,vFld,vort3,r_hFacZ,uCf,myThid)          CALL MOM_VI_CORIOLIS(bi,bj,k,uFld,vFld,hFacZ,r_hFacZ,
619        DO j=jMin,jMax       &                       uCf,vCf,myThid)
620         DO i=iMin,iMax         ENDIF
621          gU(i,j,k,bi,bj) = (gU(i,j,k,bi,bj)+uCf(i,j))         DO j=jMin,jMax
622       &                    *_maskW(i,j,k,bi,bj)          DO i=iMin,iMax
623             gU(i,j,k,bi,bj) = uCf(i,j)
624             gV(i,j,k,bi,bj) = vCf(i,j)
625            ENDDO
626         ENDDO         ENDDO
627        ENDDO         IF ( writeDiag ) THEN
628  c     CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,r_hFacZ,vCf,myThid)           IF (snapshot_mdsio) THEN
629        CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,vort3,r_hFacZ,vCf,myThid)             CALL WRITE_LOCAL_RL('fV','I10',1,uCf,bi,bj,k,myIter,myThid)
630  c     CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,vort3,r_hFacZ,vCf,myThid)             CALL WRITE_LOCAL_RL('fU','I10',1,vCf,bi,bj,k,myIter,myThid)
631        DO j=jMin,jMax           ENDIF
632         DO i=iMin,iMax  #ifdef ALLOW_MNC
633          gV(i,j,k,bi,bj) = (gV(i,j,k,bi,bj)+vCf(i,j))           IF (useMNC .AND. snapshot_mnc) THEN
634       &                    *_maskS(i,j,k,bi,bj)             CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'fV', uCf,
635         &          offsets, myThid)
636               CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'fU', vCf,
637         &          offsets, myThid)
638             ENDIF
639    #endif /*  ALLOW_MNC  */
640           ENDIF
641    #ifdef ALLOW_DIAGNOSTICS
642           IF ( useDiagnostics ) THEN
643             CALL DIAGNOSTICS_FILL(uCf,'Um_Cori ',k,1,2,bi,bj,myThid)
644             CALL DIAGNOSTICS_FILL(vCf,'Vm_Cori ',k,1,2,bi,bj,myThid)
645           ENDIF
646    #endif /* ALLOW_DIAGNOSTICS */
647          ELSE
648           DO j=jMin,jMax
649            DO i=iMin,iMax
650             gU(i,j,k,bi,bj) = 0. _d 0
651             gV(i,j,k,bi,bj) = 0. _d 0
652            ENDDO
653         ENDDO         ENDDO
654        ENDDO        ENDIF
655    
656    #ifdef ALLOW_AUTODIFF_TAMC
657    CADJ STORE ucf(:,:) =
658    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
659    CADJ STORE vcf(:,:) =
660    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
661    #endif
662    
663        IF (momAdvection) THEN        IF (momAdvection) THEN
664  C--   Vertical shear terms (Coriolis)  C--   Horizontal advection of relative (or absolute) vorticity
665        CALL MOM_VI_U_VERTSHEAR(bi,bj,K,uVel,wVel,uCf,myThid)         IF ( (highOrderVorticity.OR.upwindVorticity)
666        DO j=jMin,jMax       &     .AND.useAbsVorticity ) THEN
667         DO i=iMin,iMax          CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,selectVortScheme,
668          gU(i,j,k,bi,bj) = (gU(i,j,k,bi,bj)+uCf(i,j))       &                         highOrderVorticity,upwindVorticity,
669       &                    *_maskW(i,j,k,bi,bj)       &                         vFld,omega3,r_hFacZ,
670         &                         uCf,myThid)
671           ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
672            CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,selectVortScheme,
673         &                         highOrderVorticity, upwindVorticity,
674         &                         vFld,vort3, r_hFacZ,
675         &                         uCf,myThid)
676           ELSEIF ( useAbsVorticity ) THEN
677            CALL MOM_VI_U_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
678         &                         vFld,omega3,hFacZ,r_hFacZ,
679         &                         uCf,myThid)
680           ELSE
681            CALL MOM_VI_U_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
682         &                         vFld,vort3, hFacZ,r_hFacZ,
683         &                         uCf,myThid)
684           ENDIF
685           DO j=jMin,jMax
686            DO i=iMin,iMax
687             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
688            ENDDO
689         ENDDO         ENDDO
690        ENDDO         IF ( (highOrderVorticity.OR.upwindVorticity)
691        CALL MOM_VI_V_VERTSHEAR(bi,bj,K,vVel,wVel,vCf,myThid)       &     .AND.useAbsVorticity ) THEN
692        DO j=jMin,jMax          CALL MOM_VI_V_CORIOLIS_C4(bi,bj,k,selectVortScheme,
693         DO i=iMin,iMax       &                         highOrderVorticity, upwindVorticity,
694          gV(i,j,k,bi,bj) = (gV(i,j,k,bi,bj)+vCf(i,j))       &                         uFld,omega3,r_hFacZ,
695       &                    *_maskS(i,j,k,bi,bj)       &                         vCf,myThid)
696           ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
697            CALL MOM_VI_V_CORIOLIS_C4(bi,bj,k,selectVortScheme,
698         &                         highOrderVorticity, upwindVorticity,
699         &                         uFld,vort3, r_hFacZ,
700         &                         vCf,myThid)
701           ELSEIF ( useAbsVorticity ) THEN
702            CALL MOM_VI_V_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
703         &                         uFld,omega3,hFacZ,r_hFacZ,
704         &                         vCf,myThid)
705           ELSE
706            CALL MOM_VI_V_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
707         &                         uFld,vort3, hFacZ,r_hFacZ,
708         &                         vCf,myThid)
709           ENDIF
710           DO j=jMin,jMax
711            DO i=iMin,iMax
712             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
713            ENDDO
714         ENDDO         ENDDO
715        ENDDO  
716    #ifdef ALLOW_AUTODIFF_TAMC
717    CADJ STORE ucf(:,:) =
718    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
719    CADJ STORE vcf(:,:) =
720    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
721    #endif
722    
723           IF ( writeDiag ) THEN
724             IF (snapshot_mdsio) THEN
725               CALL WRITE_LOCAL_RL('zV','I10',1,uCf,bi,bj,k,myIter,myThid)
726               CALL WRITE_LOCAL_RL('zU','I10',1,vCf,bi,bj,k,myIter,myThid)
727             ENDIF
728    #ifdef ALLOW_MNC
729             IF (useMNC .AND. snapshot_mnc) THEN
730               CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'zV', uCf,
731         &          offsets, myThid)
732               CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'zU', vCf,
733         &          offsets, myThid)
734             ENDIF
735    #endif /*  ALLOW_MNC  */
736           ENDIF
737    
738    #ifdef ALLOW_TIMEAVE
739           IF (taveFreq.GT.0.) THEN
740             CALL TIMEAVE_CUMUL_1K1T(uZetatave,vCf,deltaTClock,
741         &                           Nr, k, bi, bj, myThid)
742             CALL TIMEAVE_CUMUL_1K1T(vZetatave,uCf,deltaTClock,
743         &                           Nr, k, bi, bj, myThid)
744           ENDIF
745    #endif /* ALLOW_TIMEAVE */
746    #ifdef ALLOW_DIAGNOSTICS
747           IF ( useDiagnostics ) THEN
748             CALL DIAGNOSTICS_FILL(uCf,'Um_AdvZ3',k,1,2,bi,bj,myThid)
749             CALL DIAGNOSTICS_FILL(vCf,'Vm_AdvZ3',k,1,2,bi,bj,myThid)
750           ENDIF
751    #endif /* ALLOW_DIAGNOSTICS */
752    
753    C--   Vertical shear terms (-w*du/dr & -w*dv/dr)
754           IF ( .NOT. momImplVertAdv ) THEN
755            CALL MOM_VI_U_VERTSHEAR(bi,bj,k,uVel,wVel,uCf,myThid)
756            DO j=jMin,jMax
757             DO i=iMin,iMax
758              gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
759             ENDDO
760            ENDDO
761            CALL MOM_VI_V_VERTSHEAR(bi,bj,k,vVel,wVel,vCf,myThid)
762            DO j=jMin,jMax
763             DO i=iMin,iMax
764              gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
765             ENDDO
766            ENDDO
767    #ifdef ALLOW_DIAGNOSTICS
768            IF ( useDiagnostics ) THEN
769             CALL DIAGNOSTICS_FILL(uCf,'Um_AdvRe',k,1,2,bi,bj,myThid)
770             CALL DIAGNOSTICS_FILL(vCf,'Vm_AdvRe',k,1,2,bi,bj,myThid)
771            ENDIF
772    #endif /* ALLOW_DIAGNOSTICS */
773           ENDIF
774    
775  C--   Bernoulli term  C--   Bernoulli term
776        CALL MOM_VI_U_GRAD_KE(bi,bj,K,KE,uCf,myThid)         CALL MOM_VI_U_GRAD_KE(bi,bj,k,KE,uCf,myThid)
777        DO j=jMin,jMax         DO j=jMin,jMax
778         DO i=iMin,iMax          DO i=iMin,iMax
779          gU(i,j,k,bi,bj) = (gU(i,j,k,bi,bj)+uCf(i,j))           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
780       &                    *_maskW(i,j,k,bi,bj)          ENDDO
781         ENDDO         ENDDO
782        ENDDO         CALL MOM_VI_V_GRAD_KE(bi,bj,k,KE,vCf,myThid)
783        CALL MOM_VI_V_GRAD_KE(bi,bj,K,KE,vCf,myThid)         DO j=jMin,jMax
784            DO i=iMin,iMax
785             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
786            ENDDO
787           ENDDO
788           IF ( writeDiag ) THEN
789             IF (snapshot_mdsio) THEN
790               CALL WRITE_LOCAL_RL('KEx','I10',1,uCf,bi,bj,k,myIter,myThid)
791               CALL WRITE_LOCAL_RL('KEy','I10',1,vCf,bi,bj,k,myIter,myThid)
792             ENDIF
793    #ifdef ALLOW_MNC
794             IF (useMNC .AND. snapshot_mnc) THEN
795               CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'KEx', uCf,
796         &          offsets, myThid)
797               CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'KEy', vCf,
798         &          offsets, myThid)
799             ENDIF
800    #endif /*  ALLOW_MNC  */
801           ENDIF
802    
803    C--   end if momAdvection
804          ENDIF
805    
806    C--   3.D Coriolis term (horizontal momentum, Eastward component: -fprime*w)
807          IF ( use3dCoriolis ) THEN
808            CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,uCf,myThid)
809            DO j=jMin,jMax
810             DO i=iMin,iMax
811              gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
812             ENDDO
813            ENDDO
814           IF ( usingCurvilinearGrid ) THEN
815    C-     presently, non zero angleSinC array only supported with Curvilinear-Grid
816            CALL MOM_V_CORIOLIS_NH(bi,bj,k,wVel,vCf,myThid)
817            DO j=jMin,jMax
818             DO i=iMin,iMax
819              gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
820             ENDDO
821            ENDDO
822           ENDIF
823          ENDIF
824    
825    C--   Non-Hydrostatic (spherical) metric terms
826          IF ( useNHMTerms ) THEN
827           CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,uCf,myThid)
828           DO j=jMin,jMax
829            DO i=iMin,iMax
830             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
831            ENDDO
832           ENDDO
833           CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,vCf,myThid)
834           DO j=jMin,jMax
835            DO i=iMin,iMax
836             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
837            ENDDO
838           ENDDO
839          ENDIF
840    
841    C--   Set du/dt & dv/dt on boundaries to zero
842        DO j=jMin,jMax        DO j=jMin,jMax
843         DO i=iMin,iMax         DO i=iMin,iMax
844          gV(i,j,k,bi,bj) = (gV(i,j,k,bi,bj)+vCf(i,j))          gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
845       &                    *_maskS(i,j,k,bi,bj)          gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
846         ENDDO         ENDDO
847        ENDDO        ENDDO
848    
849    #ifdef ALLOW_DEBUG
850          IF ( debugLevel .GE. debLevC
851         &   .AND. k.EQ.4 .AND. myIter.EQ.nIter0
852         &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
853         &   .AND. useCubedSphereExchange ) THEN
854            CALL DEBUG_CS_CORNER_UV( ' uDiss,vDiss from MOM_VECINV',
855         &             guDiss,gvDiss, k, standardMessageUnit,bi,bj,myThid )
856        ENDIF        ENDIF
857    #endif /* ALLOW_DEBUG */
858    
859        IF (        IF ( writeDiag ) THEN
860       &  DIFFERENT_MULTIPLE(diagFreq,myCurrentTime,          IF (useBiharmonicVisc) THEN
861       &                     myCurrentTime-deltaTClock)           CALL WRITE_LOCAL_RL( 'del2u', 'I10', 1, del2u,
862       & ) THEN       &                        bi,bj,k, myIter, myThid )
863         CALL WRITE_LOCAL_RL('Ds','I10',1,strain,bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL( 'del2v', 'I10', 1, del2v,
864         CALL WRITE_LOCAL_RL('Dt','I10',1,tension,bi,bj,k,myIter,myThid)       &                        bi,bj,k, myIter, myThid )
865         CALL WRITE_LOCAL_RL('fV','I10',1,uCf,bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL( 'dStar', 'I10', 1, dStar,
866         CALL WRITE_LOCAL_RL('fU','I10',1,vCf,bi,bj,k,myIter,myThid)       &                        bi,bj,k, myIter, myThid )
867         CALL WRITE_LOCAL_RL('Du','I10',1,uDiss,bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL( 'zStar', 'I10', 1, zStar,
868         CALL WRITE_LOCAL_RL('Dv','I10',1,vDiss,bi,bj,k,myIter,myThid)       &                        bi,bj,k, myIter, myThid )
869         CALL WRITE_LOCAL_RL('Z3','I10',1,vort3,bi,bj,k,myIter,myThid)          ENDIF
870         CALL WRITE_LOCAL_RL('W3','I10',1,omega3,bi,bj,k,myIter,myThid)          IF (snapshot_mdsio) THEN
871         CALL WRITE_LOCAL_RL('KE','I10',1,KE,bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL('W3','I10',1,omega3, bi,bj,k,myIter,myThid)
872         CALL WRITE_LOCAL_RL('D','I10',1,hdiv,bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL('Z3','I10',1,vort3BC,bi,bj,k,myIter,myThid)
873             CALL WRITE_LOCAL_RL('KE','I10',1,KE,     bi,bj,k,myIter,myThid)
874             CALL WRITE_LOCAL_RL('D', 'I10',1,hDiv,   bi,bj,k,myIter,myThid)
875             CALL WRITE_LOCAL_RL('Dt','I10',1,tension,bi,bj,k,myIter,myThid)
876             CALL WRITE_LOCAL_RL( 'Ds', 'I10', 1, strainBC,
877         &                        bi,bj,k, myIter, myThid )
878             CALL WRITE_LOCAL_RL('Du','I10',1,guDiss, bi,bj,k,myIter,myThid)
879             CALL WRITE_LOCAL_RL('Dv','I10',1,gvDiss, bi,bj,k,myIter,myThid)
880            ENDIF
881    #ifdef ALLOW_MNC
882            IF (useMNC .AND. snapshot_mnc) THEN
883              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'W3',omega3,
884         &          offsets, myThid)
885              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Z3',vort3BC,
886         &          offsets, myThid)
887              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'KE',KE,
888         &          offsets, myThid)
889              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'D', hDiv,
890         &          offsets, myThid)
891              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dt',tension,
892         &          offsets, myThid)
893              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Ds',strainBC,
894         &          offsets, myThid)
895              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Du',guDiss,
896         &          offsets, myThid)
897              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dv',gvDiss,
898         &          offsets, myThid)
899            ENDIF
900    #endif /*  ALLOW_MNC  */
901        ENDIF        ENDIF
902    
903    #ifdef ALLOW_DIAGNOSTICS
904          IF ( useDiagnostics ) THEN
905            CALL DIAGNOSTICS_FILL(vort3BC,'momVort3',k,1,2,bi,bj,myThid)
906            CALL DIAGNOSTICS_FILL(KE,     'momKE   ',k,1,2,bi,bj,myThid)
907           IF (momViscosity) THEN
908            CALL DIAGNOSTICS_FILL(hDiv,   'momHDiv ',k,1,2,bi,bj,myThid)
909           ENDIF
910           IF ( useVariableVisc .OR. useStrainTensionVisc ) THEN
911            CALL DIAGNOSTICS_FILL(tension, 'Tension ',k,1,2,bi,bj,myThid)
912            CALL DIAGNOSTICS_FILL(strainBC,'Strain  ',k,1,2,bi,bj,myThid)
913           ENDIF
914            CALL DIAGNOSTICS_FILL(gU(1-OLx,1-OLy,k,bi,bj),
915         &                                'Um_Advec',k,1,2,bi,bj,myThid)
916            CALL DIAGNOSTICS_FILL(gV(1-OLx,1-OLy,k,bi,bj),
917         &                                'Vm_Advec',k,1,2,bi,bj,myThid)
918          ENDIF
919    #endif /* ALLOW_DIAGNOSTICS */
920    
921    #endif /* ALLOW_MOM_VECINV */
922    
923        RETURN        RETURN
924        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.80
  ViewVC Help
Powered by ViewVC 1.1.22