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