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

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

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

revision 1.2 by heimbach, Mon Aug 13 18:05:26 2001 UTC revision 1.74 by jmc, Tue Jul 13 16:48:48 2004 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7    #ifdef ALLOW_AUTODIFF_TAMC
8    # ifdef ALLOW_GMREDI
9    #  include "GMREDI_OPTIONS.h"
10    # endif
11    # ifdef ALLOW_KPP
12    #  include "KPP_OPTIONS.h"
13    # endif
14    #endif /* ALLOW_AUTODIFF_TAMC */
15    
16    CBOP
17    C     !ROUTINE: THERMODYNAMICS
18    C     !INTERFACE:
19        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
20  C     /==========================================================\  C     !DESCRIPTION: \bv
21  C     | SUBROUTINE THERMODYNAMICS                                |  C     *==========================================================*
22  C     | o Controlling routine for the prognostic part of the     |  C     | SUBROUTINE THERMODYNAMICS                                
23  C     |   thermo-dynamics.                                       |  C     | o Controlling routine for the prognostic part of the      
24  C     |==========================================================|  C     |   thermo-dynamics.                                        
25  C     \==========================================================/  C     *===========================================================
26        IMPLICIT NONE  C     | The algorithm...
27    C     |
28    C     | "Correction Step"
29    C     | =================
30    C     | Here we update the horizontal velocities with the surface
31    C     | pressure such that the resulting flow is either consistent
32    C     | with the free-surface evolution or the rigid-lid:
33    C     |   U[n] = U* + dt x d/dx P
34    C     |   V[n] = V* + dt x d/dy P
35    C     |
36    C     | "Calculation of Gs"
37    C     | ===================
38    C     | This is where all the accelerations and tendencies (ie.
39    C     | physics, parameterizations etc...) are calculated
40    C     |   rho = rho ( theta[n], salt[n] )
41    C     |   b   = b(rho, theta)
42    C     |   K31 = K31 ( rho )
43    C     |   Gu[n] = Gu( u[n], v[n], wVel, b, ... )
44    C     |   Gv[n] = Gv( u[n], v[n], wVel, b, ... )
45    C     |   Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
46    C     |   Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
47    C     |
48    C     | "Time-stepping" or "Prediction"
49    C     | ================================
50    C     | The models variables are stepped forward with the appropriate
51    C     | time-stepping scheme (currently we use Adams-Bashforth II)
52    C     | - For momentum, the result is always *only* a "prediction"
53    C     | in that the flow may be divergent and will be "corrected"
54    C     | later with a surface pressure gradient.
55    C     | - Normally for tracers the result is the new field at time
56    C     | level [n+1} *BUT* in the case of implicit diffusion the result
57    C     | is also *only* a prediction.
58    C     | - We denote "predictors" with an asterisk (*).
59    C     |   U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
60    C     |   V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
61    C     |   theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
62    C     |   salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
63    C     | With implicit diffusion:
64    C     |   theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
65    C     |   salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
66    C     |   (1 + dt * K * d_zz) theta[n] = theta*
67    C     |   (1 + dt * K * d_zz) salt[n] = salt*
68    C     |
69    C     *==========================================================*
70    C     \ev
71    
72    C     !USES:
73          IMPLICIT NONE
74  C     == Global variables ===  C     == Global variables ===
75  #include "SIZE.h"  #include "SIZE.h"
76  #include "EEPARAMS.h"  #include "EEPARAMS.h"
77  #include "PARAMS.h"  #include "PARAMS.h"
78  #include "DYNVARS.h"  #include "DYNVARS.h"
79  #include "GRID.h"  #include "GRID.h"
80    #include "GAD.h"
81  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
82  #include "TR1.h"  #include "TR1.h"
83  #endif  #endif
84    #ifdef ALLOW_PTRACERS
85    #include "PTRACERS_SIZE.h"
86    #include "PTRACERS.h"
87    #endif
88    #ifdef ALLOW_TIMEAVE
89    #include "TIMEAVE_STATV.h"
90    #endif
91    
92  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
93  # include "tamc.h"  # include "tamc.h"
94  # include "tamc_keys.h"  # include "tamc_keys.h"
95  # include "FFIELDS.h"  # include "FFIELDS.h"
96    # include "EOS.h"
97  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
98  #  include "KPP.h"  #  include "KPP.h"
99  # endif  # endif
100  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
101  #  include "GMREDI.h"  #  include "GMREDI.h"
102  # endif  # endif
103    # ifdef ALLOW_EBM
104    #  include "EBM.h"
105    # endif
106  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
107    
108  #ifdef ALLOW_TIMEAVE  C     !INPUT/OUTPUT PARAMETERS:
 #include "TIMEAVE_STATV.h"  
 #endif  
   
109  C     == Routine arguments ==  C     == Routine arguments ==
110  C     myTime - Current time in simulation  C     myTime - Current time in simulation
111  C     myIter - Current iteration number in simulation  C     myIter - Current iteration number in simulation
# Line 46  C     myThid - Thread number for this in Line 114  C     myThid - Thread number for this in
114        INTEGER myIter        INTEGER myIter
115        INTEGER myThid        INTEGER myThid
116    
117    C     !LOCAL VARIABLES:
118  C     == Local variables  C     == Local variables
119  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
120  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C     uTrans, vTrans, rTrans - Per block temporaries holding flow
# Line 53  C                              transport Line 122  C                              transport
122  C                              o uTrans: Zonal transport  C                              o uTrans: Zonal transport
123  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
124  C                              o rTrans: Vertical transport  C                              o rTrans: Vertical transport
125    C     rTransKp1                o vertical volume transp. at interface k+1
126  C     maskUp                   o maskUp: land/water mask for W points  C     maskUp                   o maskUp: land/water mask for W points
127  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
128  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
129  C                                      so we need an fVer for each  C                                      so we need an fVer for each
130  C                                      variable.  C                                      variable.
 C     rhoK, rhoKM1   - Density at current level, and level above  
 C     phiHyd         - Hydrostatic part of the potential phiHydi.  
 C                      In z coords phiHydiHyd is the hydrostatic  
 C                      Potential (=pressure/rho0) anomaly  
 C                      In p coords phiHydiHyd is the geopotential  
 C                      surface height anomaly.  
 C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  
 C     phiSurfY             or geopotentiel (atmos) in X and Y direction  
131  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
132  C     KappaRS          (background + spatially varying, isopycnal term).  C     KappaRS          (background + spatially varying, isopycnal term).
133    C     useVariableK   = T when vertical diffusion is not constant
134  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
135  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
136  C     bi, bj  C     bi, bj
137  C     k, kup,        - Index for layer above and below. kup and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
138  C     kDown, km1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
139  C                      index into fVerTerm.  C                      index into fVerTerm.
 C     tauAB - Adams-Bashforth timestepping weight: 0=forward ; 1/2=Adams-Bashf.  
140        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
141        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
142        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
143        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
144        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
148        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
149    #ifdef ALLOW_PASSIVE_TRACER
150        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
151        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  #endif
152        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #ifdef ALLOW_PTRACERS
153        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerP   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
154        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #endif
       _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
155        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
156        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
157        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
158        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
159        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
160        _RL tauAB        _RL kp1Msk
161          LOGICAL useVariableK
 C This is currently used by IVDC and Diagnostics  
       _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
   
162        INTEGER iMin, iMax        INTEGER iMin, iMax
163        INTEGER jMin, jMax        INTEGER jMin, jMax
164        INTEGER bi, bj        INTEGER bi, bj
165        INTEGER i, j        INTEGER i, j
166        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
167          INTEGER iTracer, ip
168    
169  Cjmc : add for phiHyd output <- but not working if multi tile per CPU  CEOP
 c     CHARACTER*(MAX_LEN_MBUF) suff  
 c     LOGICAL  DIFFERENT_MULTIPLE  
 c     EXTERNAL DIFFERENT_MULTIPLE  
 Cjmc(end)  
   
 C---    The algorithm...  
 C  
 C       "Correction Step"  
 C       =================  
 C       Here we update the horizontal velocities with the surface  
 C       pressure such that the resulting flow is either consistent  
 C       with the free-surface evolution or the rigid-lid:  
 C         U[n] = U* + dt x d/dx P  
 C         V[n] = V* + dt x d/dy P  
 C  
 C       "Calculation of Gs"  
 C       ===================  
 C       This is where all the accelerations and tendencies (ie.  
 C       physics, parameterizations etc...) are calculated  
 C         rho = rho ( theta[n], salt[n] )  
 C         b   = b(rho, theta)  
 C         K31 = K31 ( rho )  
 C         Gu[n] = Gu( u[n], v[n], wVel, b, ... )  
 C         Gv[n] = Gv( u[n], v[n], wVel, b, ... )  
 C         Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )  
 C         Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )  
 C  
 C       "Time-stepping" or "Prediction"  
 C       ================================  
 C       The models variables are stepped forward with the appropriate  
 C       time-stepping scheme (currently we use Adams-Bashforth II)  
 C       - For momentum, the result is always *only* a "prediction"  
 C       in that the flow may be divergent and will be "corrected"  
 C       later with a surface pressure gradient.  
 C       - Normally for tracers the result is the new field at time  
 C       level [n+1} *BUT* in the case of implicit diffusion the result  
 C       is also *only* a prediction.  
 C       - We denote "predictors" with an asterisk (*).  
 C         U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )  
 C         V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )  
 C         theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )  
 C         salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )  
 C       With implicit diffusion:  
 C         theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )  
 C         salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )  
 C         (1 + dt * K * d_zz) theta[n] = theta*  
 C         (1 + dt * K * d_zz) salt[n] = salt*  
 C---  
170    
171    #ifdef ALLOW_DEBUG
172             IF ( debugLevel .GE. debLevB )
173         &    CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
174    #endif
175    
176  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
177  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
178        ikey = 1        ikey = 1
179          itdkey = 1
180  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
181    
 C--   Set up work arrays with valid (i.e. not NaN) values  
 C     These inital values do not alter the numerical results. They  
 C     just ensure that all memory references are to valid floating  
 C     point numbers. This prevents spurious hardware signals due to  
 C     uninitialised but inert locations.  
       DO j=1-OLy,sNy+OLy  
        DO i=1-OLx,sNx+OLx  
         xA(i,j)      = 0. _d 0  
         yA(i,j)      = 0. _d 0  
         uTrans(i,j)  = 0. _d 0  
         vTrans(i,j)  = 0. _d 0  
         DO k=1,Nr  
          phiHyd(i,j,k)  = 0. _d 0  
          sigmaX(i,j,k) = 0. _d 0  
          sigmaY(i,j,k) = 0. _d 0  
          sigmaR(i,j,k) = 0. _d 0  
         ENDDO  
         rhoKM1 (i,j) = 0. _d 0  
         rhok   (i,j) = 0. _d 0  
         phiSurfX(i,j) = 0. _d 0  
         phiSurfY(i,j) = 0. _d 0  
        ENDDO  
       ENDDO  
   
   
182  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
183  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
184  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
# Line 195  CHPF$ INDEPENDENT Line 189  CHPF$ INDEPENDENT
189  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
190  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
191  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
192  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
193  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,KappaRT,KappaRS
194  CHPF$&                  )  CHPF$&                  )
195  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 205  CHPF$&                  ) Line 199  CHPF$&                  )
199  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
200            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
201            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
202            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
203            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
204            act3 = myThid - 1            act3 = myThid - 1
205            max3 = nTx*nTy            max3 = nTx*nTy
   
206            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
207              itdkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
208       &                      + act3*max1*max2       &                      + act3*max1*max2
209       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
210  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
211    
212  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
213    C     These inital values do not alter the numerical results. They
214    C     just ensure that all memory references are to valid floating
215    C     point numbers. This prevents spurious hardware signals due to
216    C     uninitialised but inert locations.
217    
218          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
219           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
220              xA(i,j)        = 0. _d 0
221              yA(i,j)        = 0. _d 0
222              uTrans(i,j)    = 0. _d 0
223              vTrans(i,j)    = 0. _d 0
224            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
225              rTransKp1(i,j) = 0. _d 0
226            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
227            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
228            fVerS  (i,j,1) = 0. _d 0            fVerS  (i,j,1) = 0. _d 0
229            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
230    #ifdef ALLOW_PASSIVE_TRACER
231            fVerTr1(i,j,1) = 0. _d 0            fVerTr1(i,j,1) = 0. _d 0
232            fVerTr1(i,j,2) = 0. _d 0            fVerTr1(i,j,2) = 0. _d 0
233    #endif
234    #ifdef ALLOW_PTRACERS
235              DO ip=1,PTRACERS_num
236               fVerP  (i,j,1,ip) = 0. _d 0
237               fVerP  (i,j,2,ip) = 0. _d 0
238              ENDDO
239    #endif
240           ENDDO           ENDDO
241          ENDDO          ENDDO
242    
# Line 236  C--     Set up work arrays that need val Line 244  C--     Set up work arrays that need val
244           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
245            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
246  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
247             ConvectCount(i,j,k) = 0.             KappaRT(i,j,k)    = 0. _d 0
248             KappaRT(i,j,k) = 0. _d 0             KappaRS(i,j,k)    = 0. _d 0
249             KappaRS(i,j,k) = 0. _d 0  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
250               gT(i,j,k,bi,bj)   = 0. _d 0
251               gS(i,j,k,bi,bj)   = 0. _d 0
252    # ifdef ALLOW_PASSIVE_TRACER
253    ceh3 needs an IF ( use PASSIVE_TRACER) THEN
254               gTr1(i,j,k,bi,bj) = 0. _d 0
255    # endif
256    # ifdef ALLOW_PTRACERS
257    ceh3 this should have an   IF ( usePTRACERS ) THEN
258               DO iTracer=1,PTRACERS_numInUse
259                gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
260               ENDDO
261    # endif
262            ENDDO            ENDDO
263           ENDDO           ENDDO
264          ENDDO          ENDDO
265    
266          iMin = 1-OLx+1  c       iMin = 1-OLx
267          iMax = sNx+OLx  c       iMax = sNx+OLx
268          jMin = 1-OLy+1  c       jMin = 1-OLy
269          jMax = sNy+OLy  c       jMax = sNy+OLy
   
270    
271  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
272  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
273  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
274    CADJ STORE totphihyd
275    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
276    #ifdef ALLOW_KPP
277    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
278    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
279    #endif
280  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
281    
 C--     Start of diagnostic loop  
         DO k=Nr,1,-1  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick, is this formula correct now that we change the loop range?  
 C? Do we still need this?  
 cph kkey formula corrected.  
 cph Needed for rhok, rhokm1, in the case useGMREDI.  
          kkey = (ikey-1)*Nr + k  
 CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--       Integrate continuity vertically for vertical velocity  
           CALL INTEGRATE_FOR_W(  
      I                         bi, bj, k, uVel, vVel,  
      O                         wVel,  
      I                         myThid )  
   
 #ifdef    ALLOW_OBCS  
 #ifdef    ALLOW_NONHYDROSTATIC  
 C--       Apply OBC to W if in N-H mode  
           IF (useOBCS.AND.nonHydrostatic) THEN  
             CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  
           ENDIF  
 #endif    /* ALLOW_NONHYDROSTATIC */  
 #endif    /* ALLOW_OBCS */  
   
 C--       Calculate gradients of potential density for isoneutral  
 C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  
 c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  
           IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  
      I        theta, salt,  
      O        rhoK,  
      I        myThid )  
             IF (k.GT.1) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
              CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,  
      I        theta, salt,  
      O        rhoKm1,  
      I        myThid )  
             ENDIF  
             CALL GRAD_SIGMA(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             rhoK, rhoKm1, rhoK,  
      O             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDIF  
   
 C--       Implicit Vertical Diffusion for Convection  
 c ==> should use sigmaR !!!  
           IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN  
             CALL CALC_IVDC(  
      I        bi, bj, iMin, iMax, jMin, jMax, k,  
      I        rhoKm1, rhoK,  
      U        ConvectCount, KappaRT, KappaRS,  
      I        myTime, myIter, myThid)  
           ENDIF  
   
 C--     end of diagnostic k loop (Nr:1)  
         ENDDO  
   
282  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
283  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
284  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
285  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
286    
287  #ifdef  ALLOW_OBCS  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
288  C--     Calculate future values on open boundaries  C--     MOST of THERMODYNAMICS will be disabled
289          IF (useOBCS) THEN  #ifndef SINGLE_LAYER_MODE
           CALL OBCS_CALC( bi, bj, myTime+deltaT,  
      I            uVel, vVel, wVel, theta, salt,  
      I            myThid )  
         ENDIF  
 #endif  /* ALLOW_OBCS */  
   
 C--     Determines forcing terms based on external fields  
 C       relaxation terms, etc.  
         CALL EXTERNAL_FORCING_SURF(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph needed for KPP  
 CADJ STORE surfacetendencyU(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyV(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyS(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyT(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef  ALLOW_GMREDI  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  
 CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  
 CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
 C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  
         IF (useGMRedi) THEN  
           DO k=1,Nr  
             CALL GMREDI_CALC_TENSOR(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDDO  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           DO k=1, Nr  
             CALL GMREDI_CALC_TENSOR_DUMMY(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDDO  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         ENDIF  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  /* ALLOW_GMREDI */  
   
 #ifdef  ALLOW_KPP  
 C--     Compute KPP mixing coefficients  
         IF (useKPP) THEN  
           CALL KPP_CALC(  
      I                  bi, bj, myTime, myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           CALL KPP_CALC_DUMMY(  
      I                  bi, bj, myTime, myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         ENDIF  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KPPghat   (:,:,:,bi,bj)  
 CADJ &   , KPPviscAz (:,:,:,bi,bj)  
 CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  
 CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  
 CADJ &   , KPPfrac   (:,:  ,bi,bj)  
 CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  /* ALLOW_KPP */  
290    
291  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
292  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
293  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
294  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
295  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
296  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
297  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
298    #endif
299    #ifdef ALLOW_PTRACERS
300    cph-- moved to forward_step to avoid key computation
301    cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
302    cphCADJ &                              key=itdkey, byte=isbyte
303  #endif  #endif
304  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
305    
306  #ifdef ALLOW_AIM  #ifndef DISABLE_MULTIDIM_ADVECTION
307  C       AIM - atmospheric intermediate model, physics package code.  C--     Some advection schemes are better calculated using a multi-dimensional
308  C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics  C       method in the absence of any other terms and, if used, is done here.
309          IF ( useAIM ) THEN  C
310           CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
311           CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )  C The default is to use multi-dimensinal advection for non-linear advection
312           CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  C schemes. However, for the sake of efficiency of the adjoint it is necessary
313    C to be able to exclude this scheme to avoid excessive storage and
314    C recomputation. It *is* differentiable, if you need it.
315    C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
316    C disable this section of code.
317            IF (tempMultiDimAdvec) THEN
318    #ifdef ALLOW_DEBUG
319              IF ( debugLevel .GE. debLevB )
320         &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
321    #endif
322              CALL GAD_ADVECTION(
323         I             tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
324         I             GAD_TEMPERATURE,
325         I             uVel, vVel, wVel, theta,
326         O             gT,
327         I             bi,bj,myTime,myIter,myThid)
328            ENDIF
329            IF (saltMultiDimAdvec) THEN
330    #ifdef ALLOW_DEBUG
331              IF ( debugLevel .GE. debLevB )
332         &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
333    #endif
334              CALL GAD_ADVECTION(
335         I             saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
336         I             GAD_SALINITY,
337         I             uVel, vVel, wVel, salt,
338         O             gS,
339         I             bi,bj,myTime,myIter,myThid)
340            ENDIF
341    C Since passive tracers are configurable separately from T,S we
342    C call the multi-dimensional method for PTRACERS regardless
343    C of whether multiDimAdvection is set or not.
344    #ifdef ALLOW_PTRACERS
345            IF ( usePTRACERS ) THEN
346    #ifdef ALLOW_DEBUG
347              IF ( debugLevel .GE. debLevB )
348         &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
349    #endif
350             CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
351          ENDIF          ENDIF
352  #endif /* ALLOW_AIM */  #endif /* ALLOW_PTRACERS */
353    #endif /* DISABLE_MULTIDIM_ADVECTION */
354    
355    #ifdef ALLOW_DEBUG
356            IF ( debugLevel .GE. debLevB )
357         &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
358    #endif
359    
360  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
361          DO k=Nr,1,-1          DO k=Nr,1,-1
# Line 443  C--     Start of thermodynamics loop Line 363  C--     Start of thermodynamics loop
363  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
364  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
365  cph Also, the KappaR? need the index and subscript k!  cph Also, the KappaR? need the index and subscript k!
366           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
367  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
368    
369  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
# Line 459  C--       kDown  Cycles through 2,1 to p Line 379  C--       kDown  Cycles through 2,1 to p
379            jMin = 1-OLy            jMin = 1-OLy
380            jMax = sNy+OLy            jMax = sNy+OLy
381    
382              kp1Msk=1.
383              IF (k.EQ.Nr) kp1Msk=0.
384               DO j=1-Oly,sNy+Oly
385                DO i=1-Olx,sNx+Olx
386                 rTransKp1(i,j) = kp1Msk*rTrans(i,j)
387                ENDDO
388               ENDDO
389    #ifdef ALLOW_AUTODIFF_TAMC
390    CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
391    #endif
392    
393  C--       Get temporary terms used by tendency routines  C--       Get temporary terms used by tendency routines
394            CALL CALC_COMMON_FACTORS (            CALL CALC_COMMON_FACTORS (
395       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
396       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
397       I         myThid)       I         myThid)
398    
399  #ifdef ALLOW_AUTODIFF_TAMC            IF (k.EQ.1) THEN
400  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  C- Surface interface :
401  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte             DO j=1-Oly,sNy+Oly
402                DO i=1-Olx,sNx+Olx
403                 rTrans(i,j) = 0.
404                ENDDO
405               ENDDO
406              ELSE
407    C- Interior interface :
408               DO j=1-Oly,sNy+Oly
409                DO i=1-Olx,sNx+Olx
410                 rTrans(i,j) = rTrans(i,j)*maskC(i,j,k-1,bi,bj)
411                ENDDO
412               ENDDO
413              ENDIF
414    
415    #ifdef ALLOW_GMREDI
416    
417    C--   Residual transp = Bolus transp + Eulerian transp
418              IF (useGMRedi) THEN
419                CALL GMREDI_CALC_UVFLOW(
420         &            uTrans, vTrans, bi, bj, k, myThid)
421                IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
422         &                    rTrans, bi, bj, k, myThid)
423              ENDIF
424    
425    #ifdef ALLOW_AUTODIFF_TAMC
426    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
427    #ifdef GM_BOLUS_ADVEC
428    CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
429    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
430    #endif
431  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
432    
433    #endif /* ALLOW_GMREDI */
434    
435  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
436  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
437           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
# Line 477  C--      Calculate the total vertical di Line 439  C--      Calculate the total vertical di
439       I        maskUp,       I        maskUp,
440       O        KappaRT,KappaRS,       O        KappaRT,KappaRS,
441       I        myThid)       I        myThid)
442    # ifdef ALLOW_AUTODIFF_TAMC
443    CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
444    CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
445    # endif /* ALLOW_AUTODIFF_TAMC */
446  #endif  #endif
447    
448            iMin = 1-OLx+2            iMin = 1-OLx+2
# Line 489  C        and step forward storing result Line 455  C        and step forward storing result
455           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
456             CALL CALC_GT(             CALL CALC_GT(
457       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
458       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
459       I         KappaRT,       I         KappaRT,
460       U         fVerT,       U         fVerT,
461       I         myTime, myThid)       I         myTime,myIter,myThid)
            tauAB = 0.5d0 + abEps  
462             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
463       I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
464       I         theta, gT,       I         theta, gT,
      U         gTnm1,  
465       I         myIter, myThid)       I         myIter, myThid)
466           ENDIF           ENDIF
467    
468           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
469             CALL CALC_GS(             CALL CALC_GS(
470       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
471       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
472       I         KappaRS,       I         KappaRS,
473       U         fVerS,       U         fVerS,
474       I         myTime, myThid)       I         myTime,myIter,myThid)
            tauAB = 0.5d0 + abEps  
475             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
476       I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
477       I         salt, gS,       I         salt, gS,
      U         gSnm1,  
478       I         myIter, myThid)       I         myIter, myThid)
479           ENDIF           ENDIF
480  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
481    ceh3 needs an IF ( usePASSIVE_TRACER ) THEN
482           IF ( tr1Stepping ) THEN           IF ( tr1Stepping ) THEN
483             CALL CALC_GTR1(             CALL CALC_GTR1(
484       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
485       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
486       I         KappaRT,       I         KappaRT,
487       U         fVerTr1,       U         fVerTr1,
488       I         myTime, myThid)       I         myTime,myIter,myThid)
            tauAB = 0.5d0 + abEps  
489             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
490       I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
491       I         Tr1, gTr1,       I         Tr1, gTr1,
492       U         gTr1NM1,       I         myIter,myThid)
      I         myIter, myThid)  
493           ENDIF           ENDIF
494  #endif  #endif
495    #ifdef ALLOW_PTRACERS
496             IF ( usePTRACERS ) THEN
497               CALL PTRACERS_INTEGRATE(
498         I         bi,bj,k,
499         I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
500         X         fVerP, KappaRS,
501         I         myIter,myTime,myThid)
502             ENDIF
503    #endif /* ALLOW_PTRACERS */
504    
505  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
506  C--      Apply open boundary conditions  C--      Apply open boundary conditions
507           IF (useOBCS) THEN           IF (useOBCS) THEN
508             CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )             CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
509           END IF           END IF
510  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
511    
512  C--      Freeze water  C--      Freeze water
513           IF (allowFreezing) THEN  C  this bit of code is left here for backward compatibility.
514    C  freezing at surface level has been moved to FORWARD_STEP
515             IF ( useOldFreezing .AND. .NOT. useSEAICE
516         &       .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN
517  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
518  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
519  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
520  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
521              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
522           END IF           ENDIF
523    
524  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
525          ENDDO          ENDDO
526    
527    
528    C--     Implicit vertical advection & diffusion
529    #ifdef INCLUDE_IMPLVERTADV_CODE
530            IF ( tempImplVertAdv ) THEN
531              CALL GAD_IMPLICIT_R(
532         I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
533         I         kappaRT, wVel, theta,
534         U         gT,
535         I         bi, bj, myTime, myIter, myThid )
536            ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
537    #else /* INCLUDE_IMPLVERTADV_CODE */
538            IF     ( tempStepping .AND. implicitDiffusion ) THEN
539    #endif /* INCLUDE_IMPLVERTADV_CODE */
540  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
541  C? Patrick? What about this one?  CADJ STORE KappaRT(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
542  cph Keys iikey and idkey don't seem to be needed  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
 cph since storing occurs on different tape for each  
 cph impldiff call anyways.  
 cph Thus, common block comlev1_impl isn't needed either.  
 cph Storing below needed in the case useGMREDI.  
         iikey = (ikey-1)*maximpl  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--     Implicit diffusion  
         IF (implicitDiffusion) THEN  
   
          IF (tempStepping) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
             idkey = iikey + 1  
 CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
543  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
544              CALL IMPLDIFF(            CALL IMPLDIFF(
545       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
546       I         deltaTtracer, KappaRT, recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
547       U         gTNm1,       U         gT,
548       I         myThid )       I         myThid )
549           ENDIF          ENDIF
550    
551           IF (saltStepping) THEN  #ifdef INCLUDE_IMPLVERTADV_CODE
552            IF ( saltImplVertAdv ) THEN
553              CALL GAD_IMPLICIT_R(
554         I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
555         I         kappaRS, wVel, salt,
556         U         gS,
557         I         bi, bj, myTime, myIter, myThid )
558            ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
559    #else /* INCLUDE_IMPLVERTADV_CODE */
560            IF     ( saltStepping .AND. implicitDiffusion ) THEN
561    #endif /* INCLUDE_IMPLVERTADV_CODE */
562  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
563           idkey = iikey + 2  CADJ STORE KappaRS(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
564  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
565  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
566              CALL IMPLDIFF(            CALL IMPLDIFF(
567       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
568       I         deltaTtracer, KappaRS, recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
569       U         gSNm1,       U         gS,
570       I         myThid )       I         myThid )
571           ENDIF          ENDIF
572    
573  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
574           IF (tr1Stepping) THEN          IF ( tr1Stepping .AND. implicitDiffusion ) THEN
575  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
576  CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
577  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
578            CALL IMPLDIFF(            CALL IMPLDIFF(
579       I      bi, bj, iMin, iMax, jMin, jMax,       I      bi, bj, iMin, iMax, jMin, jMax,
580       I      deltaTtracer, KappaRT, recip_HFacC,       I      deltaTtracer, KappaRT, recip_HFacC,
581       U      gTr1Nm1,       U      gTr1,
582       I      myThid )       I      myThid )
583           ENDIF          ENDIF
584  #endif  #endif
585    
586    #ifdef ALLOW_PTRACERS
587    c #ifdef INCLUDE_IMPLVERTADV_CODE
588    c       IF ( usePTRACERS .AND. ptracerImplVertAdv ) THEN
589    c       ELSEIF ( usePTRACERS .AND. implicitDiffusion ) THEN
590    c #else
591            IF     ( usePTRACERS .AND. implicitDiffusion ) THEN
592    C--     Vertical diffusion (implicit) for passive tracers
593               CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
594            ENDIF
595    #endif /* ALLOW_PTRACERS */
596    
597  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
598  C--      Apply open boundary conditions  C--      Apply open boundary conditions
599           IF (useOBCS) THEN          IF ( ( implicitDiffusion
600         &    .OR. tempImplVertAdv
601         &    .OR. saltImplVertAdv
602         &       ) .AND. useOBCS     ) THEN
603             DO K=1,Nr             DO K=1,Nr
604               CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
605             ENDDO             ENDDO
606           END IF          ENDIF
607  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
608    
609  C--     End If implicitDiffusion  #ifdef ALLOW_TIMEAVE
610            IF ( taveFreq.GT. 0. _d 0 .AND.
611         &       buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN
612              CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
613            ENDIF
614    #ifndef HRCUBE
615            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
616              CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
617         I                           Nr, deltaTclock, bi, bj, myThid)
618          ENDIF          ENDIF
619            useVariableK = useKPP .OR. useGMredi .OR. ivdc_kappa.NE.0.
620            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
621             IF (implicitDiffusion) THEN
622              CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRT,
623         I                        Nr, 3, deltaTclock, bi, bj, myThid)
624             ELSE
625              CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
626         I                        Nr, 3, deltaTclock, bi, bj, myThid)
627             ENDIF
628            ENDIF
629    #endif /* ndef HRCUBE */
630    #endif /* ALLOW_TIMEAVE */
631    
632    #endif /* SINGLE_LAYER_MODE */
633    
634  Ccs-  C--   end bi,bj loops.
635         ENDDO         ENDDO
636        ENDDO        ENDDO
637    
638  #ifdef ALLOW_AIM  #ifdef ALLOW_DEBUG
639        IF ( useAIM ) THEN        If (debugMode) THEN
640         CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
641           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
642           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
643           CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
644           CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
645           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
646           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
647           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
648           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
649    #ifdef ALLOW_PTRACERS
650           IF ( usePTRACERS ) THEN
651             CALL PTRACERS_DEBUG(myThid)
652           ENDIF
653    #endif /* ALLOW_PTRACERS */
654        ENDIF        ENDIF
655         _EXCH_XYZ_R8(gTnm1,myThid)  #endif
656         _EXCH_XYZ_R8(gSnm1,myThid)  
657  #else  #ifdef ALLOW_DEBUG
658        IF (staggerTimeStep.AND.useCubedSphereExchange) THEN           IF ( debugLevel .GE. debLevB )
659         _EXCH_XYZ_R8(gTnm1,myThid)       &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
660         _EXCH_XYZ_R8(gSnm1,myThid)  #endif
       ENDIF  
 #endif /* ALLOW_AIM */  
661    
662        RETURN        RETURN
663        END        END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.74
  ViewVC Help
Powered by ViewVC 1.1.22