/[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.64 by jmc, Mon Apr 25 02:03:13 2011 UTC revision 1.80 by mlosch, Fri Apr 28 17:17:14 2017 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "MOM_VECINV_OPTIONS.h"  #include "MOM_VECINV_OPTIONS.h"
5    #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        KappaRU, KappaRV,       I        kappaRU, kappaRV,
15       U        fVerU, fVerV,       I        fVerUkm, fVerVkm,
16         O        fVerUkp, fVerVkp,
17       O        guDiss, gvDiss,       O        guDiss, gvDiss,
18       I        myTime, myIter, myThid)       I        myTime, myIter, myThid )
19  C     /==========================================================\  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"
 #ifdef ALLOW_MNC  
 #include "MNC_PARAMS.h"  
 #endif  
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  #endif
52  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
53  # include "tamc.h"  # include "tamc.h"
# Line 44  C     == Global variables == Line 55  C     == Global variables ==
55  #endif  #endif
56    
57  C     == Routine arguments ==  C     == Routine arguments ==
58  C     fVerU  :: Flux of momentum in the vertical direction, out of the upper  C     bi,bj   :: current tile indices
59  C     fVerV  :: face of a cell K ( flux into the cell above ).  C     k       :: current vertical level
60  C     guDiss :: dissipation tendency (all explicit terms), u component  C     iMin,iMax,jMin,jMax :: loop ranges
61  C     gvDiss :: dissipation tendency (all explicit terms), v component  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 :: my Thread Id number  C     fVerUkp :: vertical viscous flux of U, interface below (k+1/2)
66        _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     fVerVkp :: vertical viscous flux of V, interface below (k+1/2)
67        _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
68        _RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  C     guDiss  :: dissipation tendency (all explicit terms), u component
69        _RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  C     gvDiss  :: dissipation tendency (all explicit terms), v component
70    C     myTime  :: current time
71    C     myIter  :: current time-step number
72    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)        _RL guDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82        _RL gvDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL gvDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       INTEGER kUp,kDown  
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 71  C     == Functions == Line 91  C     == Functions ==
91        EXTERNAL DIFFERENT_MULTIPLE        EXTERNAL DIFFERENT_MULTIPLE
92    
93  C     == Local variables ==  C     == Local variables ==
94    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        _RS hFacZ   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS hFacZ   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101          _RS h0FacZ  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RS r_hFacZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS r_hFacZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103        _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104        _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 85  C     == Local variables == Line 108  C     == Local variables ==
108        _RL zStar   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL zStar   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109        _RL tension (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL tension (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110        _RL strain  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL strain  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111          _RL strainBC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112        _RL KE      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KE      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113        _RL omega3  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL omega3  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114        _RL vort3   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vort3   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115          _RL vort3BC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116        _RL hDiv    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL hDiv    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117        _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118        _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119        _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120        _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121  C     i,j,k  :: Loop counters  C     i,j    :: Loop counters
122        INTEGER i,j,k        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.
124        _RL  ArDudrFac        _RL  ArDudrFac
125        _RL  ArDvdrFac        _RL  ArDvdrFac
126        _RL  sideMaskFac        _RL  sideMaskFac
127        LOGICAL bottomDragTerms        LOGICAL bottomDragTerms
128        LOGICAL writeDiag        LOGICAL writeDiag
       LOGICAL harmonic,biharmonic,useVariableViscosity  
129  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
130        INTEGER imomkey        INTEGER imomkey
131  #endif  #endif
# Line 111  C     xxxFac - On-off tracer parameters Line 135  C     xxxFac - On-off tracer parameters
135        CHARACTER*(1) pf        CHARACTER*(1) pf
136  #endif  #endif
137    
138  #ifdef ALLOW_AUTODIFF_TAMC  #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  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 178  C--   Initialise intermediate terms Line 202  C--   Initialise intermediate terms
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)  C-    need to initialise hDiv for MOM_VI_DEL2UV(call FILL_CS_CORNER_TR_RL)
204          hDiv(i,j)  = 0.          hDiv(i,j)  = 0.
205          viscAh_Z(i,j) = 0.  c       viscAh_Z(i,j) = 0.
206          viscAh_D(i,j) = 0.  c       viscAh_D(i,j) = 0.
207          viscA4_Z(i,j) = 0.  c       viscA4_Z(i,j) = 0.
208          viscA4_D(i,j) = 0.  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_TAMC  #ifdef ALLOW_AUTODIFF
213          hFacZ(i,j)   = 0. _d 0          hFacZ(i,j)   = 0. _d 0
214  #endif  #endif
215         ENDDO         ENDDO
# Line 205  C       vorticity at a no-slip boundary Line 229  C       vorticity at a no-slip boundary
229          sideMaskFac = 0. _d 0          sideMaskFac = 0. _d 0
230        ENDIF        ENDIF
231    
232        IF (     no_slip_bottom        IF ( selectImplicitDrag.EQ.0 .AND.
233       &    .OR. bottomDragQuadratic.NE.0.       &      (  no_slip_bottom
234       &    .OR. bottomDragLinear.NE.0.) THEN       &    .OR. selectBotDragQuadr.GE.0
235         &    .OR. bottomDragLinear.NE.0. ) ) THEN
236         bottomDragTerms=.TRUE.         bottomDragTerms=.TRUE.
237        ELSE        ELSE
238         bottomDragTerms=.FALSE.         bottomDragTerms=.FALSE.
# Line 224  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)
# Line 232  c     CALL MOM_VI_HFACZ_DISS(bi,bj,k,hFa Line 272  c     CALL MOM_VI_HFACZ_DISS(bi,bj,k,hFa
272    
273        CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)        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
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  C--    For viscous term, compute horizontal divergence, tension & strain
304  C      and mask relative vorticity (free-slip case):  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  #ifdef ALLOW_AUTODIFF_TAMC
324  CADJ STORE vort3(:,:) =  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  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
328  #endif  #endif
329    
330         CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)         CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
331    
332         CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid)         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_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,strain,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-     account for no-slip / free-slip BC:  C--    Calculate Lateral Viscosities
359         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
360          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
361            IF ( hFacZ(i,j).EQ.0. ) THEN           viscAh_D(i,j) = viscAhD
362              vort3(i,j)  = sideMaskFac*vort3(i,j)           viscAh_Z(i,j) = viscAhZ
363              strain(i,j) = sideMaskFac*strain(i,j)           viscA4_D(i,j) = viscA4D
364            ENDIF           viscA4_Z(i,j) = viscA4Z
365          ENDDO          ENDDO
366         ENDDO         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  C--    Calculate Viscosities  #ifdef ALLOW_AUTODIFF_TAMC
376         CALL MOM_CALC_VISC(  CADJ STORE viscAh_Z(:,:) =
377       I        bi,bj,k,  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
378       O        viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,  CADJ STORE viscAh_D(:,:) =
379       O        harmonic,biharmonic,useVariableViscosity,  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
380       I        hDiv,vort3,tension,strain,KE,hfacZ,  CADJ STORE viscA4_Z(:,:) =
381       I        myThid)  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 (biharmonic) 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    #ifdef ALLOW_AUTODIFF_TAMC
401    CADJ STORE del2u(:,:) =
402    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
403    CADJ STORE del2v(:,:) =
404    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
405    #endif
406           CALL MOM_CALC_HDIV(bi,bj,k,2,del2u,del2v,dStar,myThid)           CALL MOM_CALC_HDIV(bi,bj,k,2,del2u,del2v,dStar,myThid)
407           CALL MOM_CALC_RELVORT3(bi,bj,k,           CALL MOM_CALC_RELVORT3(bi,bj,k,
408       &                          del2u,del2v,hFacZ,zStar,myThid)       &                          del2u,del2v,hFacZ,zStar,myThid)
          IF ( writeDiag ) THEN  
            CALL WRITE_LOCAL_RL( 'del2u', 'I10', 1, del2u,  
      &                           bi,bj,k, myIter, myThid )  
            CALL WRITE_LOCAL_RL( 'del2v', 'I10', 1, del2v,  
      &                           bi,bj,k, myIter, myThid )  
            CALL WRITE_LOCAL_RL( 'dStar', 'I10', 1, dStar,  
      &                           bi,bj,k, myIter, myThid )  
            CALL WRITE_LOCAL_RL( 'zStar', 'I10', 1, zStar,  
      &                           bi,bj,k, myIter, myThid )  
          ENDIF  
409         ENDIF         ENDIF
410    
411  C-    Strain diagnostics:  #ifdef ALLOW_AUTODIFF_TAMC
412         IF ( writeDiag ) THEN  CADJ STORE del2u(:,:) =
413          IF (snapshot_mdsio) THEN  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
414            CALL WRITE_LOCAL_RL('Ds','I10',1,strain,bi,bj,k,myIter,myThid)  CADJ STORE del2v(:,:) =
415          ENDIF  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
416  #ifdef ALLOW_MNC  CADJ STORE dStar(:,:) =
417          IF (useMNC .AND. snapshot_mnc) THEN  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
418            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Ds',strain,  CADJ STORE zStar(:,:) =
419       &          offsets, myThid)  CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
420          ENDIF  #endif
 #endif /*  ALLOW_MNC  */  
        ENDIF  
 #ifdef ALLOW_DIAGNOSTICS  
        IF ( useDiagnostics ) THEN  
         CALL DIAGNOSTICS_FILL(strain, 'Strain  ',k,1,2,bi,bj,myThid)  
        ENDIF  
 #endif /* ALLOW_DIAGNOSTICS */  
421    
422  C---   Calculate dissipation terms for U and V equations  C---   Calculate dissipation terms for U and V equations
423    
424  C      in terms of tension and strain  C-     in terms of tension and strain
425         IF (useStrainTensionVisc) THEN         IF (useStrainTensionVisc) THEN
426  C        mask strain as if free-slip since side-drag is computed separately  C      use masked strain as if free-slip since side-drag is computed separately
427           DO j=1-Oly,sNy+Oly           CALL MOM_HDISSIP( bi, bj, k,
428            DO i=1-Olx,sNx+Olx       I            tension, strain, hFacZ,
429              IF ( hFacZ(i,j).EQ.0. ) strain(i,j) = 0. _d 0       I            viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
430            ENDDO       I            useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
431           ENDDO       O            guDiss, gvDiss,
432           CALL MOM_HDISSIP(bi,bj,k,hDiv,vort3,tension,strain,KE,       I            myThid )
      I                    hFacZ,  
      I                    viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,  
      I                    harmonic,biharmonic,useVariableViscosity,  
      O                    guDiss,gvDiss,  
      I                    myThid)  
433         ELSE         ELSE
434  C      in terms of vorticity and divergence  C-     in terms of vorticity and divergence
435           CALL MOM_VI_HDISSIP(bi,bj,k,hDiv,vort3,tension,strain,KE,           CALL MOM_VI_HDISSIP( bi, bj, k,
436       I                       hFacZ,dStar,zStar,       I            hDiv, vort3, dStar, zStar, hFacZ,
437       I                       viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,       I            viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
438       I                       harmonic,biharmonic,useVariableViscosity,       I            useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
439       O                       guDiss,gvDiss,       O            guDiss, gvDiss,
440       &                       myThid)       I            myThid )
441         ENDIF         ENDIF
 C--   if (momViscosity) end of block.  
       ENDIF  
   
 C-    Return to standard hfacZ (min-4) and mask vort3 accordingly:  
 c     CALL MOM_VI_MASK_VORT3(bi,bj,k,hFacZ,r_hFacZ,vort3,myThid)  
442    
443  C---  Other dissipation terms in Zonal momentum equation  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) THEN        IF ( .NOT.implicitViscosity ) THEN
448         CALL MOM_U_RVISCFLUX(bi,bj,k+1,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          ENDDO
454         ENDDO         ENDDO
455    
456    #ifdef ALLOW_AUTODIFF_TAMC
457    CADJ STORE fVerUkp(:,:) =
458    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
459    #endif
460    
461  C--   Tendency is minus divergence of the fluxes  C--   Tendency is minus divergence of the fluxes
462         DO j=2-Oly,sNy+Oly-1  C     vert.visc.flx is scaled by deepFac2F (deep-atmos) and rhoFac (anelastic)
463          DO i=2-Olx,sNx+Olx-1         DO j=jMin,jMax
464            DO i=iMin,iMax
465           guDiss(i,j) = guDiss(i,j)           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,kDown) - fVerU(i,j,kUp)       &   *recip_deepFac2C(k)*recip_rhoFacC(k)
      &   )*rkSign  
470          ENDDO          ENDDO
471         ENDDO         ENDDO
472        ENDIF        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(         CALL MOM_U_SIDEDRAG( bi, bj, k,
478       I        bi,bj,k,       I          uFld, del2u, h0FacZ,
479       I        uFld, del2u, hFacZ,       I          viscAh_Z, viscA4_Z,
480       I        viscAh_Z,viscA4_Z,       I          useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
481       I        harmonic,biharmonic,useVariableViscosity,       O          vF,
482       O        vF,       I          myThid )
      I        myThid)  
483         DO j=jMin,jMax         DO j=jMin,jMax
484          DO i=iMin,iMax          DO i=iMin,iMax
485           guDiss(i,j) = guDiss(i,j)+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           guDiss(i,j) = guDiss(i,j)+vF(i,j)           guDiss(i,j) = guDiss(i,j)+vF(i,j)
# Line 388  C-    No-slip BCs impose a drag at botto Line 500  C-    No-slip BCs impose a drag at botto
500         ENDDO         ENDDO
501        ENDIF        ENDIF
502  #ifdef ALLOW_SHELFICE  #ifdef ALLOW_SHELFICE
503        IF (useShelfIce.AND.momViscosity.AND.bottomDragTerms) THEN        IF ( useShelfIce ) THEN
504         CALL SHELFICE_U_DRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)         CALL SHELFICE_U_DRAG( bi, bj, k,
505         I            uFld, vFld, KE, kappaRU,
506         O            vF,
507         I            myThid )
508         DO j=jMin,jMax         DO j=jMin,jMax
509          DO i=iMin,iMax          DO i=iMin,iMax
510           guDiss(i,j) = guDiss(i,j) + vF(i,j)           guDiss(i,j) = guDiss(i,j) + vF(i,j)
# Line 398  C-    No-slip BCs impose a drag at botto Line 513  C-    No-slip BCs impose a drag at botto
513        ENDIF        ENDIF
514  #endif /* ALLOW_SHELFICE */  #endif /* ALLOW_SHELFICE */
515    
   
516  C---  Other dissipation terms in Meridional momentum equation  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) THEN        IF ( .NOT.implicitViscosity ) THEN
521         CALL MOM_V_RVISCFLUX(bi,bj,k+1,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          ENDDO
527         ENDDO         ENDDO
528    #ifdef ALLOW_AUTODIFF_TAMC
529    CADJ STORE fVerVkp(:,:) =
530    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
531    #endif
532  C--   Tendency is minus divergence of the fluxes  C--   Tendency is minus divergence of the fluxes
533    C     vert.visc.flx is scaled by deepFac2F (deep-atmos) and rhoFac (anelastic)
534         DO j=jMin,jMax         DO j=jMin,jMax
535          DO i=iMin,iMax          DO i=iMin,iMax
536           gvDiss(i,j) = gvDiss(i,j)           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,kDown) - fVerV(i,j,kUp)       &   *recip_deepFac2C(k)*recip_rhoFacC(k)
      &   )*rkSign  
541          ENDDO          ENDDO
542         ENDDO         ENDDO
543        ENDIF        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(         CALL MOM_V_SIDEDRAG( bi, bj, k,
549       I        bi,bj,k,       I          vFld, del2v, h0FacZ,
550       I        vFld, del2v, hFacZ,       I          viscAh_Z, viscA4_Z,
551       I        viscAh_Z,viscA4_Z,       I          useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
552       I        harmonic,biharmonic,useVariableViscosity,       O          vF,
553       O        vF,       I          myThid )
      I        myThid)  
554         DO j=jMin,jMax         DO j=jMin,jMax
555          DO i=iMin,iMax          DO i=iMin,iMax
556           gvDiss(i,j) = gvDiss(i,j)+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           gvDiss(i,j) = gvDiss(i,j)+vF(i,j)           gvDiss(i,j) = gvDiss(i,j)+vF(i,j)
# Line 453  C-    No-slip BCs impose a drag at botto Line 571  C-    No-slip BCs impose a drag at botto
571         ENDDO         ENDDO
572        ENDIF        ENDIF
573  #ifdef ALLOW_SHELFICE  #ifdef ALLOW_SHELFICE
574        IF (useShelfIce.AND.momViscosity.AND.bottomDragTerms) THEN        IF ( useShelfIce ) THEN
575           CALL SHELFICE_V_DRAG(bi,bj,k,vFld,KE,KappaRU,vF,myThid)         CALL SHELFICE_V_DRAG( bi, bj, k,
576           DO j=jMin,jMax       I            uFld, vFld, KE, kappaRV,
577            DO i=iMin,iMax       O            vF,
578             gvDiss(i,j) = gvDiss(i,j) + vF(i,j)       I            myThid )
579            ENDDO         DO j=jMin,jMax
580           ENDDO          DO i=iMin,iMax
581          ENDIF           gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
582            ENDDO
583           ENDDO
584          ENDIF
585  #endif /* ALLOW_SHELFICE */  #endif /* ALLOW_SHELFICE */
586    
587    C--   if (momViscosity) end of block.
 C-    Vorticity diagnostics:  
       IF ( writeDiag ) THEN  
         IF (snapshot_mdsio) THEN  
           CALL WRITE_LOCAL_RL('Z3','I10',1,vort3, bi,bj,k,myIter,myThid)  
         ENDIF  
 #ifdef ALLOW_MNC  
         IF (useMNC .AND. snapshot_mnc) THEN  
           CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Z3',vort3,  
      &          offsets, myThid)  
         ENDIF  
 #endif /*  ALLOW_MNC  */  
       ENDIF  
 #ifdef ALLOW_DIAGNOSTICS  
       IF ( useDiagnostics ) THEN  
         CALL DIAGNOSTICS_FILL(vort3,  'momVort3',k,1,2,bi,bj,myThid)  
588        ENDIF        ENDIF
589  #endif /* ALLOW_DIAGNOSTICS */  
590    C-    Return to standard hfacZ (min-4) and mask vort3 accordingly:
591    c     CALL MOM_VI_MASK_VORT3(bi,bj,k,hFacZ,r_hFacZ,vort3,myThid)
592    
593  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
594    
595  C---  Prepare for Advection & Coriolis terms:  C---  Prepare for Advection & Coriolis terms:
596  C-    Mask relative vorticity and calculate absolute vorticity  C-    calculate absolute vorticity
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
          IF ( hFacZ(i,j).EQ.0. ) vort3(i,j) = 0.  
        ENDDO  
       ENDDO  
597        IF (useAbsVorticity)        IF (useAbsVorticity)
598       &  CALL MOM_CALC_ABSVORT3(bi,bj,k,vort3,omega3,myThid)       &  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  c     IF (useCoriolis .AND. .NOT.useCDscheme  c     IF (useCoriolis .AND. .NOT.useCDscheme
607  c    &    .AND. .NOT. useAbsVorticity) THEN  c    &    .AND. .NOT. useAbsVorticity) THEN
# Line 502  C- jmc: change it to keep the Coriolis t Line 610  C- jmc: change it to keep the Coriolis t
610       &     .NOT.( useCDscheme .OR. useAbsVorticity.AND.momAdvection )       &     .NOT.( useCDscheme .OR. useAbsVorticity.AND.momAdvection )
611       &   ) THEN       &   ) THEN
612         IF (useAbsVorticity) THEN         IF (useAbsVorticity) THEN
613          CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,hFacZ,r_hFacZ,          CALL MOM_VI_U_CORIOLIS(bi,bj,k,vFld,omega3,hFacZ,r_hFacZ,
614       &                         uCf,myThid)       &                         uCf,myThid)
615          CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,hFacZ,r_hFacZ,          CALL MOM_VI_V_CORIOLIS(bi,bj,k,uFld,omega3,hFacZ,r_hFacZ,
616       &                         vCf,myThid)       &                         vCf,myThid)
617         ELSE         ELSE
618          CALL MOM_VI_CORIOLIS(bi,bj,k,uFld,vFld,hFacZ,r_hFacZ,          CALL MOM_VI_CORIOLIS(bi,bj,k,uFld,vFld,hFacZ,r_hFacZ,
# Line 545  C- jmc: change it to keep the Coriolis t Line 653  C- jmc: change it to keep the Coriolis t
653         ENDDO         ENDDO
654        ENDIF        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--   Horizontal advection of relative (or absolute) vorticity  C--   Horizontal advection of relative (or absolute) vorticity
665         IF ( (highOrderVorticity.OR.upwindVorticity)         IF ( (highOrderVorticity.OR.upwindVorticity)
666       &     .AND.useAbsVorticity ) THEN       &     .AND.useAbsVorticity ) THEN
667          CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,vFld,omega3,r_hFacZ,          CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,selectVortScheme,
668         &                         highOrderVorticity,upwindVorticity,
669         &                         vFld,omega3,r_hFacZ,
670       &                         uCf,myThid)       &                         uCf,myThid)
671         ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN         ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
672          CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,vFld,vort3, r_hFacZ,          CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,selectVortScheme,
673         &                         highOrderVorticity, upwindVorticity,
674         &                         vFld,vort3, r_hFacZ,
675       &                         uCf,myThid)       &                         uCf,myThid)
676         ELSEIF ( useAbsVorticity ) THEN         ELSEIF ( useAbsVorticity ) THEN
677          CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,hFacZ,r_hFacZ,          CALL MOM_VI_U_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
678         &                         vFld,omega3,hFacZ,r_hFacZ,
679       &                         uCf,myThid)       &                         uCf,myThid)
680         ELSE         ELSE
681          CALL MOM_VI_U_CORIOLIS(bi,bj,k,vFld,vort3, hFacZ,r_hFacZ,          CALL MOM_VI_U_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
682         &                         vFld,vort3, hFacZ,r_hFacZ,
683       &                         uCf,myThid)       &                         uCf,myThid)
684         ENDIF         ENDIF
685         DO j=jMin,jMax         DO j=jMin,jMax
# Line 568  C--   Horizontal advection of relative ( Line 689  C--   Horizontal advection of relative (
689         ENDDO         ENDDO
690         IF ( (highOrderVorticity.OR.upwindVorticity)         IF ( (highOrderVorticity.OR.upwindVorticity)
691       &     .AND.useAbsVorticity ) THEN       &     .AND.useAbsVorticity ) THEN
692          CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,omega3,r_hFacZ,          CALL MOM_VI_V_CORIOLIS_C4(bi,bj,k,selectVortScheme,
693         &                         highOrderVorticity, upwindVorticity,
694         &                         uFld,omega3,r_hFacZ,
695       &                         vCf,myThid)       &                         vCf,myThid)
696         ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN         ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
697          CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,vort3, r_hFacZ,          CALL MOM_VI_V_CORIOLIS_C4(bi,bj,k,selectVortScheme,
698         &                         highOrderVorticity, upwindVorticity,
699         &                         uFld,vort3, r_hFacZ,
700       &                         vCf,myThid)       &                         vCf,myThid)
701         ELSEIF ( useAbsVorticity ) THEN         ELSEIF ( useAbsVorticity ) THEN
702          CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,hFacZ,r_hFacZ,          CALL MOM_VI_V_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
703         &                         uFld,omega3,hFacZ,r_hFacZ,
704       &                         vCf,myThid)       &                         vCf,myThid)
705         ELSE         ELSE
706          CALL MOM_VI_V_CORIOLIS(bi,bj,k,uFld,vort3, hFacZ,r_hFacZ,          CALL MOM_VI_V_CORIOLIS(bi,bj,k,selectVortScheme,useJamartMomAdv,
707         &                         uFld,vort3, hFacZ,r_hFacZ,
708       &                         vCf,myThid)       &                         vCf,myThid)
709         ENDIF         ENDIF
710         DO j=jMin,jMax         DO j=jMin,jMax
# Line 586  C--   Horizontal advection of relative ( Line 713  C--   Horizontal advection of relative (
713          ENDDO          ENDDO
714         ENDDO         ENDDO
715    
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         IF ( writeDiag ) THEN
724           IF (snapshot_mdsio) THEN           IF (snapshot_mdsio) THEN
725             CALL WRITE_LOCAL_RL('zV','I10',1,uCf,bi,bj,k,myIter,myThid)             CALL WRITE_LOCAL_RL('zV','I10',1,uCf,bi,bj,k,myIter,myThid)
# Line 618  C--   Horizontal advection of relative ( Line 752  C--   Horizontal advection of relative (
752    
753  C--   Vertical shear terms (-w*du/dr & -w*dv/dr)  C--   Vertical shear terms (-w*du/dr & -w*dv/dr)
754         IF ( .NOT. momImplVertAdv ) THEN         IF ( .NOT. momImplVertAdv ) THEN
755          CALL MOM_VI_U_VERTSHEAR(bi,bj,K,uVel,wVel,uCf,myThid)          CALL MOM_VI_U_VERTSHEAR(bi,bj,k,uVel,wVel,uCf,myThid)
756          DO j=jMin,jMax          DO j=jMin,jMax
757           DO i=iMin,iMax           DO i=iMin,iMax
758            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)
759           ENDDO           ENDDO
760          ENDDO          ENDDO
761          CALL MOM_VI_V_VERTSHEAR(bi,bj,K,vVel,wVel,vCf,myThid)          CALL MOM_VI_V_VERTSHEAR(bi,bj,k,vVel,wVel,vCf,myThid)
762          DO j=jMin,jMax          DO j=jMin,jMax
763           DO i=iMin,iMax           DO i=iMin,iMax
764            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)
# Line 639  C--   Vertical shear terms (-w*du/dr & - Line 773  C--   Vertical shear terms (-w*du/dr & -
773         ENDIF         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          ENDDO          ENDDO
781         ENDDO         ENDDO
782         CALL MOM_VI_V_GRAD_KE(bi,bj,K,KE,vCf,myThid)         CALL MOM_VI_V_GRAD_KE(bi,bj,k,KE,vCf,myThid)
783         DO j=jMin,jMax         DO j=jMin,jMax
784          DO i=iMin,iMax          DO i=iMin,iMax
785           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)
# Line 713  C--   Set du/dt & dv/dt on boundaries to Line 847  C--   Set du/dt & dv/dt on boundaries to
847        ENDDO        ENDDO
848    
849  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
850        IF ( debugLevel .GE. debLevB        IF ( debugLevel .GE. debLevC
851       &   .AND. k.EQ.4 .AND. myIter.EQ.nIter0       &   .AND. k.EQ.4 .AND. myIter.EQ.nIter0
852       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
853       &   .AND. useCubedSphereExchange ) THEN       &   .AND. useCubedSphereExchange ) THEN
# Line 723  C--   Set du/dt & dv/dt on boundaries to Line 857  C--   Set du/dt & dv/dt on boundaries to
857  #endif /* ALLOW_DEBUG */  #endif /* ALLOW_DEBUG */
858    
859        IF ( writeDiag ) THEN        IF ( writeDiag ) THEN
860            IF (useBiharmonicVisc) THEN
861             CALL WRITE_LOCAL_RL( 'del2u', 'I10', 1, del2u,
862         &                        bi,bj,k, myIter, myThid )
863             CALL WRITE_LOCAL_RL( 'del2v', 'I10', 1, del2v,
864         &                        bi,bj,k, myIter, myThid )
865             CALL WRITE_LOCAL_RL( 'dStar', 'I10', 1, dStar,
866         &                        bi,bj,k, myIter, myThid )
867             CALL WRITE_LOCAL_RL( 'zStar', 'I10', 1, zStar,
868         &                        bi,bj,k, myIter, myThid )
869            ENDIF
870          IF (snapshot_mdsio) THEN          IF (snapshot_mdsio) THEN
871           CALL WRITE_LOCAL_RL('W3','I10',1,omega3, bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL('W3','I10',1,omega3, bi,bj,k,myIter,myThid)
872             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)           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)           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)           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)           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)           CALL WRITE_LOCAL_RL('Dv','I10',1,gvDiss, bi,bj,k,myIter,myThid)
880          ENDIF          ENDIF
# Line 735  C--   Set du/dt & dv/dt on boundaries to Line 882  C--   Set du/dt & dv/dt on boundaries to
882          IF (useMNC .AND. snapshot_mnc) THEN          IF (useMNC .AND. snapshot_mnc) THEN
883            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'W3',omega3,            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'W3',omega3,
884       &          offsets, myThid)       &          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,            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'KE',KE,
888       &          offsets, myThid)       &          offsets, myThid)
889            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'D', hDiv,            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'D', hDiv,
890       &          offsets, myThid)       &          offsets, myThid)
891            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dt',tension,            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dt',tension,
892       &          offsets, myThid)       &          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,            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Du',guDiss,
896       &          offsets, myThid)       &          offsets, myThid)
897            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dv',gvDiss,            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dv',gvDiss,
# Line 751  C--   Set du/dt & dv/dt on boundaries to Line 902  C--   Set du/dt & dv/dt on boundaries to
902    
903  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
904        IF ( useDiagnostics ) THEN        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)          CALL DIAGNOSTICS_FILL(KE,     'momKE   ',k,1,2,bi,bj,myThid)
907         IF (momViscosity) THEN         IF (momViscosity) THEN
908          CALL DIAGNOSTICS_FILL(hDiv,   'momHDiv ',k,1,2,bi,bj,myThid)          CALL DIAGNOSTICS_FILL(hDiv,   'momHDiv ',k,1,2,bi,bj,myThid)
         CALL DIAGNOSTICS_FILL(tension,'Tension ',k,1,2,bi,bj,myThid)  
         CALL DIAGNOSTICS_FILL(guDiss, 'Um_Diss ',k,1,2,bi,bj,myThid)  
         CALL DIAGNOSTICS_FILL(gvDiss, 'Vm_Diss ',k,1,2,bi,bj,myThid)  
909         ENDIF         ENDIF
910          CALL DIAGNOSTICS_FILL(gU(1-Olx,1-Oly,k,bi,bj),         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)       &                                'Um_Advec',k,1,2,bi,bj,myThid)
916          CALL DIAGNOSTICS_FILL(gV(1-Olx,1-Oly,k,bi,bj),          CALL DIAGNOSTICS_FILL(gV(1-OLx,1-OLy,k,bi,bj),
917       &                                'Vm_Advec',k,1,2,bi,bj,myThid)       &                                'Vm_Advec',k,1,2,bi,bj,myThid)
918        ENDIF        ENDIF
919  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
Legend:
Removed from v.1.64  
changed lines
  Added in v.1.80
  ViewVC Help
Powered by ViewVC 1.1.22