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

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

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

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