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