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