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