/[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.7 by cnh, Mon May 25 16:17:36 1998 UTC revision 1.53 by heimbach, Mon Sep 11 23:07:29 2000 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    
3  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
4    
5        SUBROUTINE DYNAMICS(myThid)        SUBROUTINE DYNAMICS(myTime, myIter, myThid)
6  C     /==========================================================\  C     /==========================================================\
7  C     | SUBROUTINE DYNAMICS                                      |  C     | SUBROUTINE DYNAMICS                                      |
8  C     | o Controlling routine for the explicit part of the model |  C     | o Controlling routine for the explicit part of the model |
# 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
43    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.
45          _RL myTime
46          INTEGER myIter
47        INTEGER myThid        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                              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 rVel:   Vertical velocity at upper and
57    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 53  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     iMin, iMax - Ranges and sub-block indices on which calculations  C     rhoK, rhoKM1   - Density at current level, level above and level
74  C     jMin, jMax   are applied.  C                      below.
75    C     rhoKP1                                                                  
76    C     buoyK, buoyKM1 - Buoyancy at current level and level above.
77    C     phiHyd         - Hydrostatic part of the potential phiHydi.
78    C                      In z coords phiHydiHyd is the hydrostatic
79    C                      pressure anomaly
80    C                      In p coords phiHydiHyd is the geopotential
81    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 wTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rVel    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
99        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RL aTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL xTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL aTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL cTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL xTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103        _RL mTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL cTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104        _RL pTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL mTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL pTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fMer    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108        _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
109        _RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
110        _RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
111        _RL pH    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
112        _RL rhokm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
113        _RL rhokp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114        _RL pSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115        _RL pSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116        _RL K13   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)        _RL buoyKM1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117        _RL K23   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)        _RL buoyK   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118        _RL K33   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)        _RL rhotmp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119        _RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL etaSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120          _RL etaSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121          _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
122          _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
123          _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
124          _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
125          _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
126          _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
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...
152    C
153    C       "Correction Step"
154    C       =================
155    C       Here we update the horizontal velocities with the surface
156    C       pressure such that the resulting flow is either consistent
157    C       with the free-surface evolution or the rigid-lid:
158    C         U[n] = U* + dt x d/dx P
159    C         V[n] = V* + dt x d/dy P
160    C
161    C       "Calculation of Gs"
162    C       ===================
163    C       This is where all the accelerations and tendencies (ie.
164    C       physics, parameterizations etc...) are calculated
165    C         rVel = sum_r ( div. u[n] )
166    C         rho = rho ( theta[n], salt[n] )
167    C         b   = b(rho, theta)
168    C         K31 = K31 ( rho )
169    C         Gu[n] = Gu( u[n], v[n], rVel, b, ... )
170    C         Gv[n] = Gv( u[n], v[n], rVel, b, ... )
171    C         Gt[n] = Gt( theta[n], u[n], v[n], rVel, K31, ... )
172    C         Gs[n] = Gs( salt[n], u[n], v[n], rVel, K31, ... )
173    C
174    C       "Time-stepping" or "Prediction"
175    C       ================================
176    C       The models variables are stepped forward with the appropriate
177    C       time-stepping scheme (currently we use Adams-Bashforth II)
178    C       - For momentum, the result is always *only* a "prediction"
179    C       in that the flow may be divergent and will be "corrected"
180    C       later with a surface pressure gradient.
181    C       - Normally for tracers the result is the new field at time
182    C       level [n+1} *BUT* in the case of implicit diffusion the result
183    C       is also *only* a prediction.
184    C       - We denote "predictors" with an asterisk (*).
185    C         U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
186    C         V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
187    C         theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
188    C         salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
189    C       With implicit diffusion:
190    C         theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
191    C         salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
192    C         (1 + dt * K * d_zz) theta[n] = theta*
193    C         (1 + dt * K * d_zz) salt[n] = salt*
194    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
# Line 110  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             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          rhokp1(i,j)  = 0. _d 0          rhok   (i,j) = 0. _d 0
229            rhoKP1 (i,j) = 0. _d 0
230            rhoTMP (i,j) = 0. _d 0
231            buoyKM1(i,j) = 0. _d 0
232            buoyK  (i,j) = 0. _d 0
233            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  C--     Boundary condition on hydrostatic pressure is pH(z=0)=0  #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
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            pH(i,j,1) = 0. _d 0            rTrans(i,j)   = 0. _d 0
276            K13(i,j,1) = 0. _d 0            rVel  (i,j,1) = 0. _d 0
277            K23(i,j,1) = 0. _d 0            rVel  (i,j,2) = 0. _d 0
278            K33(i,j,1) = 0. _d 0            fVerT (i,j,1) = 0. _d 0
279            KapGM(i,j) = 0. _d 0            fVerT (i,j,2) = 0. _d 0
280              fVerS (i,j,1) = 0. _d 0
281              fVerS (i,j,2) = 0. _d 0
282              fVerU (i,j,1) = 0. _d 0
283              fVerU (i,j,2) = 0. _d 0
284              fVerV (i,j,1) = 0. _d 0
285              fVerV (i,j,2) = 0. _d 0
286              phiHyd(i,j,1) = 0. _d 0
287           ENDDO           ENDDO
288          ENDDO          ENDDO
289    
290  C--     Set up work arrays that need valid initial values          DO k=1,Nr
291          DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
292           DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
293            wTrans(i,j)  = 0. _d 0  #ifdef INCLUDE_CONVECT_CALL
294            fVerT(i,j,1) = 0. _d 0             ConvectCount(i,j,k) = 0.
295            fVerT(i,j,2) = 0. _d 0  #endif
296            fVerS(i,j,1) = 0. _d 0             KappaRT(i,j,k) = 0. _d 0
297            fVerS(i,j,2) = 0. _d 0             KappaRS(i,j,k) = 0. _d 0
298            fVerU(i,j,1) = 0. _d 0            ENDDO
           fVerU(i,j,2) = 0. _d 0  
           fVerV(i,j,1) = 0. _d 0  
           fVerV(i,j,2) = 0. _d 0  
299           ENDDO           ENDDO
300          ENDDO          ENDDO
301    
# Line 155  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    
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 GRAD_PSURF(          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 TIMESTEP(          CALL CORRECTION_STEP(
319       I       bi,bj,iMin,iMax,jMin,jMax,1,pSurfX,pSurfY,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
335    C--      Update fields in layer below according to tendency terms
336             CALL CORRECTION_STEP(
337         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
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, 1, 1, 'LINEAR',       I     bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
360       O     rhoKm1,       O     rhoKm1,
361       I     myThid )       I     myThid )
362  C--     Integrate hydrostatic balance for pH with BC of pH(z=0)=0  #endif
363          CALL CALC_PH(  
364       I      bi,bj,iMin,iMax,jMin,jMax,1,rhoKm1,rhoKm1,          IF (.NOT. BOTTOM_LAYER) THEN
365       U      pH,  
366    C--      Check static stability with layer below
367    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(
376         I      bi, bj, iMin, iMax, jMin, jMax, k+1, k, eosType,
377         O      rhoKp1,
378       I      myThid )       I      myThid )
379    #endif
380    
381          DO K=2,Nz  #ifdef  INCLUDE_CONVECT_CALL
382  C--     Update fields in Kth level according to tendency terms  
383          CALL TIMESTEP(  #ifdef ALLOW_AUTODIFF_TAMC
384       I       bi,bj,iMin,iMax,jMin,jMax,K,pSurfX,pSurfY,myThid)  CADJ STORE rhoKm1(:,:)  = comlev1_bibj, key = ikey, byte = isbyte
385  C--     Density of K-1 level (above W(K)) reference to K level  CADJ STORE rhoKp1(:,:)  = comlev1_bibj, key = ikey, byte = isbyte
386          CALL FIND_RHO(  #endif /* ALLOW_AUTODIFF_TAMC */
387       I     bi, bj, iMin, iMax, jMin, jMax,  K-1, K, 'LINEAR',           CALL CONVECT(
388       O     rhoKm1,       I       bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,
389       I     myThid )       U       ConvectCount,
390  C--     Density of K level (below W(K)) reference to K level       I       myTime,myIter,myThid)
391          CALL FIND_RHO(  
392       I     bi, bj, iMin, iMax, jMin, jMax,  K, K, 'LINEAR',  #ifdef ALLOW_AUTODIFF_TAMC
393       O     rhoKp1,  CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)
394       I     myThid )  CADJ &     = comlev1_bibj, key = ikey, byte = isbyte
395  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)
396          CALL CALC_ISOSLOPES(  CADJ &     = comlev1_bibj, key = ikey, byte = isbyte
397       I            bi, bj, iMin, iMax, jMin, jMax, K,  #endif /* ALLOW_AUTODIFF_TAMC */
398       I            rhoKm1, rhoKp1,  
399       O            K13, K23, K33, KapGM,  #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
410    #ifdef  INCLUDE_FIND_RHO_CALL
411             CALL FIND_RHO(
412         I      bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
413         O      rhoKm1,
414         I      myThid )
415    #endif
416            ENDIF
417    
418    C--     Calculate buoyancy
419            CALL CALC_BUOYANCY(
420         I      bi,bj,iMin,iMax,jMin,jMax,k,rhoKm1,
421         O      buoyKm1,
422         I      myThid )
423    
424    C--     Integrate hydrostatic balance for phiHyd with BC of
425    C--     phiHyd(z=0)=0
426            CALL CALC_PHI_HYD(
427         I      bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyKm1,
428         U      phiHyd,
429         I      myThid )
430    
431    #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 )       I            myThid )
437  C--     Calculate static stability and mix where convectively unstable  #endif
438          CALL CONVECT(  
439       I      bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,rhoKp1,myThid)  C--     Start of downward loop
440  C--     Density of K-1 level (above W(K)) reference to K-1 level          DO k=2,Nr
441          CALL FIND_RHO(  
442       I     bi, bj, iMin, iMax, jMin, jMax,  K-1, K-1, 'LINEAR',  #ifdef ALLOW_AUTODIFF_TAMC
443       O     rhoKm1,           kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1
444       I     myThid )  #endif /* ALLOW_AUTODIFF_TAMC */
445  C--     Density of K level (below W(K)) referenced to K level  
446          CALL FIND_RHO(           BOTTOM_LAYER = k .EQ. Nr
447       I     bi, bj, iMin, iMax, jMin, jMax,  K, K, 'LINEAR',  
448       O     rhoKp1,  #ifdef DO_PIPELINED_CORRECTION_STEP
449       I     myThid )           IF ( .NOT. BOTTOM_LAYER ) THEN
450  C--     Integrate hydrostatic balance for pH with BC of pH(z=0)=0  C--       Update fields in layer below according to tendency terms
451          CALL CALC_PH(            CALL CORRECTION_STEP(
452       I      bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,rhoKp1,       I         bi,bj,iMin,iMax,jMin,jMax,k+1,
453       U      pH,       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
470    #endif /* DO_PIPELINED_CORRECTION_STEP */
471    
472    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(
481         I      bi, bj, iMin, iMax, jMin, jMax,  k, k, eosType,
482         O      rhoK,
483       I      myThid )       I      myThid )
484    
485    #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.
493    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(
502         I       bi, bj, iMin, iMax, jMin, jMax,  k+1, k, eosType,
503         O       rhoKp1,
504         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(
512         I        bi,bj,iMin,iMax,jMin,jMax,k+1,rhoK,rhoKp1,
513         U        ConvectCount,
514         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
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(
538         I       bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
539         O       rhoK,
540         I       myThid )
541    #endif
542    
543    C--            IF (.NOT. BOTTOM_LAYER) ends here
544             ENDIF
545    
546    C--      Calculate buoyancy
547             CALL CALC_BUOYANCY(
548         I       bi,bj,iMin,iMax,jMin,jMax,k,rhoK,
549         O       buoyK,
550         I       myThid )
551    
552    C--      Integrate hydrostatic balance for phiHyd with BC of
553    C--      phiHyd(z=0)=0
554             CALL CALC_PHI_HYD(
555         I        bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyK,
556         U        phiHyd,
557         I        myThid )
558    
559    #ifdef  INCLUDE_FIND_RHO_CALL
560    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(
570         I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
571         O        rhoTmp,
572         I        myThid )
573    #endif
574    
575    
576    #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 )
582    #endif
583    
584             DO J=jMin,jMax
585              DO I=iMin,iMax
586    #ifdef  INCLUDE_FIND_RHO_CALL
587               rhoKm1 (I,J) = rhoK(I,J)
588    #endif
589               buoyKm1(I,J) = buoyK(I,J)
590              ENDDO
591             ENDDO
592    
593    C--     end of k loop
594          ENDDO          ENDDO
595    
596          DO K = Nz, 1, -1  #ifdef ALLOW_GMREDI
597           kM1  =max(1,k-1)   ! Points to level above k (=k-1)          IF (useGMRedi) THEN
598           kUp  =1+MOD(k+1,2) ! Cycles through 1,2 to point to layer above            DO k=1, Nr
599           kDown=1+MOD(k,2)   ! Cycles through 2,1 to point to current layer              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    
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,maskC,maskUp,       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,
692       I        myThid)       I        myThid)
693    
694  C--      Calculate accelerations in the momentum equations  #ifdef ALLOW_OBCS
695           CALL CALC_MOM_RHS(          IF (openBoundaries) THEN
696       I        bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,           CALL APPLY_OBCS3( bi, bj, k, kup, rTrans, rVel, myThid )
697       I        xA,yA,uTrans,vTrans,wTrans,maskC,          ENDIF
698       I        pH,  #endif
699       U        aTerm,xTerm,cTerm,mTerm,pTerm,  
700       U        fZon, fMer, fVerU, fVerV,  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
701    C--      Calculate the total vertical diffusivity
702             CALL CALC_DIFFUSIVITY(
703         I        bi,bj,iMin,iMax,jMin,jMax,k,
704         I        maskC,maskUp,
705         O        KappaRT,KappaRS,KappaRU,KappaRV,
706       I        myThid)       I        myThid)
707    #endif
708    C--      Calculate accelerations in the momentum equations
709             IF ( momStepping ) THEN
710              CALL CALC_MOM_RHS(
711         I         bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
712         I         xA,yA,uTrans,vTrans,rTrans,rVel,maskC,
713         I         phiHyd,KappaRU,KappaRV,
714         U         aTerm,xTerm,cTerm,mTerm,pTerm,
715         U         fZon, fMer, fVerU, fVerV,
716         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
729  C--      Calculate active tracer tendencies  C--      Calculate active tracer tendencies
730           CALL CALC_GT(           IF ( tempStepping ) THEN
731       I        bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,            CALL CALC_GT(
732       I        xA,yA,uTrans,vTrans,wTrans,maskUp,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
733       I        K13,K23,K33,KapGM,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
734       U        aTerm,xTerm,fZon,fMer,fVerT,       I         KappaRT,
735       I        myThid)       U         aTerm,xTerm,fZon,fMer,fVerT,
736  Cdbg     CALL CALC_GS(       I         myTime, myThid)
737  Cdbg I        bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,           ENDIF
738  Cdbg I        xA,yA,uTrans,vTrans,wTrans,maskUp,           IF ( saltStepping ) THEN
739  Cdbg I        K13,K23,K33,KapGM,            CALL CALC_GS(
740  Cdbg U        aTerm,xTerm,fZon,fMer,fVerS,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
741  Cdbg I        myThid)       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
742         I         KappaRS,
743         U         aTerm,xTerm,fZon,fMer,fVerS,
744         I         myTime, myThid)
745             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)
754             CALL TIMESTEP(
755         I       bi,bj,iMin,iMax,jMin,jMax,k,
756         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
782             CALL CALC_DIV_GHAT(
783         I       bi,bj,iMin,iMax,jMin,jMax,k,
784         I       xA,yA,
785         I       myThid)
786    #endif /* DIVG_IN_DYNAMICS */
787    
788    C--      Cumulative diagnostic calculations (ie. time-averaging)
789    #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE
790             IF (taveFreq.GT.0.) THEN
791              CALL DO_TIME_AVERAGES(
792         I                           myTime, myIter, bi, bj, k, kup, kDown,
793         I                           rVel, ConvectCount,
794         I                           myThid )
795             ENDIF
796    #endif
797    
798    
799    C--     k loop
800          ENDDO          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
808            IF (implicitDiffusion) THEN
809    
810             IF (tempStepping) THEN
811    #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
882    
883         ENDDO         ENDDO
884        ENDDO        ENDDO
885    
 !dbg  write(0,*) 'dynamics: pS',minval(cg2d_x),maxval(cg2d_x)  
 !dbg  write(0,*) 'dynamics: U',minval(uVel(1:sNx,1:sNy,:,:,:)),  
 !dbg &                         maxval(uVel(1:sNx,1:sNy,:,:,:))  
 !dbg  write(0,*) 'dynamics: V',minval(vVel(1:sNx,1:sNy,:,:,:)),  
 !dbg &                         maxval(vVel(1:sNx,1:sNy,:,:,:))  
 !dbg  write(0,*) 'dynamics: gT',minval(gT(1:sNx,1:sNy,:,:,:)),  
 !dbg &                         maxval(gT(1:sNx,1:sNy,:,:,:))  
 !dbg  write(0,*) 'dynamics: T',minval(Theta(1:sNx,1:sNy,:,:,:)),  
 !dbg &                         maxval(Theta(1:sNx,1:sNy,:,:,:))  
 !dbg  write(0,*) 'dynamics: pH',minval(pH/(Gravity*Rhonil)),  
 !dbg &                          maxval(pH/(Gravity*Rhonil))  
886    
887        RETURN        RETURN
888        END        END
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.53
  ViewVC Help
Powered by ViewVC 1.1.22