/[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.49 by heimbach, Fri Jun 9 02:45:04 2000 UTC revision 1.54.2.5 by adcroft, Tue Jan 9 16:21:18 2001 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     \==========================================================/
 c  
 c     changed: Patrick Heimbach heimbach@mit.edu 6-Jun-2000  
 c              - computation of ikey wrong for nTx,nTy > 1  
 c                and/or nsx,nsy > 1: act1 and act2 were  
 c                mixed up.  
   
23        IMPLICIT NONE        IMPLICIT NONE
24    
25  C     == Global variables ===  C     == Global variables ===
# Line 36  C     == Global variables === Line 30  C     == Global variables ===
30  #include "DYNVARS.h"  #include "DYNVARS.h"
31  #include "GRID.h"  #include "GRID.h"
32    
 #ifdef ALLOW_KPP  
 #include "KPPMIX.h"  
 #endif  
   
33  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
34  #include "tamc.h"  # include "tamc.h"
35  #include "tamc_keys.h"  # include "tamc_keys.h"
36    #endif /* ALLOW_AUTODIFF_TAMC */
37    
38    #ifdef ALLOW_KPP
39    # include "KPP.h"
40  #endif  #endif
41    
42  C     == Routine arguments ==  C     == Routine arguments ==
# Line 64  C                              o rVel: Line 58  C                              o rVel:
58  C                                        lower cell faces.  C                                        lower cell faces.
59  C     maskC,maskUp             o maskC: land/water mask for tracer cells  C     maskC,maskUp             o maskC: land/water mask for tracer cells
60  C                              o maskUp: land/water mask for W points  C                              o maskUp: land/water mask for W points
61  C     aTerm, xTerm, cTerm    - Work arrays for holding separate terms in  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
 C     mTerm, pTerm,            tendency equations.  
 C     fZon, fMer, fVer[STUV]   o aTerm: Advection term  
 C                              o xTerm: Mixing term  
 C                              o cTerm: Coriolis term  
 C                              o mTerm: Metric term  
 C                              o pTerm: Pressure term  
 C                              o fZon: Zonal flux term  
 C                              o fMer: Meridional flux term  
 C                              o fVer: Vertical flux term - note fVer  
62  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
63  C                                      so we need an fVer for each  C                                      so we need an fVer for each
64  C                                      variable.  C                                      variable.
# Line 89  C                      surface height Line 74  C                      surface height
74  C                      anomaly.  C                      anomaly.
75  C     etaSurfX,      - Holds surface elevation gradient in X and Y.  C     etaSurfX,      - Holds surface elevation gradient in X and Y.
76  C     etaSurfY  C     etaSurfY
 C     K13, K23, K33  - Non-zero elements of small-angle approximation  
 C                      diffusion tensor.  
 C     KapGM          - Spatially varying Visbeck et. al mixing coeff.  
77  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
78  C     KappaRS          (background + spatially varying, isopycnal term).  C     KappaRS          (background + spatially varying, isopycnal term).
79  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
80  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
81  C     bi, bj  C     bi, bj
82  C     k, kUp,        - Index for layer above and below. kUp and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
83  C     kDown, kM1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
84  C                      index into fVerTerm.  C                      index into fVerTerm.
85        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 108  C                      index into fVerTe Line 90  C                      index into fVerTe
90        _RL rVel    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL rVel    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
91        _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL aTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL xTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL cTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL mTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL pTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL fZon    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL fMer    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
93        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
94        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
95        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
# Line 126  C                      index into fVerTe Line 101  C                      index into fVerTe
101        _RL buoyKM1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL buoyKM1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL buoyK   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL buoyK   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103        _RL rhotmp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhotmp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL etaSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL etaSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL K13     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL K23     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL K33     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL KapGM   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
104        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
105        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
106        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
107        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
108          _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
109          _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
110          _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
111    
112  #ifdef INCLUDE_CONVECT_CALL  C This is currently also used by IVDC and Diagnostics
113    C #ifdef INCLUDE_CONVECT_CALL
114        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
115  #endif  C #endif
116    
117        INTEGER iMin, iMax        INTEGER iMin, iMax
118        INTEGER jMin, jMax        INTEGER jMin, jMax
119        INTEGER bi, bj        INTEGER bi, bj
120        INTEGER i, j        INTEGER i, j
121        INTEGER k, kM1, kUp, kDown        INTEGER k, km1, kup, kDown
       LOGICAL BOTTOM_LAYER  
122    
123  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
124        INTEGER    isbyte        INTEGER    isbyte
# Line 154  C                      index into fVerTe Line 126  C                      index into fVerTe
126    
127        INTEGER act1, act2, act3, act4        INTEGER act1, act2, act3, act4
128        INTEGER max1, max2, max3        INTEGER max1, max2, max3
129        INTEGER ikact, iikey,kkey        INTEGER iikey, kkey
130        INTEGER maximpl        INTEGER maximpl
131  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
132    
133  C---    The algorithm...  C---    The algorithm...
134  C  C
# Line 171  C Line 143  C
143  C       "Calculation of Gs"  C       "Calculation of Gs"
144  C       ===================  C       ===================
145  C       This is where all the accelerations and tendencies (ie.  C       This is where all the accelerations and tendencies (ie.
146  C       phiHydysics, parameterizations etc...) are calculated  C       physics, parameterizations etc...) are calculated
147  C         rVel = sum_r ( div. u[n] )  C         rVel = sum_r ( div. u[n] )
148  C         rho = rho ( theta[n], salt[n] )  C         rho = rho ( theta[n], salt[n] )
149  C         b   = b(rho, theta)  C         b   = b(rho, theta)
# Line 206  C--- Line 178  C---
178  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
179  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
180        ikey = 1        ikey = 1
181  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
182    
183  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
184  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
# Line 219  C     uninitialised but inert locations. Line 191  C     uninitialised but inert locations.
191          yA(i,j)      = 0. _d 0          yA(i,j)      = 0. _d 0
192          uTrans(i,j)  = 0. _d 0          uTrans(i,j)  = 0. _d 0
193          vTrans(i,j)  = 0. _d 0          vTrans(i,j)  = 0. _d 0
194          aTerm(i,j)   = 0. _d 0          DO k=1,Nr
         xTerm(i,j)   = 0. _d 0  
         cTerm(i,j)   = 0. _d 0  
         mTerm(i,j)   = 0. _d 0  
         pTerm(i,j)   = 0. _d 0  
         fZon(i,j)    = 0. _d 0  
         fMer(i,j)    = 0. _d 0  
         DO K=1,Nr  
195           phiHyd (i,j,k)  = 0. _d 0           phiHyd (i,j,k)  = 0. _d 0
          K13(i,j,k)  = 0. _d 0  
          K23(i,j,k)  = 0. _d 0  
          K33(i,j,k)  = 0. _d 0  
196           KappaRU(i,j,k) = 0. _d 0           KappaRU(i,j,k) = 0. _d 0
197           KappaRV(i,j,k) = 0. _d 0           KappaRV(i,j,k) = 0. _d 0
198             sigmaX(i,j,k) = 0. _d 0
199             sigmaY(i,j,k) = 0. _d 0
200             sigmaR(i,j,k) = 0. _d 0
201          ENDDO          ENDDO
202          rhoKM1 (i,j) = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
203          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
# Line 247  C     uninitialised but inert locations. Line 212  C     uninitialised but inert locations.
212    
213  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
214  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
215  !HPF$ INDEPENDENT  CHPF$ INDEPENDENT
216  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
217    
218        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
219    
220  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
221  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
222  !HPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV
223  !HPF$&                  ,phiHyd,K13,K23,K33,KapGM  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA
224  !HPF$&                  ,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV
225  !HPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  )
226  !HPF$&                  )  #endif /* ALLOW_AUTODIFF_TAMC */
 #endif  
227    
228         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
229    
# Line 278  C--    HPF directive to help TAMC Line 242  C--    HPF directive to help TAMC
242            ikey = (act1 + 1) + act2*max1            ikey = (act1 + 1) + act2*max1
243       &                      + act3*max1*max2       &                      + act3*max1*max2
244       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
245  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
246    
247  C--     Set up work arrays that need valid initial values  C--     Set up work arrays that need valid initial values
248          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
# Line 295  C--     Set up work arrays that need val Line 259  C--     Set up work arrays that need val
259            fVerV (i,j,1) = 0. _d 0            fVerV (i,j,1) = 0. _d 0
260            fVerV (i,j,2) = 0. _d 0            fVerV (i,j,2) = 0. _d 0
261            phiHyd(i,j,1) = 0. _d 0            phiHyd(i,j,1) = 0. _d 0
           K13   (i,j,1) = 0. _d 0  
           K23   (i,j,1) = 0. _d 0  
           K33   (i,j,1) = 0. _d 0  
           KapGM (i,j)   = GMkbackground  
262           ENDDO           ENDDO
263          ENDDO          ENDDO
264    
# Line 320  C--     Set up work arrays that need val Line 280  C--     Set up work arrays that need val
280          jMax = sNy+OLy          jMax = sNy+OLy
281    
282    
283          K = 1  C--     Start of diagnostic loop
284          BOTTOM_LAYER = K .EQ. Nr          DO k=Nr,1,-1
   
 #ifdef DO_PIPELINED_CORRECTION_STEP  
 C--     Calculate gradient of surface pressure  
         CALL CALC_GRAD_ETA_SURF(  
      I       bi,bj,iMin,iMax,jMin,jMax,  
      O       etaSurfX,etaSurfY,  
      I       myThid)  
 C--     Update fields in top level according to tendency terms  
         CALL CORRECTION_STEP(  
      I       bi,bj,iMin,iMax,jMin,jMax,K,  
      I       etaSurfX,etaSurfY,myTime,myThid)  
285    
 #ifdef ALLOW_OBCS  
         IF (openBoundaries) THEN  
286  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
287  CADJ STORE uvel (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  C? Patrick, is this formula correct now that we change the loop range?
288  CADJ STORE vvel (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  C? Do we still need this?
289  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte           kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1
290  CADJ STORE salt(:,:,k,bi,bj)   = comlev1_2d, key = ikey, byte = isbyte  #endif /* ALLOW_AUTODIFF_TAMC */
 #endif  
            CALL APPLY_OBCS1( bi, bj, K, myThid )  
         END IF  
 #endif  
291    
         IF ( .NOT. BOTTOM_LAYER ) THEN  
 C--      Update fields in layer below according to tendency terms  
          CALL CORRECTION_STEP(  
      I        bi,bj,iMin,iMax,jMin,jMax,K+1,  
      I        etaSurfX,etaSurfY,myTime,myThid)  
292  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
293           IF (openBoundaries) THEN  C--       Calculate future values on open boundaries
294  #ifdef ALLOW_AUTODIFF_TAMC            IF (openBoundaries) THEN
295  CADJ STORE uvel (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  Caja        CALL CYCLE_OBCS( k, bi, bj, myThid )
296  CADJ STORE vvel (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  c new args! CALL SET_OBCS( k, bi, bj, myTime, myThid )
297  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  c +deltaT?
298  CADJ STORE salt(:,:,k,bi,bj)   = comlev1_2d, key = ikey, byte = isbyte            ENDIF
 #endif  
             CALL APPLY_OBCS1( bi, bj, K+1, myThid )  
          END IF  
 #endif  
         ENDIF  
 #endif  
 C--     Density of 1st level (below W(1)) reference to level 1  
 #ifdef  INCLUDE_FIND_RHO_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
         CALL FIND_RHO(  
      I     bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,  
      O     rhoKm1,  
      I     myThid )  
299  #endif  #endif
300    
301          IF (       (.NOT. BOTTOM_LAYER)  C--       Integrate continuity vertically for vertical velocity
302  #ifdef ALLOW_KPP            CALL INTEGRATE_FOR_W(
303       &       .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing       I                         bi, bj, k, uVel, vVel,
304  #endif       O                         wVel,
305       &     ) THEN       I                         myThid )
306  C--      Check static stability with layer below  #ifdef ALLOW_OBCS
307  C--      and mix as needed.            IF (openBoundaries) THEN
308  #ifdef  INCLUDE_FIND_RHO_CALL  c new subr  CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
309  #ifdef ALLOW_AUTODIFF_TAMC            ENDIF
 CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
          CALL FIND_RHO(  
      I      bi, bj, iMin, iMax, jMin, jMax, K+1, K, eosType,  
      O      rhoKp1,  
      I      myThid )  
310  #endif  #endif
311    
312  #ifdef  INCLUDE_CONVECT_CALL  C--       Calculate gradients of potential density for isoneutral
313    C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
314  #ifdef ALLOW_AUTODIFF_TAMC            IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
315  CADJ STORE rhoKm1(:,:)  = comlev1_2d, key = ikey, byte = isbyte              CALL FIND_RHO(
316  CADJ STORE rhoKp1(:,:)  = comlev1_2d, key = ikey, byte = isbyte       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
317  #endif       O        rhoK,
318           CALL CONVECT(       I        myThid )
319       I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,              CALL FIND_RHO(
320       U       ConvectCount,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
321       I       myTime,myIter,myThid)       O        rhoKm1,
322  #ifdef ALLOW_AUTODIFF_TAMC       I        myThid )
323  CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)              CALL GRAD_SIGMA(
324  CADJ &     = comlev1_2d, key = ikey, byte = isbyte       I             bi, bj, iMin, iMax, jMin, jMax, k,
325  CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)       I             rhoK, rhoKm1, rhoK,
326  CADJ &     = comlev1_2d, key = ikey, byte = isbyte       O             sigmaX, sigmaY, sigmaR,
327  #endif       I             myThid )
328              ENDIF
329    
330  #endif  C--       Implicit Vertical Diffusion for Convection
331              IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
332                CALL CALC_IVDC(
333         I        bi, bj, iMin, iMax, jMin, jMax, k,
334         I        rhoKm1, rhoK,
335    c should use sigmaR !!!
336         U        ConvectCount, KappaRT, KappaRS,
337         I        myTime, myIter, myThid)
338              END IF
339    
340  C--      Implicit Vertical Diffusion for Convection  C--     end of diagnostic k loop (Nr:1)
341           IF (ivdc_kappa.NE.0.) CALL CALC_IVDC(          ENDDO
      I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,  
      U       ConvectCount, KappaRT, KappaRS,  
      I       myTime,myIter,myThid)  
 CRG: do we need do store STORE KappaRT, KappaRS ?  
342    
343  C--      Recompute density after mixing  #ifdef  ALLOW_GMREDI
344  #ifdef  INCLUDE_FIND_RHO_CALL  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
345           CALL FIND_RHO(          IF (useGMRedi) THEN
346       I      bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,            DO k=1,Nr
347       O      rhoKm1,              CALL GMREDI_CALC_TENSOR(
348       I      myThid )       I             bi, bj, iMin, iMax, jMin, jMax, k,
349  #endif       I             sigmaX, sigmaY, sigmaR,
350         I             myThid )
351              ENDDO
352          ENDIF          ENDIF
353  C--     Calculate buoyancy  #endif  /* ALLOW_GMREDI */
         CALL CALC_BUOYANCY(  
      I      bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,  
      O      buoyKm1,  
      I      myThid )  
 C--     Integrate hydrostatic balance for phiHyd with BC of  
 C--     phiHyd(z=0)=0  
         CALL CALC_PHI_HYD(  
      I      bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyKm1,  
      U      phiHyd,  
      I      myThid )  
   
 C----------------------------------------------  
 C--     start of downward loop  
 C----------------------------------------------  
         DO K=2,Nr  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
          kkey = ikact*(Nr-2+1) + (k-2) + 1  
 #endif  
354    
355           BOTTOM_LAYER = K .EQ. Nr  #ifdef  ALLOW_KPP
356    C--     Compute KPP mixing coefficients
357            IF (useKPP) THEN
358              CALL KPP_CALC(
359         I                  bi, bj, myTime, myThid )
360            ENDIF
361    #endif  /* ALLOW_KPP */
362    
363  #ifdef DO_PIPELINED_CORRECTION_STEP  C--     Determines forcing terms based on external fields
364           IF ( .NOT. BOTTOM_LAYER ) THEN  C       relaxation terms, etc.
365  C--       Update fields in layer below according to tendency terms          CALL EXTERNAL_FORCING_SURF(
366            CALL CORRECTION_STEP(       I             bi, bj, iMin, iMax, jMin, jMax,
367       I         bi,bj,iMin,iMax,jMin,jMax,K+1,       I             myThid )
      I         etaSurfX,etaSurfY,myTime,myThid)  
 #ifdef ALLOW_OBCS  
           IF (openBoundaries) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE uvel (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE vvel (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE salt(:,:,k,bi,bj)   = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
              CALL APPLY_OBCS1( bi, bj, K+1, myThid )  
           END IF  
 #endif  
          ENDIF  
 #endif  
368    
 C--      Density of K level (below W(K)) reference to K level  
 #ifdef  INCLUDE_FIND_RHO_CALL  
369  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
370  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
371  CADJ STORE salt (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
372  #endif  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
373           CALL FIND_RHO(  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
374       I      bi, bj, iMin, iMax, jMin, jMax,  K, K, eosType,  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
375       O      rhoK,  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
376       I      myThid )  #endif /* ALLOW_AUTODIFF_TAMC */
 #endif  
          IF (       (.NOT. BOTTOM_LAYER)  
 #ifdef ALLOW_KPP  
      &       .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing  
 #endif  
      &      ) THEN  
 C--       Check static stability with layer below and mix as needed.  
 C--       Density of K+1 level (below W(K+1)) reference to K level.  
 #ifdef  INCLUDE_FIND_RHO_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
           CALL FIND_RHO(  
      I       bi, bj, iMin, iMax, jMin, jMax,  K+1, K, eosType,  
      O       rhoKp1,  
      I       myThid )  
 #endif  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhok  (:,:)   = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE rhoKm1(:,:)   = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE rhoKp1(:,:)   = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
   
 #ifdef  INCLUDE_CONVECT_CALL  
           CALL CONVECT(  
      I        bi,bj,iMin,iMax,jMin,jMax,K+1,rhoK,rhoKp1,  
      U        ConvectCount,  
      I        myTime,myIter,myThid)  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)  
 CADJ &     = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)  
 CADJ &     = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
 #endif  
   
 C--      Implicit Vertical Diffusion for Convection  
          IF (ivdc_kappa.NE.0.) THEN  
             CALL CALC_IVDC(  
      I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,  
      U       ConvectCount, KappaRT, KappaRS,  
      I       myTime,myIter,myThid)  
 CRG: do we need do store STORE KappaRT, KappaRS ?  
          END IF  
377    
 C--       Recompute density after mixing  
 #ifdef  INCLUDE_FIND_RHO_CALL  
           CALL FIND_RHO(  
      I       bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,  
      O       rhoK,  
      I       myThid )  
 #endif  
          ENDIF  
 C--      Calculate buoyancy  
          CALL CALC_BUOYANCY(  
      I       bi,bj,iMin,iMax,jMin,jMax,K,rhoK,  
      O       buoyK,  
      I       myThid )  
 C--      Integrate hydrostatic balance for phiHyd with BC of  
 C--      phiHyd(z=0)=0  
          CALL CALC_PHI_HYD(  
      I        bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyK,  
      U        phiHyd,  
      I        myThid )  
 C--      Calculate iso-neutral slopes for the GM/Redi parameterisation  
 #ifdef  INCLUDE_FIND_RHO_CALL  
          CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, K-1, K, eosType,  
      O        rhoTmp,  
      I        myThid )  
 #endif  
378    
 #ifdef  INCLUDE_CALC_ISOSLOPES_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhoTmp(:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE rhok  (:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE rhoKm1(:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE kapgm (:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
          CALL CALC_ISOSLOPES(  
      I        bi, bj, iMin, iMax, jMin, jMax, K,  
      I        rhoKm1, rhoK, rhotmp,  
      O        K13, K23, K33, KapGM,  
      I        myThid )  
 #endif  
379    
380           DO J=jMin,jMax  C--     Start of thermodynamics loop
381            DO I=iMin,iMax          DO k=Nr,1,-1
 #ifdef  INCLUDE_FIND_RHO_CALL  
            rhoKm1 (I,J) = rhoK(I,J)  
 #endif  
            buoyKm1(I,J) = buoyK(I,J)  
           ENDDO  
          ENDDO  
         ENDDO  
 C--     end of k loop  
382    
383  #ifdef ALLOW_AUTODIFF_TAMC  C--       km1    Points to level above k (=k-1)
384  CADJ STORE theta(:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  C--       kup    Cycles through 1,2 to point to layer above
385  CADJ STORE salt (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  C--       kDown  Cycles through 2,1 to point to current layer
 CADJ STORE uvel (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
386    
387  #ifdef ALLOW_KPP            km1  = MAX(1,k-1)
388  C----------------------------------------------            kup  = 1+MOD(k+1,2)
389  C--     Compute KPP mixing coefficients            kDown= 1+MOD(k,2)
 C----------------------------------------------  
         IF (usingKPPmixing) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE fu  (:,:  ,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE fv  (:,:  ,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
          CALL TIMER_START('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]'  
      I          , myThid)  
          CALL KVMIX(  
      I               bi, bj, myTime, myThid )  
          CALL TIMER_STOP ('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]'  
      I        , myThid)  
         ENDIF  
 #endif  
390    
391  C----------------------------------------------            iMin = 1-OLx+2
392  C--     start of upward loop            iMax = sNx+OLx-1
393  C----------------------------------------------            jMin = 1-OLy+2
394          DO K = Nr, 1, -1            jMax = sNy+OLy-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  
   
          iMin = 1-OLx+2  
          iMax = sNx+OLx-1  
          jMin = 1-OLy+2  
          jMax = sNy+OLy-1  
395    
396  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
397           kkey = ikact*(Nr-1+1) + (k-1) + 1  CPatrick Is this formula correct?
398  #endif           kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1
399    CADJ STORE rvel  (:,:,kDown) = comlev1_bibj_k, key = kkey, byte = isbyte
400  #ifdef ALLOW_AUTODIFF_TAMC  CADJ STORE rTrans(:,:)       = comlev1_bibj_k, key = kkey, byte = isbyte
401  CADJ STORE rvel  (:,:,kDown)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRT(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
402  CADJ STORE rTrans(:,:)        = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRS(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
403  CADJ STORE KappaRT(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  #endif /* ALLOW_AUTODIFF_TAMC */
 CADJ STORE KappaRS(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
404    
405  C--      Get temporary terms used by tendency routines  C--      Get temporary terms used by tendency routines
406           CALL CALC_COMMON_FACTORS (           CALL CALC_COMMON_FACTORS (
407       I        bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
408       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,
409       I        myThid)       I        myThid)
410    
 #ifdef ALLOW_OBCS  
         IF (openBoundaries) THEN  
          CALL APPLY_OBCS3( bi, bj, K, Kup, rTrans, rVel, myThid )  
         ENDIF  
 #endif  
   
411  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
412  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
413           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
414       I        bi,bj,iMin,iMax,jMin,jMax,K,       I        bi,bj,iMin,iMax,jMin,jMax,k,
415       I        maskC,maskUp,KapGM,K33,       I        maskC,maskup,
416       O        KappaRT,KappaRS,KappaRU,KappaRV,       O        KappaRT,KappaRS,KappaRU,KappaRV,
417       I        myThid)       I        myThid)
418  #endif  #endif
419  C--      Calculate accelerations in the momentum equations  
420           IF ( momStepping ) THEN  C--      Calculate active tracer tendencies (gT,gS,...)
421            CALL CALC_MOM_RHS(  C        and step forward storing result in gTnm1, gSnm1, etc.
      I         bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,rVel,maskC,  
      I         phiHyd,KappaRU,KappaRV,  
      U         aTerm,xTerm,cTerm,mTerm,pTerm,  
      U         fZon, fMer, fVerU, fVerV,  
      I         myTime, myThid)  
 #ifdef ALLOW_AUTODIFF_TAMC  
 #ifdef INCLUDE_CD_CODE  
          ELSE  
             DO j=1-OLy,sNy+OLy  
                DO i=1-OLx,sNx+OLx  
                   guCD(i,j,k,bi,bj) = 0.0  
                   gvCD(i,j,k,bi,bj) = 0.0  
                END DO  
             END DO  
 #endif  
 #endif  
          ENDIF  
 C--      Calculate active tracer tendencies  
422           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
423            CALL CALC_GT(             CALL CALC_GT(
424       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
425       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
426       I         K13,K23,KappaRT,KapGM,       I         KappaRT,
427       U         aTerm,xTerm,fZon,fMer,fVerT,       U         fVerT,
428       I         myTime, myThid)       I         myTime, myThid)
429               CALL TIMESTEP_TRACER(
430         I         bi,bj,iMin,iMax,jMin,jMax,k,
431         I         theta, gT,
432         U         gTnm1,
433         I         myIter, myThid)
434           ENDIF           ENDIF
435           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
436            CALL CALC_GS(             CALL CALC_GS(
437       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
438       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
439       I         K13,K23,KappaRS,KapGM,       I         KappaRS,
440       U         aTerm,xTerm,fZon,fMer,fVerS,       U         fVerS,
441       I         myTime, myThid)       I         myTime, myThid)
442               CALL TIMESTEP_TRACER(
443         I         bi,bj,iMin,iMax,jMin,jMax,k,
444         I         salt, gS,
445         U         gSnm1,
446         I         myIter, myThid)
447           ENDIF           ENDIF
448  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
 C--      Calculate future values on open boundaries  
          IF (openBoundaries) THEN  
 Caja      CALL CYCLE_OBCS( K, bi, bj, myThid )  
           CALL SET_OBCS( K, bi, bj, myTime+deltaTclock, myThid )  
          ENDIF  
 #endif  
 C--      Prediction step (step forward all model variables)  
          CALL TIMESTEP(  
      I       bi,bj,iMin,iMax,jMin,jMax,K,  
      I       myIter, myThid)  
 #ifdef ALLOW_OBCS  
449  C--      Apply open boundary conditions  C--      Apply open boundary conditions
450           IF (openBoundaries) THEN           IF (openBoundaries) THEN
451  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
452  CADJ STORE gunm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE gwnm1(:,:,k,bi,bj) = comlev1_bibj_k
453  CADJ STORE gvnm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
454  CADJ STORE gwnm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  #endif /* ALLOW_AUTODIFF_TAMC */
455  #endif  c new subr  CALL OBCS_APPLY_TS( bi, bj, k, myThid )
             CALL APPLY_OBCS2( bi, bj, K, myThid )  
456           END IF           END IF
457  #endif  #endif
458  C--      Freeze water  C--      Freeze water
459           IF (allowFreezing) THEN           IF (allowFreezing) THEN
460  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
461  CADJ STORE gTNm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
462  #endif  CADJ &   , key = kkey, byte = isbyte
463              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, K, myThid )  #endif /* ALLOW_AUTODIFF_TAMC */
464                CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
465           END IF           END IF
466    
467  #ifdef DIVG_IN_DYNAMICS  C--     end of thermodynamic k loop (Nr:1)
468  C--      Diagnose barotropic divergence of predicted fields          ENDDO
          CALL CALC_DIV_GHAT(  
      I       bi,bj,iMin,iMax,jMin,jMax,K,  
      I       xA,yA,  
      I       myThid)  
 #endif /* DIVG_IN_DYNAMICS */  
   
 C--      Cumulative diagnostic calculations (ie. time-averaging)  
 #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE  
          IF (taveFreq.GT.0.) THEN  
           CALL DO_TIME_AVERAGES(  
      I                           myTime, myIter, bi, bj, K, kUp, kDown,  
      I                           K13, K23, rVel, KapGM, ConvectCount,  
      I                           myThid )  
          ENDIF  
 #endif  
469    
470    
471          ENDDO ! K  #ifdef ALLOW_AUTODIFF_TAMC
472    CPatrick? What about this one?
473               maximpl = 6
474               iikey = (ikey-1)*maximpl
475    #endif /* ALLOW_AUTODIFF_TAMC */
476    
477  C--     Implicit diffusion  C--     Implicit diffusion
478          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
479    
480  #ifdef ALLOW_AUTODIFF_TAMC            IF (tempStepping) THEN
            maximpl = 6  
            iikey = ikact*maximpl  
 #endif  
   
          IF (tempStepping) THEN  
481  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
482              idkey = iikey + 1              idkey = iikey + 1
483  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
484              CALL IMPLDIFF(              CALL IMPLDIFF(
485       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
486       I         deltaTtracer, KappaRT,recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
487       U         gTNm1,       U         gTNm1,
488       I         myThid )       I         myThid )
489           END IF           ENDIF
490    
491           IF (saltStepping) THEN           IF (saltStepping) THEN
492  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
493           idkey = iikey + 2           idkey = iikey + 2
494  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
495              CALL IMPLDIFF(              CALL IMPLDIFF(
496       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
497       I         deltaTtracer, KappaRS,recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
498       U         gSNm1,       U         gSNm1,
499       I         myThid )       I         myThid )
500           END IF           ENDIF
501    
502          ENDIF ! implicitDiffusion  C--     End If implicitDiffusion
503            ENDIF
504    
 C--     Implicit viscosity  
         IF (implicitViscosity) THEN  
505    
506           IF (momStepping) THEN  
507  #ifdef ALLOW_AUTODIFF_TAMC  C--     Start of dynamics loop
508           idkey = iikey + 3          DO k=1,Nr
509  #endif  
510    C--       km1    Points to level above k (=k-1)
511    C--       kup    Cycles through 1,2 to point to layer above
512    C--       kDown  Cycles through 2,1 to point to current layer
513    
514              km1  = MAX(1,k-1)
515              kup  = 1+MOD(k+1,2)
516              kDown= 1+MOD(k,2)
517    
518              iMin = 1-OLx+2
519              iMax = sNx+OLx-1
520              jMin = 1-OLy+2
521              jMax = sNy+OLy-1
522    
523    C--      Calculate buoyancy
524             CALL FIND_RHO(
525         I        bi, bj, iMin, iMax, jMin, jMax, km1, km1, eosType,
526         O        rhoKm1,
527         I        myThid )
528             CALL CALC_BUOYANCY(
529         I        bi,bj,iMin,iMax,jMin,jMax,k,rhoKm1,
530         O        buoyKm1,
531         I        myThid )
532             CALL FIND_RHO(
533         I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
534         O        rhoK,
535         I        myThid )
536             CALL CALC_BUOYANCY(
537         I        bi,bj,iMin,iMax,jMin,jMax,k,rhoK,
538         O        buoyK,
539         I        myThid )
540    
541    C--      Integrate hydrostatic balance for phiHyd with BC of
542    C--      phiHyd(z=0)=0
543             CALL CALC_PHI_HYD(
544         I        bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyK,
545         U        phiHyd,
546         I        myThid )
547    
548    C--      Calculate accelerations in the momentum equations (gU, gV, ...)
549    C        and step forward storing the result in gUnm1, gVnm1, etc...
550             IF ( momStepping ) THEN
551               CALL CALC_MOM_RHS(
552         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
553         I         phiHyd,KappaRU,KappaRV,
554         U         fVerU, fVerV,
555         I         myTime, myThid)
556               CALL TIMESTEP(
557         I         bi,bj,iMin,iMax,jMin,jMax,k,
558         I         myIter, myThid)
559    #ifdef   ALLOW_AUTODIFF_TAMC
560    #ifdef   INCLUDE_CD_CODE
561             ELSE
562               DO j=1-OLy,sNy+OLy
563                 DO i=1-OLx,sNx+OLx
564                   guCD(i,j,k,bi,bj) = 0.0
565                   gvCD(i,j,k,bi,bj) = 0.0
566                 END DO
567               END DO
568    #endif   /* INCLUDE_CD_CODE */
569    #endif   /* ALLOW_AUTODIFF_TAMC */
570             ENDIF
571    
572    
573    C--     end of dynamics k loop (1:Nr)
574            ENDDO
575    
576    
577    
578    C--     Implicit viscosity
579            IF (implicitViscosity.AND.momStepping) THEN
580    #ifdef    ALLOW_AUTODIFF_TAMC
581              idkey = iikey + 3
582    #endif    /* ALLOW_AUTODIFF_TAMC */
583            CALL IMPLDIFF(            CALL IMPLDIFF(
584       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
585       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
586       U         gUNm1,       U         gUNm1,
587       I         myThid )       I         myThid )
588  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
589           idkey = iikey + 4            idkey = iikey + 4
590  #endif  #endif    /* ALLOW_AUTODIFF_TAMC */
591            CALL IMPLDIFF(            CALL IMPLDIFF(
592       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
593       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
594       U         gVNm1,       U         gVNm1,
595       I         myThid )       I         myThid )
596    
597  #ifdef INCLUDE_CD_CODE  #ifdef    INCLUDE_CD_CODE
598    #ifdef    ALLOW_AUTODIFF_TAMC
599  #ifdef ALLOW_AUTODIFF_TAMC            idkey = iikey + 5
600           idkey = iikey + 5  #endif    /* ALLOW_AUTODIFF_TAMC */
 #endif  
601            CALL IMPLDIFF(            CALL IMPLDIFF(
602       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
603       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
604       U         vVelD,       U         vVelD,
605       I         myThid )       I         myThid )
606  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
607          idkey = iikey + 6            idkey = iikey + 6
608  #endif  #endif    /* ALLOW_AUTODIFF_TAMC */
609            CALL IMPLDIFF(            CALL IMPLDIFF(
610       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
611       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
612       U         uVelD,       U         uVelD,
613       I         myThid )       I         myThid )
614    #endif    /* INCLUDE_CD_CODE */
615  #endif  C--     End If implicitViscosity.AND.momStepping
616            ENDIF
          ENDIF ! momStepping  
         ENDIF ! implicitViscosity  
617    
618         ENDDO         ENDDO
619        ENDDO        ENDDO
620    
 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: rVel(1) ',  
 C    &            minval(rVel(1:sNx,1:sNy,1),mask=rVel(1:sNx,1:sNy,1).NE.0.),  
 C    &            maxval(rVel(1:sNx,1:sNy,1),mask=rVel(1:sNx,1:sNy,1).NE.0.)  
 C     write(0,*) 'dynamics: rVel(2) ',  
 C    &            minval(rVel(1:sNx,1:sNy,2),mask=rVel(1:sNx,1:sNy,2).NE.0.),  
 C    &            maxval(rVel(1:sNx,1:sNy,2),mask=rVel(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: phiHyd ',minval(phiHyd/(Gravity*Rhonil),mask=phiHyd.NE.0.),  
 C    &                           maxval(phiHyd/(Gravity*Rhonil))  
 C     CALL PLOT_FIELD_XYZRL( gU, ' GU exiting dyanmics ' ,  
 C    &Nr, 1, myThid )  
 C     CALL PLOT_FIELD_XYZRL( gV, ' GV exiting dyanmics ' ,  
 C    &Nr, 1, myThid )  
 C     CALL PLOT_FIELD_XYZRL( gS, ' GS exiting dyanmics ' ,  
 C    &Nr, 1, myThid )  
 C     CALL PLOT_FIELD_XYZRL( gT, ' GT exiting dyanmics ' ,  
 C    &Nr, 1, myThid )  
 C     CALL PLOT_FIELD_XYZRL( phiHyd, ' phiHyd exiting dyanmics ' ,  
 C    &Nr, 1, myThid )  
   
   
621        RETURN        RETURN
622        END        END
623    
624    
625    C--      Cumulative diagnostic calculations (ie. time-averaging)
626    #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE
627    c        IF (taveFreq.GT.0.) THEN
628    c         CALL DO_TIME_AVERAGES(
629    c    I                           myTime, myIter, bi, bj, k, kup, kDown,
630    c    I                           ConvectCount,
631    c    I                           myThid )
632    c        ENDIF
633    #endif
634    
Legend:
Removed from v.1.49  
changed lines
  Added in v.1.54.2.5
  ViewVC Help
Powered by ViewVC 1.1.22