/[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.56 by heimbach, Mon Jan 29 20:05:46 2001 UTC revision 1.71 by cnh, Mon Jun 18 17:39:58 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
# Line 25  C     \================================= Line 26  C     \=================================
26  C     == Global variables ===  C     == Global variables ===
27  #include "SIZE.h"  #include "SIZE.h"
28  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #include "CG2D.h"  
29  #include "PARAMS.h"  #include "PARAMS.h"
30  #include "DYNVARS.h"  #include "DYNVARS.h"
31  #include "GRID.h"  #include "GRID.h"
# Line 42  C     == Global variables === Line 42  C     == Global variables ===
42  # endif  # endif
43  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
44    
45    #ifdef ALLOW_TIMEAVE
46    #include "TIMEAVE_STATV.h"
47    #endif
48    
49  C     == Routine arguments ==  C     == Routine arguments ==
50  C     myTime - Current time in simulation  C     myTime - Current time in simulation
51  C     myIter - Current iteration number in simulation  C     myIter - Current iteration number in simulation
# Line 54  C     == Local variables Line 58  C     == Local variables
58  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
59  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C     uTrans, vTrans, rTrans - Per block temporaries holding flow
60  C                              transport  C                              transport
61  C     rVel                     o uTrans: Zonal transport  C                              o uTrans: Zonal transport
62  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
63  C                              o rTrans: Vertical transport  C                              o rTrans: Vertical transport
64  C                              o rVel:   Vertical velocity at upper and  C     maskUp                   o maskUp: land/water mask for W points
65  C                                        lower cell faces.  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
 C     maskC,maskUp             o maskC: land/water mask for tracer cells  
 C                              o maskUp: land/water mask for W points  
 C     aTerm, xTerm, cTerm    - Work arrays for holding separate terms in  
 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  
66  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
67  C                                      so we need an fVer for each  C                                      so we need an fVer for each
68  C                                      variable.  C                                      variable.
69  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.  
70  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     phiHyd         - Hydrostatic part of the potential phiHydi.
71  C                      In z coords phiHydiHyd is the hydrostatic  C                      In z coords phiHydiHyd is the hydrostatic
72  C                      pressure anomaly  C                      Potential (=pressure/rho0) anomaly
73  C                      In p coords phiHydiHyd is the geopotential  C                      In p coords phiHydiHyd is the geopotential
74  C                      surface height  C                      surface height anomaly.
75  C                      anomaly.  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
76  C     etaSurfX,      - Holds surface elevation gradient in X and Y.  C     phiSurfY             or geopotentiel (atmos) in X and Y direction
 C     etaSurfY  
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
# Line 94  C     bi, bj Line 82  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    C     tauAB - Adams-Bashforth timestepping weight: 0=forward ; 1/2=Adams-Bashf.
86        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rVel    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
       _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
91        _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)  
92        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
93        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
94        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
95        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
96        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97        _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)  
98        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99        _RL buoyKM1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RL buoyK   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfY(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)  
101        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
102        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
103        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
# Line 129  C                      index into fVerTe Line 105  C                      index into fVerTe
105        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
106        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
107        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
108          _RL tauAB
109    
110  C This is currently also used by IVDC and Diagnostics  C This is currently used by IVDC and Diagnostics
 C #ifdef INCLUDE_CONVECT_CALL  
111        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
 C #endif  
112    
113        INTEGER iMin, iMax        INTEGER iMin, iMax
114        INTEGER jMin, jMax        INTEGER jMin, jMax
115        INTEGER bi, bj        INTEGER bi, bj
116        INTEGER i, j        INTEGER i, j
117        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
       LOGICAL BOTTOM_LAYER  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
       INTEGER    isbyte  
       PARAMETER( isbyte = 4 )  
   
       INTEGER act1, act2, act3, act4  
       INTEGER max1, max2, max3  
       INTEGER iikey, kkey  
       INTEGER maximpl  
 #endif /* ALLOW_AUTODIFF_TAMC */  
118    
119    Cjmc : add for phiHyd output <- but not working if multi tile per CPU
120    c     CHARACTER*(MAX_LEN_MBUF) suff
121    c     LOGICAL  DIFFERENT_MULTIPLE
122    c     EXTERNAL DIFFERENT_MULTIPLE
123    Cjmc(end)
124    
125  C---    The algorithm...  C---    The algorithm...
126  C  C
127  C       "Correction Step"  C       "Correction Step"
# Line 166  C       "Calculation of Gs" Line 136  C       "Calculation of Gs"
136  C       ===================  C       ===================
137  C       This is where all the accelerations and tendencies (ie.  C       This is where all the accelerations and tendencies (ie.
138  C       physics, parameterizations etc...) are calculated  C       physics, parameterizations etc...) are calculated
 C         rVel = sum_r ( div. u[n] )  
139  C         rho = rho ( theta[n], salt[n] )  C         rho = rho ( theta[n], salt[n] )
140  C         b   = b(rho, theta)  C         b   = b(rho, theta)
141  C         K31 = K31 ( rho )  C         K31 = K31 ( rho )
142  C         Gu[n] = Gu( u[n], v[n], rVel, b, ... )  C         Gu[n] = Gu( u[n], v[n], wVel, b, ... )
143  C         Gv[n] = Gv( u[n], v[n], rVel, b, ... )  C         Gv[n] = Gv( u[n], v[n], wVel, b, ... )
144  C         Gt[n] = Gt( theta[n], u[n], v[n], rVel, K31, ... )  C         Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
145  C         Gs[n] = Gs( salt[n], u[n], v[n], rVel, K31, ... )  C         Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
146  C  C
147  C       "Time-stepping" or "Prediction"  C       "Time-stepping" or "Prediction"
148  C       ================================  C       ================================
# Line 202  C--   dummy statement to end declaration Line 171  C--   dummy statement to end declaration
171        ikey = 1        ikey = 1
172  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
173    
   
174  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
175  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
176  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
# Line 214  C     uninitialised but inert locations. Line 182  C     uninitialised but inert locations.
182          yA(i,j)      = 0. _d 0          yA(i,j)      = 0. _d 0
183          uTrans(i,j)  = 0. _d 0          uTrans(i,j)  = 0. _d 0
184          vTrans(i,j)  = 0. _d 0          vTrans(i,j)  = 0. _d 0
         aTerm(i,j)   = 0. _d 0  
         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  
185          DO k=1,Nr          DO k=1,Nr
186           phiHyd (i,j,k)  = 0. _d 0           phiHyd(i,j,k)  = 0. _d 0
187           KappaRU(i,j,k) = 0. _d 0           KappaRU(i,j,k) = 0. _d 0
188           KappaRV(i,j,k) = 0. _d 0           KappaRV(i,j,k) = 0. _d 0
189           sigmaX(i,j,k) = 0. _d 0           sigmaX(i,j,k) = 0. _d 0
# Line 231  C     uninitialised but inert locations. Line 192  C     uninitialised but inert locations.
192          ENDDO          ENDDO
193          rhoKM1 (i,j) = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
194          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
195          rhoKP1 (i,j) = 0. _d 0          phiSurfX(i,j) = 0. _d 0
196          rhoTMP (i,j) = 0. _d 0          phiSurfY(i,j) = 0. _d 0
         buoyKM1(i,j) = 0. _d 0  
         buoyK  (i,j) = 0. _d 0  
         maskC  (i,j) = 0. _d 0  
197         ENDDO         ENDDO
198        ENDDO        ENDDO
199    
# Line 249  CHPF$ INDEPENDENT Line 207  CHPF$ INDEPENDENT
207    
208  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
209  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
210  CHPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV
211  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA
212  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV
213  CHPF$&                  )  CHPF$&                  )
214  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 278  C--     Set up work arrays that need val Line 236  C--     Set up work arrays that need val
236          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
237           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
238            rTrans(i,j)   = 0. _d 0            rTrans(i,j)   = 0. _d 0
           rVel  (i,j,1) = 0. _d 0  
           rVel  (i,j,2) = 0. _d 0  
239            fVerT (i,j,1) = 0. _d 0            fVerT (i,j,1) = 0. _d 0
240            fVerT (i,j,2) = 0. _d 0            fVerT (i,j,2) = 0. _d 0
241            fVerS (i,j,1) = 0. _d 0            fVerS (i,j,1) = 0. _d 0
# Line 288  C--     Set up work arrays that need val Line 244  C--     Set up work arrays that need val
244            fVerU (i,j,2) = 0. _d 0            fVerU (i,j,2) = 0. _d 0
245            fVerV (i,j,1) = 0. _d 0            fVerV (i,j,1) = 0. _d 0
246            fVerV (i,j,2) = 0. _d 0            fVerV (i,j,2) = 0. _d 0
           phiHyd(i,j,1) = 0. _d 0  
247           ENDDO           ENDDO
248          ENDDO          ENDDO
249    
250          DO k=1,Nr          DO k=1,Nr
251           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
252            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
253  #ifdef INCLUDE_CONVECT_CALL  C This is currently also used by IVDC and Diagnostics
254             ConvectCount(i,j,k) = 0.             ConvectCount(i,j,k) = 0.
 #endif  
255             KappaRT(i,j,k) = 0. _d 0             KappaRT(i,j,k) = 0. _d 0
256             KappaRS(i,j,k) = 0. _d 0             KappaRS(i,j,k) = 0. _d 0
257            ENDDO            ENDDO
# Line 309  C--     Set up work arrays that need val Line 263  C--     Set up work arrays that need val
263          jMin = 1-OLy+1          jMin = 1-OLy+1
264          jMax = sNy+OLy          jMax = sNy+OLy
265    
         k = 1  
         BOTTOM_LAYER = k .EQ. Nr  
   
 #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_bibj, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
            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  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL APPLY_OBCS1( bi, bj, k+1, myThid )  
          END IF  
 #endif  
         ENDIF  
 #endif  
   
 C--     Density of 1st level (below W(1)) reference to level 1  
 #ifdef  INCLUDE_FIND_RHO_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         CALL FIND_RHO(  
      I     bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  
      O     rhoKm1,  
      I     myThid )  
 #endif  
   
         IF (.NOT. BOTTOM_LAYER) 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+1,bi,bj) = comlev1_bibj  
 CADJ &   , key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,k+1,bi,bj) = comlev1_bibj  
 CADJ &   , key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
          CALL FIND_RHO(  
      I      bi, bj, iMin, iMax, jMin, jMax, k+1, k, eosType,  
      O      rhoKp1,  
      I      myThid )  
 #endif  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhoKm1(:,:)  = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE rhoKp1(:,:)  = comlev1_bibj, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef  INCLUDE_CONVECT_CALL  
   
          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_bibj, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)  
 CADJ &     = comlev1_bibj, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  
   
 C--      Implicit Vertical Diffusion for Convection  
          IF (ivdc_kappa.NE.0.) THEN  
             CALL CALC_IVDC(  
      I       bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,  
      U       ConvectCount, KappaRT, KappaRS,  
      I       myTime,myIter,myThid)  
          ENDIF  
   
 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 )  
   
 #ifdef ALLOW_GMREDI  
         IF ( useGMRedi ) THEN  
         CALL GRAD_SIGMA(  
      I            bi, bj, iMin, iMax, jMin, jMax, k,  
      I            rhoKm1, rhoKm1, rhoKm1,  
      O            sigmaX, sigmaY, sigmaR,  
      I            myThid )  
         ELSE  
          DO j=1-OLy,sNy+OLy  
           DO i=1-OLx,sNx+OLx  
            sigmaX(i,j,k) = 0. _d 0  
            sigmaY(i,j,k) = 0. _d 0  
            sigmaR(i,j,k) = 0. _d 0  
           ENDDO  
          ENDDO  
         ENDIF  
 #endif  
   
 C--     Start of downward loop  
         DO k=2,Nr  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
          kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
          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_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
              CALL APPLY_OBCS1( bi, bj, k+1, myThid )  
           END IF  
 #endif  
          ENDIF  
 #endif /* DO_PIPELINED_CORRECTION_STEP */  
   
 C--      Density of k level (below W(k)) reference to k level  
 #ifdef  INCLUDE_FIND_RHO_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
          CALL FIND_RHO(  
      I      bi, bj, iMin, iMax, jMin, jMax,  k, k, eosType,  
      O      rhoK,  
      I      myThid )  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph(   storing not necessary  
 cphCADJ STORE rhoK(:,:)  = comlev1_bibj_k, key = kkey, byte = isbyte  
 cph)  
 #endif /* ALLOW_AUTODIFF_TAMC */  
 #endif  
   
          IF (.NOT. BOTTOM_LAYER) 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+1,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 CADJ STORE salt (:,:,k+1,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
           CALL FIND_RHO(  
      I       bi, bj, iMin, iMax, jMin, jMax,  k+1, k, eosType,  
      O       rhoKp1,  
      I       myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhoKp1(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
 #endif  
   
 #ifdef  INCLUDE_CONVECT_CALL  
           CALL CONVECT(  
      I        bi,bj,iMin,iMax,jMin,jMax,k+1,rhoK,rhoKp1,  
      U        ConvectCount,  
      I        myTime,myIter,myThid)  
   
 #endif  
266    
 C--      Implicit Vertical Diffusion for Convection  
          IF (ivdc_kappa.NE.0.) THEN  
267  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
268  CADJ STORE rhoKm1(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
269  #endif /* ALLOW_AUTODIFF_TAMC */  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
270              CALL CALC_IVDC(  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
271       I       bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
272       U       ConvectCount, KappaRT, KappaRS,  #endif /* ALLOW_AUTODIFF_TAMC */
273       I       myTime,myIter,myThid)  
274           END IF  C--     Start of diagnostic loop
275            DO k=Nr,1,-1
276  C--       Recompute density after mixing  
277  #ifdef  INCLUDE_FIND_RHO_CALL  #ifdef ALLOW_AUTODIFF_TAMC
278  #ifdef ALLOW_AUTODIFF_TAMC  C? Patrick, is this formula correct now that we change the loop range?
279  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k  C? Do we still need this?
280  CADJ &   , key = kkey, byte = isbyte  cph kkey formula corrected.
281  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k  cph Needed for rhok, rhokm1, in the case useGMREDI.
282  CADJ &   , key = kkey, byte = isbyte           kkey = (ikey-1)*Nr + k
283  #endif /* ALLOW_AUTODIFF_TAMC */  CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
284            CALL FIND_RHO(  CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
285       I       bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  #endif /* ALLOW_AUTODIFF_TAMC */
286       O       rhoK,  
287       I       myThid )  C--       Integrate continuity vertically for vertical velocity
288  #endif            CALL INTEGRATE_FOR_W(
289         I                         bi, bj, k, uVel, vVel,
290  C--            IF (.NOT. BOTTOM_LAYER) ends here       O                         wVel,
291           ENDIF       I                         myThid )
292    
293  C--      Calculate buoyancy  #ifdef    ALLOW_OBCS
294           CALL CALC_BUOYANCY(  #ifdef    ALLOW_NONHYDROSTATIC
295       I       bi,bj,iMin,iMax,jMin,jMax,k,rhoK,  C--       Apply OBC to W if in N-H mode
296       O       buoyK,            IF (useOBCS.AND.nonHydrostatic) THEN
297       I       myThid )              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
298              ENDIF
299  C--      Integrate hydrostatic balance for phiHyd with BC of  #endif    /* ALLOW_NONHYDROSTATIC */
300  C--      phiHyd(z=0)=0  #endif    /* ALLOW_OBCS */
301           CALL CALC_PHI_HYD(  
302       I        bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyK,  C--       Calculate gradients of potential density for isoneutral
303       U        phiHyd,  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
304    c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
305              IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
306    #ifdef ALLOW_AUTODIFF_TAMC
307    CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
308    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
309    #endif /* ALLOW_AUTODIFF_TAMC */
310                CALL FIND_RHO(
311         I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
312         I        theta, salt,
313         O        rhoK,
314       I        myThid )       I        myThid )
315                IF (k.GT.1) THEN
 #ifdef  INCLUDE_FIND_RHO_CALL  
 C--      Calculate iso-neutral slopes for the GM/Redi parameterisation  
   
316  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
317  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
318  CADJ &   , key = kkey, byte = isbyte  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
319  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
320                 CALL FIND_RHO(
          CALL FIND_RHO(  
321       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
322       O        rhoTmp,       I        theta, salt,
323         O        rhoKm1,
324       I        myThid )       I        myThid )
325  #endif              ENDIF
326                CALL GRAD_SIGMA(
   
 #ifdef ALLOW_GMREDI  
          IF ( useGMRedi ) THEN  
          CALL GRAD_SIGMA(  
327       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
328       I             rhoK, rhotmp, rhoK,       I             rhoK, rhoKm1, rhoK,
329       O             sigmaX, sigmaY, sigmaR,       O             sigmaX, sigmaY, sigmaR,
330       I             myThid )       I             myThid )
331           ELSE            ENDIF
           DO j=1-OLy,sNy+OLy  
            DO i=1-OLx,sNx+OLx  
             sigmaX(i,j,k) = 0. _d 0  
             sigmaY(i,j,k) = 0. _d 0  
             sigmaR(i,j,k) = 0. _d 0  
            ENDDO  
           ENDDO  
          ENDIF  
 #endif  
332    
333           DO J=jMin,jMax  C--       Implicit Vertical Diffusion for Convection
334            DO I=iMin,iMax  c ==> should use sigmaR !!!
335  #ifdef  INCLUDE_FIND_RHO_CALL            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
336             rhoKm1 (I,J) = rhoK(I,J)              CALL CALC_IVDC(
337  #endif       I        bi, bj, iMin, iMax, jMin, jMax, k,
338             buoyKm1(I,J) = buoyK(I,J)       I        rhoKm1, rhoK,
339            ENDDO       U        ConvectCount, KappaRT, KappaRS,
340           ENDDO       I        myTime, myIter, myThid)
341              ENDIF
342    
343  C--     end of k loop  C--     end of diagnostic k loop (Nr:1)
344          ENDDO          ENDDO
345    
346  C     Determines forcing terms based on external fields  #ifdef ALLOW_AUTODIFF_TAMC
347  C     relaxation terms, etc.  cph avoids recomputation of integrate_for_w
348        CALL EXTERNAL_FORCING_SURF(  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
349    #endif /* ALLOW_AUTODIFF_TAMC */
350    
351    #ifdef  ALLOW_OBCS
352    C--     Calculate future values on open boundaries
353            IF (useOBCS) THEN
354              CALL OBCS_CALC( bi, bj, myTime+deltaT,
355         I            uVel, vVel, wVel, theta, salt,
356         I            myThid )
357            ENDIF
358    #endif  /* ALLOW_OBCS */
359    
360    C--     Determines forcing terms based on external fields
361    C       relaxation terms, etc.
362            CALL EXTERNAL_FORCING_SURF(
363       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
364       I             myThid )       I             myThid )
   
365  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
366    cph needed for KPP
367    CADJ STORE surfacetendencyU(:,:,bi,bj)
368    CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
369    CADJ STORE surfacetendencyV(:,:,bi,bj)
370    CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
371    CADJ STORE surfacetendencyS(:,:,bi,bj)
372    CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
373    CADJ STORE surfacetendencyT(:,:,bi,bj)
374    CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
375    #endif /* ALLOW_AUTODIFF_TAMC */
376    
377  CADJ STORE surfacetendencyu(:,:,bi,bj)  #ifdef  ALLOW_GMREDI
 CADJ &   , surfacetendencyv(:,:,bi,bj)  
 CADJ &   , surfacetendencys(:,:,bi,bj)  
 CADJ &   , surfacetendencyt(:,:,bi,bj)  
 CADJ &                        = comlev1_bibj, key=ikey, byte=isbyte  
378    
379  # ifdef ALLOW_GMREDI  #ifdef ALLOW_AUTODIFF_TAMC
380  CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte
381  CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte
382  CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte
 # endif /* ALLOW_GMREDI */  
   
383  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
384    C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
 #ifdef ALLOW_GMREDI  
385          IF (useGMRedi) THEN          IF (useGMRedi) THEN
386            DO k=1, Nr            DO k=1,Nr
387              CALL GMREDI_CALC_TENSOR(              CALL GMREDI_CALC_TENSOR(
388       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
389       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
# Line 661  CADJ STORE sigmaR(:,:,:) = comlev1, key= Line 399  CADJ STORE sigmaR(:,:,:) = comlev1, key=
399            ENDDO            ENDDO
400  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
401          ENDIF          ENDIF
 #endif  
402    
403  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
404  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte
405  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte
406    CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte
 #ifdef ALLOW_GMREDI  
 C-- R.G. We need to define a new tape since Kw use mythid instead of bi,bj  
 CADJ STORE Kwx(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwy(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwz(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte  
 #endif  
   
 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 =======  
 C-- R.G. We need to define a new tape since Kw use mythid instead of bi,bj  
 CADJ STORE Kwx(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwy(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwz(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte  
   
 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
   
 C--     dummy initialization to break data flow because  
 C--     calc_div_ghat has a condition for initialization  
         DO J=jMin,jMax  
            DO I=iMin,iMax  
               cg2d_b(i,j,bi,bj) = 0.0  
            ENDDO  
         ENDDO  
407  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
408    
409  #ifdef ALLOW_KPP  #endif  /* ALLOW_GMREDI */
 C--   Compute KPP mixing coefficients  
         IF (useKPP) THEN  
410    
411            CALL TIMER_START('KPP_CALC               [DYNAMICS]', myThid)  #ifdef  ALLOW_KPP
412    C--     Compute KPP mixing coefficients
413            IF (useKPP) THEN
414            CALL KPP_CALC(            CALL KPP_CALC(
415       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myThid )
           CALL TIMER_STOP ('KPP_CALC               [DYNAMICS]', myThid)  
   
416  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
417          ELSE          ELSE
418            DO j=1-OLy,sNy+OLy            CALL KPP_CALC_DUMMY(
419              DO i=1-OLx,sNx+OLx       I                  bi, bj, myTime, myThid )
               KPPhbl (i,j,bi,bj) = 1.0  
               KPPfrac(i,j,bi,bj) = 0.0  
               DO k = 1,Nr  
                  KPPghat   (i,j,k,bi,bj) = 0.0  
                  KPPviscAz (i,j,k,bi,bj) = viscAz  
                  KPPdiffKzT(i,j,k,bi,bj) = diffKzT  
                  KPPdiffKzS(i,j,k,bi,bj) = diffKzS  
               ENDDO  
             ENDDO  
           ENDDO  
420  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
421          ENDIF          ENDIF
422    
# Line 733  CADJ &   , KPPfrac   (:,:  ,bi,bj) Line 429  CADJ &   , KPPfrac   (:,:  ,bi,bj)
429  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte
430  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
431    
432  #endif /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
433    
434  C--     Start of upward loop  #ifdef ALLOW_AUTODIFF_TAMC
435          DO k = Nr, 1, -1  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
436    CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
437    CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
438    CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
439    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
440    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
441    #endif /* ALLOW_AUTODIFF_TAMC */
442    
443    #ifdef ALLOW_AIM
444    C       AIM - atmospheric intermediate model, physics package code.
445    C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics
446            IF ( useAIM ) THEN
447             CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
448             CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )
449             CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
450            ENDIF
451    #endif /* ALLOW_AIM */
452    
 C--      km1    Points to level above k (=k-1)  
 C--      kup    Cycles through 1,2 to point to layer above  
 C--      kDown  Cycles through 2,1 to point to current layer  
   
          km1  =max(1,k-1)  
          kup  =1+MOD(k+1,2)  
          kDown=1+MOD(k,2)  
   
          iMin = 1-OLx+2  
          iMax = sNx+OLx-1  
          jMin = 1-OLy+2  
          jMax = sNy+OLy-1  
453    
454    C--     Start of thermodynamics loop
455            DO k=Nr,1,-1
456  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
457           kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1  C? Patrick Is this formula correct?
458  CADJ STORE rvel  (:,:,kdown)  = comlev1_bibj_k, key=kkey, byte=isbyte  cph Yes, but I rewrote it.
459  CADJ STORE rTrans(:,:)        = comlev1_bibj_k, key=kkey, byte=isbyte  cph Also, the KappaR? need the index and subscript k!
460  CADJ STORE KappaRT(:,:,k)     = comlev1_bibj_k, key=kkey, byte=isbyte           kkey = (ikey-1)*Nr + k
461  CADJ STORE KappaRS(:,:,k)     = comlev1_bibj_k, key=kkey, byte=isbyte  #endif /* ALLOW_AUTODIFF_TAMC */
462  #endif /* ALLOW_AUTODIFF_TAMC */  
463    C--       km1    Points to level above k (=k-1)
464    C--       kup    Cycles through 1,2 to point to layer above
465    C--       kDown  Cycles through 2,1 to point to current layer
466    
467              km1  = MAX(1,k-1)
468              kup  = 1+MOD(k+1,2)
469              kDown= 1+MOD(k,2)
470    
471              iMin = 1-OLx+2
472              iMax = sNx+OLx-1
473              jMin = 1-OLy+2
474              jMax = sNy+OLy-1
475    
476  C--      Get temporary terms used by tendency routines  C--      Get temporary terms used by tendency routines
477           CALL CALC_COMMON_FACTORS (           CALL CALC_COMMON_FACTORS (
478       I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,       I        bi,bj,iMin,iMax,jMin,jMax,k,
479       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,       O        xA,yA,uTrans,vTrans,rTrans,maskUp,
480       I        myThid)       I        myThid)
481    
482  #ifdef ALLOW_OBCS  #ifdef ALLOW_AUTODIFF_TAMC
483          IF (openBoundaries) THEN  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
484           CALL APPLY_OBCS3( bi, bj, k, kup, rTrans, rVel, myThid )  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
485          ENDIF  #endif /* ALLOW_AUTODIFF_TAMC */
 #endif  
486    
487  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
488  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
489           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
490       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
491       I        maskC,maskUp,       I        maskUp,
492       O        KappaRT,KappaRS,KappaRU,KappaRV,       O        KappaRT,KappaRS,KappaRU,KappaRV,
493       I        myThid)       I        myThid)
494  #endif  #endif
495  C--      Calculate accelerations in the momentum equations  
496           IF ( momStepping ) THEN  C--      Calculate active tracer tendencies (gT,gS,...)
497            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 /* ALLOW_AUTODIFF_TAMC */  
          ENDIF  
 C--      Calculate active tracer tendencies  
498           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
499            CALL CALC_GT(             CALL CALC_GT(
500       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
501       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
502       I         KappaRT,       I         KappaRT,
503       U         aTerm,xTerm,fZon,fMer,fVerT,       U         fVerT,
504       I         myTime, myThid)       I         myTime, myThid)
505               tauAB = 0.5d0 + abEps
506               CALL TIMESTEP_TRACER(
507         I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
508         I         theta, gT,
509         U         gTnm1,
510         I         myIter, myThid)
511           ENDIF           ENDIF
512           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
513            CALL CALC_GS(             CALL CALC_GS(
514       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
515       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
516       I         KappaRS,       I         KappaRS,
517       U         aTerm,xTerm,fZon,fMer,fVerS,       U         fVerS,
518       I         myTime, myThid)       I         myTime, myThid)
519               tauAB = 0.5d0 + abEps
520               CALL TIMESTEP_TRACER(
521         I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
522         I         salt, gS,
523         U         gSnm1,
524         I         myIter, myThid)
525           ENDIF           ENDIF
 #ifdef ALLOW_OBCS  
 C--      Calculate future values on open boundaries  
          IF (openBoundaries) THEN  
 Caja      CALL CYCLE_OBCS( k, bi, bj, myThid )  
           CALL SET_OBCS( k, bi, bj, myTime+deltaTclock, myThid )  
          ENDIF  
 #endif  
 C--      Prediction step (step forward all model variables)  
          CALL TIMESTEP(  
      I       bi,bj,iMin,iMax,jMin,jMax,k,  
      I       myIter, myThid)  
 #ifdef ALLOW_OBCS  
 C--      Apply open boundary conditions  
          IF (openBoundaries) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE gunm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE gvnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE gwnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
526    
527              CALL APPLY_OBCS2( bi, bj, k, myThid )  #ifdef   ALLOW_OBCS
528    C--      Apply open boundary conditions
529             IF (useOBCS) THEN
530               CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
531           END IF           END IF
532  #endif  #endif   /* ALLOW_OBCS */
533    
534  C--      Freeze water  C--      Freeze water
535           IF (allowFreezing) THEN           IF (allowFreezing) THEN
536  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
537  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
538    CADJ &   , key = kkey, byte = isbyte
539  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
540              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
541           END IF           END IF
542    
543  #ifdef DIVG_IN_DYNAMICS  C--     end of thermodynamic k loop (Nr:1)
 C--      Diagnose barotropic divergence of predicted fields  
          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                           rVel, ConvectCount,  
      I                           myThid )  
          ENDIF  
 #endif  
   
   
 C--     k loop  
544          ENDDO          ENDDO
545    
546    
547  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
548             maximpl = 6  C? Patrick? What about this one?
549             iikey = (ikey-1)*maximpl  cph Keys iikey and idkey don't seem to be needed
550    cph since storing occurs on different tape for each
551    cph impldiff call anyways.
552    cph Thus, common block comlev1_impl isn't needed either.
553    cph Storing below needed in the case useGMREDI.
554            iikey = (ikey-1)*maximpl
555  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
556    
557  C--     Implicit diffusion  C--     Implicit diffusion
# Line 885  CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_ Line 564  CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_
564  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
565              CALL IMPLDIFF(              CALL IMPLDIFF(
566       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
567       I         deltaTtracer, KappaRT,recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
568       U         gTNm1,       U         gTNm1,
569       I         myThid )       I         myThid )
570           END IF           ENDIF
571    
572           IF (saltStepping) THEN           IF (saltStepping) THEN
573  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 897  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_ Line 576  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_
576  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
577              CALL IMPLDIFF(              CALL IMPLDIFF(
578       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
579       I         deltaTtracer, KappaRS,recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
580       U         gSNm1,       U         gSNm1,
581       I         myThid )       I         myThid )
582             ENDIF
583    
584    #ifdef   ALLOW_OBCS
585    C--      Apply open boundary conditions
586             IF (useOBCS) THEN
587               DO K=1,Nr
588                 CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
589               ENDDO
590           END IF           END IF
591    #endif   /* ALLOW_OBCS */
592    
593  C--     implicitDiffusion  C--     End If implicitDiffusion
594          ENDIF          ENDIF
595    
596  C--     Implicit viscosity  C--     Start computation of dynamics
597          IF (implicitViscosity) THEN          iMin = 1-OLx+2
598            iMax = sNx+OLx-1
599            jMin = 1-OLy+2
600            jMax = sNy+OLy-1
601    
602    C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
603    C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
604            IF (implicSurfPress.NE.1.) THEN
605              CALL CALC_GRAD_PHI_SURF(
606         I         bi,bj,iMin,iMax,jMin,jMax,
607         I         etaN,
608         O         phiSurfX,phiSurfY,
609         I         myThid )                        
610            ENDIF
611    
612           IF (momStepping) THEN  C--     Start of dynamics loop
613  #ifdef ALLOW_AUTODIFF_TAMC          DO k=1,Nr
614           idkey = iikey + 3  
615    C--       km1    Points to level above k (=k-1)
616    C--       kup    Cycles through 1,2 to point to layer above
617    C--       kDown  Cycles through 2,1 to point to current layer
618    
619              km1  = MAX(1,k-1)
620              kup  = 1+MOD(k+1,2)
621              kDown= 1+MOD(k,2)
622    
623    C--      Integrate hydrostatic balance for phiHyd with BC of
624    C        phiHyd(z=0)=0
625    C        distinguishe between Stagger and Non Stagger time stepping
626             IF (staggerTimeStep) THEN
627               CALL CALC_PHI_HYD(
628         I        bi,bj,iMin,iMax,jMin,jMax,k,
629         I        gTnm1, gSnm1,
630         U        phiHyd,
631         I        myThid )
632             ELSE
633               CALL CALC_PHI_HYD(
634         I        bi,bj,iMin,iMax,jMin,jMax,k,
635         I        theta, salt,
636         U        phiHyd,
637         I        myThid )
638             ENDIF
639    
640    C--      Calculate accelerations in the momentum equations (gU, gV, ...)
641    C        and step forward storing the result in gUnm1, gVnm1, etc...
642             IF ( momStepping ) THEN
643               CALL CALC_MOM_RHS(
644         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
645         I         phiHyd,KappaRU,KappaRV,
646         U         fVerU, fVerV,
647         I         myTime, myThid)
648               CALL TIMESTEP(
649         I         bi,bj,iMin,iMax,jMin,jMax,k,
650         I         phiHyd, phiSurfX, phiSurfY,
651         I         myIter, myThid)
652    
653    #ifdef   ALLOW_OBCS
654    C--      Apply open boundary conditions
655             IF (useOBCS) THEN
656               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
657             END IF
658    #endif   /* ALLOW_OBCS */
659    
660    #ifdef   ALLOW_AUTODIFF_TAMC
661    #ifdef   INCLUDE_CD_CODE
662             ELSE
663               DO j=1-OLy,sNy+OLy
664                 DO i=1-OLx,sNx+OLx
665                   guCD(i,j,k,bi,bj) = 0.0
666                   gvCD(i,j,k,bi,bj) = 0.0
667                 END DO
668               END DO
669    #endif   /* INCLUDE_CD_CODE */
670    #endif   /* ALLOW_AUTODIFF_TAMC */
671             ENDIF
672    
673    
674    C--     end of dynamics k loop (1:Nr)
675            ENDDO
676    
677    
678    
679    C--     Implicit viscosity
680            IF (implicitViscosity.AND.momStepping) THEN
681    #ifdef    ALLOW_AUTODIFF_TAMC
682              idkey = iikey + 3
683  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
684  #endif /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
685            CALL IMPLDIFF(            CALL IMPLDIFF(
686       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
687       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
688       U         gUNm1,       U         gUNm1,
689       I         myThid )       I         myThid )
690  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
691           idkey = iikey + 4            idkey = iikey + 4
692  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
693  #endif /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
694            CALL IMPLDIFF(            CALL IMPLDIFF(
695       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
696       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
697       U         gVNm1,       U         gVNm1,
698       I         myThid )       I         myThid )
699    
700  #ifdef INCLUDE_CD_CODE  #ifdef   ALLOW_OBCS
701    C--      Apply open boundary conditions
702             IF (useOBCS) THEN
703               DO K=1,Nr
704                 CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
705               ENDDO
706             END IF
707    #endif   /* ALLOW_OBCS */
708    
709  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    INCLUDE_CD_CODE
710           idkey = iikey + 5  #ifdef    ALLOW_AUTODIFF_TAMC
711              idkey = iikey + 5
712  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
713  #endif /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
714            CALL IMPLDIFF(            CALL IMPLDIFF(
715       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
716       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
717       U         vVelD,       U         vVelD,
718       I         myThid )       I         myThid )
719  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
720          idkey = iikey + 6            idkey = iikey + 6
721  CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
722  #endif /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
723            CALL IMPLDIFF(            CALL IMPLDIFF(
724       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
725       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
726       U         uVelD,       U         uVelD,
727       I         myThid )       I         myThid )
728    #endif    /* INCLUDE_CD_CODE */
729  #endif  C--     End If implicitViscosity.AND.momStepping
   
 C--      momStepping  
          ENDIF  
   
 C--     implicitViscosity  
730          ENDIF          ENDIF
731    
732    Cjmc : add for phiHyd output <- but not working if multi tile per CPU
733    c       IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)
734    c    &  .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
735    c         WRITE(suff,'(I10.10)') myIter+1
736    c         CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)
737    c       ENDIF
738    Cjmc(end)
739    
740    #ifdef ALLOW_TIMEAVE
741            IF (taveFreq.GT.0.) THEN
742              CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr,
743         I                              deltaTclock, bi, bj, myThid)
744              IF (ivdc_kappa.NE.0.) THEN
745                CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
746         I                              deltaTclock, bi, bj, myThid)
747              ENDIF
748            ENDIF
749    #endif /* ALLOW_TIMEAVE */
750    
751         ENDDO         ENDDO
752        ENDDO        ENDDO
753    
754    #ifndef EXCLUDE_DEBUGMODE
755          If (debugMode) THEN
756           CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
757           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
758           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
759           CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
760           CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
761           CALL DEBUG_STATS_RL(Nr,Gu,'Gu (DYNAMICS)',myThid)
762           CALL DEBUG_STATS_RL(Nr,Gv,'Gv (DYNAMICS)',myThid)
763           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (DYNAMICS)',myThid)
764           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (DYNAMICS)',myThid)
765           CALL DEBUG_STATS_RL(Nr,GuNm1,'GuNm1 (DYNAMICS)',myThid)
766           CALL DEBUG_STATS_RL(Nr,GvNm1,'GvNm1 (DYNAMICS)',myThid)
767           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (DYNAMICS)',myThid)
768           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (DYNAMICS)',myThid)
769          ENDIF
770    #endif
771    
772        RETURN        RETURN
773        END        END
Legend:
Removed from v.1.56  
changed lines
  Added in v.1.71
  ViewVC Help
Powered by ViewVC 1.1.22