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