/[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.64 by jmc, Tue Mar 6 16:59:44 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 20  C     | ===== Line 21  C     | =====
21  C     | C*P* comments indicating place holders for which code is |  C     | C*P* comments indicating place holders for which code is |
22  C     |      presently being developed.                          |  C     |      presently being developed.                          |
23  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.  
   
24        IMPLICIT NONE        IMPLICIT NONE
25    
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"
32    
33    #ifdef ALLOW_AUTODIFF_TAMC
34    # include "tamc.h"
35    # include "tamc_keys.h"
36    #endif /* ALLOW_AUTODIFF_TAMC */
37    
38  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
39  #include "KPPMIX.h"  # include "KPP.h"
40  #endif  #endif
41    
42  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_TIMEAVE
43  #include "tamc.h"  #include "TIMEAVE_STATV.h"
 #include "tamc_keys.h"  
44  #endif  #endif
45    
46  C     == Routine arguments ==  C     == Routine arguments ==
# Line 57  C     == Local variables Line 55  C     == Local variables
55  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
56  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C     uTrans, vTrans, rTrans - Per block temporaries holding flow
57  C                              transport  C                              transport
58  C     rVel                     o uTrans: Zonal transport  C                              o uTrans: Zonal transport
59  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
60  C                              o rTrans: Vertical transport  C                              o rTrans: Vertical transport
 C                              o rVel:   Vertical velocity at upper and  
 C                                        lower cell faces.  
61  C     maskC,maskUp             o maskC: land/water mask for tracer cells  C     maskC,maskUp             o maskC: land/water mask for tracer cells
62  C                              o maskUp: land/water mask for W points  C                              o maskUp: land/water mask for W points
63  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  
64  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
65  C                                      so we need an fVer for each  C                                      so we need an fVer for each
66  C                                      variable.  C                                      variable.
67  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.  
68  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     phiHyd         - Hydrostatic part of the potential phiHydi.
69  C                      In z coords phiHydiHyd is the hydrostatic  C                      In z coords phiHydiHyd is the hydrostatic
70  C                      pressure anomaly  C                      pressure anomaly
71  C                      In p coords phiHydiHyd is the geopotential  C                      In p coords phiHydiHyd is the geopotential
72  C                      surface height  C                      surface height
73  C                      anomaly.  C                      anomaly.
74  C     etaSurfX,      - Holds surface elevation gradient in X and Y.  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
75  C     etaSurfY  C     phiSurfY             or geopotentiel (atmos) in X and Y direction
 C     K13, K23, K33  - Non-zero elements of small-angle approximation  
 C                      diffusion tensor.  
 C     KapGM          - Spatially varying Visbeck et. al mixing coeff.  
76  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
77  C     KappaRS          (background + spatially varying, isopycnal term).  C     KappaRS          (background + spatially varying, isopycnal term).
78  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
79  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
80  C     bi, bj  C     bi, bj
81  C     k, kUp,        - Index for layer above and below. kUp and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
82  C     kDown, kM1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
83  C                      index into fVerTerm.  C                      index into fVerTerm.
84        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88        _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)  
89        _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90        _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)  
91        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
92        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
93        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
94        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
95        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
96        _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)  
97        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RL buoyKM1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99        _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)  
       _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)  
100        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
101        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
102        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
103        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
104          _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
105          _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
106          _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
107    
108  #ifdef INCLUDE_CONVECT_CALL  C This is currently used by IVDC and Diagnostics
109        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
 #endif  
110    
111        INTEGER iMin, iMax        INTEGER iMin, iMax
112        INTEGER jMin, jMax        INTEGER jMin, jMax
113        INTEGER bi, bj        INTEGER bi, bj
114        INTEGER i, j        INTEGER i, j
115        INTEGER k, kM1, kUp, kDown        INTEGER k, km1, kup, kDown
       LOGICAL BOTTOM_LAYER  
116    
117    Cjmc : add for phiHyd output <- but not working if multi tile per CPU
118    c     CHARACTER*(MAX_LEN_MBUF) suff
119    c     LOGICAL  DIFFERENT_MULTIPLE
120    c     EXTERNAL DIFFERENT_MULTIPLE
121    Cjmc(end)
122    
123  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
124        INTEGER    isbyte        INTEGER    isbyte
125        PARAMETER( isbyte = 4 )        PARAMETER( isbyte = 4 )
126    
127        INTEGER act1, act2, act3, act4        INTEGER act1, act2, act3, act4
128        INTEGER max1, max2, max3        INTEGER max1, max2, max3
129        INTEGER ikact, iikey,kkey        INTEGER iikey, kkey
130        INTEGER maximpl        INTEGER maximpl
131  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
132    
133  C---    The algorithm...  C---    The algorithm...
134  C  C
# Line 171  C Line 143  C
143  C       "Calculation of Gs"  C       "Calculation of Gs"
144  C       ===================  C       ===================
145  C       This is where all the accelerations and tendencies (ie.  C       This is where all the accelerations and tendencies (ie.
146  C       phiHydysics, parameterizations etc...) are calculated  C       physics, parameterizations etc...) are calculated
 C         rVel = sum_r ( div. u[n] )  
147  C         rho = rho ( theta[n], salt[n] )  C         rho = rho ( theta[n], salt[n] )
148  C         b   = b(rho, theta)  C         b   = b(rho, theta)
149  C         K31 = K31 ( rho )  C         K31 = K31 ( rho )
150  C         Gu[n] = Gu( u[n], v[n], rVel, b, ... )  C         Gu[n] = Gu( u[n], v[n], wVel, b, ... )
151  C         Gv[n] = Gv( u[n], v[n], rVel, b, ... )  C         Gv[n] = Gv( u[n], v[n], wVel, b, ... )
152  C         Gt[n] = Gt( theta[n], u[n], v[n], rVel, K31, ... )  C         Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
153  C         Gs[n] = Gs( salt[n], u[n], v[n], rVel, K31, ... )  C         Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
154  C  C
155  C       "Time-stepping" or "Prediction"  C       "Time-stepping" or "Prediction"
156  C       ================================  C       ================================
# Line 206  C--- Line 177  C---
177  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
178  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
179        ikey = 1        ikey = 1
180  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
181    
182  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
183  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 190  C     uninitialised but inert locations.
190          yA(i,j)      = 0. _d 0          yA(i,j)      = 0. _d 0
191          uTrans(i,j)  = 0. _d 0          uTrans(i,j)  = 0. _d 0
192          vTrans(i,j)  = 0. _d 0          vTrans(i,j)  = 0. _d 0
193          aTerm(i,j)   = 0. _d 0          DO k=1,Nr
194          xTerm(i,j)   = 0. _d 0           phiHyd(i,j,k)  = 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  
195           KappaRU(i,j,k) = 0. _d 0           KappaRU(i,j,k) = 0. _d 0
196           KappaRV(i,j,k) = 0. _d 0           KappaRV(i,j,k) = 0. _d 0
197             sigmaX(i,j,k) = 0. _d 0
198             sigmaY(i,j,k) = 0. _d 0
199             sigmaR(i,j,k) = 0. _d 0
200          ENDDO          ENDDO
201          rhoKM1 (i,j) = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
202          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  
203          maskC  (i,j) = 0. _d 0          maskC  (i,j) = 0. _d 0
204            phiSurfX(i,j) = 0. _d 0
205            phiSurfY(i,j) = 0. _d 0
206         ENDDO         ENDDO
207        ENDDO        ENDDO
208    
209    
210  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
211  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
212  !HPF$ INDEPENDENT  CHPF$ INDEPENDENT
213  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
214    
215        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
216    
217  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
218  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
219  !HPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV
220  !HPF$&                  ,phiHyd,K13,K23,K33,KapGM  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA
221  !HPF$&                  ,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV
222  !HPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  )
223  !HPF$&                  )  #endif /* ALLOW_AUTODIFF_TAMC */
 #endif  
224    
225         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
226    
# Line 278  C--    HPF directive to help TAMC Line 239  C--    HPF directive to help TAMC
239            ikey = (act1 + 1) + act2*max1            ikey = (act1 + 1) + act2*max1
240       &                      + act3*max1*max2       &                      + act3*max1*max2
241       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
242  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
243    
244  C--     Set up work arrays that need valid initial values  C--     Set up work arrays that need valid initial values
245          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
246           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
247            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  
248            fVerT (i,j,1) = 0. _d 0            fVerT (i,j,1) = 0. _d 0
249            fVerT (i,j,2) = 0. _d 0            fVerT (i,j,2) = 0. _d 0
250            fVerS (i,j,1) = 0. _d 0            fVerS (i,j,1) = 0. _d 0
# Line 294  C--     Set up work arrays that need val Line 253  C--     Set up work arrays that need val
253            fVerU (i,j,2) = 0. _d 0            fVerU (i,j,2) = 0. _d 0
254            fVerV (i,j,1) = 0. _d 0            fVerV (i,j,1) = 0. _d 0
255            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  
256           ENDDO           ENDDO
257          ENDDO          ENDDO
258    
259          DO k=1,Nr          DO k=1,Nr
260           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
261            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
262  #ifdef INCLUDE_CONVECT_CALL  C This is currently also used by IVDC and Diagnostics
263             ConvectCount(i,j,k) = 0.             ConvectCount(i,j,k) = 0.
 #endif  
264             KappaRT(i,j,k) = 0. _d 0             KappaRT(i,j,k) = 0. _d 0
265             KappaRS(i,j,k) = 0. _d 0             KappaRS(i,j,k) = 0. _d 0
266            ENDDO            ENDDO
# Line 320  C--     Set up work arrays that need val Line 273  C--     Set up work arrays that need val
273          jMax = sNy+OLy          jMax = sNy+OLy
274    
275    
276          K = 1  C--     Start of diagnostic loop
277          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  
   
         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  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
         CALL FIND_RHO(  
      I     bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,  
      O     rhoKm1,  
      I     myThid )  
 #endif  
   
         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  
 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  
278    
279  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
280  CADJ STORE rhoKm1(:,:)  = comlev1_2d, key = ikey, byte = isbyte  C? Patrick, is this formula correct now that we change the loop range?
281  CADJ STORE rhoKp1(:,:)  = comlev1_2d, key = ikey, byte = isbyte  C? Do we still need this?
282  #endif           kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1
283           CALL CONVECT(  #endif /* ALLOW_AUTODIFF_TAMC */
284       I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,  
285       U       ConvectCount,  C--       Integrate continuity vertically for vertical velocity
286       I       myTime,myIter,myThid)            CALL INTEGRATE_FOR_W(
287  #ifdef ALLOW_AUTODIFF_TAMC       I                         bi, bj, k, uVel, vVel,
288  CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)       O                         wVel,
289  CADJ &     = comlev1_2d, key = ikey, byte = isbyte       I                         myThid )
290  CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)  
291  CADJ &     = comlev1_2d, key = ikey, byte = isbyte  #ifdef    ALLOW_OBCS
292  #endif  #ifdef    ALLOW_NONHYDROSTATIC
293    C--       Apply OBC to W if in N-H mode
294              IF (useOBCS.AND.nonHydrostatic) THEN
295                CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
296              ENDIF
297    #endif    /* ALLOW_NONHYDROSTATIC */
298    #endif    /* ALLOW_OBCS */
299    
300    C--       Calculate gradients of potential density for isoneutral
301    C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
302    c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
303              IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
304                CALL FIND_RHO(
305         I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
306         I        theta, salt,
307         O        rhoK,
308         I        myThid )
309                IF (k.GT.1) CALL FIND_RHO(
310         I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
311         I        theta, salt,
312         O        rhoKm1,
313         I        myThid )
314                CALL GRAD_SIGMA(
315         I             bi, bj, iMin, iMax, jMin, jMax, k,
316         I             rhoK, rhoKm1, rhoK,
317         O             sigmaX, sigmaY, sigmaR,
318         I             myThid )
319              ENDIF
320    
321    C--       Implicit Vertical Diffusion for Convection
322    c ==> should use sigmaR !!!
323              IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
324                CALL CALC_IVDC(
325         I        bi, bj, iMin, iMax, jMin, jMax, k,
326         I        rhoKm1, rhoK,
327         U        ConvectCount, KappaRT, KappaRS,
328         I        myTime, myIter, myThid)
329              ENDIF
330    
331  #endif  C--     end of diagnostic k loop (Nr:1)
332            ENDDO
333    
334  C--      Implicit Vertical Diffusion for Convection  #ifdef  ALLOW_OBCS
335           IF (ivdc_kappa.NE.0.) CALL CALC_IVDC(  C--     Calculate future values on open boundaries
336       I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,          IF (useOBCS) THEN
337       U       ConvectCount, KappaRT, KappaRS,            CALL OBCS_CALC( bi, bj, myTime+deltaT,
338       I       myTime,myIter,myThid)       I            uVel, vVel, wVel, theta, salt,
339  CRG: do we need do store STORE KappaRT, KappaRS ?       I            myThid )
   
 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  
340          ENDIF          ENDIF
341  C--     Calculate buoyancy  #endif  /* ALLOW_OBCS */
         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  
342    
343    C--     Determines forcing terms based on external fields
344    C       relaxation terms, etc.
345            CALL EXTERNAL_FORCING_SURF(
346         I             bi, bj, iMin, iMax, jMin, jMax,
347         I             myThid )
348    
349    #ifdef  ALLOW_GMREDI
350    C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
351            IF (useGMRedi) THEN
352              DO k=1,Nr
353                CALL GMREDI_CALC_TENSOR(
354         I             bi, bj, iMin, iMax, jMin, jMax, k,
355         I             sigmaX, sigmaY, sigmaR,
356         I             myThid )
357              ENDDO
358  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
359           kkey = ikact*(Nr-2+1) + (k-2) + 1          ELSE
360  #endif            DO k=1, Nr
361                CALL GMREDI_CALC_TENSOR_DUMMY(
362           BOTTOM_LAYER = K .EQ. Nr       I             bi, bj, iMin, iMax, jMin, jMax, k,
363         I             sigmaX, sigmaY, sigmaR,
364  #ifdef DO_PIPELINED_CORRECTION_STEP       I             myThid )
365           IF ( .NOT. BOTTOM_LAYER ) THEN            ENDDO
366  C--       Update fields in layer below according to tendency terms  #endif /* ALLOW_AUTODIFF_TAMC */
367            CALL CORRECTION_STEP(          ENDIF
368       I         bi,bj,iMin,iMax,jMin,jMax,K+1,  #endif  /* ALLOW_GMREDI */
      I         etaSurfX,etaSurfY,myTime,myThid)  
 #ifdef ALLOW_OBCS  
           IF (openBoundaries) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE uvel (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE vvel (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE salt(:,:,k,bi,bj)   = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
              CALL APPLY_OBCS1( bi, bj, K+1, myThid )  
           END IF  
 #endif  
          ENDIF  
 #endif  
369    
370  C--      Density of K level (below W(K)) reference to K level  #ifdef  ALLOW_KPP
371  #ifdef  INCLUDE_FIND_RHO_CALL  C--     Compute KPP mixing coefficients
372  #ifdef ALLOW_AUTODIFF_TAMC          IF (useKPP) THEN
373  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte            CALL KPP_CALC(
374  CADJ STORE salt (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte       I                  bi, bj, myTime, myThid )
375  #endif          ENDIF
376           CALL FIND_RHO(  #endif  /* ALLOW_KPP */
      I      bi, bj, iMin, iMax, jMin, jMax,  K, K, eosType,  
      O      rhoK,  
      I      myThid )  
 #endif  
          IF (       (.NOT. BOTTOM_LAYER)  
 #ifdef ALLOW_KPP  
      &       .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing  
 #endif  
      &      ) THEN  
 C--       Check static stability with layer below and mix as needed.  
 C--       Density of K+1 level (below W(K+1)) reference to K level.  
 #ifdef  INCLUDE_FIND_RHO_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
           CALL FIND_RHO(  
      I       bi, bj, iMin, iMax, jMin, jMax,  K+1, K, eosType,  
      O       rhoKp1,  
      I       myThid )  
 #endif  
377    
378  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
379  CADJ STORE rhok  (:,:)   = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
380  CADJ STORE rhoKm1(:,:)   = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
381  CADJ STORE rhoKp1(:,:)   = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
382  #endif  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
383    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
384    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
385    #endif /* ALLOW_AUTODIFF_TAMC */
386    
387    #ifdef ALLOW_AIM
388    C       AIM - atmospheric intermediate model, physics package code.
389    C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics
390            IF ( useAIM ) THEN
391             CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
392             CALL AIM_DO_ATMOS_PHYSICS( phiHyd, myTime, myThid )
393             CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
394            ENDIF
395    #endif /* ALLOW_AIM */
396    
 #ifdef  INCLUDE_CONVECT_CALL  
           CALL CONVECT(  
      I        bi,bj,iMin,iMax,jMin,jMax,K+1,rhoK,rhoKp1,  
      U        ConvectCount,  
      I        myTime,myIter,myThid)  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)  
 CADJ &     = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)  
 CADJ &     = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
 #endif  
397    
398  C--      Implicit Vertical Diffusion for Convection  C--     Start of thermodynamics loop
399           IF (ivdc_kappa.NE.0.) THEN          DO k=Nr,1,-1
             CALL CALC_IVDC(  
      I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,  
      U       ConvectCount, KappaRT, KappaRS,  
      I       myTime,myIter,myThid)  
 CRG: do we need do store STORE KappaRT, KappaRS ?  
          END IF  
400    
401  C--       Recompute density after mixing  C--       km1    Points to level above k (=k-1)
402  #ifdef  INCLUDE_FIND_RHO_CALL  C--       kup    Cycles through 1,2 to point to layer above
403            CALL FIND_RHO(  C--       kDown  Cycles through 2,1 to point to current layer
      I       bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,  
      O       rhoK,  
      I       myThid )  
 #endif  
          ENDIF  
 C--      Calculate buoyancy  
          CALL CALC_BUOYANCY(  
      I       bi,bj,iMin,iMax,jMin,jMax,K,rhoK,  
      O       buoyK,  
      I       myThid )  
 C--      Integrate hydrostatic balance for phiHyd with BC of  
 C--      phiHyd(z=0)=0  
          CALL CALC_PHI_HYD(  
      I        bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyK,  
      U        phiHyd,  
      I        myThid )  
 C--      Calculate iso-neutral slopes for the GM/Redi parameterisation  
 #ifdef  INCLUDE_FIND_RHO_CALL  
          CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, K-1, K, eosType,  
      O        rhoTmp,  
      I        myThid )  
 #endif  
   
 #ifdef  INCLUDE_CALC_ISOSLOPES_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhoTmp(:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE rhok  (:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE rhoKm1(:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE kapgm (:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
          CALL CALC_ISOSLOPES(  
      I        bi, bj, iMin, iMax, jMin, jMax, K,  
      I        rhoKm1, rhoK, rhotmp,  
      O        K13, K23, K33, KapGM,  
      I        myThid )  
 #endif  
   
          DO J=jMin,jMax  
           DO I=iMin,iMax  
 #ifdef  INCLUDE_FIND_RHO_CALL  
            rhoKm1 (I,J) = rhoK(I,J)  
 #endif  
            buoyKm1(I,J) = buoyK(I,J)  
           ENDDO  
          ENDDO  
         ENDDO  
 C--     end of k loop  
404    
405  #ifdef ALLOW_AUTODIFF_TAMC            km1  = MAX(1,k-1)
406  CADJ STORE theta(:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte            kup  = 1+MOD(k+1,2)
407  CADJ STORE salt (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte            kDown= 1+MOD(k,2)
 CADJ STORE uvel (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
408    
409  #ifdef ALLOW_KPP            iMin = 1-OLx+2
410  C----------------------------------------------            iMax = sNx+OLx-1
411  C--     Compute KPP mixing coefficients            jMin = 1-OLy+2
412  C----------------------------------------------            jMax = sNy+OLy-1
         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  
   
 C----------------------------------------------  
 C--     start of upward loop  
 C----------------------------------------------  
         DO K = Nr, 1, -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  
413    
414  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
415           kkey = ikact*(Nr-1+1) + (k-1) + 1  CPatrick Is this formula correct?
416  #endif           kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1
417    CADJ STORE rTrans(:,:)       = comlev1_bibj_k, key = kkey, byte = isbyte
418  #ifdef ALLOW_AUTODIFF_TAMC  CADJ STORE KappaRT(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
419  CADJ STORE rvel  (:,:,kDown)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRS(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
420  CADJ STORE rTrans(:,:)        = comlev1_3d, key = kkey, byte = isbyte  #endif /* ALLOW_AUTODIFF_TAMC */
 CADJ STORE KappaRT(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE KappaRS(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
421    
422  C--      Get temporary terms used by tendency routines  C--      Get temporary terms used by tendency routines
423           CALL CALC_COMMON_FACTORS (           CALL CALC_COMMON_FACTORS (
424       I        bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
425       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,       O        xA,yA,uTrans,vTrans,rTrans,maskC,maskUp,
426       I        myThid)       I        myThid)
427    
 #ifdef ALLOW_OBCS  
         IF (openBoundaries) THEN  
          CALL APPLY_OBCS3( bi, bj, K, Kup, rTrans, rVel, myThid )  
         ENDIF  
 #endif  
   
428  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
429  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
430           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
431       I        bi,bj,iMin,iMax,jMin,jMax,K,       I        bi,bj,iMin,iMax,jMin,jMax,k,
432       I        maskC,maskUp,KapGM,K33,       I        maskC,maskup,
433       O        KappaRT,KappaRS,KappaRU,KappaRV,       O        KappaRT,KappaRS,KappaRU,KappaRV,
434       I        myThid)       I        myThid)
435  #endif  #endif
436  C--      Calculate accelerations in the momentum equations  
437           IF ( momStepping ) THEN  C--      Calculate active tracer tendencies (gT,gS,...)
438            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  
439           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
440            CALL CALC_GT(             CALL CALC_GT(
441       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
442       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
443       I         K13,K23,KappaRT,KapGM,       I         KappaRT,
444       U         aTerm,xTerm,fZon,fMer,fVerT,       U         fVerT,
445       I         myTime, myThid)       I         myTime, myThid)
446               CALL TIMESTEP_TRACER(
447         I         bi,bj,iMin,iMax,jMin,jMax,k,
448         I         theta, gT,
449         U         gTnm1,
450         I         myIter, myThid)
451           ENDIF           ENDIF
452           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
453            CALL CALC_GS(             CALL CALC_GS(
454       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
455       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
456       I         K13,K23,KappaRS,KapGM,       I         KappaRS,
457       U         aTerm,xTerm,fZon,fMer,fVerS,       U         fVerS,
458       I         myTime, myThid)       I         myTime, myThid)
459               CALL TIMESTEP_TRACER(
460         I         bi,bj,iMin,iMax,jMin,jMax,k,
461         I         salt, gS,
462         U         gSnm1,
463         I         myIter, myThid)
464           ENDIF           ENDIF
465  #ifdef ALLOW_OBCS  
466  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  
467  C--      Apply open boundary conditions  C--      Apply open boundary conditions
468           IF (openBoundaries) THEN           IF (useOBCS) THEN
469  #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 )  
470           END IF           END IF
471  #endif  #endif   /* ALLOW_OBCS */
472    
473  C--      Freeze water  C--      Freeze water
474           IF (allowFreezing) THEN           IF (allowFreezing) THEN
475  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
476  CADJ STORE gTNm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
477  #endif  CADJ &   , key = kkey, byte = isbyte
478              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, K, myThid )  #endif /* ALLOW_AUTODIFF_TAMC */
479                CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
480           END IF           END IF
481    
482  #ifdef DIVG_IN_DYNAMICS  C--     end of thermodynamic k loop (Nr:1)
483  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  
484    
485    
486          ENDDO ! K  #ifdef ALLOW_AUTODIFF_TAMC
487    CPatrick? What about this one?
488               maximpl = 6
489               iikey = (ikey-1)*maximpl
490    #endif /* ALLOW_AUTODIFF_TAMC */
491    
492  C--     Implicit diffusion  C--     Implicit diffusion
493          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
494    
 #ifdef ALLOW_AUTODIFF_TAMC  
            maximpl = 6  
            iikey = ikact*maximpl  
 #endif  
   
495           IF (tempStepping) THEN           IF (tempStepping) THEN
496  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
497              idkey = iikey + 1              idkey = iikey + 1
498  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
499              CALL IMPLDIFF(              CALL IMPLDIFF(
500       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
501       I         deltaTtracer, KappaRT,recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
502       U         gTNm1,       U         gTNm1,
503       I         myThid )       I         myThid )
504           END IF           ENDIF
505    
506           IF (saltStepping) THEN           IF (saltStepping) THEN
507  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
508           idkey = iikey + 2           idkey = iikey + 2
509  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
510              CALL IMPLDIFF(              CALL IMPLDIFF(
511       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
512       I         deltaTtracer, KappaRS,recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
513       U         gSNm1,       U         gSNm1,
514       I         myThid )       I         myThid )
515             ENDIF
516    
517    #ifdef   ALLOW_OBCS
518    C--      Apply open boundary conditions
519             IF (useOBCS) THEN
520               DO K=1,Nr
521                 CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
522               ENDDO
523           END IF           END IF
524    #endif   /* ALLOW_OBCS */
525    
526          ENDIF ! implicitDiffusion  C--     End If implicitDiffusion
527            ENDIF
528    
529  C--     Implicit viscosity  C--     Start computation of dynamics
530          IF (implicitViscosity) THEN          iMin = 1-OLx+2
531            iMax = sNx+OLx-1
532            jMin = 1-OLy+2
533            jMax = sNy+OLy-1
534    
535    C--     Explicit part of the Surface Pressure Gradient (add in TIMESTEP)
536    C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
537            IF (implicSurfPress.NE.1.) THEN
538              DO j=jMin,jMax
539                DO i=iMin,iMax
540                  phiSurfX(i,j) = _recip_dxC(i,j,bi,bj)*gBaro
541         &           *(etaN(i,j,bi,bj)-etaN(i-1,j,bi,bj))
542                  phiSurfY(i,j) = _recip_dyC(i,j,bi,bj)*gBaro
543         &           *(etaN(i,j,bi,bj)-etaN(i,j-1,bi,bj))
544                ENDDO
545              ENDDO
546            ENDIF
547    
548           IF (momStepping) THEN  C--     Start of dynamics loop
549  #ifdef ALLOW_AUTODIFF_TAMC          DO k=1,Nr
550           idkey = iikey + 3  
551  #endif  C--       km1    Points to level above k (=k-1)
552    C--       kup    Cycles through 1,2 to point to layer above
553    C--       kDown  Cycles through 2,1 to point to current layer
554    
555              km1  = MAX(1,k-1)
556              kup  = 1+MOD(k+1,2)
557              kDown= 1+MOD(k,2)
558    
559    C--      Integrate hydrostatic balance for phiHyd with BC of
560    C        phiHyd(z=0)=0
561    C        distinguishe between Stagger and Non Stagger time stepping
562             IF (staggerTimeStep) THEN
563               CALL CALC_PHI_HYD(
564         I        bi,bj,iMin,iMax,jMin,jMax,k,
565         I        gTnm1, gSnm1,
566         U        phiHyd,
567         I        myThid )
568             ELSE
569               CALL CALC_PHI_HYD(
570         I        bi,bj,iMin,iMax,jMin,jMax,k,
571         I        theta, salt,
572         U        phiHyd,
573         I        myThid )
574             ENDIF
575    
576    C--      Calculate accelerations in the momentum equations (gU, gV, ...)
577    C        and step forward storing the result in gUnm1, gVnm1, etc...
578             IF ( momStepping ) THEN
579               CALL CALC_MOM_RHS(
580         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
581         I         phiHyd,KappaRU,KappaRV,
582         U         fVerU, fVerV,
583         I         myTime, myThid)
584               CALL TIMESTEP(
585         I         bi,bj,iMin,iMax,jMin,jMax,k,
586         I         phiHyd, phiSurfX, phiSurfY,
587         I         myIter, myThid)
588    
589    #ifdef   ALLOW_OBCS
590    C--      Apply open boundary conditions
591             IF (useOBCS) THEN
592               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
593             END IF
594    #endif   /* ALLOW_OBCS */
595    
596    #ifdef   ALLOW_AUTODIFF_TAMC
597    #ifdef   INCLUDE_CD_CODE
598             ELSE
599               DO j=1-OLy,sNy+OLy
600                 DO i=1-OLx,sNx+OLx
601                   guCD(i,j,k,bi,bj) = 0.0
602                   gvCD(i,j,k,bi,bj) = 0.0
603                 END DO
604               END DO
605    #endif   /* INCLUDE_CD_CODE */
606    #endif   /* ALLOW_AUTODIFF_TAMC */
607             ENDIF
608    
609    
610    C--     end of dynamics k loop (1:Nr)
611            ENDDO
612    
613    
614    
615    C--     Implicit viscosity
616            IF (implicitViscosity.AND.momStepping) THEN
617    #ifdef    ALLOW_AUTODIFF_TAMC
618              idkey = iikey + 3
619    #endif    /* ALLOW_AUTODIFF_TAMC */
620            CALL IMPLDIFF(            CALL IMPLDIFF(
621       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
622       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
623       U         gUNm1,       U         gUNm1,
624       I         myThid )       I         myThid )
625  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
626           idkey = iikey + 4            idkey = iikey + 4
627  #endif  #endif    /* ALLOW_AUTODIFF_TAMC */
628            CALL IMPLDIFF(            CALL IMPLDIFF(
629       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
630       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
631       U         gVNm1,       U         gVNm1,
632       I         myThid )       I         myThid )
633    
634  #ifdef INCLUDE_CD_CODE  #ifdef   ALLOW_OBCS
635    C--      Apply open boundary conditions
636             IF (useOBCS) THEN
637               DO K=1,Nr
638                 CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
639               ENDDO
640             END IF
641    #endif   /* ALLOW_OBCS */
642    
643  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    INCLUDE_CD_CODE
644           idkey = iikey + 5  #ifdef    ALLOW_AUTODIFF_TAMC
645  #endif            idkey = iikey + 5
646    #endif    /* ALLOW_AUTODIFF_TAMC */
647            CALL IMPLDIFF(            CALL IMPLDIFF(
648       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
649       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
650       U         vVelD,       U         vVelD,
651       I         myThid )       I         myThid )
652  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
653          idkey = iikey + 6            idkey = iikey + 6
654  #endif  #endif    /* ALLOW_AUTODIFF_TAMC */
655            CALL IMPLDIFF(            CALL IMPLDIFF(
656       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
657       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
658       U         uVelD,       U         uVelD,
659       I         myThid )       I         myThid )
660    #endif    /* INCLUDE_CD_CODE */
661  #endif  C--     End If implicitViscosity.AND.momStepping
662            ENDIF
          ENDIF ! momStepping  
         ENDIF ! implicitViscosity  
663    
664    Cjmc : add for phiHyd output <- but not working if multi tile per CPU
665    c       IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)
666    c    &  .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
667    c         WRITE(suff,'(I10.10)') myIter+1
668    c         CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)
669    c       ENDIF
670    Cjmc(end)
671    
672    #ifdef ALLOW_TIMEAVE
673            IF (taveFreq.GT.0.) THEN
674              CALL TIMEAVE_CUMULATE(phiHydtave, phiHyd, Nr,
675         I                              deltaTclock, bi, bj, myThid)
676              IF (ivdc_kappa.NE.0.) THEN
677                CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
678         I                              deltaTclock, bi, bj, myThid)
679              ENDIF
680            ENDIF
681    #endif /* ALLOW_TIMEAVE */
682    
683         ENDDO         ENDDO
684        ENDDO        ENDDO
685    
 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 )  
   
   
686        RETURN        RETURN
687        END        END
Legend:
Removed from v.1.49  
changed lines
  Added in v.1.64
  ViewVC Help
Powered by ViewVC 1.1.22