/[MITgcm]/MITgcm/model/src/dynamics.F
ViewVC logotype

Diff of /MITgcm/model/src/dynamics.F

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

revision 1.28 by cnh, Thu Aug 20 19:25:05 1998 UTC revision 1.53 by heimbach, Mon Sep 11 23:07:29 2000 UTC
# Line 20  C     | ===== Line 20  C     | =====
20  C     | C*P* comments indicating place holders for which code is |  C     | C*P* comments indicating place holders for which code is |
21  C     |      presently being developed.                          |  C     |      presently being developed.                          |
22  C     \==========================================================/  C     \==========================================================/
23          IMPLICIT NONE
24    
25  C     == Global variables ===  C     == Global variables ===
26  #include "SIZE.h"  #include "SIZE.h"
# Line 27  C     == Global variables === Line 28  C     == Global variables ===
28  #include "CG2D.h"  #include "CG2D.h"
29  #include "PARAMS.h"  #include "PARAMS.h"
30  #include "DYNVARS.h"  #include "DYNVARS.h"
31    #include "GRID.h"
32    
33    #ifdef ALLOW_AUTODIFF_TAMC
34    # include "tamc.h"
35    # include "tamc_keys.h"
36    # ifdef ALLOW_KPP
37    #  include "KPP.h"
38    # endif
39    #endif /* ALLOW_AUTODIFF_TAMC */
40    
41  C     == Routine arguments ==  C     == Routine arguments ==
42  C     myTime - Current time in simulation  C     myTime - Current time in simulation
43  C     myIter - Current iteration number in simulation  C     myIter - Current iteration number in simulation
44  C     myThid - Thread number for this instance of the routine.  C     myThid - Thread number for this instance of the routine.
       INTEGER myThid  
45        _RL myTime        _RL myTime
46        INTEGER myIter        INTEGER myIter
47          INTEGER myThid
48    
49  C     == Local variables  C     == Local variables
50  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
51  C     uTrans, vTrans, wTrans - Per block temporaries holding flow transport  C     uTrans, vTrans, rTrans - Per block temporaries holding flow
52  C     wVel                     o uTrans: Zonal transport  C                              transport
53    C     rVel                     o uTrans: Zonal transport
54  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
55  C                              o wTrans: Vertical transport  C                              o rTrans: Vertical transport
56  C                              o wVel:   Vertical velocity at upper and lower  C                              o rVel:   Vertical velocity at upper and
57  C                                        cell faces.  C                                        lower cell faces.
58  C     maskC,maskUp             o maskC: land/water mask for tracer cells  C     maskC,maskUp             o maskC: land/water mask for tracer cells
59  C                              o maskUp: land/water mask for W points  C                              o maskUp: land/water mask for W points
60  C     aTerm, xTerm, cTerm    - Work arrays for holding separate terms in  C     aTerm, xTerm, cTerm    - Work arrays for holding separate terms in
# Line 59  C                              o fVer: V Line 70  C                              o fVer: V
70  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
71  C                                      so we need an fVer for each  C                                      so we need an fVer for each
72  C                                      variable.  C                                      variable.
73  C     rhoK, rhoKM1   - Density at current level, level above and level below.  C     rhoK, rhoKM1   - Density at current level, level above and level
74    C                      below.
75  C     rhoKP1                                                                    C     rhoKP1                                                                  
76  C     buoyK, buoyKM1 - Buoyancy at current level and level above.  C     buoyK, buoyKM1 - Buoyancy at current level and level above.
77  C     phiHyd         - Hydrostatic part of the potential phi.  C     phiHyd         - Hydrostatic part of the potential phiHydi.
78  C                      In z coords phiHyd is the hydrostatic pressure anomaly  C                      In z coords phiHydiHyd is the hydrostatic
79  C                      In p coords phiHyd is the geopotential surface height anomaly.  C                      pressure anomaly
80  C     iMin, iMax - Ranges and sub-block indices on which calculations  C                      In p coords phiHydiHyd is the geopotential
81  C     jMin, jMax   are applied.  C                      surface height
82    C                      anomaly.
83    C     etaSurfX,      - Holds surface elevation gradient in X and Y.
84    C     etaSurfY
85    C     KappaRT,       - Total diffusion in vertical for T and S.
86    C     KappaRS          (background + spatially varying, isopycnal term).
87    C     iMin, iMax     - Ranges and sub-block indices on which calculations
88    C     jMin, jMax       are applied.
89  C     bi, bj  C     bi, bj
90  C     k, kUp, kDown, kM1 - Index for layer above and below. kUp and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
91  C                          are switched with layer to be the appropriate index  C     kDown, km1       are switched with layer to be the appropriate
92  C                          into fVerTerm  C                      index into fVerTerm.
93        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RL rVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL rVel    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
99        _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL aTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL aTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL xTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL xTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103        _RL cTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL cTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104        _RL mTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL mTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105        _RL pTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL pTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fMer    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
109        _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
110        _RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
111        _RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
112        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
113        _RL rhokm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114        _RL rhokp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115        _RL rhok  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116        _RL buoyKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL buoyKM1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117        _RL buoyK (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL buoyK   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118        _RL rhotmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhotmp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119        _RL pSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL etaSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120        _RL pSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL etaSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121        _RL K13   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
122        _RL K23   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
123        _RL K33   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
124        _RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
125        _RL KappaZT(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nz)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
126        _RL KappaZS(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nz)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
127          _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
128    
129    C This is currently also used by IVDC and Diagnostics
130    C #ifdef INCLUDE_CONVECT_CALL
131          _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
132    C #endif
133    
134        INTEGER iMin, iMax        INTEGER iMin, iMax
135        INTEGER jMin, jMax        INTEGER jMin, jMax
136        INTEGER bi, bj        INTEGER bi, bj
137        INTEGER i, j        INTEGER i, j
138        INTEGER k, kM1, kUp, kDown        INTEGER k, km1, kup, kDown
139        LOGICAL BOTTOM_LAYER        LOGICAL BOTTOM_LAYER
140    
141    #ifdef ALLOW_AUTODIFF_TAMC
142          INTEGER    isbyte
143          PARAMETER( isbyte = 4 )
144    
145          INTEGER act1, act2, act3, act4
146          INTEGER max1, max2, max3
147          INTEGER iikey, kkey
148          INTEGER maximpl
149    #endif /* ALLOW_AUTODIFF_TAMC */
150    
151  C---    The algorithm...  C---    The algorithm...
152  C  C
153  C       "Correction Step"  C       "Correction Step"
# Line 158  C         (1 + dt * K * d_zz) theta[n] = Line 193  C         (1 + dt * K * d_zz) theta[n] =
193  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
194  C---  C---
195    
196    #ifdef ALLOW_AUTODIFF_TAMC
197    C--   dummy statement to end declaration part
198          ikey = 1
199    #endif /* ALLOW_AUTODIFF_TAMC */
200    
201  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
202  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
203  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
# Line 176  C     uninitialised but inert locations. Line 216  C     uninitialised but inert locations.
216          pTerm(i,j)   = 0. _d 0          pTerm(i,j)   = 0. _d 0
217          fZon(i,j)    = 0. _d 0          fZon(i,j)    = 0. _d 0
218          fMer(i,j)    = 0. _d 0          fMer(i,j)    = 0. _d 0
219          DO K=1,nZ          DO k=1,Nr
220           pH (i,j,k)  = 0. _d 0           phiHyd (i,j,k)  = 0. _d 0
221           K13(i,j,k) = 0. _d 0           KappaRU(i,j,k) = 0. _d 0
222           K23(i,j,k) = 0. _d 0           KappaRV(i,j,k) = 0. _d 0
223           K33(i,j,k) = 0. _d 0           sigmaX(i,j,k) = 0. _d 0
224           KappaZT(i,j,k) = 0. _d 0           sigmaY(i,j,k) = 0. _d 0
225             sigmaR(i,j,k) = 0. _d 0
226          ENDDO          ENDDO
227          rhokm1(i,j)  = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
228          rhok  (i,j)  = 0. _d 0          rhok   (i,j) = 0. _d 0
229          rhokp1(i,j)  = 0. _d 0          rhoKP1 (i,j) = 0. _d 0
230          rhotmp(i,j)  = 0. _d 0          rhoTMP (i,j) = 0. _d 0
231          buoyKM1(i,j) = 0. _d 0          buoyKM1(i,j) = 0. _d 0
232          buoyK  (i,j) = 0. _d 0          buoyK  (i,j) = 0. _d 0
233          maskC (i,j)  = 0. _d 0          maskC  (i,j) = 0. _d 0
234         ENDDO         ENDDO
235        ENDDO        ENDDO
236    
237    
238    #ifdef ALLOW_AUTODIFF_TAMC
239    C--   HPF directive to help TAMC
240    CHPF$ INDEPENDENT
241    #endif /* ALLOW_AUTODIFF_TAMC */
242    
243        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
244    
245    #ifdef ALLOW_AUTODIFF_TAMC
246    C--    HPF directive to help TAMC
247    CHPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV
248    CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA
249    CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV
250    CHPF$&                  )
251    #endif /* ALLOW_AUTODIFF_TAMC */
252    
253         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
254    
255    #ifdef ALLOW_AUTODIFF_TAMC
256              act1 = bi - myBxLo(myThid)
257              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
258    
259              act2 = bj - myByLo(myThid)
260              max2 = myByHi(myThid) - myByLo(myThid) + 1
261    
262              act3 = myThid - 1
263              max3 = nTx*nTy
264    
265              act4 = ikey_dynamics - 1
266    
267              ikey = (act1 + 1) + act2*max1
268         &                      + act3*max1*max2
269         &                      + act4*max1*max2*max3
270    #endif /* ALLOW_AUTODIFF_TAMC */
271    
272  C--     Set up work arrays that need valid initial values  C--     Set up work arrays that need valid initial values
273          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
274           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
275            rTrans(i,j)   = 0. _d 0            rTrans(i,j)   = 0. _d 0
276            rVel  (i,j,1) = 0. _d 0            rVel  (i,j,1) = 0. _d 0
277            rVel  (i,j,2) = 0. _d 0            rVel  (i,j,2) = 0. _d 0
278            fVerT(i,j,1)  = 0. _d 0            fVerT (i,j,1) = 0. _d 0
279            fVerT(i,j,2)  = 0. _d 0            fVerT (i,j,2) = 0. _d 0
280            fVerS(i,j,1)  = 0. _d 0            fVerS (i,j,1) = 0. _d 0
281            fVerS(i,j,2)  = 0. _d 0            fVerS (i,j,2) = 0. _d 0
282            fVerU(i,j,1)  = 0. _d 0            fVerU (i,j,1) = 0. _d 0
283            fVerU(i,j,2)  = 0. _d 0            fVerU (i,j,2) = 0. _d 0
284            fVerV(i,j,1)  = 0. _d 0            fVerV (i,j,1) = 0. _d 0
285            fVerV(i,j,2)  = 0. _d 0            fVerV (i,j,2) = 0. _d 0
286            phiHyd(i,j,1) = 0. _d 0            phiHyd(i,j,1) = 0. _d 0
287            K13(i,j,1)    = 0. _d 0           ENDDO
288            K23(i,j,1)    = 0. _d 0          ENDDO
289            K33(i,j,1)    = 0. _d 0  
290            KapGM(i,j)    = GMkbackground          DO k=1,Nr
291             DO j=1-OLy,sNy+OLy
292              DO i=1-OLx,sNx+OLx
293    #ifdef INCLUDE_CONVECT_CALL
294               ConvectCount(i,j,k) = 0.
295    #endif
296               KappaRT(i,j,k) = 0. _d 0
297               KappaRS(i,j,k) = 0. _d 0
298              ENDDO
299           ENDDO           ENDDO
300          ENDDO          ENDDO
301    
# Line 223  C--     Set up work arrays that need val Line 304  C--     Set up work arrays that need val
304          jMin = 1-OLy+1          jMin = 1-OLy+1
305          jMax = sNy+OLy          jMax = sNy+OLy
306    
         K = 1  
         BOTTOM_LAYER = K .EQ. Nz  
307    
308            k = 1
309            BOTTOM_LAYER = k .EQ. Nr
310    
311    #ifdef DO_PIPELINED_CORRECTION_STEP
312  C--     Calculate gradient of surface pressure  C--     Calculate gradient of surface pressure
313          CALL CALC_GRAD_ETA_SURF(          CALL CALC_GRAD_ETA_SURF(
314       I       bi,bj,iMin,iMax,jMin,jMax,       I       bi,bj,iMin,iMax,jMin,jMax,
315       O       pSurfX,pSurfY,       O       etaSurfX,etaSurfY,
316       I       myThid)       I       myThid)
   
317  C--     Update fields in top level according to tendency terms  C--     Update fields in top level according to tendency terms
318          CALL CORRECTION_STEP(          CALL CORRECTION_STEP(
319       I       bi,bj,iMin,iMax,jMin,jMax,K,pSurfX,pSurfY,myTime,myThid)       I       bi,bj,iMin,iMax,jMin,jMax,k,
320         I       etaSurfX,etaSurfY,myTime,myThid)
321    
322    #ifdef ALLOW_OBCS
323            IF (openBoundaries) THEN
324    #ifdef ALLOW_AUTODIFF_TAMC
325    CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
326    CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
327    CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
328    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
329    #endif /* ALLOW_AUTODIFF_TAMC */
330               CALL APPLY_OBCS1( bi, bj, k, myThid )
331            END IF
332    #endif
333    
334          IF ( .NOT. BOTTOM_LAYER ) THEN          IF ( .NOT. BOTTOM_LAYER ) THEN
335  C--      Update fields in layer below according to tendency terms  C--      Update fields in layer below according to tendency terms
336           CALL CORRECTION_STEP(           CALL CORRECTION_STEP(
337       I        bi,bj,iMin,iMax,jMin,jMax,K+1,pSurfX,pSurfY,myTime,myThid)       I        bi,bj,iMin,iMax,jMin,jMax,k+1,
338         I        etaSurfX,etaSurfY,myTime,myThid)
339    #ifdef ALLOW_OBCS
340             IF (openBoundaries) THEN
341    #ifdef ALLOW_AUTODIFF_TAMC
342    CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
343    CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
344    CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
345    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
346    #endif /* ALLOW_AUTODIFF_TAMC */
347                CALL APPLY_OBCS1( bi, bj, k+1, myThid )
348             END IF
349    #endif
350          ENDIF          ENDIF
351    #endif
352  C--     Density of 1st level (below W(1)) reference to level 1  C--     Density of 1st level (below W(1)) reference to level 1
353    #ifdef  INCLUDE_FIND_RHO_CALL
354    #ifdef ALLOW_AUTODIFF_TAMC
355    CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
356    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
357    #endif /* ALLOW_AUTODIFF_TAMC */
358          CALL FIND_RHO(          CALL FIND_RHO(
359       I     bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,       I     bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
360       O     rhoKm1,       O     rhoKm1,
361       I     myThid )       I     myThid )
362    #endif
363    
364          IF ( .NOT. BOTTOM_LAYER ) THEN          IF (.NOT. BOTTOM_LAYER) THEN
365    
366  C--      Check static stability with layer below  C--      Check static stability with layer below
367  C        and mix as needed.  C--      and mix as needed.
368    #ifdef  INCLUDE_FIND_RHO_CALL
369    #ifdef ALLOW_AUTODIFF_TAMC
370    CADJ STORE theta(:,:,k+1,bi,bj) = comlev1_bibj
371    CADJ &   , key = ikey, byte = isbyte
372    CADJ STORE salt (:,:,k+1,bi,bj) = comlev1_bibj
373    CADJ &   , key = ikey, byte = isbyte
374    #endif /* ALLOW_AUTODIFF_TAMC */
375           CALL FIND_RHO(           CALL FIND_RHO(
376       I      bi, bj, iMin, iMax, jMin, jMax, K+1, K, eosType,       I      bi, bj, iMin, iMax, jMin, jMax, k+1, k, eosType,
377       O      rhoKp1,       O      rhoKp1,
378       I      myThid )       I      myThid )
379    #endif
380    
381    #ifdef  INCLUDE_CONVECT_CALL
382    
383    #ifdef ALLOW_AUTODIFF_TAMC
384    CADJ STORE rhoKm1(:,:)  = comlev1_bibj, key = ikey, byte = isbyte
385    CADJ STORE rhoKp1(:,:)  = comlev1_bibj, key = ikey, byte = isbyte
386    #endif /* ALLOW_AUTODIFF_TAMC */
387           CALL CONVECT(           CALL CONVECT(
388       I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,       I       bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,
389         U       ConvectCount,
390       I       myTime,myIter,myThid)       I       myTime,myIter,myThid)
391    
392    #ifdef ALLOW_AUTODIFF_TAMC
393    CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)
394    CADJ &     = comlev1_bibj, key = ikey, byte = isbyte
395    CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)
396    CADJ &     = comlev1_bibj, key = ikey, byte = isbyte
397    #endif /* ALLOW_AUTODIFF_TAMC */
398    
399    #endif
400    
401    C--      Implicit Vertical Diffusion for Convection
402             IF (ivdc_kappa.NE.0.) THEN
403                CALL CALC_IVDC(
404         I       bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,
405         U       ConvectCount, KappaRT, KappaRS,
406         I       myTime,myIter,myThid)
407             ENDIF
408    
409  C--      Recompute density after mixing  C--      Recompute density after mixing
410    #ifdef  INCLUDE_FIND_RHO_CALL
411           CALL FIND_RHO(           CALL FIND_RHO(
412       I      bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,       I      bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
413       O      rhoKm1,       O      rhoKm1,
414       I      myThid )       I      myThid )
415    #endif
416          ENDIF          ENDIF
417    
418  C--     Calculate buoyancy  C--     Calculate buoyancy
419          CALL CALC_BUOY(          CALL CALC_BUOYANCY(
420       I      bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,       I      bi,bj,iMin,iMax,jMin,jMax,k,rhoKm1,
421       O      buoyKm1,       O      buoyKm1,
422       I      myThid )       I      myThid )
423    
424  C--     Integrate hydrostatic balance for pH with BC of pH(z=0)=0  C--     Integrate hydrostatic balance for phiHyd with BC of
425    C--     phiHyd(z=0)=0
426          CALL CALC_PHI_HYD(          CALL CALC_PHI_HYD(
427       I      bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyKm1,       I      bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyKm1,
428       U      phiHyd,       U      phiHyd,
429       I      myThid )       I      myThid )
430    
431          DO K=2,Nz  #ifdef ALLOW_GMREDI
432            CALL GRAD_SIGMA(
433         I            bi, bj, iMin, iMax, jMin, jMax, k,
434         I            rhoKm1, rhoKm1, rhoKm1,
435         O            sigmaX, sigmaY, sigmaR,
436         I            myThid )
437    #endif
438    
439    C--     Start of downward loop
440            DO k=2,Nr
441    
442           BOTTOM_LAYER = K .EQ. Nz  #ifdef ALLOW_AUTODIFF_TAMC
443             kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1
444    #endif /* ALLOW_AUTODIFF_TAMC */
445    
446             BOTTOM_LAYER = k .EQ. Nr
447    
448    #ifdef DO_PIPELINED_CORRECTION_STEP
449           IF ( .NOT. BOTTOM_LAYER ) THEN           IF ( .NOT. BOTTOM_LAYER ) THEN
450  C--       Update fields in layer below according to tendency terms  C--       Update fields in layer below according to tendency terms
451            CALL CORRECTION_STEP(            CALL CORRECTION_STEP(
452       I         bi,bj,iMin,iMax,jMin,jMax,K+1,pSurfX,pSurfY,myTime,myThid)       I         bi,bj,iMin,iMax,jMin,jMax,k+1,
453         I         etaSurfX,etaSurfY,myTime,myThid)
454    #ifdef ALLOW_OBCS
455              IF (openBoundaries) THEN
456    #ifdef ALLOW_AUTODIFF_TAMC
457    CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj_k
458    CADJ &   , key = kkey, byte = isbyte
459    CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj_k
460    CADJ &   , key = kkey, byte = isbyte
461    CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
462    CADJ &   , key = kkey, byte = isbyte
463    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
464    CADJ &   , key = kkey, byte = isbyte
465    #endif /* ALLOW_AUTODIFF_TAMC */
466                 CALL APPLY_OBCS1( bi, bj, k+1, myThid )
467              END IF
468    #endif
469           ENDIF           ENDIF
470    #endif /* DO_PIPELINED_CORRECTION_STEP */
471    
472  C--      Density of K level (below W(K)) reference to K level  C--      Density of k level (below W(k)) reference to k level
473    #ifdef  INCLUDE_FIND_RHO_CALL
474    #ifdef ALLOW_AUTODIFF_TAMC
475    CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
476    CADJ &   , key = kkey, byte = isbyte
477    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
478    CADJ &   , key = kkey, byte = isbyte
479    #endif /* ALLOW_AUTODIFF_TAMC */
480           CALL FIND_RHO(           CALL FIND_RHO(
481       I      bi, bj, iMin, iMax, jMin, jMax,  K, K, eosType,       I      bi, bj, iMin, iMax, jMin, jMax,  k, k, eosType,
482       O      rhoK,       O      rhoK,
483       I      myThid )       I      myThid )
484    
485           IF ( .NOT. BOTTOM_LAYER ) THEN  #ifdef ALLOW_AUTODIFF_TAMC
486    CADJ STORE rhoK(:,:)  = comlev1_bibj_k, key = kkey, byte = isbyte
487    #endif /* ALLOW_AUTODIFF_TAMC */
488    #endif
489    
490             IF (.NOT. BOTTOM_LAYER) THEN
491    
492  C--       Check static stability with layer below and mix as needed.  C--       Check static stability with layer below and mix as needed.
493  C--       Density of K+1 level (below W(K+1)) reference to K level.  C--       Density of k+1 level (below W(k+1)) reference to k level.
494    #ifdef  INCLUDE_FIND_RHO_CALL
495    #ifdef ALLOW_AUTODIFF_TAMC
496    CADJ STORE theta(:,:,k+1,bi,bj) = comlev1_bibj_k
497    CADJ &   , key = kkey, byte = isbyte
498    CADJ STORE salt (:,:,k+1,bi,bj) = comlev1_bibj_k
499    CADJ &   , key = kkey, byte = isbyte
500    #endif /* ALLOW_AUTODIFF_TAMC */
501            CALL FIND_RHO(            CALL FIND_RHO(
502       I       bi, bj, iMin, iMax, jMin, jMax,  K+1, K, eosType,       I       bi, bj, iMin, iMax, jMin, jMax,  k+1, k, eosType,
503       O       rhoKp1,       O       rhoKp1,
504       I       myThid )       I       myThid )
505    #ifdef ALLOW_AUTODIFF_TAMC
506    CADJ STORE rhoKp1(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
507    #endif /* ALLOW_AUTODIFF_TAMC */
508    #endif
509    
510    #ifdef  INCLUDE_CONVECT_CALL
511            CALL CONVECT(            CALL CONVECT(
512       I        bi,bj,iMin,iMax,jMin,jMax,K+1,rhoK,rhoKp1,       I        bi,bj,iMin,iMax,jMin,jMax,k+1,rhoK,rhoKp1,
513         U        ConvectCount,
514       I        myTime,myIter,myThid)       I        myTime,myIter,myThid)
515    
516    #endif
517    
518    C--      Implicit Vertical Diffusion for Convection
519             IF (ivdc_kappa.NE.0.) THEN
520    #ifdef ALLOW_AUTODIFF_TAMC
521    CADJ STORE rhoKm1(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
522    #endif /* ALLOW_AUTODIFF_TAMC */
523                CALL CALC_IVDC(
524         I       bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,
525         U       ConvectCount, KappaRT, KappaRS,
526         I       myTime,myIter,myThid)
527             END IF
528    
529  C--       Recompute density after mixing  C--       Recompute density after mixing
530    #ifdef  INCLUDE_FIND_RHO_CALL
531    #ifdef ALLOW_AUTODIFF_TAMC
532    CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
533    CADJ &   , key = kkey, byte = isbyte
534    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
535    CADJ &   , key = kkey, byte = isbyte
536    #endif /* ALLOW_AUTODIFF_TAMC */
537            CALL FIND_RHO(            CALL FIND_RHO(
538       I       bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,       I       bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
539       O       rhoK,       O       rhoK,
540       I       myThid )       I       myThid )
541    #endif
542    
543    C--            IF (.NOT. BOTTOM_LAYER) ends here
544           ENDIF           ENDIF
545    
546  C--      Calculate buoyancy  C--      Calculate buoyancy
547           CALL CALC_BUOY(           CALL CALC_BUOYANCY(
548       I       bi,bj,iMin,iMax,jMin,jMax,K,rhoK,       I       bi,bj,iMin,iMax,jMin,jMax,k,rhoK,
549       O       buoyK,       O       buoyK,
550       I       myThid )       I       myThid )
551    
552  C--      Integrate hydrostatic balance for pH with BC of pH(z=0)=0  C--      Integrate hydrostatic balance for phiHyd with BC of
553    C--      phiHyd(z=0)=0
554           CALL CALC_PHI_HYD(           CALL CALC_PHI_HYD(
555       I       bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyK,       I        bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyK,
556       U       phiHyd,       U        phiHyd,
557       I       myThid )       I        myThid )
558    
559    #ifdef  INCLUDE_FIND_RHO_CALL
560  C--      Calculate iso-neutral slopes for the GM/Redi parameterisation  C--      Calculate iso-neutral slopes for the GM/Redi parameterisation
561    
562    #ifdef ALLOW_AUTODIFF_TAMC
563    CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k
564    CADJ &   , key = kkey, byte = isbyte
565    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k
566    CADJ &   , key = kkey, byte = isbyte
567    #endif /* ALLOW_AUTODIFF_TAMC */
568    
569           CALL FIND_RHO(           CALL FIND_RHO(
570       I      bi, bj, iMin, iMax, jMin, jMax, K-1, K, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
571       O      rhoTmp,       O        rhoTmp,
572       I      myThid )       I        myThid )
573           CALL CALC_ISOSLOPES(  #endif
574       I             bi, bj, iMin, iMax, jMin, jMax, K,  
575       I             rhoKm1, rhoK, rhotmp,  
576       O             K13, K23, K33, KapGM,  #ifdef ALLOW_GMREDI
577             CALL GRAD_SIGMA(
578         I             bi, bj, iMin, iMax, jMin, jMax, k,
579         I             rhoK, rhotmp, rhoK,
580         O             sigmaX, sigmaY, sigmaR,
581       I             myThid )       I             myThid )
582    #endif
583    
584           DO J=jMin,jMax           DO J=jMin,jMax
585            DO I=iMin,iMax            DO I=iMin,iMax
586    #ifdef  INCLUDE_FIND_RHO_CALL
587             rhoKm1 (I,J) = rhoK(I,J)             rhoKm1 (I,J) = rhoK(I,J)
588    #endif
589             buoyKm1(I,J) = buoyK(I,J)             buoyKm1(I,J) = buoyK(I,J)
590            ENDDO            ENDDO
591           ENDDO           ENDDO
592    
593          ENDDO ! K  C--     end of k loop
594            ENDDO
595    
596    #ifdef ALLOW_GMREDI
597            IF (useGMRedi) THEN
598              DO k=1, Nr
599                CALL GMREDI_CALC_TENSOR(
600         I             bi, bj, iMin, iMax, jMin, jMax, k,
601         I             sigmaX, sigmaY, sigmaR,
602         I             myThid )
603              ENDDO
604            ENDIF
605    #endif
606    
607    #ifdef ALLOW_AUTODIFF_TAMC
608    CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
609    CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
610    
611    CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
612    CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
613    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
614    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
615    
616    C--     dummy initialization to break data flow because
617    C--     calc_div_ghat has a condition for initialization
618            DO J=jMin,jMax
619               DO I=iMin,iMax
620                  cg2d_b(i,j,bi,bj) = 0.0
621               ENDDO
622            ENDDO
623    #endif /* ALLOW_AUTODIFF_TAMC */
624    
625    
626    #ifdef ALLOW_KPP
627    C--   Compute KPP mixing coefficients
628            IF (useKPP) THEN
629    
630              CALL TIMER_START('KPP_CALC               [DYNAMICS]', myThid)
631              CALL KPP_CALC(
632         I                 bi, bj, myTime, myThid )
633              CALL TIMER_STOP ('KPP_CALC               [DYNAMICS]', myThid)
634    
635    #ifdef ALLOW_AUTODIFF_TAMC
636            ELSE
637              DO j=1-OLy,sNy+OLy
638                DO i=1-OLx,sNx+OLx
639                  KPPhbl (i,j,bi,bj) = 1.0
640                  KPPfrac(i,j,bi,bj) = 0.0
641                  DO k = 1,Nr
642                     KPPghat   (i,j,k,bi,bj) = 0.0
643                     KPPviscAz (i,j,k,bi,bj) = viscAz
644                     KPPdiffKzT(i,j,k,bi,bj) = diffKzT
645                     KPPdiffKzS(i,j,k,bi,bj) = diffKzS
646                  ENDDO
647                ENDDO
648              ENDDO
649    #endif /* ALLOW_AUTODIFF_TAMC */
650            ENDIF
651    
652    #ifdef ALLOW_AUTODIFF_TAMC
653    CADJ STORE KPPghat   (:,:,:,bi,bj)
654    CADJ &   , KPPviscAz (:,:,:,bi,bj)
655    CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
656    CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
657    CADJ &   , KPPfrac   (:,:  ,bi,bj)
658    CADJ &                 = comlev1_bibj, key = ikey, byte = isbyte
659    #endif /* ALLOW_AUTODIFF_TAMC */
660    
661    #endif /* ALLOW_KPP */
662    
663    C--     Start of upward loop
664            DO k = Nr, 1, -1
665    
666    C--      km1    Points to level above k (=k-1)
667    C--      kup    Cycles through 1,2 to point to layer above
668    C--      kDown  Cycles through 2,1 to point to current layer
669    
670             km1  =max(1,k-1)
671             kup  =1+MOD(k+1,2)
672             kDown=1+MOD(k,2)
673    
         DO K = Nz, 1, -1  
          kM1  =max(1,k-1)   ! Points to level above k (=k-1)  
          kUp  =1+MOD(k+1,2) ! Cycles through 1,2 to point to layer above  
          kDown=1+MOD(k,2)   ! Cycles through 2,1 to point to current layer  
674           iMin = 1-OLx+2           iMin = 1-OLx+2
675           iMax = sNx+OLx-1           iMax = sNx+OLx-1
676           jMin = 1-OLy+2           jMin = 1-OLy+2
677           jMax = sNy+OLy-1           jMax = sNy+OLy-1
678    
679    #ifdef ALLOW_AUTODIFF_TAMC
680             kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1
681    
682    CADJ STORE rvel  (:,:,kDown) = comlev1_bibj_k, key = kkey, byte = isbyte
683    CADJ STORE rTrans(:,:)       = comlev1_bibj_k, key = kkey, byte = isbyte
684    CADJ STORE KappaRT(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
685    CADJ STORE KappaRS(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
686    #endif /* ALLOW_AUTODIFF_TAMC */
687    
688  C--      Get temporary terms used by tendency routines  C--      Get temporary terms used by tendency routines
689           CALL CALC_COMMON_FACTORS (           CALL CALC_COMMON_FACTORS (
690       I        bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
691       O        xA,yA,uTrans,vTrans,wTrans,wVel,maskC,maskUp,       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,
692       I        myThid)       I        myThid)
693    
694    #ifdef ALLOW_OBCS
695            IF (openBoundaries) THEN
696             CALL APPLY_OBCS3( bi, bj, k, kup, rTrans, rVel, myThid )
697            ENDIF
698    #endif
699    
700    #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
701  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
702           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
703       I        bi,bj,iMin,iMax,jMin,jMax,K,       I        bi,bj,iMin,iMax,jMin,jMax,k,
704       I        maskC,maskUp,KapGM,K33,       I        maskC,maskUp,
705       O        KappaZT,KappaZS,       O        KappaRT,KappaRS,KappaRU,KappaRV,
706       I        myThid)       I        myThid)
707    #endif
708  C--      Calculate accelerations in the momentum equations  C--      Calculate accelerations in the momentum equations
709           IF ( momStepping ) THEN           IF ( momStepping ) THEN
710            CALL CALC_MOM_RHS(            CALL CALC_MOM_RHS(
711       I         bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
712       I         xA,yA,uTrans,vTrans,wTrans,wVel,maskC,       I         xA,yA,uTrans,vTrans,rTrans,rVel,maskC,
713       I         phiHyd,       I         phiHyd,KappaRU,KappaRV,
714       U         aTerm,xTerm,cTerm,mTerm,pTerm,       U         aTerm,xTerm,cTerm,mTerm,pTerm,
715       U         fZon, fMer, fVerU, fVerV,       U         fZon, fMer, fVerU, fVerV,
716       I         myThid)       I         myTime, myThid)
717    #ifdef ALLOW_AUTODIFF_TAMC
718    #ifdef INCLUDE_CD_CODE
719             ELSE
720                DO j=1-OLy,sNy+OLy
721                   DO i=1-OLx,sNx+OLx
722                      guCD(i,j,k,bi,bj) = 0.0
723                      gvCD(i,j,k,bi,bj) = 0.0
724                   END DO
725                END DO
726    #endif
727    #endif /* ALLOW_AUTODIFF_TAMC */
728           ENDIF           ENDIF
   
729  C--      Calculate active tracer tendencies  C--      Calculate active tracer tendencies
730           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
731            CALL CALC_GT(            CALL CALC_GT(
732       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
733       I         xA,yA,uTrans,vTrans,wTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
734       I         K13,K23,KappaZT,KapGM,       I         KappaRT,
735       U         aTerm,xTerm,fZon,fMer,fVerT,       U         aTerm,xTerm,fZon,fMer,fVerT,
736       I         myThid)       I         myTime, myThid)
737           ENDIF           ENDIF
738           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
739            CALL CALC_GS(            CALL CALC_GS(
740       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
741       I         xA,yA,uTrans,vTrans,wTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
742       I         K13,K23,KappaZS,KapGM,       I         KappaRS,
743       U         aTerm,xTerm,fZon,fMer,fVerS,       U         aTerm,xTerm,fZon,fMer,fVerS,
744       I         myThid)       I         myTime, myThid)
745           ENDIF           ENDIF
746    #ifdef ALLOW_OBCS
747    C--      Calculate future values on open boundaries
748             IF (openBoundaries) THEN
749    Caja      CALL CYCLE_OBCS( k, bi, bj, myThid )
750              CALL SET_OBCS( k, bi, bj, myTime+deltaTclock, myThid )
751             ENDIF
752    #endif
753  C--      Prediction step (step forward all model variables)  C--      Prediction step (step forward all model variables)
754           CALL TIMESTEP(           CALL TIMESTEP(
755       I       bi,bj,iMin,iMax,jMin,jMax,K,       I       bi,bj,iMin,iMax,jMin,jMax,k,
756       I       myThid)       I       myIter, myThid)
757    #ifdef ALLOW_OBCS
758    C--      Apply open boundary conditions
759             IF (openBoundaries) THEN
760    #ifdef ALLOW_AUTODIFF_TAMC
761    CADJ STORE gunm1(:,:,k,bi,bj) = comlev1_bibj_k
762    CADJ &   , key = kkey, byte = isbyte
763    CADJ STORE gvnm1(:,:,k,bi,bj) = comlev1_bibj_k
764    CADJ &   , key = kkey, byte = isbyte
765    CADJ STORE gwnm1(:,:,k,bi,bj) = comlev1_bibj_k, key = kkey, byte = isbyte
766    #endif /* ALLOW_AUTODIFF_TAMC */
767    CADJ &  
768                CALL APPLY_OBCS2( bi, bj, k, myThid )
769             END IF
770    #endif
771    C--      Freeze water
772             IF (allowFreezing) THEN
773    #ifdef ALLOW_AUTODIFF_TAMC
774    CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
775    CADJ &   , key = kkey, byte = isbyte
776    #endif /* ALLOW_AUTODIFF_TAMC */
777                CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
778             END IF
779    
780    #ifdef DIVG_IN_DYNAMICS
781  C--      Diagnose barotropic divergence of predicted fields  C--      Diagnose barotropic divergence of predicted fields
782           CALL DIV_G(           CALL CALC_DIV_GHAT(
783       I       bi,bj,iMin,iMax,jMin,jMax,K,       I       bi,bj,iMin,iMax,jMin,jMax,k,
784       I       xA,yA,       I       xA,yA,
785       I       myThid)       I       myThid)
786    #endif /* DIVG_IN_DYNAMICS */
787    
788  C--      Cumulative diagnostic calculations (ie. time-averaging)  C--      Cumulative diagnostic calculations (ie. time-averaging)
789  #ifdef ALLOW_DIAGNOSTICS  #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE
790           IF (taveFreq.GT.0.) THEN           IF (taveFreq.GT.0.) THEN
791            CALL DO_TIME_AVERAGES(            CALL DO_TIME_AVERAGES(
792       I                           myTime, myIter, bi, bj, K, kUp, kDown,       I                           myTime, myIter, bi, bj, k, kup, kDown,
793       I                           K13, K23, wVel, KapGM,       I                           rVel, ConvectCount,
794       I                           myThid )       I                           myThid )
795           ENDIF           ENDIF
796  #endif  #endif
797    
798          ENDDO ! K  
799    C--     k loop
800            ENDDO
801    
802    #ifdef ALLOW_AUTODIFF_TAMC
803               maximpl = 6
804               iikey = (ikey-1)*maximpl
805    #endif /* ALLOW_AUTODIFF_TAMC */
806    
807  C--     Implicit diffusion  C--     Implicit diffusion
808          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
809           CALL IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,  
810       I                  KappaZT,KappaZS,           IF (tempStepping) THEN
811       I                  myThid )  #ifdef ALLOW_AUTODIFF_TAMC
812                idkey = iikey + 1
813    #endif /* ALLOW_AUTODIFF_TAMC */
814                CALL IMPLDIFF(
815         I         bi, bj, iMin, iMax, jMin, jMax,
816         I         deltaTtracer, KappaRT,recip_HFacC,
817         U         gTNm1,
818         I         myThid )
819             END IF
820    
821             IF (saltStepping) THEN
822    #ifdef ALLOW_AUTODIFF_TAMC
823             idkey = iikey + 2
824    #endif /* ALLOW_AUTODIFF_TAMC */
825                CALL IMPLDIFF(
826         I         bi, bj, iMin, iMax, jMin, jMax,
827         I         deltaTtracer, KappaRS,recip_HFacC,
828         U         gSNm1,
829         I         myThid )
830             END IF
831    
832    C--     implicitDiffusion
833            ENDIF
834    
835    C--     Implicit viscosity
836            IF (implicitViscosity) THEN
837    
838             IF (momStepping) THEN
839    #ifdef ALLOW_AUTODIFF_TAMC
840             idkey = iikey + 3
841    #endif /* ALLOW_AUTODIFF_TAMC */
842              CALL IMPLDIFF(
843         I         bi, bj, iMin, iMax, jMin, jMax,
844         I         deltaTmom, KappaRU,recip_HFacW,
845         U         gUNm1,
846         I         myThid )
847    #ifdef ALLOW_AUTODIFF_TAMC
848             idkey = iikey + 4
849    #endif /* ALLOW_AUTODIFF_TAMC */
850              CALL IMPLDIFF(
851         I         bi, bj, iMin, iMax, jMin, jMax,
852         I         deltaTmom, KappaRV,recip_HFacS,
853         U         gVNm1,
854         I         myThid )
855    
856    #ifdef INCLUDE_CD_CODE
857    
858    #ifdef ALLOW_AUTODIFF_TAMC
859             idkey = iikey + 5
860    #endif /* ALLOW_AUTODIFF_TAMC */
861              CALL IMPLDIFF(
862         I         bi, bj, iMin, iMax, jMin, jMax,
863         I         deltaTmom, KappaRU,recip_HFacW,
864         U         vVelD,
865         I         myThid )
866    #ifdef ALLOW_AUTODIFF_TAMC
867            idkey = iikey + 6
868    #endif /* ALLOW_AUTODIFF_TAMC */
869              CALL IMPLDIFF(
870         I         bi, bj, iMin, iMax, jMin, jMax,
871         I         deltaTmom, KappaRV,recip_HFacS,
872         U         uVelD,
873         I         myThid )
874    
875    #endif
876    
877    C--      momStepping
878             ENDIF
879    
880    C--     implicitViscosity
881          ENDIF          ENDIF
882    
883         ENDDO         ENDDO
884        ENDDO        ENDDO
885    
 C     write(0,*) 'dynamics: pS ',minval(cg2d_x(1:sNx,1:sNy,:,:)),  
 C    &                           maxval(cg2d_x(1:sNx,1:sNy,:,:))  
 C     write(0,*) 'dynamics: U  ',minval(uVel(1:sNx,1:sNy,1,:,:),mask=uVel(1:sNx,1:sNy,1,:,:).NE.0.),  
 C    &                           maxval(uVel(1:sNx,1:sNy,1,:,:),mask=uVel(1:sNx,1:sNy,1,:,:).NE.0.)  
 C     write(0,*) 'dynamics: V  ',minval(vVel(1:sNx,1:sNy,1,:,:),mask=vVel(1:sNx,1:sNy,1,:,:).NE.0.),  
 C    &                           maxval(vVel(1:sNx,1:sNy,1,:,:),mask=vVel(1:sNx,1:sNy,1,:,:).NE.0.)  
 C     write(0,*) 'dynamics: wVel(1) ',  
 C    &            minval(wVel(1:sNx,1:sNy,1),mask=wVel(1:sNx,1:sNy,1).NE.0.),  
 C    &            maxval(wVel(1:sNx,1:sNy,1),mask=wVel(1:sNx,1:sNy,1).NE.0.)  
 C     write(0,*) 'dynamics: wVel(2) ',  
 C    &            minval(wVel(1:sNx,1:sNy,2),mask=wVel(1:sNx,1:sNy,2).NE.0.),  
 C    &            maxval(wVel(1:sNx,1:sNy,2),mask=wVel(1:sNx,1:sNy,2).NE.0.)  
 cblk  write(0,*) 'dynamics: K13',minval(K13(1:sNx,1:sNy,:)),  
 cblk &                           maxval(K13(1:sNx,1:sNy,:))  
 cblk  write(0,*) 'dynamics: K23',minval(K23(1:sNx,1:sNy,:)),  
 cblk &                           maxval(K23(1:sNx,1:sNy,:))  
 cblk  write(0,*) 'dynamics: K33',minval(K33(1:sNx,1:sNy,:)),  
 cblk &                           maxval(K33(1:sNx,1:sNy,:))  
 C     write(0,*) 'dynamics: gT ',minval(gT(1:sNx,1:sNy,:,:,:)),  
 C    &                           maxval(gT(1:sNx,1:sNy,:,:,:))  
 C     write(0,*) 'dynamics: T  ',minval(Theta(1:sNx,1:sNy,:,:,:)),  
 C    &                           maxval(Theta(1:sNx,1:sNy,:,:,:))  
 C     write(0,*) 'dynamics: gS ',minval(gS(1:sNx,1:sNy,:,:,:)),  
 C    &                           maxval(gS(1:sNx,1:sNy,:,:,:))  
 C     write(0,*) 'dynamics: S  ',minval(salt(1:sNx,1:sNy,:,:,:)),  
 C    &                           maxval(salt(1:sNx,1:sNy,:,:,:))  
 C     write(0,*) 'dynamics: pH ',minval(pH/(Gravity*Rhonil),mask=ph.NE.0.),  
 C    &                           maxval(pH/(Gravity*Rhonil))  
886    
887        RETURN        RETURN
888        END        END
Legend:
Removed from v.1.28  
changed lines
  Added in v.1.53
  ViewVC Help
Powered by ViewVC 1.1.22