/[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.29 by heimbach, Tue Nov 12 20:45:41 2002 UTC revision 1.65 by dimitri, Sun Jan 25 00:31:52 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    #ifdef ALLOW_PTRACERS
7    # include "PTRACERS_OPTIONS.h"
8    #endif
9    
10  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
11  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
12  #  include "GMREDI_OPTIONS.h"  #  include "GMREDI_OPTIONS.h"
# Line 79  C     == Global variables === Line 84  C     == Global variables ===
84  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
85  #include "TR1.h"  #include "TR1.h"
86  #endif  #endif
87    #ifdef ALLOW_PTRACERS
88    #include "PTRACERS.h"
89    #endif
90    #ifdef ALLOW_TIMEAVE
91    #include "TIMEAVE_STATV.h"
92    #endif
93    
94  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
95  # include "tamc.h"  # include "tamc.h"
96  # include "tamc_keys.h"  # include "tamc_keys.h"
97  # include "FFIELDS.h"  # include "FFIELDS.h"
98    # include "EOS.h"
99  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
100  #  include "KPP.h"  #  include "KPP.h"
101  # endif  # endif
# Line 90  C     == Global variables === Line 103  C     == Global variables ===
103  #  include "GMREDI.h"  #  include "GMREDI.h"
104  # endif  # endif
105  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 #ifdef ALLOW_TIMEAVE  
 #include "TIMEAVE_STATV.h"  
 #endif  
106    
107  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
108  C     == Routine arguments ==  C     == Routine arguments ==
# Line 111  C                              transport Line 121  C                              transport
121  C                              o uTrans: Zonal transport  C                              o uTrans: Zonal transport
122  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
123  C                              o rTrans: Vertical transport  C                              o rTrans: Vertical transport
124    C     rTransKp1                o vertical volume transp. at interface k+1
125  C     maskUp                   o maskUp: land/water mask for W points  C     maskUp                   o maskUp: land/water mask for W points
126  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
127  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
128  C                                      so we need an fVer for each  C                                      so we need an fVer for each
129  C                                      variable.  C                                      variable.
130  C     rhoK, rhoKM1   - Density at current level, and level above  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.  
131  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
132  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  C     phiSurfY             or geopotentiel (atmos) in X and Y direction
133  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
134  C     KappaRS          (background + spatially varying, isopycnal term).  C     KappaRS          (background + spatially varying, isopycnal term).
135    C     useVariableK   = T when vertical diffusion is not constant
136  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
137  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
138  C     bi, bj  C     bi, bj
# Line 137  C                      index into fVerTe Line 144  C                      index into fVerTe
144        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
150        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
151    #ifdef ALLOW_PASSIVE_TRACER
152        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
153        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  #endif
154    #ifdef ALLOW_PTRACERS
155          _RL fVerP   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
156    #endif
157        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
158        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
159        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 153  C                      index into fVerTe Line 165  C                      index into fVerTe
165        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
166  C     This is currently used by IVDC and Diagnostics  C     This is currently used by IVDC and Diagnostics
167        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
168          _RL kp1Msk
169          LOGICAL useVariableK
170        INTEGER iMin, iMax        INTEGER iMin, iMax
171        INTEGER jMin, jMax        INTEGER jMin, jMax
172        INTEGER bi, bj        INTEGER bi, bj
173        INTEGER i, j        INTEGER i, j
174        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
175          INTEGER iTracer, ip
176    
177  CEOP  CEOP
178    
179    #ifdef ALLOW_DEBUG
180             IF ( debugLevel .GE. debLevB )
181         &    CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
182    #endif
183    
184  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
185  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
186        ikey = 1        ikey = 1
187          itdkey = 1
188  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
189    
 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  
         rhok   (i,j) = 0. _d 0  
         phiSurfX(i,j) = 0. _d 0  
         phiSurfY(i,j) = 0. _d 0  
        ENDDO  
       ENDDO  
   
   
190  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
191  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
192  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
# Line 194  CHPF$ INDEPENDENT Line 197  CHPF$ INDEPENDENT
197  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
198  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
199  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
200  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
201  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,KappaRT,KappaRS
202  CHPF$&                  )  CHPF$&                  )
203  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 209  CHPF$&                  ) Line 212  CHPF$&                  )
212            act3 = myThid - 1            act3 = myThid - 1
213            max3 = nTx*nTy            max3 = nTx*nTy
214            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
215            ikey = (act1 + 1) + act2*max1            itdkey = (act1 + 1) + act2*max1
216       &                      + act3*max1*max2       &                      + act3*max1*max2
217       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
218  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
219    
220  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
221    C     These inital values do not alter the numerical results. They
222    C     just ensure that all memory references are to valid floating
223    C     point numbers. This prevents spurious hardware signals due to
224    C     uninitialised but inert locations.
225    
226          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
227           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
228              xA(i,j)        = 0. _d 0
229              yA(i,j)        = 0. _d 0
230              uTrans(i,j)    = 0. _d 0
231              vTrans(i,j)    = 0. _d 0
232              rhok   (i,j)   = 0. _d 0
233              rhoKM1 (i,j)   = 0. _d 0
234              phiSurfX(i,j)  = 0. _d 0
235              phiSurfY(i,j)  = 0. _d 0
236            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
237              rTransKp1(i,j) = 0. _d 0
238            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
239            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
240            fVerS  (i,j,1) = 0. _d 0            fVerS  (i,j,1) = 0. _d 0
241            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
242    #ifdef ALLOW_PASSIVE_TRACER
243            fVerTr1(i,j,1) = 0. _d 0            fVerTr1(i,j,1) = 0. _d 0
244            fVerTr1(i,j,2) = 0. _d 0            fVerTr1(i,j,2) = 0. _d 0
245            rhoKM1 (i,j)   = 0. _d 0  #endif
246    #ifdef ALLOW_PTRACERS
247              DO ip=1,PTRACERS_num
248               fVerP  (i,j,1,ip) = 0. _d 0
249               fVerP  (i,j,2,ip) = 0. _d 0
250              ENDDO
251    #endif
252           ENDDO           ENDDO
253          ENDDO          ENDDO
254    
# Line 232  C--     Set up work arrays that need val Line 256  C--     Set up work arrays that need val
256           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
257            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
258  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
            phiHyd(i,j,k)  = 0. _d 0  
259             sigmaX(i,j,k) = 0. _d 0             sigmaX(i,j,k) = 0. _d 0
260             sigmaY(i,j,k) = 0. _d 0             sigmaY(i,j,k) = 0. _d 0
261             sigmaR(i,j,k) = 0. _d 0             sigmaR(i,j,k) = 0. _d 0
262             ConvectCount(i,j,k) = 0.             ConvectCount(i,j,k) = 0.
263             KappaRT(i,j,k) = 0. _d 0             KappaRT(i,j,k)    = 0. _d 0
264             KappaRS(i,j,k) = 0. _d 0             KappaRS(i,j,k)    = 0. _d 0
265  #ifdef ALLOW_AUTODIFF_TAMC  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
266             gT(i,j,k,bi,bj) = 0. _d 0             gT(i,j,k,bi,bj)   = 0. _d 0
267             gS(i,j,k,bi,bj) = 0. _d 0             gS(i,j,k,bi,bj)   = 0. _d 0
268  #ifdef ALLOW_PASSIVE_TRACER  # ifdef ALLOW_PASSIVE_TRACER
269    ceh3 needs an IF ( use PASSIVE_TRACER) THEN
270             gTr1(i,j,k,bi,bj) = 0. _d 0             gTr1(i,j,k,bi,bj) = 0. _d 0
271  #endif  # endif
272  #ifdef ALLOW_GMREDI  # ifdef ALLOW_PTRACERS
273             Kwx(i,j,k,bi,bj)    = 0. _d 0  ceh3 this should have an   IF ( usePTRACERS ) THEN
274             Kwy(i,j,k,bi,bj)    = 0. _d 0             DO iTracer=1,PTRACERS_numInUse
275             Kwz(i,j,k,bi,bj)    = 0. _d 0              gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
276  #ifdef GM_NON_UNITY_DIAGONAL             ENDDO
277             Kux(i,j,k,bi,bj)    = 0. _d 0  # endif
278             Kvy(i,j,k,bi,bj)    = 0. _d 0  #ifdef ALLOW_AUTODIFF_TAMC
279  #endif  cph all the following init. are necessary for TAF
280  #endif /* ALLOW_GMREDI */  cph although some of these are re-initialised later.
281  #endif  # ifdef ALLOW_GMREDI
282               Kwx(i,j,k,bi,bj)  = 0. _d 0
283               Kwy(i,j,k,bi,bj)  = 0. _d 0
284               Kwz(i,j,k,bi,bj)  = 0. _d 0
285    #  ifdef GM_NON_UNITY_DIAGONAL
286               Kux(i,j,k,bi,bj)  = 0. _d 0
287               Kvy(i,j,k,bi,bj)  = 0. _d 0
288    #  endif
289    #  ifdef GM_EXTRA_DIAGONAL
290               Kuz(i,j,k,bi,bj)  = 0. _d 0
291               Kvz(i,j,k,bi,bj)  = 0. _d 0
292    #  endif
293    #  ifdef GM_BOLUS_ADVEC
294               GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
295               GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
296    #  endif
297    #  ifdef GM_VISBECK_VARIABLE_K
298               VisbeckK(i,j,bi,bj)   = 0. _d 0
299    #  endif
300    # endif /* ALLOW_GMREDI */
301    #endif /* ALLOW_AUTODIFF_TAMC */
302            ENDDO            ENDDO
303           ENDDO           ENDDO
304          ENDDO          ENDDO
# Line 264  C This is currently also used by IVDC an Line 308  C This is currently also used by IVDC an
308          jMin = 1-OLy          jMin = 1-OLy
309          jMax = sNy+OLy          jMax = sNy+OLy
310    
   
311  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
312  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
313  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
314    CADJ STORE totphihyd
315    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
316  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
317  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
318  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
319  #endif  #endif
320  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
321    
322    #ifdef ALLOW_DEBUG
323            IF ( debugLevel .GE. debLevB )
324         &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
325    #endif
326    
327  C--     Start of diagnostic loop  C--     Start of diagnostic loop
328          DO k=Nr,1,-1          DO k=Nr,1,-1
329    
# Line 282  C? Patrick, is this formula correct now Line 332  C? Patrick, is this formula correct now
332  C? Do we still need this?  C? Do we still need this?
333  cph kkey formula corrected.  cph kkey formula corrected.
334  cph Needed for rhok, rhokm1, in the case useGMREDI.  cph Needed for rhok, rhokm1, in the case useGMREDI.
335           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
 CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
336  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
337    
338  C--       Integrate continuity vertically for vertical velocity  C--       Integrate continuity vertically for vertical velocity
# Line 310  C--       Calculate gradients of potenti Line 358  C--       Calculate gradients of potenti
358  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
359  c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
360            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
361    #ifdef ALLOW_DEBUG
362                IF ( debugLevel .GE. debLevB )
363         &       CALL DEBUG_CALL('FIND_RHO',myThid)
364    #endif
365  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
366  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
367  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
# Line 319  CADJ STORE salt (:,:,k,bi,bj) = comlev1_ Line 371  CADJ STORE salt (:,:,k,bi,bj) = comlev1_
371       I        theta, salt,       I        theta, salt,
372       O        rhoK,       O        rhoK,
373       I        myThid )       I        myThid )
374    
375              IF (k.GT.1) THEN              IF (k.GT.1) THEN
376  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
377  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
# Line 330  CADJ STORE salt (:,:,k-1,bi,bj) = comlev Line 383  CADJ STORE salt (:,:,k-1,bi,bj) = comlev
383       O        rhoKm1,       O        rhoKm1,
384       I        myThid )       I        myThid )
385              ENDIF              ENDIF
386    #ifdef ALLOW_DEBUG
387                IF ( debugLevel .GE. debLevB )
388         &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)
389    #endif
390              CALL GRAD_SIGMA(              CALL GRAD_SIGMA(
391       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
392       I             rhoK, rhoKm1, rhoK,       I             rhoK, rhoKm1, rhoK,
# Line 337  CADJ STORE salt (:,:,k-1,bi,bj) = comlev Line 394  CADJ STORE salt (:,:,k-1,bi,bj) = comlev
394       I             myThid )       I             myThid )
395            ENDIF            ENDIF
396    
397    #ifdef ALLOW_AUTODIFF_TAMC
398    CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
399    CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
400    #endif /* ALLOW_AUTODIFF_TAMC */
401  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
402  c ==> should use sigmaR !!!  c ==> should use sigmaR !!!
403            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
404    #ifdef ALLOW_DEBUG
405                IF ( debugLevel .GE. debLevB )
406         &       CALL DEBUG_CALL('CALC_IVDC',myThid)
407    #endif
408              CALL CALC_IVDC(              CALL CALC_IVDC(
409       I        bi, bj, iMin, iMax, jMin, jMax, k,       I        bi, bj, iMin, iMax, jMin, jMax, k,
410       I        rhoKm1, rhoK,       I        rhoKm1, rhoK,
# Line 354  C--     end of diagnostic k loop (Nr:1) Line 419  C--     end of diagnostic k loop (Nr:1)
419    
420  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
421  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
422  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
423  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
424    
425  #ifdef  ALLOW_OBCS  #ifdef  ALLOW_OBCS
426  C--     Calculate future values on open boundaries  C--     Calculate future values on open boundaries
427          IF (useOBCS) THEN          IF (useOBCS) THEN
428    #ifdef ALLOW_DEBUG
429              IF ( debugLevel .GE. debLevB )
430         &     CALL DEBUG_CALL('OBCS_CALC',myThid)
431    #endif
432            CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,            CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
433       I            uVel, vVel, wVel, theta, salt,       I            uVel, vVel, wVel, theta, salt,
434       I            myThid )       I            myThid )
435          ENDIF          ENDIF
436  #endif  /* ALLOW_OBCS */  #endif  /* ALLOW_OBCS */
437    
438            IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN
439  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
440  C       relaxation terms, etc.  C       relaxation terms, etc.
441          CALL EXTERNAL_FORCING_SURF(  #ifdef ALLOW_DEBUG
442            IF ( debugLevel .GE. debLevB )
443         &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
444    #endif
445            CALL EXTERNAL_FORCING_SURF(
446       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
447       I             myThid )       I             myTime, myIter, myThid )
448    
449  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
450  cph needed for KPP  cph needed for KPP
451  CADJ STORE surfacetendencyU(:,:,bi,bj)  CADJ STORE surfacetendencyU(:,:,bi,bj)
452  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
453  CADJ STORE surfacetendencyV(:,:,bi,bj)  CADJ STORE surfacetendencyV(:,:,bi,bj)
454  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
455  CADJ STORE surfacetendencyS(:,:,bi,bj)  CADJ STORE surfacetendencyS(:,:,bi,bj)
456  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
457  CADJ STORE surfacetendencyT(:,:,bi,bj)  CADJ STORE surfacetendencyT(:,:,bi,bj)
458  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
459    # ifdef ALLOW_SEAICE
460    CADJ STORE surfacetendencyTice(:,:,bi,bj)
461    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
462    # endif
463  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
464            ENDIF
465    
466  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
467  C--     MOST of THERMODYNAMICS will be disabled  C--     MOST of THERMODYNAMICS will be disabled
# Line 390  C--     MOST of THERMODYNAMICS will be d Line 470  C--     MOST of THERMODYNAMICS will be d
470  #ifdef  ALLOW_GMREDI  #ifdef  ALLOW_GMREDI
471    
472  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
473  CADJ STORE sigmaX(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte  cph storing here is needed only for one GMREDI_OPTIONS:
474  CADJ STORE sigmaY(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte  cph define GM_BOLUS_ADVEC
475  CADJ STORE sigmaR(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte  cph but I've avoided the #ifdef for now, in case more things change
476    CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
477    CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
478    CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
479  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
480    
481  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
482          IF (useGMRedi) THEN          IF (useGMRedi) THEN
483    #ifdef ALLOW_DEBUG
484              IF ( debugLevel .GE. debLevB )
485         &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
486    #endif
487            CALL GMREDI_CALC_TENSOR(            CALL GMREDI_CALC_TENSOR(
488       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
489       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
# Line 410  C--     Calculate iso-neutral slopes for Line 498  C--     Calculate iso-neutral slopes for
498          ENDIF          ENDIF
499    
500  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
501  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
502  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
503  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
504  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
505    
506  #endif  /* ALLOW_GMREDI */  #endif  /* ALLOW_GMREDI */
# Line 420  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_ Line 508  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_
508  #ifdef  ALLOW_KPP  #ifdef  ALLOW_KPP
509  C--     Compute KPP mixing coefficients  C--     Compute KPP mixing coefficients
510          IF (useKPP) THEN          IF (useKPP) THEN
511    #ifdef ALLOW_DEBUG
512              IF ( debugLevel .GE. debLevB )
513         &     CALL DEBUG_CALL('KPP_CALC',myThid)
514    #endif
515            CALL KPP_CALC(            CALL KPP_CALC(
516       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myThid )
517  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 434  CADJ STORE KPPghat   (:,:,:,bi,bj) Line 526  CADJ STORE KPPghat   (:,:,:,bi,bj)
526  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
527  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
528  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &   , KPPfrac   (:,:  ,bi,bj)
529  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte
530  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
531    
532  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
533    
534  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
535  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
536  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
537  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
538  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  
539  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
540  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
541    #endif
542    #ifdef ALLOW_PTRACERS
543    cph-- moved to forward_step to avoid key computation
544    cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
545    cphCADJ &                              key=itdkey, byte=isbyte
546  #endif  #endif
547  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
548    
549  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
550  C       AIM - atmospheric intermediate model, physics package code.  C       AIM - atmospheric intermediate model, physics package code.
 C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics  
551          IF ( useAIM ) THEN          IF ( useAIM ) THEN
552           CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  #ifdef ALLOW_DEBUG
553           CALL AIM_DO_ATMOS_PHYSICS( bi, bj, myTime, myThid )            IF ( debugLevel .GE. debLevB )
554           CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)       &     CALL DEBUG_CALL('AIM_DO_PHYSICS',myThid)
555    #endif
556             CALL TIMER_START('AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
557             CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
558             CALL TIMER_STOP( 'AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
559          ENDIF          ENDIF
560  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
561    
 #ifdef ALLOW_TIMEAVE  
         IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN  
           CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,  
      I                               deltaTclock, bi, bj, myThid)  
         ENDIF  
 #endif /* ALLOW_TIMEAVE */  
   
562  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
563  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
564  C       method in the absence of any other terms and, if used, is done here.  C       method in the absence of any other terms and, if used, is done here.
# Line 480  C recomputation. It *is* differentiable, Line 571  C recomputation. It *is* differentiable,
571  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
572  C disable this section of code.  C disable this section of code.
573          IF (tempMultiDimAdvec) THEN          IF (tempMultiDimAdvec) THEN
574            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,  #ifdef ALLOW_DEBUG
575       U                      theta,gT,            IF ( debugLevel .GE. debLevB )
576       I                      myTime,myIter,myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
577    #endif
578              CALL GAD_ADVECTION(
579         I             tempImplVertAdv,tempAdvScheme,GAD_TEMPERATURE,
580         I             uVel, vVel, wVel, theta,
581         O             gT,
582         I             bi,bj,myTime,myIter,myThid)
583          ENDIF          ENDIF
584          IF (saltMultiDimAdvec) THEN          IF (saltMultiDimAdvec) THEN
585            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,  #ifdef ALLOW_DEBUG
586       U                      salt,gS,            IF ( debugLevel .GE. debLevB )
587       I                      myTime,myIter,myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
588    #endif
589              CALL GAD_ADVECTION(
590         I             saltImplVertAdv,saltAdvScheme,GAD_SALINITY,
591         I             uVel, vVel, wVel, salt,
592         O             gS,
593         I             bi,bj,myTime,myIter,myThid)
594          ENDIF          ENDIF
595  C Since passive tracers are configurable separately from T,S we  C Since passive tracers are configurable separately from T,S we
596  C call the multi-dimensional method for PTRACERS regardless  C call the multi-dimensional method for PTRACERS regardless
597  C of whether multiDimAdvection is set or not.  C of whether multiDimAdvection is set or not.
598  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
599          IF ( usePTRACERS ) THEN          IF ( usePTRACERS ) THEN
600    #ifdef ALLOW_DEBUG
601              IF ( debugLevel .GE. debLevB )
602         &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
603    #endif
604           CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )           CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
605          ENDIF          ENDIF
606  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
607  #endif /* DISABLE_MULTIDIM_ADVECTION */  #endif /* DISABLE_MULTIDIM_ADVECTION */
608    
609    #ifdef ALLOW_DEBUG
610            IF ( debugLevel .GE. debLevB )
611         &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
612    #endif
613    
614  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
615          DO k=Nr,1,-1          DO k=Nr,1,-1
616  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
617  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
618  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
619  cph Also, the KappaR? need the index and subscript k!  cph Also, the KappaR? need the index and subscript k!
620           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
621  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
622    
623  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
# Line 521  C--       kDown  Cycles through 2,1 to p Line 633  C--       kDown  Cycles through 2,1 to p
633            jMin = 1-OLy            jMin = 1-OLy
634            jMax = sNy+OLy            jMax = sNy+OLy
635    
636              kp1Msk=1.
637              IF (k.EQ.Nr) kp1Msk=0.
638               DO j=1-Oly,sNy+Oly
639                DO i=1-Olx,sNx+Olx
640                 rTransKp1(i,j) = kp1Msk*rTrans(i,j)
641                ENDDO
642               ENDDO
643    
644  C--       Get temporary terms used by tendency routines  C--       Get temporary terms used by tendency routines
645            CALL CALC_COMMON_FACTORS (            CALL CALC_COMMON_FACTORS (
646       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
647       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
648       I         myThid)       I         myThid)
649    
650              IF (k.EQ.1) THEN
651    C- Surface interface :
652               DO j=1-Oly,sNy+Oly
653                DO i=1-Olx,sNx+Olx
654                 rTrans(i,j) = 0.
655                ENDDO
656               ENDDO
657              ELSE
658    C- Interior interface :
659               DO j=1-Oly,sNy+Oly
660                DO i=1-Olx,sNx+Olx
661                 rTrans(i,j) = rTrans(i,j)*maskC(i,j,k-1,bi,bj)
662                ENDDO
663               ENDDO
664              ENDIF
665    
666  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
667    
668  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
669            IF (useGMRedi) THEN            IF (useGMRedi) THEN
670              CALL GMREDI_CALC_UVFLOW(              CALL GMREDI_CALC_UVFLOW(
# Line 535  C--   Residual transp = Bolus transp + E Line 672  C--   Residual transp = Bolus transp + E
672              IF (K.GE.2) CALL GMREDI_CALC_WFLOW(              IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
673       &                    rTrans, bi, bj, k, myThid)       &                    rTrans, bi, bj, k, myThid)
674            ENDIF            ENDIF
 #endif /* ALLOW_GMREDI */  
675    
676  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
677  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  #ifdef GM_BOLUS_ADVEC
678  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
679    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
680    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
681    #endif
682  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
683    
684    #endif /* ALLOW_GMREDI */
685    
686  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
687  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
688           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
# Line 549  C--      Calculate the total vertical di Line 690  C--      Calculate the total vertical di
690       I        maskUp,       I        maskUp,
691       O        KappaRT,KappaRS,       O        KappaRT,KappaRS,
692       I        myThid)       I        myThid)
693    # ifdef ALLOW_AUTODIFF_TAMC
694    CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
695    CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
696    # endif /* ALLOW_AUTODIFF_TAMC */
697  #endif  #endif
698    
699            iMin = 1-OLx+2            iMin = 1-OLx+2
# Line 561  C        and step forward storing result Line 706  C        and step forward storing result
706           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
707             CALL CALC_GT(             CALL CALC_GT(
708       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
709       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
710       I         KappaRT,       I         KappaRT,
711       U         fVerT,       U         fVerT,
712       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
# Line 570  C        and step forward storing result Line 715  C        and step forward storing result
715       I         theta, gT,       I         theta, gT,
716       I         myIter, myThid)       I         myIter, myThid)
717           ENDIF           ENDIF
718    
719           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
720             CALL CALC_GS(             CALL CALC_GS(
721       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
722       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
723       I         KappaRS,       I         KappaRS,
724       U         fVerS,       U         fVerS,
725       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
# Line 583  C        and step forward storing result Line 729  C        and step forward storing result
729       I         myIter, myThid)       I         myIter, myThid)
730           ENDIF           ENDIF
731  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
732    ceh3 needs an IF ( usePASSIVE_TRACER ) THEN
733           IF ( tr1Stepping ) THEN           IF ( tr1Stepping ) THEN
734             CALL CALC_GTR1(             CALL CALC_GTR1(
735       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
736       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
737       I         KappaRT,       I         KappaRT,
738       U         fVerTr1,       U         fVerTr1,
739       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
# Line 598  C        and step forward storing result Line 745  C        and step forward storing result
745  #endif  #endif
746  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
747           IF ( usePTRACERS ) THEN           IF ( usePTRACERS ) THEN
748             CALL PTRACERS_INTEGERATE(             CALL PTRACERS_INTEGRATE(
749       I         bi,bj,k,       I         bi,bj,k,
750       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
751       X         KappaRS,       X         fVerP, KappaRS,
752       I         myIter,myTime,myThid)       I         myIter,myTime,myThid)
753           ENDIF           ENDIF
754  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
# Line 614  C--      Apply open boundary conditions Line 761  C--      Apply open boundary conditions
761  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
762    
763  C--      Freeze water  C--      Freeze water
764           IF ( allowFreezing .AND. .NOT. useSEAICE ) THEN  C  this bit of code is left here for backward compatibility.
765    C  freezing at surface level has been moved to FORWARD_STEP
766             IF ( useOldFreezing .AND. .NOT. useSEAICE
767         &       .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN
768  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
769  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
770  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
771  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
772              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
773           END IF           ENDIF
774    
775  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
776          ENDDO          ENDDO
777    
 C--     Implicit diffusion  
         IF (implicitDiffusion) THEN  
778    
779           IF (tempStepping) THEN  C--     Implicit vertical advection & diffusion
780    #ifdef INCLUDE_IMPLVERTADV_CODE
781            IF ( tempImplVertAdv ) THEN
782              CALL GAD_IMPLICIT_R(
783         I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
784         I         kappaRT, wVel, theta,
785         U         gT,
786         I         bi, bj, myTime, myIter, myThid )
787            ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
788    #else /* INCLUDE_IMPLVERTADV_CODE */
789            IF     ( tempStepping .AND. implicitDiffusion ) THEN
790    #endif /* INCLUDE_IMPLVERTADV_CODE */
791  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
792  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE KappaRT(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
793    CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
794  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
795              CALL IMPLDIFF(            CALL IMPLDIFF(
796       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
797       I         deltaTtracer, KappaRT, recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
798       U         gT,       U         gT,
799       I         myThid )       I         myThid )
800           ENDIF          ENDIF
801    
802           IF (saltStepping) THEN  #ifdef INCLUDE_IMPLVERTADV_CODE
803            IF ( saltImplVertAdv ) THEN
804              CALL GAD_IMPLICIT_R(
805         I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
806         I         kappaRS, wVel, salt,
807         U         gS,
808         I         bi, bj, myTime, myIter, myThid )
809            ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
810    #else /* INCLUDE_IMPLVERTADV_CODE */
811            IF     ( saltStepping .AND. implicitDiffusion ) THEN
812    #endif /* INCLUDE_IMPLVERTADV_CODE */
813  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
814  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE KappaRS(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
815    CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
816  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
817              CALL IMPLDIFF(            CALL IMPLDIFF(
818       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
819       I         deltaTtracer, KappaRS, recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
820       U         gS,       U         gS,
821       I         myThid )       I         myThid )
822           ENDIF          ENDIF
823    
824  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
825           IF (tr1Stepping) THEN          IF ( tr1Stepping .AND. implicitDiffusion ) THEN
826  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
827  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
828  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
829            CALL IMPLDIFF(            CALL IMPLDIFF(
830       I      bi, bj, iMin, iMax, jMin, jMax,       I      bi, bj, iMin, iMax, jMin, jMax,
831       I      deltaTtracer, KappaRT, recip_HFacC,       I      deltaTtracer, KappaRT, recip_HFacC,
832       U      gTr1,       U      gTr1,
833       I      myThid )       I      myThid )
834           ENDIF          ENDIF
835  #endif  #endif
836    
837  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
838  C Vertical diffusion (implicit) for passive tracers  c #ifdef INCLUDE_IMPLVERTADV_CODE
839           IF ( usePTRACERS ) THEN  c       IF ( usePTRACERS .AND. ptracerImplVertAdv ) THEN
840    c       ELSEIF ( usePTRACERS .AND. implicitDiffusion ) THEN
841    c #else
842            IF     ( usePTRACERS .AND. implicitDiffusion ) THEN
843    C--     Vertical diffusion (implicit) for passive tracers
844             CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )             CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
845           ENDIF          ENDIF
846  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
847    
848  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
849  C--      Apply open boundary conditions  C--      Apply open boundary conditions
850           IF (useOBCS) THEN          IF ( ( implicitDiffusion
851         &    .OR. tempImplVertAdv
852         &    .OR. saltImplVertAdv
853         &       ) .AND. useOBCS     ) THEN
854             DO K=1,Nr             DO K=1,Nr
855               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
856             ENDDO             ENDDO
857           END IF          ENDIF
858  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
859    
860  C--     End If implicitDiffusion  #ifdef ALLOW_TIMEAVE
861    #ifndef HRCUBE
862            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
863              CALL TIMEAVE_CUMUL_1T(ConvectCountTave, ConvectCount,
864         I                           Nr, deltaTclock, bi, bj, myThid)
865          ENDIF          ENDIF
866            useVariableK = useKPP .OR. useGMredi .OR. ivdc_kappa.NE.0.
867            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
868             IF (implicitDiffusion) THEN
869              CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRT,
870         I                        Nr, 3, deltaTclock, bi, bj, myThid)
871             ELSE
872              CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
873         I                        Nr, 3, deltaTclock, bi, bj, myThid)
874             ENDIF
875            ENDIF
876    #endif /* ndef HRCUBE */
877    #endif /* ALLOW_TIMEAVE */
878    
879  #endif /* SINGLE_LAYER_MODE */  #endif /* SINGLE_LAYER_MODE */
880    
881  Ccs-  C--   end bi,bj loops.
882         ENDDO         ENDDO
883        ENDDO        ENDDO
884    
885  #ifdef ALLOW_AIM  #ifdef ALLOW_DEBUG
 c     IF ( useAIM ) THEN  
 c      CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )  
 c     ENDIF  
 #endif /* ALLOW_AIM */  
 c     IF ( staggerTimeStep ) THEN  
 c      IF ( useAIM .OR. useCubedSphereExchange ) THEN  
 c        IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)  
 c        IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)  
 c      ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN  
 cc       .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN  
 c        IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)  
 c        IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)  
 c      ENDIF    
 c     ENDIF    
   
 #ifndef DISABLE_DEBUGMODE  
886        If (debugMode) THEN        If (debugMode) THEN
887         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
888         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
# Line 723  c     ENDIF Line 901  c     ENDIF
901        ENDIF        ENDIF
902  #endif  #endif
903    
904    #ifdef ALLOW_DEBUG
905             IF ( debugLevel .GE. debLevB )
906         &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
907    #endif
908    
909        RETURN        RETURN
910        END        END
Legend:
Removed from v.1.29  
changed lines
  Added in v.1.65
  ViewVC Help
Powered by ViewVC 1.1.22