/[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.61 by jmc, Wed Feb 7 21:48:02 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  #endif
42    
43  C     == Routine arguments ==  C     == Routine arguments ==
# Line 53  C     == Local variables Line 52  C     == Local variables
52  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
53  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C     uTrans, vTrans, rTrans - Per block temporaries holding flow
54  C                              transport  C                              transport
55  C     rVel                     o uTrans: Zonal transport  C                              o uTrans: Zonal transport
56  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
57  C                              o rTrans: Vertical transport  C                              o rTrans: Vertical transport
 C                              o rVel:   Vertical velocity at upper and  
 C                                        lower cell faces.  
58  C     maskC,maskUp             o maskC: land/water mask for tracer cells  C     maskC,maskUp             o maskC: land/water mask for tracer cells
59  C                              o maskUp: land/water mask for W points  C                              o maskUp: land/water mask for W points
60  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  
61  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
62  C                                      so we need an fVer for each  C                                      so we need an fVer for each
63  C                                      variable.  C                                      variable.
64  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.  
65  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     phiHyd         - Hydrostatic part of the potential phiHydi.
66  C                      In z coords phiHydiHyd is the hydrostatic  C                      In z coords phiHydiHyd is the hydrostatic
67  C                      pressure anomaly  C                      pressure anomaly
# Line 90  C     KappaRS          (background + spa Line 75  C     KappaRS          (background + spa
75  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
76  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
77  C     bi, bj  C     bi, bj
78  C     k, kUp,        - Index for layer above and below. kUp and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
79  C     kDown, kM1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
80  C                      index into fVerTerm.  C                      index into fVerTerm.
81        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85        _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)  
86        _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87        _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)  
88        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
89        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
90        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
91        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
92        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
93        _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)  
94        _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)  
95        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
96        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
97        _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 100  C                      index into fVerTe
100        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
102    
103  #ifdef INCLUDE_CONVECT_CALL  C This is currently also used by IVDC and Diagnostics
104    C #ifdef INCLUDE_CONVECT_CALL
105        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
106  #endif  C #endif
107    
108        INTEGER iMin, iMax        INTEGER iMin, iMax
109        INTEGER jMin, jMax        INTEGER jMin, jMax
110        INTEGER bi, bj        INTEGER bi, bj
111        INTEGER i, j        INTEGER i, j
112        INTEGER k, kM1, kUp, kDown        INTEGER k, km1, kup, kDown
       LOGICAL BOTTOM_LAYER  
113    
114  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
115        INTEGER    isbyte        INTEGER    isbyte
# Line 146  C                      index into fVerTe Line 117  C                      index into fVerTe
117    
118        INTEGER act1, act2, act3, act4        INTEGER act1, act2, act3, act4
119        INTEGER max1, max2, max3        INTEGER max1, max2, max3
120        INTEGER ikact, iikey,kkey        INTEGER iikey, kkey
121        INTEGER maximpl        INTEGER maximpl
122  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
123    
124  C---    The algorithm...  C---    The algorithm...
125  C  C
# Line 163  C Line 134  C
134  C       "Calculation of Gs"  C       "Calculation of Gs"
135  C       ===================  C       ===================
136  C       This is where all the accelerations and tendencies (ie.  C       This is where all the accelerations and tendencies (ie.
137  C       phiHydysics, parameterizations etc...) are calculated  C       physics, parameterizations etc...) are calculated
 C         rVel = sum_r ( div. u[n] )  
138  C         rho = rho ( theta[n], salt[n] )  C         rho = rho ( theta[n], salt[n] )
139  C         b   = b(rho, theta)  C         b   = b(rho, theta)
140  C         K31 = K31 ( rho )  C         K31 = K31 ( rho )
141  C         Gu[n] = Gu( u[n], v[n], rVel, b, ... )  C         Gu[n] = Gu( u[n], v[n], wVel, b, ... )
142  C         Gv[n] = Gv( u[n], v[n], rVel, b, ... )  C         Gv[n] = Gv( u[n], v[n], wVel, b, ... )
143  C         Gt[n] = Gt( theta[n], u[n], v[n], rVel, K31, ... )  C         Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
144  C         Gs[n] = Gs( salt[n], u[n], v[n], rVel, K31, ... )  C         Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
145  C  C
146  C       "Time-stepping" or "Prediction"  C       "Time-stepping" or "Prediction"
147  C       ================================  C       ================================
# Line 198  C--- Line 168  C---
168  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
169  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
170        ikey = 1        ikey = 1
171  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
172    
173  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
174  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 181  C     uninitialised but inert locations.
181          yA(i,j)      = 0. _d 0          yA(i,j)      = 0. _d 0
182          uTrans(i,j)  = 0. _d 0          uTrans(i,j)  = 0. _d 0
183          vTrans(i,j)  = 0. _d 0          vTrans(i,j)  = 0. _d 0
184          aTerm(i,j)   = 0. _d 0          DO k=1,Nr
185          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  
186           KappaRU(i,j,k) = 0. _d 0           KappaRU(i,j,k) = 0. _d 0
187           KappaRV(i,j,k) = 0. _d 0           KappaRV(i,j,k) = 0. _d 0
188           sigmaX(i,j,k) = 0. _d 0           sigmaX(i,j,k) = 0. _d 0
# Line 228  C     uninitialised but inert locations. Line 191  C     uninitialised but inert locations.
191          ENDDO          ENDDO
192          rhoKM1 (i,j) = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
193          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  
194          maskC  (i,j) = 0. _d 0          maskC  (i,j) = 0. _d 0
195         ENDDO         ENDDO
196        ENDDO        ENDDO
# Line 239  C     uninitialised but inert locations. Line 198  C     uninitialised but inert locations.
198    
199  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
200  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
201  !HPF$ INDEPENDENT  CHPF$ INDEPENDENT
202  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
203    
204        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
205    
206  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
207  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
208  !HPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV
209  !HPF$&                  ,phiHyd,  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA
210  !HPF$&                  ,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV
211  !HPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  )
212  !HPF$&                  )  #endif /* ALLOW_AUTODIFF_TAMC */
 #endif  
213    
214         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
215    
# Line 270  C--    HPF directive to help TAMC Line 228  C--    HPF directive to help TAMC
228            ikey = (act1 + 1) + act2*max1            ikey = (act1 + 1) + act2*max1
229       &                      + act3*max1*max2       &                      + act3*max1*max2
230       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
231  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
232    
233  C--     Set up work arrays that need valid initial values  C--     Set up work arrays that need valid initial values
234          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
235           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
236            rTrans(i,j)   = 0. _d 0            rTrans(i,j)   = 0. _d 0
           rVel  (i,j,1) = 0. _d 0  
           rVel  (i,j,2) = 0. _d 0  
237            fVerT (i,j,1) = 0. _d 0            fVerT (i,j,1) = 0. _d 0
238            fVerT (i,j,2) = 0. _d 0            fVerT (i,j,2) = 0. _d 0
239            fVerS (i,j,1) = 0. _d 0            fVerS (i,j,1) = 0. _d 0
# Line 286  C--     Set up work arrays that need val Line 242  C--     Set up work arrays that need val
242            fVerU (i,j,2) = 0. _d 0            fVerU (i,j,2) = 0. _d 0
243            fVerV (i,j,1) = 0. _d 0            fVerV (i,j,1) = 0. _d 0
244            fVerV (i,j,2) = 0. _d 0            fVerV (i,j,2) = 0. _d 0
           phiHyd(i,j,1) = 0. _d 0  
245           ENDDO           ENDDO
246          ENDDO          ENDDO
247    
# Line 308  C--     Set up work arrays that need val Line 263  C--     Set up work arrays that need val
263          jMax = sNy+OLy          jMax = sNy+OLy
264    
265    
266          K = 1  C--     Start of diagnostic loop
267          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  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
          kkey = ikact*(Nr-2+1) + (k-2) + 1  
 #endif  
   
          BOTTOM_LAYER = K .EQ. Nr  
   
 #ifdef DO_PIPELINED_CORRECTION_STEP  
          IF ( .NOT. BOTTOM_LAYER ) THEN  
 C--       Update fields in layer below according to tendency terms  
           CALL CORRECTION_STEP(  
      I         bi,bj,iMin,iMax,jMin,jMax,K+1,  
      I         etaSurfX,etaSurfY,myTime,myThid)  
 #ifdef ALLOW_OBCS  
           IF (openBoundaries) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE uvel (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE vvel (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE salt(:,:,k,bi,bj)   = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
              CALL APPLY_OBCS1( bi, bj, K+1, myThid )  
           END IF  
 #endif  
          ENDIF  
 #endif  
   
 C--      Density of K level (below W(K)) reference to K level  
 #ifdef  INCLUDE_FIND_RHO_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj)  = comlev1_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  
268    
269  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
270  CADJ STORE rhok  (:,:)   = comlev1_3d, key = kkey, byte = isbyte  C? Patrick, is this formula correct now that we change the loop range?
271  CADJ STORE rhoKm1(:,:)   = comlev1_3d, key = kkey, byte = isbyte  C? Do we still need this?
272  CADJ STORE rhoKp1(:,:)   = comlev1_3d, key = kkey, byte = isbyte           kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1
273  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
274    
275  #ifdef  INCLUDE_CONVECT_CALL  C--       Integrate continuity vertically for vertical velocity
276            CALL CONVECT(            CALL INTEGRATE_FOR_W(
277       I        bi,bj,iMin,iMax,jMin,jMax,K+1,rhoK,rhoKp1,       I                         bi, bj, k, uVel, vVel,
278       U        ConvectCount,       O                         wVel,
279       I        myTime,myIter,myThid)       I                         myThid )
280  #ifdef ALLOW_AUTODIFF_TAMC  
281  CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)  #ifdef    ALLOW_OBCS
282  CADJ &     = comlev1_3d, key = kkey, byte = isbyte  #ifdef    ALLOW_NONHYDROSTATIC
283  CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)  C--       Apply OBC to W if in N-H mode
284  CADJ &     = comlev1_3d, key = kkey, byte = isbyte            IF (useOBCS.AND.nonHydrostatic) THEN
285  #endif              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
286  #endif            ENDIF
287    #endif    /* ALLOW_NONHYDROSTATIC */
288  C--      Implicit Vertical Diffusion for Convection  #endif    /* ALLOW_OBCS */
289           IF (ivdc_kappa.NE.0.) THEN  
290              CALL CALC_IVDC(  C--       Calculate gradients of potential density for isoneutral
291       I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
292       U       ConvectCount, KappaRT, KappaRS,  c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
293       I       myTime,myIter,myThid)            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
294  CRG: do we need do store STORE KappaRT, KappaRS ?              CALL FIND_RHO(
295           END IF       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
296         I        theta, salt,
297  C--       Recompute density after mixing       O        rhoK,
 #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,  
298       I        myThid )       I        myThid )
299  C--      Calculate iso-neutral slopes for the GM/Redi parameterisation              IF (k.GT.1) CALL FIND_RHO(
300  #ifdef  INCLUDE_FIND_RHO_CALL       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
301           CALL FIND_RHO(       I        theta, salt,
302       I        bi, bj, iMin, iMax, jMin, jMax, K-1, K, eosType,       O        rhoKm1,
      O        rhoTmp,  
303       I        myThid )       I        myThid )
304  #endif              CALL GRAD_SIGMA(
305           CALL GRAD_SIGMA(       I             bi, bj, iMin, iMax, jMin, jMax, k,
306       I             bi, bj, iMin, iMax, jMin, jMax, K,       I             rhoK, rhoKm1, rhoK,
      I             rhoK, rhotmp, rhoK,  
307       O             sigmaX, sigmaY, sigmaR,       O             sigmaX, sigmaY, sigmaR,
308       I             myThid )       I             myThid )
309              ENDIF
310    
311    C--       Implicit Vertical Diffusion for Convection
312    c ==> should use sigmaR !!!
313              IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
314                CALL CALC_IVDC(
315         I        bi, bj, iMin, iMax, jMin, jMax, k,
316         I        rhoKm1, rhoK,
317         U        ConvectCount, KappaRT, KappaRS,
318         I        myTime, myIter, myThid)
319              END IF
320    
321           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  
322          ENDDO          ENDDO
 C--     end of k loop  
323    
324  #ifdef ALLOW_GMREDI  #ifdef  ALLOW_OBCS
325    C--     Calculate future values on open boundaries
326            IF (useOBCS) THEN
327              CALL OBCS_CALC( bi, bj, myTime+deltaT,
328         I            uVel, vVel, wVel, theta, salt,
329         I            myThid )
330            ENDIF
331    #endif  /* ALLOW_OBCS */
332    
333    C--     Determines forcing terms based on external fields
334    C       relaxation terms, etc.
335            CALL EXTERNAL_FORCING_SURF(
336         I             bi, bj, iMin, iMax, jMin, jMax,
337         I             myThid )
338    
339    #ifdef  ALLOW_GMREDI
340    C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
341            IF (useGMRedi) THEN
342              DO k=1,Nr
343                CALL GMREDI_CALC_TENSOR(
344         I             bi, bj, iMin, iMax, jMin, jMax, k,
345         I             sigmaX, sigmaY, sigmaR,
346         I             myThid )
347              ENDDO
348  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
349  CADJ STORE rhoTmp(:,:)  = comlev1_3d, key = kkey, byte = isbyte          ELSE
350  CADJ STORE rhok  (:,:)  = comlev1_3d, key = kkey, byte = isbyte            DO k=1, Nr
351  CADJ STORE rhoKm1(:,:)  = comlev1_3d, key = kkey, byte = isbyte              CALL GMREDI_CALC_TENSOR_DUMMY(
352  #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,  
353       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
354       I             myThid )       I             myThid )
355          ENDDO            ENDDO
356  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
357            ENDIF
358    #endif  /* ALLOW_GMREDI */
359    
360    #ifdef  ALLOW_KPP
361    C--     Compute KPP mixing coefficients
362            IF (useKPP) THEN
363              CALL KPP_CALC(
364         I                  bi, bj, myTime, myThid )
365            ENDIF
366    #endif  /* ALLOW_KPP */
367    
368  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
369  CADJ STORE theta(:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
370  CADJ STORE salt (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
371  CADJ STORE uvel (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
372  CADJ STORE vvel (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
373  #endif  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
374    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
375    #endif /* ALLOW_AUTODIFF_TAMC */
376    
377    #ifdef ALLOW_AIM
378    C       AIM - atmospheric intermediate model, physics package code.
379    C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics
380            IF ( useAIM ) THEN
381             CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
382             CALL AIM_DO_ATMOS_PHYSICS( phiHyd, myTime, myThid )
383             CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
384            ENDIF
385    #endif /* ALLOW_AIM */
386    
 #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  
387    
388  C--     Start of upward loop  C--     Start of thermodynamics loop
389          DO K = Nr, 1, -1          DO k=Nr,1,-1
390    
391           kM1  =max(1,k-1)   ! Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
392           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
393           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  
394    
395  #ifdef ALLOW_AUTODIFF_TAMC            km1  = MAX(1,k-1)
396           kkey = ikact*(Nr-1+1) + (k-1) + 1            kup  = 1+MOD(k+1,2)
397  #endif            kDown= 1+MOD(k,2)
398    
399  #ifdef ALLOW_AUTODIFF_TAMC            iMin = 1-OLx+2
400  CADJ STORE rvel  (:,:,kDown)  = comlev1_3d, key = kkey, byte = isbyte            iMax = sNx+OLx-1
401  CADJ STORE rTrans(:,:)        = comlev1_3d, key = kkey, byte = isbyte            jMin = 1-OLy+2
402  CADJ STORE KappaRT(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte            jMax = sNy+OLy-1
403  CADJ STORE KappaRS(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  
404  #endif  #ifdef ALLOW_AUTODIFF_TAMC
405    CPatrick Is this formula correct?
406             kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1
407    CADJ STORE rTrans(:,:)       = comlev1_bibj_k, key = kkey, byte = isbyte
408    CADJ STORE KappaRT(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
409    CADJ STORE KappaRS(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
410    #endif /* ALLOW_AUTODIFF_TAMC */
411    
412  C--      Get temporary terms used by tendency routines  C--      Get temporary terms used by tendency routines
413           CALL CALC_COMMON_FACTORS (           CALL CALC_COMMON_FACTORS (
414       I        bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
415       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,       O        xA,yA,uTrans,vTrans,rTrans,maskC,maskUp,
416       I        myThid)       I        myThid)
417    
 #ifdef ALLOW_OBCS  
         IF (openBoundaries) THEN  
          CALL APPLY_OBCS3( bi, bj, K, Kup, rTrans, rVel, myThid )  
         ENDIF  
 #endif  
   
418  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
419  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
420           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
421       I        bi,bj,iMin,iMax,jMin,jMax,K,       I        bi,bj,iMin,iMax,jMin,jMax,k,
422       I        maskC,maskUp,       I        maskC,maskup,
423       O        KappaRT,KappaRS,KappaRU,KappaRV,       O        KappaRT,KappaRS,KappaRU,KappaRV,
424       I        myThid)       I        myThid)
425  #endif  #endif
426  C--      Calculate accelerations in the momentum equations  
427           IF ( momStepping ) THEN  C--      Calculate active tracer tendencies (gT,gS,...)
428            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  
429           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
430            CALL CALC_GT(             CALL CALC_GT(
431       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
432       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
433       I         KappaRT,       I         KappaRT,
434       U         aTerm,xTerm,fZon,fMer,fVerT,       U         fVerT,
435       I         myTime, myThid)       I         myTime, myThid)
436               CALL TIMESTEP_TRACER(
437         I         bi,bj,iMin,iMax,jMin,jMax,k,
438         I         theta, gT,
439         U         gTnm1,
440         I         myIter, myThid)
441           ENDIF           ENDIF
442           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
443            CALL CALC_GS(             CALL CALC_GS(
444       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
445       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
446       I         KappaRS,       I         KappaRS,
447       U         aTerm,xTerm,fZon,fMer,fVerS,       U         fVerS,
448       I         myTime, myThid)       I         myTime, myThid)
449               CALL TIMESTEP_TRACER(
450         I         bi,bj,iMin,iMax,jMin,jMax,k,
451         I         salt, gS,
452         U         gSnm1,
453         I         myIter, myThid)
454           ENDIF           ENDIF
455  #ifdef ALLOW_OBCS  
456  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  
457  C--      Apply open boundary conditions  C--      Apply open boundary conditions
458           IF (openBoundaries) THEN           IF (useOBCS) THEN
459  #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 )  
460           END IF           END IF
461  #endif  #endif   /* ALLOW_OBCS */
462    
463  C--      Freeze water  C--      Freeze water
464           IF (allowFreezing) THEN           IF (allowFreezing) THEN
465  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
466  CADJ STORE gTNm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
467  #endif  CADJ &   , key = kkey, byte = isbyte
468              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, K, myThid )  #endif /* ALLOW_AUTODIFF_TAMC */
469                CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
470           END IF           END IF
471    
472  #ifdef DIVG_IN_DYNAMICS  C--     end of thermodynamic k loop (Nr:1)
473  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  
474    
475    
476          ENDDO ! K  #ifdef ALLOW_AUTODIFF_TAMC
477    CPatrick? What about this one?
478               maximpl = 6
479               iikey = (ikey-1)*maximpl
480    #endif /* ALLOW_AUTODIFF_TAMC */
481    
482  C--     Implicit diffusion  C--     Implicit diffusion
483          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
484    
485  #ifdef ALLOW_AUTODIFF_TAMC            IF (tempStepping) THEN
            maximpl = 6  
            iikey = ikact*maximpl  
 #endif  
   
          IF (tempStepping) THEN  
486  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
487              idkey = iikey + 1              idkey = iikey + 1
488  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
489              CALL IMPLDIFF(              CALL IMPLDIFF(
490       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
491       I         deltaTtracer, KappaRT,recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
492       U         gTNm1,       U         gTNm1,
493       I         myThid )       I         myThid )
494           END IF           ENDIF
495    
496           IF (saltStepping) THEN           IF (saltStepping) THEN
497  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
498           idkey = iikey + 2           idkey = iikey + 2
499  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
500              CALL IMPLDIFF(              CALL IMPLDIFF(
501       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
502       I         deltaTtracer, KappaRS,recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
503       U         gSNm1,       U         gSNm1,
504       I         myThid )       I         myThid )
505             ENDIF
506    
507    #ifdef   ALLOW_OBCS
508    C--      Apply open boundary conditions
509             IF (useOBCS) THEN
510               DO K=1,Nr
511                 CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
512               ENDDO
513           END IF           END IF
514    #endif   /* ALLOW_OBCS */
515    
516          ENDIF ! implicitDiffusion  C--     End If implicitDiffusion
517            ENDIF
518    
 C--     Implicit viscosity  
         IF (implicitViscosity) THEN  
519    
520           IF (momStepping) THEN  
521  #ifdef ALLOW_AUTODIFF_TAMC  C--     Start of dynamics loop
522           idkey = iikey + 3          DO k=1,Nr
523  #endif  
524    C--       km1    Points to level above k (=k-1)
525    C--       kup    Cycles through 1,2 to point to layer above
526    C--       kDown  Cycles through 2,1 to point to current layer
527    
528              km1  = MAX(1,k-1)
529              kup  = 1+MOD(k+1,2)
530              kDown= 1+MOD(k,2)
531    
532              iMin = 1-OLx+2
533              iMax = sNx+OLx-1
534              jMin = 1-OLy+2
535              jMax = sNy+OLy-1
536    
537    C--      Integrate hydrostatic balance for phiHyd with BC of
538    C        phiHyd(z=0)=0
539    C        distinguishe between Stagger and Non Stagger time stepping
540             IF (staggerTimeStep) THEN
541               CALL CALC_PHI_HYD(
542         I        bi,bj,iMin,iMax,jMin,jMax,k,
543         I        gTnm1, gSnm1,
544         U        phiHyd,
545         I        myThid )
546             ELSE
547               CALL CALC_PHI_HYD(
548         I        bi,bj,iMin,iMax,jMin,jMax,k,
549         I        theta, salt,
550         U        phiHyd,
551         I        myThid )
552             ENDIF
553    
554    C--      Calculate accelerations in the momentum equations (gU, gV, ...)
555    C        and step forward storing the result in gUnm1, gVnm1, etc...
556             IF ( momStepping ) THEN
557               CALL CALC_MOM_RHS(
558         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
559         I         phiHyd,KappaRU,KappaRV,
560         U         fVerU, fVerV,
561         I         myTime, myThid)
562               CALL TIMESTEP(
563         I         bi,bj,iMin,iMax,jMin,jMax,k,phiHyd,
564         I         myIter, myThid)
565    
566    #ifdef   ALLOW_OBCS
567    C--      Apply open boundary conditions
568             IF (useOBCS) THEN
569               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
570             END IF
571    #endif   /* ALLOW_OBCS */
572    
573    #ifdef   ALLOW_AUTODIFF_TAMC
574    #ifdef   INCLUDE_CD_CODE
575             ELSE
576               DO j=1-OLy,sNy+OLy
577                 DO i=1-OLx,sNx+OLx
578                   guCD(i,j,k,bi,bj) = 0.0
579                   gvCD(i,j,k,bi,bj) = 0.0
580                 END DO
581               END DO
582    #endif   /* INCLUDE_CD_CODE */
583    #endif   /* ALLOW_AUTODIFF_TAMC */
584             ENDIF
585    
586    
587    C--     end of dynamics k loop (1:Nr)
588            ENDDO
589    
590    
591    
592    C--     Implicit viscosity
593            IF (implicitViscosity.AND.momStepping) THEN
594    #ifdef    ALLOW_AUTODIFF_TAMC
595              idkey = iikey + 3
596    #endif    /* ALLOW_AUTODIFF_TAMC */
597            CALL IMPLDIFF(            CALL IMPLDIFF(
598       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
599       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
600       U         gUNm1,       U         gUNm1,
601       I         myThid )       I         myThid )
602  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
603           idkey = iikey + 4            idkey = iikey + 4
604  #endif  #endif    /* ALLOW_AUTODIFF_TAMC */
605            CALL IMPLDIFF(            CALL IMPLDIFF(
606       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
607       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
608       U         gVNm1,       U         gVNm1,
609       I         myThid )       I         myThid )
610    
611  #ifdef INCLUDE_CD_CODE  #ifdef   ALLOW_OBCS
612    C--      Apply open boundary conditions
613             IF (useOBCS) THEN
614               DO K=1,Nr
615                 CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
616               ENDDO
617             END IF
618    #endif   /* ALLOW_OBCS */
619    
620  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    INCLUDE_CD_CODE
621           idkey = iikey + 5  #ifdef    ALLOW_AUTODIFF_TAMC
622  #endif            idkey = iikey + 5
623    #endif    /* ALLOW_AUTODIFF_TAMC */
624            CALL IMPLDIFF(            CALL IMPLDIFF(
625       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
626       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
627       U         vVelD,       U         vVelD,
628       I         myThid )       I         myThid )
629  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
630          idkey = iikey + 6            idkey = iikey + 6
631  #endif  #endif    /* ALLOW_AUTODIFF_TAMC */
632            CALL IMPLDIFF(            CALL IMPLDIFF(
633       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
634       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
635       U         uVelD,       U         uVelD,
636       I         myThid )       I         myThid )
637    #endif    /* INCLUDE_CD_CODE */
638  #endif  C--     End If implicitViscosity.AND.momStepping
639            ENDIF
          ENDIF ! momStepping  
         ENDIF ! implicitViscosity  
640    
641         ENDDO         ENDDO
642        ENDDO        ENDDO
643    
 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 )  
   
   
644        RETURN        RETURN
645        END        END
Legend:
Removed from v.1.50  
changed lines
  Added in v.1.61
  ViewVC Help
Powered by ViewVC 1.1.22