/[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.15 by jmc, Wed Feb 25 00:56:47 2004 UTC revision 1.81 by jmc, Sat Apr 29 16:11:38 2017 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "PACKAGES_CONFIG.h"  #include "MOM_VECINV_OPTIONS.h"
5  #include "CPP_OPTIONS.h"  #ifdef ALLOW_AUTODIFF
6    # include "AUTODIFF_OPTIONS.h"
7    #endif
8    #ifdef ALLOW_MOM_COMMON
9    # include "MOM_COMMON_OPTIONS.h"
10    #endif
11    
12        SUBROUTINE MOM_VECINV(        SUBROUTINE MOM_VECINV(
13       I        bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,       I        bi,bj,k,iMin,iMax,jMin,jMax,
14       I        dPhiHydX,dPhiHydY,KappaRU,KappaRV,       I        kappaRU, kappaRV,
15       U        fVerU, fVerV,       I        fVerUkm, fVerVkm,
16       I        myTime, myIter, myThid)       O        fVerUkp, fVerVkp,
17  C     /==========================================================\       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 23  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  #ifdef ALLOW_TIMEAVE
47  #include "TIMEAVE_STATV.h"  # 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  #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          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        _RL     myTime
84        INTEGER myIter        INTEGER myIter
85        INTEGER myThid        INTEGER myThid
       INTEGER bi,bj,iMin,iMax,jMin,jMax  
86    
87  #ifdef ALLOW_MOM_VECINV  #ifdef ALLOW_MOM_VECINV
88    
# Line 64  C     == Functions == Line 91  C     == Functions ==
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 uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL strain  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111        _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL strainBC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112        _RL dStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KE      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113        _RL zStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL omega3  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114        _RL uDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vort3   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115        _RL vDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vort3BC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116  C     I,J,K - Loop counters        _RL hDiv    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117        INTEGER i,j,k        _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118  C     rVelMaskOverride - Factor for imposing special surface boundary conditions        _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119  C                        ( set according to free-surface condition ).        _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120  C     hFacROpen        - Lopped cell factos used tohold fraction of open        _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121  C     hFacRClosed        and closed cell wall.  C     i,j    :: Loop counters
122        _RL  rVelMaskOverride        INTEGER i,j
123  C     xxxFac - On-off tracer parameters used for switching terms off.  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  
       _RL wVelBottomOverride  
127        LOGICAL bottomDragTerms        LOGICAL bottomDragTerms
128        LOGICAL writeDiag        LOGICAL writeDiag
       _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL omega3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
   
129  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
130          INTEGER imomkey
131    #endif
132    
133    #ifdef ALLOW_MNC
134          INTEGER offsets(9)
135          CHARACTER*(1) pf
136    #endif
137    
138    #ifdef ALLOW_AUTODIFF
139  C--   only the kDown part of fverU/V is set in this subroutine  C--   only the kDown part of fverU/V is set in this subroutine
140  C--   the kUp is still required  C--   the kUp is still required
141  C--   In the case of mom_fluxform Kup is set as well  C--   In the case of mom_fluxform kUp is set as well
142  C--   (at least in part)  C--   (at least in part)
143        fVerU(1,1,kUp) = fVerU(1,1,kUp)        fVerUkm(1,1) = fVerUkm(1,1)
144        fVerV(1,1,kUp) = fVerV(1,1,kUp)        fVerVkm(1,1) = fVerVkm(1,1)
145  #endif  #endif
146    
147        rVelMaskOverride=1.  #ifdef ALLOW_AUTODIFF_TAMC
148        IF ( k .EQ. 1 ) rVelMaskOverride=freeSurfFac            act0 = k - 1
149        wVelBottomOverride=1.            max0 = Nr
150        IF (k.EQ.Nr) wVelBottomOverride=0.            act1 = bi - myBxLo(myThid)
151        writeDiag = DIFFERENT_MULTIPLE(diagFreq, myTime,            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
152       &                                         myTime-deltaTClock)            act2 = bj - myByLo(myThid)
153              max2 = myByHi(myThid) - myByLo(myThid) + 1
154  C     Initialise intermediate terms            act3 = myThid - 1
155        DO J=1-OLy,sNy+OLy            max3 = nTx*nTy
156         DO I=1-OLx,sNx+OLx            act4 = ikey_dynamics - 1
157          aF(i,j)   = 0.            imomkey = (act0 + 1)
158          vF(i,j)   = 0.       &                    + act1*max0
159          vrF(i,j)  = 0.       &                    + 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  #ifdef ALLOW_AUTODIFF_TAMC  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          strain(i,j)  = 0. _d 0
210            strainBC(i,j)= 0. _d 0
211          tension(i,j) = 0. _d 0          tension(i,j) = 0. _d 0
212    #ifdef ALLOW_AUTODIFF
213            hFacZ(i,j)   = 0. _d 0
214  #endif  #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 229  C     Make local copies of horizontal fl Line 249  C     Make local copies of horizontal fl
249         ENDDO         ENDDO
250        ENDDO        ENDDO
251    
252    #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  C note (jmc) : Dissipation and Vort3 advection do not necesary
268  C              use the same maskZ (and hFacZ)  => needs 2 call(s)  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)  c     CALL MOM_VI_HFACZ_DISS(bi,bj,k,hFacZ,r_hFacZ,myThid)
270    
271        CALL MOM_VI_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid)        CALL MOM_CALC_KE(bi,bj,k,selectKEscheme,uFld,vFld,KE,myThid)
272    
273        CALL MOM_VI_CALC_HDIV(bi,bj,k,uFld,vFld,hDiv,myThid)        CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
274    
275        CALL MOM_VI_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)  #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     CALL MOM_VI_CALC_ABSVORT3(bi,bj,k,vort3,omega3,myThid)  C-    mask vort3 and account for no-slip / free-slip BC in vort3BC:
285          DO j=1-OLy,sNy+OLy
286           DO i=1-OLx,sNx+OLx
287             vort3BC(i,j) = vort3(i,j)
288             IF ( hFacZ(i,j).EQ.zeroRS ) THEN
289               vort3BC(i,j) = sideMaskFac*vort3BC(i,j)
290               vort3(i,j)   = 0.
291             ENDIF
292           ENDDO
293          ENDDO
294    
295    #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        IF (momViscosity) THEN
303    C--    For viscous term, compute horizontal divergence, tension & strain
304    C      and mask relative vorticity (free-slip case):
305    
306           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    #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    #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    
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. .OR. viscA4Grid.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. .OR.           CALL MOM_CALC_HDIV(bi,bj,k,2,del2u,del2v,dStar,myThid)
407       &      viscAhGrid.NE.0. .OR. viscA4Grid.NE.0. ) THEN           CALL MOM_CALC_RELVORT3(bi,bj,k,
408           CALL MOM_VI_HDISSIP(bi,bj,k,hDiv,vort3,hFacZ,dStar,zStar,       &                          del2u,del2v,hFacZ,zStar,myThid)
      O                       uDiss,vDiss,  
      &                       myThid)  
        ENDIF  
 C      or in terms of tension and strain  
        IF (viscAstrain.NE.0. .OR. viscAtension.NE.0.) THEN  
          CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,  
      O                         tension,  
      I                         myThid)  
          CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,  
      O                        strain,  
      I                        myThid)  
          CALL MOM_HDISSIP(bi,bj,k,  
      I                    tension,strain,hFacZ,viscAtension,viscAstrain,  
      O                    uDiss,vDiss,  
      I                    myThid)  
409         ENDIF         ENDIF
       ENDIF  
410    
411  C-    Return to standard hfacZ (min-4) and mask vort3 accordingly:  #ifdef ALLOW_AUTODIFF_TAMC
412  c     CALL MOM_VI_MASK_VORT3(bi,bj,k,hFacZ,r_hFacZ,vort3,myThid)  CADJ STORE del2u(:,:) =
413    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
414    CADJ STORE del2v(:,:) =
415    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
416    CADJ STORE dStar(:,:) =
417    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
418    CADJ STORE zStar(:,:) =
419    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
420    #endif
421    
422  C---- Zonal momentum equation starts here  C---   Calculate dissipation terms for U and V equations
423    
424  C--   Vertical flux (fVer is at upper face of "u" cell)  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
442    
443  C     Eddy component of vertical flux (interior component only) -> vrF  C---  Other dissipation terms in Zonal momentum equation
       IF (momViscosity.AND..NOT.implicitViscosity)  
      & CALL MOM_U_RVISCFLUX(bi,bj,k,uVel,KappaRU,vrF,myThid)  
444    
445    C--   Vertical flux (fVer is at upper face of "u" cell)
446    C     Eddy component of vertical flux (interior component only) -> vrF
447          IF ( .NOT.implicitViscosity ) THEN
448           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          IF ( useShelfIce ) THEN
504           CALL SHELFICE_U_DRAG( bi, bj, k,
505         I            uFld, vFld, KE, kappaRU,
506         O            vF,
507         I            myThid )
508           DO j=jMin,jMax
509            DO i=iMin,iMax
510             guDiss(i,j) = guDiss(i,j) + vF(i,j)
511            ENDDO
512           ENDDO
513          ENDIF
514    #endif /* ALLOW_SHELFICE */
515    
516  C--   Metric terms for curvilinear grid systems  C---  Other dissipation terms in Meridional momentum equation
 c     IF (usingSphericalPolarMTerms) THEN  
 C      o Spherical polar grid metric terms  
 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---- Meridional momentum equation starts here  
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--   Metric terms for curvilinear grid systems  C--   if (momViscosity) end of block.
588  c     IF (usingSphericalPolarMTerms) THEN        ENDIF
589  C      o Spherical polar grid metric terms  
590  c      CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)  C-    Return to standard hfacZ (min-4) and mask vort3 accordingly:
591  c      DO j=jMin,jMax  c     CALL MOM_VI_MASK_VORT3(bi,bj,k,hFacZ,r_hFacZ,vort3,myThid)
592  c       DO i=iMin,iMax  
593  c        gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
594  c       ENDDO  
595  c      ENDDO  C---  Prepare for Advection & Coriolis terms:
596  c     ENDIF  C-    calculate absolute vorticity
597          IF (useAbsVorticity)
598         &  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        IF (useCoriolis .AND. .NOT.useCDscheme) THEN  c     IF (useCoriolis .AND. .NOT.useCDscheme
607         CALL MOM_VI_CORIOLIS(bi,bj,k,uFld,vFld,omega3,hFacZ,r_hFacZ,  c    &    .AND. .NOT. useAbsVorticity) THEN
608       &                      uCf,vCf,myThid)  C- jmc: change it to keep the Coriolis terms when useAbsVorticity=T & momAdvection=F
609          IF ( useCoriolis .AND.
610         &     .NOT.( useCDscheme .OR. useAbsVorticity.AND.momAdvection )
611         &   ) THEN
612           IF (useAbsVorticity) THEN
613            CALL MOM_VI_U_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
614         &                         vFld,omega3,hFacZ,r_hFacZ,
615         &                         uCf,myThid)
616            CALL MOM_VI_V_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
617         &                         uFld,omega3,hFacZ,r_hFacZ,
618         &                         vCf,myThid)
619           ELSE
620            CALL MOM_VI_CORIOLIS(bi,bj,k,uFld,vFld,hFacZ,r_hFacZ,
621         &                       uCf,vCf,myThid)
622           ENDIF
623         DO j=jMin,jMax         DO j=jMin,jMax
624          DO i=iMin,iMax          DO i=iMin,iMax
625           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)           gU(i,j,k,bi,bj) = uCf(i,j)
626           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)           gV(i,j,k,bi,bj) = vCf(i,j)
627          ENDDO          ENDDO
628         ENDDO         ENDDO
629         IF ( writeDiag ) THEN         IF ( writeDiag ) THEN
630          CALL WRITE_LOCAL_RL('fV','I10',1,uCf,bi,bj,k,myIter,myThid)           IF (snapshot_mdsio) THEN
631          CALL WRITE_LOCAL_RL('fU','I10',1,vCf,bi,bj,k,myIter,myThid)             CALL WRITE_LOCAL_RL('fV','I10',1,uCf,bi,bj,k,myIter,myThid)
632               CALL WRITE_LOCAL_RL('fU','I10',1,vCf,bi,bj,k,myIter,myThid)
633             ENDIF
634    #ifdef ALLOW_MNC
635             IF (useMNC .AND. snapshot_mnc) THEN
636               CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'fV', uCf,
637         &          offsets, myThid)
638               CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'fU', vCf,
639         &          offsets, myThid)
640             ENDIF
641    #endif /*  ALLOW_MNC  */
642         ENDIF         ENDIF
643    #ifdef ALLOW_DIAGNOSTICS
644           IF ( useDiagnostics ) THEN
645             CALL DIAGNOSTICS_FILL(uCf,'Um_Cori ',k,1,2,bi,bj,myThid)
646             CALL DIAGNOSTICS_FILL(vCf,'Vm_Cori ',k,1,2,bi,bj,myThid)
647           ENDIF
648    #endif /* ALLOW_DIAGNOSTICS */
649          ELSE
650           DO j=jMin,jMax
651            DO i=iMin,iMax
652             gU(i,j,k,bi,bj) = 0. _d 0
653             gV(i,j,k,bi,bj) = 0. _d 0
654            ENDDO
655           ENDDO
656        ENDIF        ENDIF
657    
658    #ifdef ALLOW_AUTODIFF_TAMC
659    CADJ STORE ucf(:,:) =
660    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
661    CADJ STORE vcf(:,:) =
662    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
663    #endif
664    
665        IF (momAdvection) THEN        IF (momAdvection) THEN
666  C--   Horizontal advection of relative vorticity  C--   Horizontal advection of relative (or absolute) vorticity
667  c      CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,r_hFacZ,uCf,myThid)         IF ( (highOrderVorticity.OR.upwindVorticity)
668         CALL MOM_VI_U_CORIOLIS(bi,bj,k,vFld,vort3,hFacZ,r_hFacZ,       &     .AND.useAbsVorticity ) THEN
669       &                        uCf,myThid)          CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,selectVortScheme,
670  c      CALL MOM_VI_U_CORIOLIS_C4(bi,bj,K,vFld,vort3,r_hFacZ,uCf,myThid)       &                         highOrderVorticity,upwindVorticity,
671         &                         vFld,omega3,r_hFacZ,
672         &                         uCf,myThid)
673           ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
674            CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,selectVortScheme,
675         &                         highOrderVorticity, upwindVorticity,
676         &                         vFld,vort3, r_hFacZ,
677         &                         uCf,myThid)
678           ELSEIF ( useAbsVorticity ) THEN
679            CALL MOM_VI_U_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
680         &                         vFld,omega3,hFacZ,r_hFacZ,
681         &                         uCf,myThid)
682           ELSE
683            CALL MOM_VI_U_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
684         &                         vFld,vort3, hFacZ,r_hFacZ,
685         &                         uCf,myThid)
686           ENDIF
687         DO j=jMin,jMax         DO j=jMin,jMax
688          DO i=iMin,iMax          DO i=iMin,iMax
689           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)
690          ENDDO          ENDDO
691         ENDDO         ENDDO
692  c      CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,r_hFacZ,vCf,myThid)         IF ( (highOrderVorticity.OR.upwindVorticity)
693         CALL MOM_VI_V_CORIOLIS(bi,bj,k,uFld,vort3,hFacZ,r_hFacZ,       &     .AND.useAbsVorticity ) THEN
694       &                        vCf,myThid)          CALL MOM_VI_V_CORIOLIS_C4(bi,bj,k,selectVortScheme,
695  c      CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,vort3,r_hFacZ,vCf,myThid)       &                         highOrderVorticity, upwindVorticity,
696         &                         uFld,omega3,r_hFacZ,
697         &                         vCf,myThid)
698           ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
699            CALL MOM_VI_V_CORIOLIS_C4(bi,bj,k,selectVortScheme,
700         &                         highOrderVorticity, upwindVorticity,
701         &                         uFld,vort3, r_hFacZ,
702         &                         vCf,myThid)
703           ELSEIF ( useAbsVorticity ) THEN
704            CALL MOM_VI_V_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
705         &                         uFld,omega3,hFacZ,r_hFacZ,
706         &                         vCf,myThid)
707           ELSE
708            CALL MOM_VI_V_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
709         &                         uFld,vort3, hFacZ,r_hFacZ,
710         &                         vCf,myThid)
711           ENDIF
712         DO j=jMin,jMax         DO j=jMin,jMax
713          DO i=iMin,iMax          DO i=iMin,iMax
714           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
715          ENDDO          ENDDO
716         ENDDO         ENDDO
717    
718    #ifdef ALLOW_AUTODIFF_TAMC
719    CADJ STORE ucf(:,:) =
720    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
721    CADJ STORE vcf(:,:) =
722    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
723    #endif
724    
725         IF ( writeDiag ) THEN         IF ( writeDiag ) THEN
726          CALL WRITE_LOCAL_RL('zV','I10',1,uCf,bi,bj,k,myIter,myThid)           IF (snapshot_mdsio) THEN
727          CALL WRITE_LOCAL_RL('zU','I10',1,vCf,bi,bj,k,myIter,myThid)             CALL WRITE_LOCAL_RL('zV','I10',1,uCf,bi,bj,k,myIter,myThid)
728               CALL WRITE_LOCAL_RL('zU','I10',1,vCf,bi,bj,k,myIter,myThid)
729             ENDIF
730    #ifdef ALLOW_MNC
731             IF (useMNC .AND. snapshot_mnc) THEN
732               CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'zV', uCf,
733         &          offsets, myThid)
734               CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'zU', vCf,
735         &          offsets, myThid)
736             ENDIF
737    #endif /*  ALLOW_MNC  */
738         ENDIF         ENDIF
739    
740  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
 #ifndef HRCUBE  
741         IF (taveFreq.GT.0.) THEN         IF (taveFreq.GT.0.) THEN
742           CALL TIMEAVE_CUMUL_1K1T(uZetatave,vCf,deltaTClock,           CALL TIMEAVE_CUMUL_1K1T(uZetatave,vCf,deltaTClock,
743       &                           Nr, k, bi, bj, myThid)       &                           Nr, k, bi, bj, myThid)
# Line 446  c      CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K Line 745  c      CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K
745       &                           Nr, k, bi, bj, myThid)       &                           Nr, k, bi, bj, myThid)
746         ENDIF         ENDIF
747  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
748  #endif /* ndef HRCUBE */  #ifdef ALLOW_DIAGNOSTICS
749           IF ( useDiagnostics ) THEN
750             CALL DIAGNOSTICS_FILL(uCf,'Um_AdvZ3',k,1,2,bi,bj,myThid)
751             CALL DIAGNOSTICS_FILL(vCf,'Vm_AdvZ3',k,1,2,bi,bj,myThid)
752           ENDIF
753    #endif /* ALLOW_DIAGNOSTICS */
754    
755  C--   Vertical shear terms (-w*du/dr & -w*dv/dr)  C--   Vertical shear terms (-w*du/dr & -w*dv/dr)
756         IF ( .NOT. momImplVertAdv ) THEN         IF ( .NOT. momImplVertAdv ) THEN
757          CALL MOM_VI_U_VERTSHEAR(bi,bj,K,uVel,wVel,uCf,myThid)          CALL MOM_VI_U_VERTSHEAR(bi,bj,k,uVel,wVel,uCf,myThid)
758          DO j=jMin,jMax          DO j=jMin,jMax
759           DO i=iMin,iMax           DO i=iMin,iMax
760            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)
761           ENDDO           ENDDO
762          ENDDO          ENDDO
763          CALL MOM_VI_V_VERTSHEAR(bi,bj,K,vVel,wVel,vCf,myThid)          CALL MOM_VI_V_VERTSHEAR(bi,bj,k,vVel,wVel,vCf,myThid)
764          DO j=jMin,jMax          DO j=jMin,jMax
765           DO i=iMin,iMax           DO i=iMin,iMax
766            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
767           ENDDO           ENDDO
768          ENDDO          ENDDO
769    #ifdef ALLOW_DIAGNOSTICS
770            IF ( useDiagnostics ) THEN
771             CALL DIAGNOSTICS_FILL(uCf,'Um_AdvRe',k,1,2,bi,bj,myThid)
772             CALL DIAGNOSTICS_FILL(vCf,'Vm_AdvRe',k,1,2,bi,bj,myThid)
773            ENDIF
774    #endif /* ALLOW_DIAGNOSTICS */
775         ENDIF         ENDIF
776    
777  C--   Bernoulli term  C--   Bernoulli term
778         CALL MOM_VI_U_GRAD_KE(bi,bj,K,KE,uCf,myThid)         CALL MOM_VI_U_GRAD_KE(bi,bj,k,KE,uCf,myThid)
779         DO j=jMin,jMax         DO j=jMin,jMax
780          DO i=iMin,iMax          DO i=iMin,iMax
781           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)
782          ENDDO          ENDDO
783         ENDDO         ENDDO
784         CALL MOM_VI_V_GRAD_KE(bi,bj,K,KE,vCf,myThid)         CALL MOM_VI_V_GRAD_KE(bi,bj,k,KE,vCf,myThid)
785         DO j=jMin,jMax         DO j=jMin,jMax
786          DO i=iMin,iMax          DO i=iMin,iMax
787           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
788          ENDDO          ENDDO
789         ENDDO         ENDDO
790         IF ( writeDiag ) THEN         IF ( writeDiag ) THEN
791          CALL WRITE_LOCAL_RL('KEx','I10',1,uCf,bi,bj,k,myIter,myThid)           IF (snapshot_mdsio) THEN
792          CALL WRITE_LOCAL_RL('KEy','I10',1,vCf,bi,bj,k,myIter,myThid)             CALL WRITE_LOCAL_RL('KEx','I10',1,uCf,bi,bj,k,myIter,myThid)
793               CALL WRITE_LOCAL_RL('KEy','I10',1,vCf,bi,bj,k,myIter,myThid)
794             ENDIF
795    #ifdef ALLOW_MNC
796             IF (useMNC .AND. snapshot_mnc) THEN
797               CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'KEx', uCf,
798         &          offsets, myThid)
799               CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'KEy', vCf,
800         &          offsets, myThid)
801             ENDIF
802    #endif /*  ALLOW_MNC  */
803         ENDIF         ENDIF
804    
805  C--   end if momAdvection  C--   end if momAdvection
806        ENDIF        ENDIF
807    
808    C--   3.D Coriolis term (horizontal momentum, Eastward component: -fprime*w)
809          IF ( use3dCoriolis ) THEN
810            CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,uCf,myThid)
811            DO j=jMin,jMax
812             DO i=iMin,iMax
813              gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
814             ENDDO
815            ENDDO
816           IF ( usingCurvilinearGrid ) THEN
817    C-     presently, non zero angleSinC array only supported with Curvilinear-Grid
818            CALL MOM_V_CORIOLIS_NH(bi,bj,k,wVel,vCf,myThid)
819            DO j=jMin,jMax
820             DO i=iMin,iMax
821              gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
822             ENDDO
823            ENDDO
824           ENDIF
825          ENDIF
826    
827    C--   Non-Hydrostatic (spherical) metric terms
828          IF ( useNHMTerms ) THEN
829           CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,uCf,myThid)
830           DO j=jMin,jMax
831            DO i=iMin,iMax
832             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
833            ENDDO
834           ENDDO
835           CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,vCf,myThid)
836           DO j=jMin,jMax
837            DO i=iMin,iMax
838             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
839            ENDDO
840           ENDDO
841          ENDIF
842    
843  C--   Set du/dt & dv/dt on boundaries to zero  C--   Set du/dt & dv/dt on boundaries to zero
844        DO j=jMin,jMax        DO j=jMin,jMax
845         DO i=iMin,iMax         DO i=iMin,iMax
# Line 493  C--   Set du/dt & dv/dt on boundaries to Line 848  C--   Set du/dt & dv/dt on boundaries to
848         ENDDO         ENDDO
849        ENDDO        ENDDO
850    
851    #ifdef ALLOW_DEBUG
852          IF ( debugLevel .GE. debLevC
853         &   .AND. k.EQ.4 .AND. myIter.EQ.nIter0
854         &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
855         &   .AND. useCubedSphereExchange ) THEN
856            CALL DEBUG_CS_CORNER_UV( ' uDiss,vDiss from MOM_VECINV',
857         &             guDiss,gvDiss, k, standardMessageUnit,bi,bj,myThid )
858          ENDIF
859    #endif /* ALLOW_DEBUG */
860    
861        IF ( writeDiag ) THEN        IF ( writeDiag ) THEN
862         CALL WRITE_LOCAL_RL('Ds','I10',1,strain,bi,bj,k,myIter,myThid)          IF (useBiharmonicVisc) THEN
863         CALL WRITE_LOCAL_RL('Dt','I10',1,tension,bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL( 'del2u', 'I10', 1, del2u,
864         CALL WRITE_LOCAL_RL('Du','I10',1,uDiss,bi,bj,k,myIter,myThid)       &                        bi,bj,k, myIter, myThid )
865         CALL WRITE_LOCAL_RL('Dv','I10',1,vDiss,bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL( 'del2v', 'I10', 1, del2v,
866         CALL WRITE_LOCAL_RL('Z3','I10',1,vort3,bi,bj,k,myIter,myThid)       &                        bi,bj,k, myIter, myThid )
867  c      CALL WRITE_LOCAL_RL('W3','I10',1,omega3,bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL( 'dStar', 'I10', 1, dStar,
868         CALL WRITE_LOCAL_RL('KE','I10',1,KE,bi,bj,k,myIter,myThid)       &                        bi,bj,k, myIter, myThid )
869         CALL WRITE_LOCAL_RL('D','I10',1,hdiv,bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL( 'zStar', 'I10', 1, zStar,
870         &                        bi,bj,k, myIter, myThid )
871            ENDIF
872            IF (snapshot_mdsio) THEN
873             CALL WRITE_LOCAL_RL('W3','I10',1,omega3, bi,bj,k,myIter,myThid)
874             CALL WRITE_LOCAL_RL('Z3','I10',1,vort3BC,bi,bj,k,myIter,myThid)
875             CALL WRITE_LOCAL_RL('KE','I10',1,KE,     bi,bj,k,myIter,myThid)
876             CALL WRITE_LOCAL_RL('D', 'I10',1,hDiv,   bi,bj,k,myIter,myThid)
877             CALL WRITE_LOCAL_RL('Dt','I10',1,tension,bi,bj,k,myIter,myThid)
878             CALL WRITE_LOCAL_RL( 'Ds', 'I10', 1, strainBC,
879         &                        bi,bj,k, myIter, myThid )
880             CALL WRITE_LOCAL_RL('Du','I10',1,guDiss, bi,bj,k,myIter,myThid)
881             CALL WRITE_LOCAL_RL('Dv','I10',1,gvDiss, bi,bj,k,myIter,myThid)
882            ENDIF
883    #ifdef ALLOW_MNC
884            IF (useMNC .AND. snapshot_mnc) THEN
885              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'W3',omega3,
886         &          offsets, myThid)
887              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Z3',vort3BC,
888         &          offsets, myThid)
889              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'KE',KE,
890         &          offsets, myThid)
891              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'D', hDiv,
892         &          offsets, myThid)
893              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dt',tension,
894         &          offsets, myThid)
895              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Ds',strainBC,
896         &          offsets, myThid)
897              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Du',guDiss,
898         &          offsets, myThid)
899              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dv',gvDiss,
900         &          offsets, myThid)
901            ENDIF
902    #endif /*  ALLOW_MNC  */
903          ENDIF
904    
905    #ifdef ALLOW_DIAGNOSTICS
906          IF ( useDiagnostics ) THEN
907            CALL DIAGNOSTICS_FILL(vort3BC,'momVort3',k,1,2,bi,bj,myThid)
908            CALL DIAGNOSTICS_FILL(KE,     'momKE   ',k,1,2,bi,bj,myThid)
909           IF (momViscosity) THEN
910            CALL DIAGNOSTICS_FILL(hDiv,   'momHDiv ',k,1,2,bi,bj,myThid)
911           ENDIF
912           IF ( useVariableVisc .OR. useStrainTensionVisc ) THEN
913            CALL DIAGNOSTICS_FILL(tension, 'Tension ',k,1,2,bi,bj,myThid)
914            CALL DIAGNOSTICS_FILL(strainBC,'Strain  ',k,1,2,bi,bj,myThid)
915           ENDIF
916            CALL DIAGNOSTICS_FILL(gU(1-OLx,1-OLy,k,bi,bj),
917         &                                'Um_Advec',k,1,2,bi,bj,myThid)
918            CALL DIAGNOSTICS_FILL(gV(1-OLx,1-OLy,k,bi,bj),
919         &                                'Vm_Advec',k,1,2,bi,bj,myThid)
920        ENDIF        ENDIF
921    #endif /* ALLOW_DIAGNOSTICS */
922    
923  #endif /* ALLOW_MOM_VECINV */  #endif /* ALLOW_MOM_VECINV */
924    
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.81
  ViewVC Help
Powered by ViewVC 1.1.22