/[MITgcm]/MITgcm/model/src/dynamics.F
ViewVC logotype

Diff of /MITgcm/model/src/dynamics.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.56 by heimbach, Mon Jan 29 20:05:46 2001 UTC revision 1.66 by heimbach, Sun Mar 25 22:33:52 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
# Line 25  C     \================================= Line 26  C     \=================================
26  C     == Global variables ===  C     == Global variables ===
27  #include "SIZE.h"  #include "SIZE.h"
28  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #include "CG2D.h"  
29  #include "PARAMS.h"  #include "PARAMS.h"
30  #include "DYNVARS.h"  #include "DYNVARS.h"
31  #include "GRID.h"  #include "GRID.h"
# Line 33  C     == Global variables === Line 33  C     == Global variables ===
33  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
34  # include "tamc.h"  # include "tamc.h"
35  # include "tamc_keys.h"  # include "tamc_keys.h"
 # include "FFIELDS.h"  
 # ifdef ALLOW_KPP  
 #  include "KPP.h"  
 # endif  
 # ifdef ALLOW_GMREDI  
 #  include "GMREDI.h"  
 # endif  
36  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
37    
38    #ifdef ALLOW_KPP
39    # include "KPP.h"
40    #endif
41    
42    #ifdef ALLOW_TIMEAVE
43    #include "TIMEAVE_STATV.h"
44    #endif
45    
46  C     == Routine arguments ==  C     == Routine arguments ==
47  C     myTime - Current time in simulation  C     myTime - Current time in simulation
48  C     myIter - Current iteration number in simulation  C     myIter - Current iteration number in simulation
# Line 54  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                      Potential (=pressure/rho0) 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 anomaly.
73  C                      anomaly.  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
74  C     etaSurfX,      - Holds surface elevation gradient in X and Y.  C     phiSurfY             or geopotentiel (atmos) in X and Y direction
 C     etaSurfY  
75  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
76  C     KappaRS          (background + spatially varying, isopycnal term).  C     KappaRS          (background + spatially varying, isopycnal term).
77  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
# Line 99  C                      index into fVerTe Line 85  C                      index into fVerTe
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)  
88        _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL aTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL xTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL cTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL mTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL pTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL fZon    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL fMer    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
90        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
91        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
92        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
93        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
94        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
95        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rhokp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
96        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL buoyKM1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _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)  
99        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
100        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
101        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
# Line 130  C                      index into fVerTe Line 104  C                      index into fVerTe
104        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
105        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
106    
107  C This is currently also used by IVDC and Diagnostics  C This is currently used by IVDC and Diagnostics
 C #ifdef INCLUDE_CONVECT_CALL  
108        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
 C #endif  
109    
110        INTEGER iMin, iMax        INTEGER iMin, iMax
111        INTEGER jMin, jMax        INTEGER jMin, jMax
112        INTEGER bi, bj        INTEGER bi, bj
113        INTEGER i, j        INTEGER i, j
114        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
       LOGICAL BOTTOM_LAYER  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
       INTEGER    isbyte  
       PARAMETER( isbyte = 4 )  
   
       INTEGER act1, act2, act3, act4  
       INTEGER max1, max2, max3  
       INTEGER iikey, kkey  
       INTEGER maximpl  
 #endif /* ALLOW_AUTODIFF_TAMC */  
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  C---    The algorithm...  C---    The algorithm...
123  C  C
124  C       "Correction Step"  C       "Correction Step"
# Line 166  C       "Calculation of Gs" Line 133  C       "Calculation of Gs"
133  C       ===================  C       ===================
134  C       This is where all the accelerations and tendencies (ie.  C       This is where all the accelerations and tendencies (ie.
135  C       physics, parameterizations etc...) are calculated  C       physics, parameterizations etc...) are calculated
 C         rVel = sum_r ( div. u[n] )  
136  C         rho = rho ( theta[n], salt[n] )  C         rho = rho ( theta[n], salt[n] )
137  C         b   = b(rho, theta)  C         b   = b(rho, theta)
138  C         K31 = K31 ( rho )  C         K31 = K31 ( rho )
139  C         Gu[n] = Gu( u[n], v[n], rVel, b, ... )  C         Gu[n] = Gu( u[n], v[n], wVel, b, ... )
140  C         Gv[n] = Gv( u[n], v[n], rVel, b, ... )  C         Gv[n] = Gv( u[n], v[n], wVel, b, ... )
141  C         Gt[n] = Gt( theta[n], u[n], v[n], rVel, K31, ... )  C         Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
142  C         Gs[n] = Gs( salt[n], u[n], v[n], rVel, K31, ... )  C         Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
143  C  C
144  C       "Time-stepping" or "Prediction"  C       "Time-stepping" or "Prediction"
145  C       ================================  C       ================================
# Line 202  C--   dummy statement to end declaration Line 168  C--   dummy statement to end declaration
168        ikey = 1        ikey = 1
169  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
170    
   
171  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
172  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
173  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
# Line 214  C     uninitialised but inert locations. Line 179  C     uninitialised but inert locations.
179          yA(i,j)      = 0. _d 0          yA(i,j)      = 0. _d 0
180          uTrans(i,j)  = 0. _d 0          uTrans(i,j)  = 0. _d 0
181          vTrans(i,j)  = 0. _d 0          vTrans(i,j)  = 0. _d 0
         aTerm(i,j)   = 0. _d 0  
         xTerm(i,j)   = 0. _d 0  
         cTerm(i,j)   = 0. _d 0  
         mTerm(i,j)   = 0. _d 0  
         pTerm(i,j)   = 0. _d 0  
         fZon(i,j)    = 0. _d 0  
         fMer(i,j)    = 0. _d 0  
182          DO k=1,Nr          DO k=1,Nr
183           phiHyd (i,j,k)  = 0. _d 0           phiHyd(i,j,k)  = 0. _d 0
184           KappaRU(i,j,k) = 0. _d 0           KappaRU(i,j,k) = 0. _d 0
185           KappaRV(i,j,k) = 0. _d 0           KappaRV(i,j,k) = 0. _d 0
186           sigmaX(i,j,k) = 0. _d 0           sigmaX(i,j,k) = 0. _d 0
# Line 231  C     uninitialised but inert locations. Line 189  C     uninitialised but inert locations.
189          ENDDO          ENDDO
190          rhoKM1 (i,j) = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
191          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  
192          maskC  (i,j) = 0. _d 0          maskC  (i,j) = 0. _d 0
193            phiSurfX(i,j) = 0. _d 0
194            phiSurfY(i,j) = 0. _d 0
195         ENDDO         ENDDO
196        ENDDO        ENDDO
197    
# Line 249  CHPF$ INDEPENDENT Line 205  CHPF$ INDEPENDENT
205    
206  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
207  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
208  CHPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV
209  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA
210  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV
211  CHPF$&                  )  CHPF$&                  )
# Line 278  C--     Set up work arrays that need val Line 234  C--     Set up work arrays that need val
234          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
235           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
236            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  
237            fVerT (i,j,1) = 0. _d 0            fVerT (i,j,1) = 0. _d 0
238            fVerT (i,j,2) = 0. _d 0            fVerT (i,j,2) = 0. _d 0
239            fVerS (i,j,1) = 0. _d 0            fVerS (i,j,1) = 0. _d 0
# Line 288  C--     Set up work arrays that need val Line 242  C--     Set up work arrays that need val
242            fVerU (i,j,2) = 0. _d 0            fVerU (i,j,2) = 0. _d 0
243            fVerV (i,j,1) = 0. _d 0            fVerV (i,j,1) = 0. _d 0
244            fVerV (i,j,2) = 0. _d 0            fVerV (i,j,2) = 0. _d 0
           phiHyd(i,j,1) = 0. _d 0  
245           ENDDO           ENDDO
246          ENDDO          ENDDO
247    
248          DO k=1,Nr          DO k=1,Nr
249           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
250            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
251  #ifdef INCLUDE_CONVECT_CALL  C This is currently also used by IVDC and Diagnostics
252             ConvectCount(i,j,k) = 0.             ConvectCount(i,j,k) = 0.
 #endif  
253             KappaRT(i,j,k) = 0. _d 0             KappaRT(i,j,k) = 0. _d 0
254             KappaRS(i,j,k) = 0. _d 0             KappaRS(i,j,k) = 0. _d 0
255            ENDDO            ENDDO
# Line 309  C--     Set up work arrays that need val Line 261  C--     Set up work arrays that need val
261          jMin = 1-OLy+1          jMin = 1-OLy+1
262          jMax = sNy+OLy          jMax = sNy+OLy
263    
         k = 1  
         BOTTOM_LAYER = k .EQ. Nr  
   
 #ifdef DO_PIPELINED_CORRECTION_STEP  
 C--     Calculate gradient of surface pressure  
         CALL CALC_GRAD_ETA_SURF(  
      I       bi,bj,iMin,iMax,jMin,jMax,  
      O       etaSurfX,etaSurfY,  
      I       myThid)  
 C--     Update fields in top level according to tendency terms  
         CALL CORRECTION_STEP(  
      I       bi,bj,iMin,iMax,jMin,jMax,k,  
      I       etaSurfX,etaSurfY,myTime,myThid)  
   
 #ifdef ALLOW_OBCS  
         IF (openBoundaries) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
            CALL APPLY_OBCS1( bi, bj, k, myThid )  
         END IF  
 #endif  
   
         IF ( .NOT. BOTTOM_LAYER ) THEN  
 C--      Update fields in layer below according to tendency terms  
          CALL CORRECTION_STEP(  
      I        bi,bj,iMin,iMax,jMin,jMax,k+1,  
      I        etaSurfX,etaSurfY,myTime,myThid)  
 #ifdef ALLOW_OBCS  
          IF (openBoundaries) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL APPLY_OBCS1( bi, bj, k+1, myThid )  
          END IF  
 #endif  
         ENDIF  
 #endif  
   
 C--     Density of 1st level (below W(1)) reference to level 1  
 #ifdef  INCLUDE_FIND_RHO_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         CALL FIND_RHO(  
      I     bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  
      O     rhoKm1,  
      I     myThid )  
 #endif  
   
         IF (.NOT. BOTTOM_LAYER) THEN  
   
 C--      Check static stability with layer below  
 C--      and mix as needed.  
 #ifdef  INCLUDE_FIND_RHO_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k+1,bi,bj) = comlev1_bibj  
 CADJ &   , key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,k+1,bi,bj) = comlev1_bibj  
 CADJ &   , key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
          CALL FIND_RHO(  
      I      bi, bj, iMin, iMax, jMin, jMax, k+1, k, eosType,  
      O      rhoKp1,  
      I      myThid )  
 #endif  
264    
265  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
266  CADJ STORE rhoKm1(:,:)  = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
267  CADJ STORE rhoKp1(:,:)  = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
268  #endif /* ALLOW_AUTODIFF_TAMC */  CADJ STORE uvel(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
269    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
270  #ifdef  INCLUDE_CONVECT_CALL  #endif /* ALLOW_AUTODIFF_TAMC */
271    
272           CALL CONVECT(  C--     Start of diagnostic loop
273       I       bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,          DO k=Nr,1,-1
274       U       ConvectCount,  
275       I       myTime,myIter,myThid)  #ifdef ALLOW_AUTODIFF_TAMC
276    C? Patrick, is this formula correct now that we change the loop range?
277  #ifdef ALLOW_AUTODIFF_TAMC  C? Do we still need this?
278  CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)  cph kkey formula corrected.
279  CADJ &     = comlev1_bibj, key = ikey, byte = isbyte  cph Needed for rhok, rhokm1, in the case useGMREDI.
280  CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)           kkey = (ikey-1)*Nr + k
281  CADJ &     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE rhokm1(:,:) = comlev1_bibj_k , key = kkey, byte = isbyte
282  #endif /* ALLOW_AUTODIFF_TAMC */  CADJ STORE rhok  (:,:) = comlev1_bibj_k , key = kkey, byte = isbyte
283    #endif /* ALLOW_AUTODIFF_TAMC */
284  #endif  
285    C--       Integrate continuity vertically for vertical velocity
286  C--      Implicit Vertical Diffusion for Convection            CALL INTEGRATE_FOR_W(
287           IF (ivdc_kappa.NE.0.) THEN       I                         bi, bj, k, uVel, vVel,
288              CALL CALC_IVDC(       O                         wVel,
289       I       bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,       I                         myThid )
290       U       ConvectCount, KappaRT, KappaRS,  
291       I       myTime,myIter,myThid)  #ifdef    ALLOW_OBCS
292           ENDIF  #ifdef    ALLOW_NONHYDROSTATIC
293    C--       Apply OBC to W if in N-H mode
294  C--      Recompute density after mixing            IF (useOBCS.AND.nonHydrostatic) THEN
295  #ifdef  INCLUDE_FIND_RHO_CALL              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
296           CALL FIND_RHO(            ENDIF
297       I      bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  #endif    /* ALLOW_NONHYDROSTATIC */
298       O      rhoKm1,  #endif    /* ALLOW_OBCS */
299       I      myThid )  
300  #endif  C--       Calculate gradients of potential density for isoneutral
301          ENDIF  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  C--     Calculate buoyancy            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
304          CALL CALC_BUOYANCY(              CALL FIND_RHO(
305       I      bi,bj,iMin,iMax,jMin,jMax,k,rhoKm1,       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
306       O      buoyKm1,       I        theta, salt,
307       I      myThid )       O        rhoK,
   
 C--     Integrate hydrostatic balance for phiHyd with BC of  
 C--     phiHyd(z=0)=0  
         CALL CALC_PHI_HYD(  
      I      bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyKm1,  
      U      phiHyd,  
      I      myThid )  
   
 #ifdef ALLOW_GMREDI  
         IF ( useGMRedi ) THEN  
         CALL GRAD_SIGMA(  
      I            bi, bj, iMin, iMax, jMin, jMax, k,  
      I            rhoKm1, rhoKm1, rhoKm1,  
      O            sigmaX, sigmaY, sigmaR,  
      I            myThid )  
         ELSE  
          DO j=1-OLy,sNy+OLy  
           DO i=1-OLx,sNx+OLx  
            sigmaX(i,j,k) = 0. _d 0  
            sigmaY(i,j,k) = 0. _d 0  
            sigmaR(i,j,k) = 0. _d 0  
           ENDDO  
          ENDDO  
         ENDIF  
 #endif  
   
 C--     Start of downward loop  
         DO k=2,Nr  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
          kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
          BOTTOM_LAYER = k .EQ. Nr  
   
 #ifdef DO_PIPELINED_CORRECTION_STEP  
          IF ( .NOT. BOTTOM_LAYER ) THEN  
 C--       Update fields in layer below according to tendency terms  
           CALL CORRECTION_STEP(  
      I         bi,bj,iMin,iMax,jMin,jMax,k+1,  
      I         etaSurfX,etaSurfY,myTime,myThid)  
 #ifdef ALLOW_OBCS  
           IF (openBoundaries) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
              CALL APPLY_OBCS1( bi, bj, k+1, myThid )  
           END IF  
 #endif  
          ENDIF  
 #endif /* DO_PIPELINED_CORRECTION_STEP */  
   
 C--      Density of k level (below W(k)) reference to k level  
 #ifdef  INCLUDE_FIND_RHO_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
          CALL FIND_RHO(  
      I      bi, bj, iMin, iMax, jMin, jMax,  k, k, eosType,  
      O      rhoK,  
      I      myThid )  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph(   storing not necessary  
 cphCADJ STORE rhoK(:,:)  = comlev1_bibj_k, key = kkey, byte = isbyte  
 cph)  
 #endif /* ALLOW_AUTODIFF_TAMC */  
 #endif  
   
          IF (.NOT. BOTTOM_LAYER) THEN  
   
 C--       Check static stability with layer below and mix as needed.  
 C--       Density of k+1 level (below W(k+1)) reference to k level.  
 #ifdef  INCLUDE_FIND_RHO_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k+1,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 CADJ STORE salt (:,:,k+1,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
           CALL FIND_RHO(  
      I       bi, bj, iMin, iMax, jMin, jMax,  k+1, k, eosType,  
      O       rhoKp1,  
      I       myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhoKp1(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
 #endif  
   
 #ifdef  INCLUDE_CONVECT_CALL  
           CALL CONVECT(  
      I        bi,bj,iMin,iMax,jMin,jMax,k+1,rhoK,rhoKp1,  
      U        ConvectCount,  
      I        myTime,myIter,myThid)  
   
 #endif  
   
 C--      Implicit Vertical Diffusion for Convection  
          IF (ivdc_kappa.NE.0.) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhoKm1(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL CALC_IVDC(  
      I       bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,  
      U       ConvectCount, KappaRT, KappaRS,  
      I       myTime,myIter,myThid)  
          END IF  
   
 C--       Recompute density after mixing  
 #ifdef  INCLUDE_FIND_RHO_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
           CALL FIND_RHO(  
      I       bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  
      O       rhoK,  
      I       myThid )  
 #endif  
   
 C--            IF (.NOT. BOTTOM_LAYER) ends here  
          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,  
308       I        myThid )       I        myThid )
309                IF (k.GT.1) CALL FIND_RHO(
 #ifdef  INCLUDE_FIND_RHO_CALL  
 C--      Calculate iso-neutral slopes for the GM/Redi parameterisation  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
          CALL FIND_RHO(  
310       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
311       O        rhoTmp,       I        theta, salt,
312         O        rhoKm1,
313       I        myThid )       I        myThid )
314  #endif              CALL GRAD_SIGMA(
   
   
 #ifdef ALLOW_GMREDI  
          IF ( useGMRedi ) THEN  
          CALL GRAD_SIGMA(  
315       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
316       I             rhoK, rhotmp, rhoK,       I             rhoK, rhoKm1, rhoK,
317       O             sigmaX, sigmaY, sigmaR,       O             sigmaX, sigmaY, sigmaR,
318       I             myThid )       I             myThid )
319           ELSE            ENDIF
           DO j=1-OLy,sNy+OLy  
            DO i=1-OLx,sNx+OLx  
             sigmaX(i,j,k) = 0. _d 0  
             sigmaY(i,j,k) = 0. _d 0  
             sigmaR(i,j,k) = 0. _d 0  
            ENDDO  
           ENDDO  
          ENDIF  
 #endif  
320    
321           DO J=jMin,jMax  C--       Implicit Vertical Diffusion for Convection
322            DO I=iMin,iMax  c ==> should use sigmaR !!!
323  #ifdef  INCLUDE_FIND_RHO_CALL            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
324             rhoKm1 (I,J) = rhoK(I,J)              CALL CALC_IVDC(
325  #endif       I        bi, bj, iMin, iMax, jMin, jMax, k,
326             buoyKm1(I,J) = buoyK(I,J)       I        rhoKm1, rhoK,
327            ENDDO       U        ConvectCount, KappaRT, KappaRS,
328           ENDDO       I        myTime, myIter, myThid)
329              ENDIF
330    
331  C--     end of k loop  C--     end of diagnostic k loop (Nr:1)
332          ENDDO          ENDDO
333    
334  C     Determines forcing terms based on external fields  #ifdef  ALLOW_OBCS
335  C     relaxation terms, etc.  C--     Calculate future values on open boundaries
336        CALL EXTERNAL_FORCING_SURF(          IF (useOBCS) THEN
337              CALL OBCS_CALC( bi, bj, myTime+deltaT,
338         I            uVel, vVel, wVel, theta, salt,
339         I            myThid )
340            ENDIF
341    #endif  /* ALLOW_OBCS */
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,       I             bi, bj, iMin, iMax, jMin, jMax,
347       I             myThid )       I             myThid )
348    
349  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef  ALLOW_GMREDI
350    C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
 CADJ STORE surfacetendencyu(:,:,bi,bj)  
 CADJ &   , surfacetendencyv(:,:,bi,bj)  
 CADJ &   , surfacetendencys(:,:,bi,bj)  
 CADJ &   , surfacetendencyt(:,:,bi,bj)  
 CADJ &                        = comlev1_bibj, key=ikey, byte=isbyte  
   
 # ifdef ALLOW_GMREDI  
 CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  
 CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  
 CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  
 # endif /* ALLOW_GMREDI */  
   
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef ALLOW_GMREDI  
351          IF (useGMRedi) THEN          IF (useGMRedi) THEN
352            DO k=1, Nr            DO k=1,Nr
353              CALL GMREDI_CALC_TENSOR(              CALL GMREDI_CALC_TENSOR(
354       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
355       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
# Line 661  CADJ STORE sigmaR(:,:,:) = comlev1, key= Line 365  CADJ STORE sigmaR(:,:,:) = comlev1, key=
365            ENDDO            ENDDO
366  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
367          ENDIF          ENDIF
368  #endif  #endif  /* ALLOW_GMREDI */
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  
   
 #ifdef ALLOW_GMREDI  
 C-- R.G. We need to define a new tape since Kw use mythid instead of bi,bj  
 CADJ STORE Kwx(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwy(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwz(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte  
 #endif  
369    
370  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  #ifdef  ALLOW_KPP
371  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  C--     Compute KPP mixing coefficients
 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 =======  
 C-- R.G. We need to define a new tape since Kw use mythid instead of bi,bj  
 CADJ STORE Kwx(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwy(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwz(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte  
   
 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
   
 C--     dummy initialization to break data flow because  
 C--     calc_div_ghat has a condition for initialization  
         DO J=jMin,jMax  
            DO I=iMin,iMax  
               cg2d_b(i,j,bi,bj) = 0.0  
            ENDDO  
         ENDDO  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef ALLOW_KPP  
 C--   Compute KPP mixing coefficients  
372          IF (useKPP) THEN          IF (useKPP) THEN
   
           CALL TIMER_START('KPP_CALC               [DYNAMICS]', myThid)  
373            CALL KPP_CALC(            CALL KPP_CALC(
374       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myThid )
           CALL TIMER_STOP ('KPP_CALC               [DYNAMICS]', myThid)  
   
375  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
376          ELSE          ELSE
377            DO j=1-OLy,sNy+OLy            DO j=1-OLy,sNy+OLy
# Line 733  CADJ &   , KPPfrac   (:,:  ,bi,bj) Line 398  CADJ &   , KPPfrac   (:,:  ,bi,bj)
398  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte
399  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
400    
401  #endif /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
402    
403  C--     Start of upward loop  #ifdef ALLOW_AUTODIFF_TAMC
404          DO k = Nr, 1, -1  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
405    CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
406    CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
407    CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
408    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
409    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
410    #endif /* ALLOW_AUTODIFF_TAMC */
411    
412    #ifdef ALLOW_AIM
413    C       AIM - atmospheric intermediate model, physics package code.
414    C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics
415            IF ( useAIM ) THEN
416             CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
417             CALL AIM_DO_ATMOS_PHYSICS( phiHyd, myTime, myThid )
418             CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
419            ENDIF
420    #endif /* ALLOW_AIM */
421    
422    
423    C--     Start of thermodynamics loop
424            DO k=Nr,1,-1
425    
426  C--      km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
427  C--      kup    Cycles through 1,2 to point to layer above  C--       kup    Cycles through 1,2 to point to layer above
428  C--      kDown  Cycles through 2,1 to point to current layer  C--       kDown  Cycles through 2,1 to point to current layer
429    
430           km1  =max(1,k-1)            km1  = MAX(1,k-1)
431           kup  =1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
432           kDown=1+MOD(k,2)            kDown= 1+MOD(k,2)
433    
434           iMin = 1-OLx+2            iMin = 1-OLx+2
435           iMax = sNx+OLx-1            iMax = sNx+OLx-1
436           jMin = 1-OLy+2            jMin = 1-OLy+2
437           jMax = sNy+OLy-1            jMax = sNy+OLy-1
438    
439  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
440           kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1  C? Patrick Is this formula correct?
441  CADJ STORE rvel  (:,:,kdown)  = comlev1_bibj_k, key=kkey, byte=isbyte  cph Yes, but I rewrote it.
442  CADJ STORE rTrans(:,:)        = comlev1_bibj_k, key=kkey, byte=isbyte  cph Also, the KappaR? need the index k!
443  CADJ STORE KappaRT(:,:,k)     = comlev1_bibj_k, key=kkey, byte=isbyte           kkey = (ikey-1)*Nr + k
444  CADJ STORE KappaRS(:,:,k)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key = kkey, byte = isbyte
445    CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key = kkey, byte = isbyte
446  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
447    
448  C--      Get temporary terms used by tendency routines  C--      Get temporary terms used by tendency routines
449           CALL CALC_COMMON_FACTORS (           CALL CALC_COMMON_FACTORS (
450       I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,       I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
451       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,       O        xA,yA,uTrans,vTrans,rTrans,maskC,maskUp,
452       I        myThid)       I        myThid)
453    
 #ifdef ALLOW_OBCS  
         IF (openBoundaries) THEN  
          CALL APPLY_OBCS3( bi, bj, k, kup, rTrans, rVel, myThid )  
         ENDIF  
 #endif  
   
454  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
455  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
456           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
457       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
458       I        maskC,maskUp,       I        maskC,maskup,
459       O        KappaRT,KappaRS,KappaRU,KappaRV,       O        KappaRT,KappaRS,KappaRU,KappaRV,
460       I        myThid)       I        myThid)
461  #endif  #endif
462  C--      Calculate accelerations in the momentum equations  
463           IF ( momStepping ) THEN  C--      Calculate active tracer tendencies (gT,gS,...)
464            CALL CALC_MOM_RHS(  C        and step forward storing result in gTnm1, gSnm1, etc.
      I         bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,rVel,maskC,  
      I         phiHyd,KappaRU,KappaRV,  
      U         aTerm,xTerm,cTerm,mTerm,pTerm,  
      U         fZon, fMer, fVerU, fVerV,  
      I         myTime, myThid)  
 #ifdef ALLOW_AUTODIFF_TAMC  
 #ifdef INCLUDE_CD_CODE  
          ELSE  
             DO j=1-OLy,sNy+OLy  
                DO i=1-OLx,sNx+OLx  
                   guCD(i,j,k,bi,bj) = 0.0  
                   gvCD(i,j,k,bi,bj) = 0.0  
                END DO  
             END DO  
 #endif  
 #endif /* ALLOW_AUTODIFF_TAMC */  
          ENDIF  
 C--      Calculate active tracer tendencies  
465           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
466            CALL CALC_GT(             CALL CALC_GT(
467       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
468       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
469       I         KappaRT,       I         KappaRT,
470       U         aTerm,xTerm,fZon,fMer,fVerT,       U         fVerT,
471       I         myTime, myThid)       I         myTime, myThid)
472               CALL TIMESTEP_TRACER(
473         I         bi,bj,iMin,iMax,jMin,jMax,k,
474         I         theta, gT,
475         U         gTnm1,
476         I         myIter, myThid)
477           ENDIF           ENDIF
478           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
479            CALL CALC_GS(             CALL CALC_GS(
480       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
481       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
482       I         KappaRS,       I         KappaRS,
483       U         aTerm,xTerm,fZon,fMer,fVerS,       U         fVerS,
484       I         myTime, myThid)       I         myTime, myThid)
485               CALL TIMESTEP_TRACER(
486         I         bi,bj,iMin,iMax,jMin,jMax,k,
487         I         salt, gS,
488         U         gSnm1,
489         I         myIter, myThid)
490           ENDIF           ENDIF
 #ifdef ALLOW_OBCS  
 C--      Calculate future values on open boundaries  
          IF (openBoundaries) THEN  
 Caja      CALL CYCLE_OBCS( k, bi, bj, myThid )  
           CALL SET_OBCS( k, bi, bj, myTime+deltaTclock, myThid )  
          ENDIF  
 #endif  
 C--      Prediction step (step forward all model variables)  
          CALL TIMESTEP(  
      I       bi,bj,iMin,iMax,jMin,jMax,k,  
      I       myIter, myThid)  
 #ifdef ALLOW_OBCS  
 C--      Apply open boundary conditions  
          IF (openBoundaries) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE gunm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE gvnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE gwnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
491    
492              CALL APPLY_OBCS2( bi, bj, k, myThid )  #ifdef   ALLOW_OBCS
493    C--      Apply open boundary conditions
494             IF (useOBCS) THEN
495               CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
496           END IF           END IF
497  #endif  #endif   /* ALLOW_OBCS */
498    
499  C--      Freeze water  C--      Freeze water
500           IF (allowFreezing) THEN           IF (allowFreezing) THEN
501  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
502  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
503    CADJ &   , key = kkey, byte = isbyte
504  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
505              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
506           END IF           END IF
507    
508  #ifdef DIVG_IN_DYNAMICS  C--     end of thermodynamic k loop (Nr:1)
 C--      Diagnose barotropic divergence of predicted fields  
          CALL CALC_DIV_GHAT(  
      I       bi,bj,iMin,iMax,jMin,jMax,k,  
      I       xA,yA,  
      I       myThid)  
 #endif /* DIVG_IN_DYNAMICS */  
   
 C--      Cumulative diagnostic calculations (ie. time-averaging)  
 #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE  
          IF (taveFreq.GT.0.) THEN  
           CALL DO_TIME_AVERAGES(  
      I                           myTime, myIter, bi, bj, k, kup, kDown,  
      I                           rVel, ConvectCount,  
      I                           myThid )  
          ENDIF  
 #endif  
   
   
 C--     k loop  
509          ENDDO          ENDDO
510    
511    
512  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
513             maximpl = 6  C? Patrick? What about this one?
514             iikey = (ikey-1)*maximpl  cph Keys iikey and idkey don't seem to be needed
515    cph since storing occurs on different tape for each
516    cph impldiff call anyways.
517    cph Thus, common block comlev1_impl isn't needed either.
518    cph Storing below needed in the case useGMREDI.
519            iikey = (ikey-1)*maximpl
520  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
521    
522  C--     Implicit diffusion  C--     Implicit diffusion
# Line 885  CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_ Line 529  CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_
529  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
530              CALL IMPLDIFF(              CALL IMPLDIFF(
531       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
532       I         deltaTtracer, KappaRT,recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
533       U         gTNm1,       U         gTNm1,
534       I         myThid )       I         myThid )
535           END IF           ENDIF
536    
537           IF (saltStepping) THEN           IF (saltStepping) THEN
538  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 897  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_ Line 541  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_
541  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
542              CALL IMPLDIFF(              CALL IMPLDIFF(
543       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
544       I         deltaTtracer, KappaRS,recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
545       U         gSNm1,       U         gSNm1,
546       I         myThid )       I         myThid )
547             ENDIF
548    
549    #ifdef   ALLOW_OBCS
550    C--      Apply open boundary conditions
551             IF (useOBCS) THEN
552               DO K=1,Nr
553                 CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
554               ENDDO
555           END IF           END IF
556    #endif   /* ALLOW_OBCS */
557    
558  C--     implicitDiffusion  C--     End If implicitDiffusion
559          ENDIF          ENDIF
560    
561  C--     Implicit viscosity  C--     Start computation of dynamics
562          IF (implicitViscosity) THEN          iMin = 1-OLx+2
563            iMax = sNx+OLx-1
564            jMin = 1-OLy+2
565            jMax = sNy+OLy-1
566    
567    C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
568    C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
569            IF (implicSurfPress.NE.1.) THEN
570              CALL CALC_GRAD_PHI_SURF(
571         I         bi,bj,iMin,iMax,jMin,jMax,
572         I         etaN,
573         O         phiSurfX,phiSurfY,
574         I         myThid )                        
575            ENDIF
576    
577           IF (momStepping) THEN  C--     Start of dynamics loop
578  #ifdef ALLOW_AUTODIFF_TAMC          DO k=1,Nr
579           idkey = iikey + 3  
580    C--       km1    Points to level above k (=k-1)
581    C--       kup    Cycles through 1,2 to point to layer above
582    C--       kDown  Cycles through 2,1 to point to current layer
583    
584              km1  = MAX(1,k-1)
585              kup  = 1+MOD(k+1,2)
586              kDown= 1+MOD(k,2)
587    
588    C--      Integrate hydrostatic balance for phiHyd with BC of
589    C        phiHyd(z=0)=0
590    C        distinguishe between Stagger and Non Stagger time stepping
591             IF (staggerTimeStep) THEN
592               CALL CALC_PHI_HYD(
593         I        bi,bj,iMin,iMax,jMin,jMax,k,
594         I        gTnm1, gSnm1,
595         U        phiHyd,
596         I        myThid )
597             ELSE
598               CALL CALC_PHI_HYD(
599         I        bi,bj,iMin,iMax,jMin,jMax,k,
600         I        theta, salt,
601         U        phiHyd,
602         I        myThid )
603             ENDIF
604    
605    C--      Calculate accelerations in the momentum equations (gU, gV, ...)
606    C        and step forward storing the result in gUnm1, gVnm1, etc...
607             IF ( momStepping ) THEN
608               CALL CALC_MOM_RHS(
609         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
610         I         phiHyd,KappaRU,KappaRV,
611         U         fVerU, fVerV,
612         I         myTime, myThid)
613               CALL TIMESTEP(
614         I         bi,bj,iMin,iMax,jMin,jMax,k,
615         I         phiHyd, phiSurfX, phiSurfY,
616         I         myIter, myThid)
617    
618    #ifdef   ALLOW_OBCS
619    C--      Apply open boundary conditions
620             IF (useOBCS) THEN
621               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
622             END IF
623    #endif   /* ALLOW_OBCS */
624    
625    #ifdef   ALLOW_AUTODIFF_TAMC
626    #ifdef   INCLUDE_CD_CODE
627             ELSE
628               DO j=1-OLy,sNy+OLy
629                 DO i=1-OLx,sNx+OLx
630                   guCD(i,j,k,bi,bj) = 0.0
631                   gvCD(i,j,k,bi,bj) = 0.0
632                 END DO
633               END DO
634    #endif   /* INCLUDE_CD_CODE */
635    #endif   /* ALLOW_AUTODIFF_TAMC */
636             ENDIF
637    
638    
639    C--     end of dynamics k loop (1:Nr)
640            ENDDO
641    
642    
643    
644    C--     Implicit viscosity
645            IF (implicitViscosity.AND.momStepping) THEN
646    #ifdef    ALLOW_AUTODIFF_TAMC
647              idkey = iikey + 3
648  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
649  #endif /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
650            CALL IMPLDIFF(            CALL IMPLDIFF(
651       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
652       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
653       U         gUNm1,       U         gUNm1,
654       I         myThid )       I         myThid )
655  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
656           idkey = iikey + 4            idkey = iikey + 4
657  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
658  #endif /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
659            CALL IMPLDIFF(            CALL IMPLDIFF(
660       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
661       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
662       U         gVNm1,       U         gVNm1,
663       I         myThid )       I         myThid )
664    
665  #ifdef INCLUDE_CD_CODE  #ifdef   ALLOW_OBCS
666    C--      Apply open boundary conditions
667             IF (useOBCS) THEN
668               DO K=1,Nr
669                 CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
670               ENDDO
671             END IF
672    #endif   /* ALLOW_OBCS */
673    
674  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    INCLUDE_CD_CODE
675           idkey = iikey + 5  #ifdef    ALLOW_AUTODIFF_TAMC
676              idkey = iikey + 5
677  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
678  #endif /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
679            CALL IMPLDIFF(            CALL IMPLDIFF(
680       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
681       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
682       U         vVelD,       U         vVelD,
683       I         myThid )       I         myThid )
684  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
685          idkey = iikey + 6            idkey = iikey + 6
686  CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
687  #endif /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
688            CALL IMPLDIFF(            CALL IMPLDIFF(
689       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
690       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
691       U         uVelD,       U         uVelD,
692       I         myThid )       I         myThid )
693    #endif    /* INCLUDE_CD_CODE */
694  #endif  C--     End If implicitViscosity.AND.momStepping
   
 C--      momStepping  
          ENDIF  
   
 C--     implicitViscosity  
695          ENDIF          ENDIF
696    
697    Cjmc : add for phiHyd output <- but not working if multi tile per CPU
698    c       IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)
699    c    &  .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
700    c         WRITE(suff,'(I10.10)') myIter+1
701    c         CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)
702    c       ENDIF
703    Cjmc(end)
704    
705    #ifdef ALLOW_TIMEAVE
706            IF (taveFreq.GT.0.) THEN
707              CALL TIMEAVE_CUMULATE(phiHydtave, phiHyd, Nr,
708         I                              deltaTclock, bi, bj, myThid)
709              IF (ivdc_kappa.NE.0.) THEN
710                CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
711         I                              deltaTclock, bi, bj, myThid)
712              ENDIF
713            ENDIF
714    #endif /* ALLOW_TIMEAVE */
715    
716         ENDDO         ENDDO
717        ENDDO        ENDDO
718    
Legend:
Removed from v.1.56  
changed lines
  Added in v.1.66
  ViewVC Help
Powered by ViewVC 1.1.22