/[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.7 by adcroft, Tue Jan 9 21:26:07 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
285    
 #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)  
   
 #ifdef ALLOW_OBCS  
         IF (openBoundaries) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE uvel (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 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 APPLY_OBCS1( bi, bj, K, myThid )  
         END IF  
 #endif  
   
         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)  
 #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 */
291  #endif  
292              CALL APPLY_OBCS1( bi, bj, K+1, myThid )  C--       Integrate continuity vertically for vertical velocity
293           END IF            CALL INTEGRATE_FOR_W(
294  #endif       I                         bi, bj, k, uVel, vVel,
295          ENDIF       O                         wVel,
296  #endif       I                         myThid )
 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 )  
 #endif  
297    
298          IF (       (.NOT. BOTTOM_LAYER)  #ifdef    ALLOW_OBCS
299  #ifdef ALLOW_KPP  C--       Calculate future values on open boundaries
300       &       .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing            IF (openBoundaries) THEN
301  #endif  #ifdef      ALLOW_NONHYDROSTATIC
302       &     ) THEN              IF (nonHydrostatic) THEN
303  C--      Check static stability with layer below                CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
304  C--      and mix as needed.              ENDIF
305  #ifdef  INCLUDE_FIND_RHO_CALL  #endif      /* ALLOW_NONHYDROSTATIC */
306  #ifdef ALLOW_AUTODIFF_TAMC              CALL OBCS_CALC( bi, bj, k, myTime+deltaT, myThid )
307  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte            ENDIF
308  CADJ STORE salt (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  #endif    /* ALLOW_OBCS */
309  #endif  
310           CALL FIND_RHO(  C--       Calculate gradients of potential density for isoneutral
311       I      bi, bj, iMin, iMax, jMin, jMax, K+1, K, eosType,  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
312       O      rhoKp1,            IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
313       I      myThid )              CALL FIND_RHO(
314  #endif       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
315         O        rhoK,
316  #ifdef  INCLUDE_CONVECT_CALL       I        myThid )
317                CALL FIND_RHO(
318  #ifdef ALLOW_AUTODIFF_TAMC       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
319  CADJ STORE rhoKm1(:,:)  = comlev1_2d, key = ikey, byte = isbyte       O        rhoKm1,
320  CADJ STORE rhoKp1(:,:)  = comlev1_2d, key = ikey, byte = isbyte       I        myThid )
321  #endif              CALL GRAD_SIGMA(
322           CALL CONVECT(       I             bi, bj, iMin, iMax, jMin, jMax, k,
323       I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,       I             rhoK, rhoKm1, rhoK,
324       U       ConvectCount,       O             sigmaX, sigmaY, sigmaR,
325       I       myTime,myIter,myThid)       I             myThid )
326  #ifdef ALLOW_AUTODIFF_TAMC            ENDIF
 CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)  
 CADJ &     = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)  
 CADJ &     = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
327    
328  #endif  C--       Implicit Vertical Diffusion for Convection
329              IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
330                CALL CALC_IVDC(
331         I        bi, bj, iMin, iMax, jMin, jMax, k,
332         I        rhoKm1, rhoK,
333    c should use sigmaR !!!
334         U        ConvectCount, KappaRT, KappaRS,
335         I        myTime, myIter, myThid)
336              END IF
337    
338  C--      Implicit Vertical Diffusion for Convection  C--     end of diagnostic k loop (Nr:1)
339           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 ?  
340    
341  C--      Recompute density after mixing  C--     Determines forcing terms based on external fields
342  #ifdef  INCLUDE_FIND_RHO_CALL  C       relaxation terms, etc.
343           CALL FIND_RHO(          CALL EXTERNAL_FORCING_SURF(
344       I      bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,       I             bi, bj, iMin, iMax, jMin, jMax,
345       O      rhoKm1,       I             myThid )
346       I      myThid )  
347  #endif  #ifdef  ALLOW_GMREDI
348    C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
349            IF (useGMRedi) THEN
350              DO k=1,Nr
351                CALL GMREDI_CALC_TENSOR(
352         I             bi, bj, iMin, iMax, jMin, jMax, k,
353         I             sigmaX, sigmaY, sigmaR,
354         I             myThid )
355              ENDDO
356          ENDIF          ENDIF
357  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  
   
          BOTTOM_LAYER = K .EQ. Nr  
   
 #ifdef DO_PIPELINED_CORRECTION_STEP  
          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)  
 #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  
358    
359  C--      Density of K level (below W(K)) reference to K level  #ifdef  ALLOW_KPP
360  #ifdef  INCLUDE_FIND_RHO_CALL  C--     Compute KPP mixing coefficients
361  #ifdef ALLOW_AUTODIFF_TAMC          IF (useKPP) THEN
362  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte            CALL KPP_CALC(
363  CADJ STORE salt (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte       I                  bi, bj, myTime, myThid )
364  #endif          ENDIF
365           CALL FIND_RHO(  #endif  /* ALLOW_KPP */
      I      bi, bj, iMin, iMax, jMin, jMax,  K, K, eosType,  
      O      rhoK,  
      I      myThid )  
 #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  
366    
367  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
368  CADJ STORE rhok  (:,:)   = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
369  CADJ STORE rhoKm1(:,:)   = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
370  CADJ STORE rhoKp1(:,:)   = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
371  #endif  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
372    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
373    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
374    #endif /* ALLOW_AUTODIFF_TAMC */
375    
 #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  
376    
 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    
378  C--       Recompute density after mixing  C--     Start of thermodynamics loop
379  #ifdef  INCLUDE_FIND_RHO_CALL          DO k=Nr,1,-1
           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  
   
 #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  
   
          DO J=jMin,jMax  
           DO I=iMin,iMax  
 #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  
380    
381  #ifdef ALLOW_AUTODIFF_TAMC  C--       km1    Points to level above k (=k-1)
382  CADJ STORE theta(:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  C--       kup    Cycles through 1,2 to point to layer above
383  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  
384    
385  #ifdef ALLOW_KPP            km1  = MAX(1,k-1)
386  C----------------------------------------------            kup  = 1+MOD(k+1,2)
387  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  
388    
389  C----------------------------------------------            iMin = 1-OLx+2
390  C--     start of upward loop            iMax = sNx+OLx-1
391  C----------------------------------------------            jMin = 1-OLy+2
392          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  
393    
394  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
395           kkey = ikact*(Nr-1+1) + (k-1) + 1  CPatrick Is this formula correct?
396  #endif           kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1
397    CADJ STORE rvel  (:,:,kDown) = comlev1_bibj_k, key = kkey, byte = isbyte
398  #ifdef ALLOW_AUTODIFF_TAMC  CADJ STORE rTrans(:,:)       = comlev1_bibj_k, key = kkey, byte = isbyte
399  CADJ STORE rvel  (:,:,kDown)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRT(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
400  CADJ STORE rTrans(:,:)        = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRS(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
401  CADJ STORE KappaRT(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  #endif /* ALLOW_AUTODIFF_TAMC */
 CADJ STORE KappaRS(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
402    
403  C--      Get temporary terms used by tendency routines  C--      Get temporary terms used by tendency routines
404           CALL CALC_COMMON_FACTORS (           CALL CALC_COMMON_FACTORS (
405       I        bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
406       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,
407       I        myThid)       I        myThid)
408    
 #ifdef ALLOW_OBCS  
         IF (openBoundaries) THEN  
          CALL APPLY_OBCS3( bi, bj, K, Kup, rTrans, rVel, myThid )  
         ENDIF  
 #endif  
   
409  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
410  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
411           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
412       I        bi,bj,iMin,iMax,jMin,jMax,K,       I        bi,bj,iMin,iMax,jMin,jMax,k,
413       I        maskC,maskUp,KapGM,K33,       I        maskC,maskup,
414       O        KappaRT,KappaRS,KappaRU,KappaRV,       O        KappaRT,KappaRS,KappaRU,KappaRV,
415       I        myThid)       I        myThid)
416  #endif  #endif
417  C--      Calculate accelerations in the momentum equations  
418           IF ( momStepping ) THEN  C--      Calculate active tracer tendencies (gT,gS,...)
419            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  
420           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
421            CALL CALC_GT(             CALL CALC_GT(
422       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
423       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
424       I         K13,K23,KappaRT,KapGM,       I         KappaRT,
425       U         aTerm,xTerm,fZon,fMer,fVerT,       U         fVerT,
426       I         myTime, myThid)       I         myTime, myThid)
427               CALL TIMESTEP_TRACER(
428         I         bi,bj,iMin,iMax,jMin,jMax,k,
429         I         theta, gT,
430         U         gTnm1,
431         I         myIter, myThid)
432           ENDIF           ENDIF
433           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
434            CALL CALC_GS(             CALL CALC_GS(
435       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
436       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
437       I         K13,K23,KappaRS,KapGM,       I         KappaRS,
438       U         aTerm,xTerm,fZon,fMer,fVerS,       U         fVerS,
439       I         myTime, myThid)       I         myTime, myThid)
440               CALL TIMESTEP_TRACER(
441         I         bi,bj,iMin,iMax,jMin,jMax,k,
442         I         salt, gS,
443         U         gSnm1,
444         I         myIter, myThid)
445           ENDIF           ENDIF
446  #ifdef ALLOW_OBCS  
447  C--      Calculate future values on open boundaries  #ifdef   ALLOW_OBCS
          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  
448  C--      Apply open boundary conditions  C--      Apply open boundary conditions
449           IF (openBoundaries) THEN           IF (openBoundaries) THEN
450  #ifdef ALLOW_AUTODIFF_TAMC             CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
 CADJ STORE gunm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE gvnm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE gwnm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
             CALL APPLY_OBCS2( bi, bj, K, myThid )  
451           END IF           END IF
452  #endif  #endif   /* ALLOW_OBCS */
453    
454  C--      Freeze water  C--      Freeze water
455           IF (allowFreezing) THEN           IF (allowFreezing) THEN
456  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
457  CADJ STORE gTNm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
458  #endif  CADJ &   , key = kkey, byte = isbyte
459              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, K, myThid )  #endif /* ALLOW_AUTODIFF_TAMC */
460                CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
461           END IF           END IF
462    
463  #ifdef DIVG_IN_DYNAMICS  C--     end of thermodynamic k loop (Nr:1)
464  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  
465    
466    
467          ENDDO ! K  #ifdef ALLOW_AUTODIFF_TAMC
468    CPatrick? What about this one?
469               maximpl = 6
470               iikey = (ikey-1)*maximpl
471    #endif /* ALLOW_AUTODIFF_TAMC */
472    
473  C--     Implicit diffusion  C--     Implicit diffusion
474          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
475    
476  #ifdef ALLOW_AUTODIFF_TAMC            IF (tempStepping) THEN
            maximpl = 6  
            iikey = ikact*maximpl  
 #endif  
   
          IF (tempStepping) THEN  
477  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
478              idkey = iikey + 1              idkey = iikey + 1
479  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
480              CALL IMPLDIFF(              CALL IMPLDIFF(
481       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
482       I         deltaTtracer, KappaRT,recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
483       U         gTNm1,       U         gTNm1,
484       I         myThid )       I         myThid )
485           END IF           ENDIF
486    
487           IF (saltStepping) THEN           IF (saltStepping) THEN
488  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
489           idkey = iikey + 2           idkey = iikey + 2
490  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
491              CALL IMPLDIFF(              CALL IMPLDIFF(
492       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
493       I         deltaTtracer, KappaRS,recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
494       U         gSNm1,       U         gSNm1,
495       I         myThid )       I         myThid )
496             ENDIF
497    
498    #ifdef   ALLOW_OBCS
499    C--      Apply open boundary conditions
500             IF (openBoundaries) THEN
501               DO K=1,Nr
502                 CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
503               ENDDO
504           END IF           END IF
505    #endif   /* ALLOW_OBCS */
506    
507          ENDIF ! implicitDiffusion  C--     End If implicitDiffusion
508            ENDIF
509    
 C--     Implicit viscosity  
         IF (implicitViscosity) THEN  
510    
511           IF (momStepping) THEN  
512  #ifdef ALLOW_AUTODIFF_TAMC  C--     Start of dynamics loop
513           idkey = iikey + 3          DO k=1,Nr
514  #endif  
515    C--       km1    Points to level above k (=k-1)
516    C--       kup    Cycles through 1,2 to point to layer above
517    C--       kDown  Cycles through 2,1 to point to current layer
518    
519              km1  = MAX(1,k-1)
520              kup  = 1+MOD(k+1,2)
521              kDown= 1+MOD(k,2)
522    
523              iMin = 1-OLx+2
524              iMax = sNx+OLx-1
525              jMin = 1-OLy+2
526              jMax = sNy+OLy-1
527    
528    C--      Calculate buoyancy
529             CALL FIND_RHO(
530         I        bi, bj, iMin, iMax, jMin, jMax, km1, km1, eosType,
531         O        rhoKm1,
532         I        myThid )
533             CALL CALC_BUOYANCY(
534         I        bi,bj,iMin,iMax,jMin,jMax,k,rhoKm1,
535         O        buoyKm1,
536         I        myThid )
537             CALL FIND_RHO(
538         I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
539         O        rhoK,
540         I        myThid )
541             CALL CALC_BUOYANCY(
542         I        bi,bj,iMin,iMax,jMin,jMax,k,rhoK,
543         O        buoyK,
544         I        myThid )
545    
546    C--      Integrate hydrostatic balance for phiHyd with BC of
547    C--      phiHyd(z=0)=0
548             CALL CALC_PHI_HYD(
549         I        bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyK,
550         U        phiHyd,
551         I        myThid )
552    
553    C--      Calculate accelerations in the momentum equations (gU, gV, ...)
554    C        and step forward storing the result in gUnm1, gVnm1, etc...
555             IF ( momStepping ) THEN
556               CALL CALC_MOM_RHS(
557         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
558         I         phiHyd,KappaRU,KappaRV,
559         U         fVerU, fVerV,
560         I         myTime, myThid)
561               CALL TIMESTEP(
562         I         bi,bj,iMin,iMax,jMin,jMax,k,
563         I         myIter, myThid)
564    
565    #ifdef   ALLOW_OBCS
566    C--      Apply open boundary conditions
567             IF (openBoundaries) THEN
568               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
569             END IF
570    #endif   /* ALLOW_OBCS */
571    
572    #ifdef   ALLOW_AUTODIFF_TAMC
573    #ifdef   INCLUDE_CD_CODE
574             ELSE
575               DO j=1-OLy,sNy+OLy
576                 DO i=1-OLx,sNx+OLx
577                   guCD(i,j,k,bi,bj) = 0.0
578                   gvCD(i,j,k,bi,bj) = 0.0
579                 END DO
580               END DO
581    #endif   /* INCLUDE_CD_CODE */
582    #endif   /* ALLOW_AUTODIFF_TAMC */
583             ENDIF
584    
585    
586    C--     end of dynamics k loop (1:Nr)
587            ENDDO
588    
589    
590    
591    C--     Implicit viscosity
592            IF (implicitViscosity.AND.momStepping) THEN
593    #ifdef    ALLOW_AUTODIFF_TAMC
594              idkey = iikey + 3
595    #endif    /* ALLOW_AUTODIFF_TAMC */
596            CALL IMPLDIFF(            CALL IMPLDIFF(
597       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
598       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
599       U         gUNm1,       U         gUNm1,
600       I         myThid )       I         myThid )
601  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
602           idkey = iikey + 4            idkey = iikey + 4
603  #endif  #endif    /* ALLOW_AUTODIFF_TAMC */
604            CALL IMPLDIFF(            CALL IMPLDIFF(
605       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
606       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
607       U         gVNm1,       U         gVNm1,
608       I         myThid )       I         myThid )
609    
610  #ifdef INCLUDE_CD_CODE  #ifdef   ALLOW_OBCS
611    C--      Apply open boundary conditions
612             IF (openBoundaries) THEN
613               DO K=1,Nr
614                 CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
615               ENDDO
616             END IF
617    #endif   /* ALLOW_OBCS */
618    
619  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    INCLUDE_CD_CODE
620           idkey = iikey + 5  #ifdef    ALLOW_AUTODIFF_TAMC
621  #endif            idkey = iikey + 5
622    #endif    /* ALLOW_AUTODIFF_TAMC */
623            CALL IMPLDIFF(            CALL IMPLDIFF(
624       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
625       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
626       U         vVelD,       U         vVelD,
627       I         myThid )       I         myThid )
628  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
629          idkey = iikey + 6            idkey = iikey + 6
630  #endif  #endif    /* ALLOW_AUTODIFF_TAMC */
631            CALL IMPLDIFF(            CALL IMPLDIFF(
632       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
633       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
634       U         uVelD,       U         uVelD,
635       I         myThid )       I         myThid )
636    #endif    /* INCLUDE_CD_CODE */
637  #endif  C--     End If implicitViscosity.AND.momStepping
638            ENDIF
          ENDIF ! momStepping  
         ENDIF ! implicitViscosity  
639    
640         ENDDO         ENDDO
641        ENDDO        ENDDO
642    
 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 )  
   
   
643        RETURN        RETURN
644        END        END
645    
646    
647    C--      Cumulative diagnostic calculations (ie. time-averaging)
648    #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE
649    c        IF (taveFreq.GT.0.) THEN
650    c         CALL DO_TIME_AVERAGES(
651    c    I                           myTime, myIter, bi, bj, k, kup, kDown,
652    c    I                           ConvectCount,
653    c    I                           myThid )
654    c        ENDIF
655    #endif
656    
Legend:
Removed from v.1.49  
changed lines
  Added in v.1.54.2.7
  ViewVC Help
Powered by ViewVC 1.1.22