/[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.11 by adcroft, Thu Jan 25 19:43:32 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.
65  C     rhoK, rhoKM1   - Density at current level, level above and level  C     rhoK, rhoKM1   - Density at current level, and level above
 C                      below.  
 C     rhoKP1                                                                    
 C     buoyK, buoyKM1 - Buoyancy at current level and level above.  
66  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     phiHyd         - Hydrostatic part of the potential phiHydi.
67  C                      In z coords phiHydiHyd is the hydrostatic  C                      In z coords phiHydiHyd is the hydrostatic
68  C                      pressure anomaly  C                      pressure anomaly
# Line 89  C                      surface height Line 71  C                      surface height
71  C                      anomaly.  C                      anomaly.
72  C     etaSurfX,      - Holds surface elevation gradient in X and Y.  C     etaSurfX,      - Holds surface elevation gradient in X and Y.
73  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.  
74  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
75  C     KappaRS          (background + spatially varying, isopycnal term).  C     KappaRS          (background + spatially varying, isopycnal term).
76  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
77  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
78  C     bi, bj  C     bi, bj
79  C     k, kUp,        - Index for layer above and below. kUp and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
80  C     kDown, kM1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
81  C                      index into fVerTerm.  C                      index into fVerTerm.
82        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83        _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 87  C                      index into fVerTe
87        _RL rVel    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL rVel    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
88        _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89        _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)  
90        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
91        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
92        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
93        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
94        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
95        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rhokp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
96        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL buoyKM1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL buoyK   (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)  
97        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
98        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
99        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
100        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
101          _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
102          _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103          _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
104    
105  #ifdef INCLUDE_CONVECT_CALL  C This is currently also used by IVDC and Diagnostics
106    C #ifdef INCLUDE_CONVECT_CALL
107        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
108  #endif  C #endif
109    
110        INTEGER iMin, iMax        INTEGER iMin, iMax
111        INTEGER jMin, jMax        INTEGER jMin, jMax
112        INTEGER bi, bj        INTEGER bi, bj
113        INTEGER i, j        INTEGER i, j
114        INTEGER k, kM1, kUp, kDown        INTEGER k, km1, kup, kDown
       LOGICAL BOTTOM_LAYER  
115    
116  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
117        INTEGER    isbyte        INTEGER    isbyte
# Line 154  C                      index into fVerTe Line 119  C                      index into fVerTe
119    
120        INTEGER act1, act2, act3, act4        INTEGER act1, act2, act3, act4
121        INTEGER max1, max2, max3        INTEGER max1, max2, max3
122        INTEGER ikact, iikey,kkey        INTEGER iikey, kkey
123        INTEGER maximpl        INTEGER maximpl
124  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
125    
126  C---    The algorithm...  C---    The algorithm...
127  C  C
# Line 171  C Line 136  C
136  C       "Calculation of Gs"  C       "Calculation of Gs"
137  C       ===================  C       ===================
138  C       This is where all the accelerations and tendencies (ie.  C       This is where all the accelerations and tendencies (ie.
139  C       phiHydysics, parameterizations etc...) are calculated  C       physics, parameterizations etc...) are calculated
140  C         rVel = sum_r ( div. u[n] )  C         rVel = sum_r ( div. u[n] )
141  C         rho = rho ( theta[n], salt[n] )  C         rho = rho ( theta[n], salt[n] )
142  C         b   = b(rho, theta)  C         b   = b(rho, theta)
# Line 206  C--- Line 171  C---
171  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
172  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
173        ikey = 1        ikey = 1
174  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
175    
176  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
177  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 184  C     uninitialised but inert locations.
184          yA(i,j)      = 0. _d 0          yA(i,j)      = 0. _d 0
185          uTrans(i,j)  = 0. _d 0          uTrans(i,j)  = 0. _d 0
186          vTrans(i,j)  = 0. _d 0          vTrans(i,j)  = 0. _d 0
187          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  
          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  
188           KappaRU(i,j,k) = 0. _d 0           KappaRU(i,j,k) = 0. _d 0
189           KappaRV(i,j,k) = 0. _d 0           KappaRV(i,j,k) = 0. _d 0
190             sigmaX(i,j,k) = 0. _d 0
191             sigmaY(i,j,k) = 0. _d 0
192             sigmaR(i,j,k) = 0. _d 0
193          ENDDO          ENDDO
194          rhoKM1 (i,j) = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
195          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
         rhoKP1 (i,j) = 0. _d 0  
         rhoTMP (i,j) = 0. _d 0  
         buoyKM1(i,j) = 0. _d 0  
         buoyK  (i,j) = 0. _d 0  
196          maskC  (i,j) = 0. _d 0          maskC  (i,j) = 0. _d 0
197         ENDDO         ENDDO
198        ENDDO        ENDDO
# Line 247  C     uninitialised but inert locations. Line 200  C     uninitialised but inert locations.
200    
201  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
202  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
203  !HPF$ INDEPENDENT  CHPF$ INDEPENDENT
204  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
205    
206        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
207    
208  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
209  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
210  !HPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV
211  !HPF$&                  ,phiHyd,K13,K23,K33,KapGM  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA
212  !HPF$&                  ,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV
213  !HPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  )
214  !HPF$&                  )  #endif /* ALLOW_AUTODIFF_TAMC */
 #endif  
215    
216         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
217    
# Line 278  C--    HPF directive to help TAMC Line 230  C--    HPF directive to help TAMC
230            ikey = (act1 + 1) + act2*max1            ikey = (act1 + 1) + act2*max1
231       &                      + act3*max1*max2       &                      + act3*max1*max2
232       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
233  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
234    
235  C--     Set up work arrays that need valid initial values  C--     Set up work arrays that need valid initial values
236          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
# Line 294  C--     Set up work arrays that need val Line 246  C--     Set up work arrays that need val
246            fVerU (i,j,2) = 0. _d 0            fVerU (i,j,2) = 0. _d 0
247            fVerV (i,j,1) = 0. _d 0            fVerV (i,j,1) = 0. _d 0
248            fVerV (i,j,2) = 0. _d 0            fVerV (i,j,2) = 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  
249           ENDDO           ENDDO
250          ENDDO          ENDDO
251    
# Line 320  C--     Set up work arrays that need val Line 267  C--     Set up work arrays that need val
267          jMax = sNy+OLy          jMax = sNy+OLy
268    
269    
270          K = 1  C--     Start of diagnostic loop
271          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)  
   
 #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  
272    
         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_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+1, myThid )  
          END IF  
 #endif  
         ENDIF  
 #endif  
 C--     Density of 1st level (below W(1)) reference to level 1  
 #ifdef  INCLUDE_FIND_RHO_CALL  
273  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
274  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  C? Patrick, is this formula correct now that we change the loop range?
275  CADJ STORE salt (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  C? Do we still need this?
276  #endif           kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1
277          CALL FIND_RHO(  #endif /* ALLOW_AUTODIFF_TAMC */
278       I     bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,  
279       O     rhoKm1,  C--       Integrate continuity vertically for vertical velocity
280       I     myThid )            CALL INTEGRATE_FOR_W(
281  #endif       I                         bi, bj, k, uVel, vVel,
282         O                         wVel,
283          IF (       (.NOT. BOTTOM_LAYER)       I                         myThid )
 #ifdef ALLOW_KPP  
      &       .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing  
 #endif  
      &     ) THEN  
 C--      Check static stability with layer below  
 C--      and mix as needed.  
 #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+1, K, eosType,  
      O      rhoKp1,  
      I      myThid )  
 #endif  
   
 #ifdef  INCLUDE_CONVECT_CALL  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhoKm1(:,:)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE rhoKp1(:,:)  = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
          CALL CONVECT(  
      I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,  
      U       ConvectCount,  
      I       myTime,myIter,myThid)  
 #ifdef ALLOW_AUTODIFF_TAMC  
 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  
284    
285  #endif  #ifdef    ALLOW_OBCS
286    C--       Calculate future values on open boundaries
 C--      Implicit Vertical Diffusion for Convection  
          IF (ivdc_kappa.NE.0.) 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 ?  
   
 C--      Recompute density after mixing  
 #ifdef  INCLUDE_FIND_RHO_CALL  
          CALL FIND_RHO(  
      I      bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,  
      O      rhoKm1,  
      I      myThid )  
 #endif  
         ENDIF  
 C--     Calculate buoyancy  
         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  
287            IF (openBoundaries) THEN            IF (openBoundaries) THEN
288  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef      ALLOW_NONHYDROSTATIC
289  CADJ STORE uvel (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte              IF (nonHydrostatic) THEN
290  CADJ STORE vvel (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte                CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
291  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte              ENDIF
292  CADJ STORE salt(:,:,k,bi,bj)   = comlev1_2d, key = ikey, byte = isbyte  #endif      /* ALLOW_NONHYDROSTATIC */
293  #endif              CALL OBCS_CALC( bi, bj, k, myTime+deltaT, myThid )
294               CALL APPLY_OBCS1( bi, bj, K+1, myThid )            ENDIF
295    #endif    /* ALLOW_OBCS */
296    
297    C--       Calculate gradients of potential density for isoneutral
298    C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
299    c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
300              IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
301                CALL FIND_RHO(
302         I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
303         I        theta, salt,
304         O        rhoK,
305         I        myThid )
306                IF (k.GT.1) CALL FIND_RHO(
307         I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
308         I        theta, salt,
309         O        rhoKm1,
310         I        myThid )
311                CALL GRAD_SIGMA(
312         I             bi, bj, iMin, iMax, jMin, jMax, k,
313         I             rhoK, rhoKm1, rhoK,
314         O             sigmaX, sigmaY, sigmaR,
315         I             myThid )
316              ENDIF
317    
318    C--       Implicit Vertical Diffusion for Convection
319    c ==> should use sigmaR !!!
320              IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
321                CALL CALC_IVDC(
322         I        bi, bj, iMin, iMax, jMin, jMax, k,
323         I        rhoKm1, rhoK,
324         U        ConvectCount, KappaRT, KappaRS,
325         I        myTime, myIter, myThid)
326            END IF            END IF
 #endif  
          ENDIF  
 #endif  
327    
328  C--      Density of K level (below W(K)) reference to K level  C--     end of diagnostic k loop (Nr:1)
329  #ifdef  INCLUDE_FIND_RHO_CALL          ENDDO
 #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, 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  
330    
331  #ifdef ALLOW_AUTODIFF_TAMC  C--     Determines forcing terms based on external fields
332  CADJ STORE rhok  (:,:)   = comlev1_3d, key = kkey, byte = isbyte  C       relaxation terms, etc.
333  CADJ STORE rhoKm1(:,:)   = comlev1_3d, key = kkey, byte = isbyte          CALL EXTERNAL_FORCING_SURF(
334  CADJ STORE rhoKp1(:,:)   = comlev1_3d, key = kkey, byte = isbyte       I             bi, bj, iMin, iMax, jMin, jMax,
335  #endif       I             myThid )
336    
337    #ifdef  ALLOW_GMREDI
338    C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
339            IF (useGMRedi) THEN
340              DO k=1,Nr
341                CALL GMREDI_CALC_TENSOR(
342         I             bi, bj, iMin, iMax, jMin, jMax, k,
343         I             sigmaX, sigmaY, sigmaR,
344         I             myThid )
345              ENDDO
346            ENDIF
347    #endif  /* ALLOW_GMREDI */
348    
349  #ifdef  INCLUDE_CONVECT_CALL  #ifdef  ALLOW_KPP
350            CALL CONVECT(  C--     Compute KPP mixing coefficients
351       I        bi,bj,iMin,iMax,jMin,jMax,K+1,rhoK,rhoKp1,          IF (useKPP) THEN
352       U        ConvectCount,            CALL KPP_CALC(
353       I        myTime,myIter,myThid)       I                  bi, bj, myTime, myThid )
354  #ifdef ALLOW_AUTODIFF_TAMC          ENDIF
355  CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)  #endif  /* ALLOW_KPP */
 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  
356    
357  C--      Implicit Vertical Diffusion for Convection  #ifdef ALLOW_AUTODIFF_TAMC
358           IF (ivdc_kappa.NE.0.) THEN  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
359              CALL CALC_IVDC(  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
360       I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
361       U       ConvectCount, KappaRT, KappaRS,  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
362       I       myTime,myIter,myThid)  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
363  CRG: do we need do store STORE KappaRT, KappaRS ?  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
364           END IF  #endif /* ALLOW_AUTODIFF_TAMC */
365    
 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  
366    
 #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  
367    
368           DO J=jMin,jMax  C--     Start of thermodynamics loop
369            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  
370    
371  #ifdef ALLOW_AUTODIFF_TAMC  C--       km1    Points to level above k (=k-1)
372  CADJ STORE theta(:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  C--       kup    Cycles through 1,2 to point to layer above
373  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  
374    
375  #ifdef ALLOW_KPP            km1  = MAX(1,k-1)
376  C----------------------------------------------            kup  = 1+MOD(k+1,2)
377  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  
378    
379  C----------------------------------------------            iMin = 1-OLx+2
380  C--     start of upward loop            iMax = sNx+OLx-1
381  C----------------------------------------------            jMin = 1-OLy+2
382          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  
383    
384  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
385           kkey = ikact*(Nr-1+1) + (k-1) + 1  CPatrick Is this formula correct?
386  #endif           kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1
387    CADJ STORE rvel  (:,:,kDown) = comlev1_bibj_k, key = kkey, byte = isbyte
388  #ifdef ALLOW_AUTODIFF_TAMC  CADJ STORE rTrans(:,:)       = comlev1_bibj_k, key = kkey, byte = isbyte
389  CADJ STORE rvel  (:,:,kDown)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRT(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
390  CADJ STORE rTrans(:,:)        = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRS(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
391  CADJ STORE KappaRT(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  #endif /* ALLOW_AUTODIFF_TAMC */
 CADJ STORE KappaRS(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
392    
393  C--      Get temporary terms used by tendency routines  C--      Get temporary terms used by tendency routines
394           CALL CALC_COMMON_FACTORS (           CALL CALC_COMMON_FACTORS (
395       I        bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
396       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,
397       I        myThid)       I        myThid)
398    
 #ifdef ALLOW_OBCS  
         IF (openBoundaries) THEN  
          CALL APPLY_OBCS3( bi, bj, K, Kup, rTrans, rVel, myThid )  
         ENDIF  
 #endif  
   
399  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
400  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
401           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
402       I        bi,bj,iMin,iMax,jMin,jMax,K,       I        bi,bj,iMin,iMax,jMin,jMax,k,
403       I        maskC,maskUp,KapGM,K33,       I        maskC,maskup,
404       O        KappaRT,KappaRS,KappaRU,KappaRV,       O        KappaRT,KappaRS,KappaRU,KappaRV,
405       I        myThid)       I        myThid)
406  #endif  #endif
407  C--      Calculate accelerations in the momentum equations  
408           IF ( momStepping ) THEN  C--      Calculate active tracer tendencies (gT,gS,...)
409            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  
410           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
411            CALL CALC_GT(             CALL CALC_GT(
412       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
413       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
414       I         K13,K23,KappaRT,KapGM,       I         KappaRT,
415       U         aTerm,xTerm,fZon,fMer,fVerT,       U         fVerT,
416       I         myTime, myThid)       I         myTime, myThid)
417               CALL TIMESTEP_TRACER(
418         I         bi,bj,iMin,iMax,jMin,jMax,k,
419         I         theta, gT,
420         U         gTnm1,
421         I         myIter, myThid)
422           ENDIF           ENDIF
423           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
424            CALL CALC_GS(             CALL CALC_GS(
425       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
426       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
427       I         K13,K23,KappaRS,KapGM,       I         KappaRS,
428       U         aTerm,xTerm,fZon,fMer,fVerS,       U         fVerS,
429       I         myTime, myThid)       I         myTime, myThid)
430               CALL TIMESTEP_TRACER(
431         I         bi,bj,iMin,iMax,jMin,jMax,k,
432         I         salt, gS,
433         U         gSnm1,
434         I         myIter, myThid)
435           ENDIF           ENDIF
436  #ifdef ALLOW_OBCS  
437  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  
438  C--      Apply open boundary conditions  C--      Apply open boundary conditions
439           IF (openBoundaries) THEN           IF (openBoundaries) THEN
440  #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 )  
441           END IF           END IF
442  #endif  #endif   /* ALLOW_OBCS */
443    
444  C--      Freeze water  C--      Freeze water
445           IF (allowFreezing) THEN           IF (allowFreezing) THEN
446  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
447  CADJ STORE gTNm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
448  #endif  CADJ &   , key = kkey, byte = isbyte
449              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, K, myThid )  #endif /* ALLOW_AUTODIFF_TAMC */
450                CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
451           END IF           END IF
452    
453  #ifdef DIVG_IN_DYNAMICS  C--     end of thermodynamic k loop (Nr:1)
454  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  
455    
456    
457          ENDDO ! K  #ifdef ALLOW_AUTODIFF_TAMC
458    CPatrick? What about this one?
459               maximpl = 6
460               iikey = (ikey-1)*maximpl
461    #endif /* ALLOW_AUTODIFF_TAMC */
462    
463  C--     Implicit diffusion  C--     Implicit diffusion
464          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
465    
466  #ifdef ALLOW_AUTODIFF_TAMC            IF (tempStepping) THEN
            maximpl = 6  
            iikey = ikact*maximpl  
 #endif  
   
          IF (tempStepping) THEN  
467  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
468              idkey = iikey + 1              idkey = iikey + 1
469  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
470              CALL IMPLDIFF(              CALL IMPLDIFF(
471       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
472       I         deltaTtracer, KappaRT,recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
473       U         gTNm1,       U         gTNm1,
474       I         myThid )       I         myThid )
475           END IF           ENDIF
476    
477           IF (saltStepping) THEN           IF (saltStepping) THEN
478  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
479           idkey = iikey + 2           idkey = iikey + 2
480  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
481              CALL IMPLDIFF(              CALL IMPLDIFF(
482       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
483       I         deltaTtracer, KappaRS,recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
484       U         gSNm1,       U         gSNm1,
485       I         myThid )       I         myThid )
486             ENDIF
487    
488    #ifdef   ALLOW_OBCS
489    C--      Apply open boundary conditions
490             IF (openBoundaries) THEN
491               DO K=1,Nr
492                 CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
493               ENDDO
494           END IF           END IF
495    #endif   /* ALLOW_OBCS */
496    
497          ENDIF ! implicitDiffusion  C--     End If implicitDiffusion
498            ENDIF
499    
 C--     Implicit viscosity  
         IF (implicitViscosity) THEN  
500    
501           IF (momStepping) THEN  
502  #ifdef ALLOW_AUTODIFF_TAMC  C--     Start of dynamics loop
503           idkey = iikey + 3          DO k=1,Nr
504  #endif  
505    C--       km1    Points to level above k (=k-1)
506    C--       kup    Cycles through 1,2 to point to layer above
507    C--       kDown  Cycles through 2,1 to point to current layer
508    
509              km1  = MAX(1,k-1)
510              kup  = 1+MOD(k+1,2)
511              kDown= 1+MOD(k,2)
512    
513              iMin = 1-OLx+2
514              iMax = sNx+OLx-1
515              jMin = 1-OLy+2
516              jMax = sNy+OLy-1
517    
518    C--      Integrate hydrostatic balance for phiHyd with BC of
519    C        phiHyd(z=0)=0
520    C        distinguishe between Stagger and Non Stagger time stepping
521             IF (staggerTimeStep) THEN
522               CALL CALC_PHI_HYD(
523         I        bi,bj,iMin,iMax,jMin,jMax,k,
524         I        gTnm1, gSnm1,
525         U        phiHyd,
526         I        myThid )
527             ELSE
528               CALL CALC_PHI_HYD(
529         I        bi,bj,iMin,iMax,jMin,jMax,k,
530         I        theta, salt,
531         U        phiHyd,
532         I        myThid )
533             ENDIF
534    
535    C--      Calculate accelerations in the momentum equations (gU, gV, ...)
536    C        and step forward storing the result in gUnm1, gVnm1, etc...
537             IF ( momStepping ) THEN
538               CALL CALC_MOM_RHS(
539         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
540         I         phiHyd,KappaRU,KappaRV,
541         U         fVerU, fVerV,
542         I         myTime, myThid)
543               CALL TIMESTEP(
544         I         bi,bj,iMin,iMax,jMin,jMax,k,phiHyd,
545         I         myIter, myThid)
546    
547    #ifdef   ALLOW_OBCS
548    C--      Apply open boundary conditions
549             IF (openBoundaries) THEN
550               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
551             END IF
552    #endif   /* ALLOW_OBCS */
553    
554    #ifdef   ALLOW_AUTODIFF_TAMC
555    #ifdef   INCLUDE_CD_CODE
556             ELSE
557               DO j=1-OLy,sNy+OLy
558                 DO i=1-OLx,sNx+OLx
559                   guCD(i,j,k,bi,bj) = 0.0
560                   gvCD(i,j,k,bi,bj) = 0.0
561                 END DO
562               END DO
563    #endif   /* INCLUDE_CD_CODE */
564    #endif   /* ALLOW_AUTODIFF_TAMC */
565             ENDIF
566    
567    
568    C--     end of dynamics k loop (1:Nr)
569            ENDDO
570    
571    
572    
573    C--     Implicit viscosity
574            IF (implicitViscosity.AND.momStepping) THEN
575    #ifdef    ALLOW_AUTODIFF_TAMC
576              idkey = iikey + 3
577    #endif    /* ALLOW_AUTODIFF_TAMC */
578            CALL IMPLDIFF(            CALL IMPLDIFF(
579       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
580       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
581       U         gUNm1,       U         gUNm1,
582       I         myThid )       I         myThid )
583  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
584           idkey = iikey + 4            idkey = iikey + 4
585  #endif  #endif    /* ALLOW_AUTODIFF_TAMC */
586            CALL IMPLDIFF(            CALL IMPLDIFF(
587       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
588       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
589       U         gVNm1,       U         gVNm1,
590       I         myThid )       I         myThid )
591    
592  #ifdef INCLUDE_CD_CODE  #ifdef   ALLOW_OBCS
593    C--      Apply open boundary conditions
594             IF (openBoundaries) THEN
595               DO K=1,Nr
596                 CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
597               ENDDO
598             END IF
599    #endif   /* ALLOW_OBCS */
600    
601  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    INCLUDE_CD_CODE
602           idkey = iikey + 5  #ifdef    ALLOW_AUTODIFF_TAMC
603  #endif            idkey = iikey + 5
604    #endif    /* ALLOW_AUTODIFF_TAMC */
605            CALL IMPLDIFF(            CALL IMPLDIFF(
606       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
607       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
608       U         vVelD,       U         vVelD,
609       I         myThid )       I         myThid )
610  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
611          idkey = iikey + 6            idkey = iikey + 6
612  #endif  #endif    /* ALLOW_AUTODIFF_TAMC */
613            CALL IMPLDIFF(            CALL IMPLDIFF(
614       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
615       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
616       U         uVelD,       U         uVelD,
617       I         myThid )       I         myThid )
618    #endif    /* INCLUDE_CD_CODE */
619  #endif  C--     End If implicitViscosity.AND.momStepping
620            ENDIF
          ENDIF ! momStepping  
         ENDIF ! implicitViscosity  
621    
622         ENDDO         ENDDO
623        ENDDO        ENDDO
624    
 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 )  
   
   
625        RETURN        RETURN
626        END        END
627    
628    
629    C--      Cumulative diagnostic calculations (ie. time-averaging)
630    #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE
631    c        IF (taveFreq.GT.0.) THEN
632    c         CALL DO_TIME_AVERAGES(
633    c    I                           myTime, myIter, bi, bj, k, kup, kDown,
634    c    I                           ConvectCount,
635    c    I                           myThid )
636    c        ENDIF
637    #endif
638    
Legend:
Removed from v.1.49  
changed lines
  Added in v.1.54.2.11
  ViewVC Help
Powered by ViewVC 1.1.22