/[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.14 by jmc, Fri Nov 16 03:25:41 2001 UTC revision 1.51 by edhill, Thu Oct 9 04:19:18 2003 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    #ifdef ALLOW_PTRACERS
15    # include "PTRACERS_OPTIONS.h"
16    #endif
17    #endif /* ALLOW_AUTODIFF_TAMC */
18    
19  CBOP  CBOP
20  C     !ROUTINE: THERMODYNAMICS  C     !ROUTINE: THERMODYNAMICS
21  C     !INTERFACE:  C     !INTERFACE:
# Line 71  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 82  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 109  C                                      i Line 127  C                                      i
127  C                                      so we need an fVer for each  C                                      so we need an fVer for each
128  C                                      variable.  C                                      variable.
129  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.  
130  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
131  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  C     phiSurfY             or geopotentiel (atmos) in X and Y direction
132  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
133  C     KappaRS          (background + spatially varying, isopycnal term).  C     KappaRS          (background + spatially varying, isopycnal term).
134    C     useVariableK   = T when vertical diffusion is not constant
135  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
136  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
137  C     bi, bj  C     bi, bj
# Line 133  C                      index into fVerTe Line 147  C                      index into fVerTe
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        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
       _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
150        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
151        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
152        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 145  C                      index into fVerTe Line 158  C                      index into fVerTe
158        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
159  C     This is currently used by IVDC and Diagnostics  C     This is currently used by IVDC and Diagnostics
160        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
161          LOGICAL useVariableK
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
168    
169  CEOP  CEOP
170    
171    #ifndef DISABLE_DEBUGMODE
172             IF ( debugLevel .GE. debLevB )
173         &    CALL DEBUG_ENTER('FORWARD_STEP',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 193  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 208  CHPF$&                  ) Line 204  CHPF$&                  )
204            act3 = myThid - 1            act3 = myThid - 1
205            max3 = nTx*nTy            max3 = nTx*nTy
206            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
207            ikey = (act1 + 1) + act2*max1            itdkey = (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              rhok   (i,j)   = 0. _d 0
225              rhoKM1 (i,j)   = 0. _d 0
226              phiSurfX(i,j)  = 0. _d 0
227              phiSurfY(i,j)  = 0. _d 0
228            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
229            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
230            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
# Line 230  C--     Set up work arrays that need val Line 239  C--     Set up work arrays that need val
239           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
240            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
241  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
242               sigmaX(i,j,k) = 0. _d 0
243               sigmaY(i,j,k) = 0. _d 0
244               sigmaR(i,j,k) = 0. _d 0
245             ConvectCount(i,j,k) = 0.             ConvectCount(i,j,k) = 0.
246             KappaRT(i,j,k) = 0. _d 0             KappaRT(i,j,k)    = 0. _d 0
247             KappaRS(i,j,k) = 0. _d 0             KappaRS(i,j,k)    = 0. _d 0
248  #ifdef ALLOW_AUTODIFF_TAMC  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
249             gT(i,j,k,bi,bj) = 0. _d 0             gT(i,j,k,bi,bj)   = 0. _d 0
250             gS(i,j,k,bi,bj) = 0. _d 0             gS(i,j,k,bi,bj)   = 0. _d 0
251  #ifdef ALLOW_PASSIVE_TRACER  # ifdef ALLOW_PASSIVE_TRACER
252    ceh3 needs an IF ( use PASSIVE_TRACER) THEN
253             gTr1(i,j,k,bi,bj) = 0. _d 0             gTr1(i,j,k,bi,bj) = 0. _d 0
254  #endif  # endif
255  #endif  # ifdef ALLOW_PTRACERS
256    ceh3 this should have an   IF ( usePTRACERS ) THEN
257               DO iTracer=1,PTRACERS_numInUse
258                gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
259               ENDDO
260    # endif
261    #ifdef ALLOW_AUTODIFF_TAMC
262    cph all the following init. are necessary for TAF
263    cph although some of these are re-initialised later.
264    # ifdef ALLOW_GMREDI
265               Kwx(i,j,k,bi,bj)  = 0. _d 0
266               Kwy(i,j,k,bi,bj)  = 0. _d 0
267               Kwz(i,j,k,bi,bj)  = 0. _d 0
268    #  ifdef GM_NON_UNITY_DIAGONAL
269               Kux(i,j,k,bi,bj)  = 0. _d 0
270               Kvy(i,j,k,bi,bj)  = 0. _d 0
271    #  endif
272    #  ifdef GM_EXTRA_DIAGONAL
273               Kuz(i,j,k,bi,bj)  = 0. _d 0
274               Kvz(i,j,k,bi,bj)  = 0. _d 0
275    #  endif
276    #  ifdef GM_BOLUS_ADVEC
277               GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
278               GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
279    #  endif
280    #  ifdef GM_VISBECK_VARIABLE_K
281               VisbeckK(i,j,bi,bj)   = 0. _d 0
282    #  endif
283    # endif /* ALLOW_GMREDI */
284    #endif /* ALLOW_AUTODIFF_TAMC */
285            ENDDO            ENDDO
286           ENDDO           ENDDO
287          ENDDO          ENDDO
288    
289          iMin = 1-OLx+1          iMin = 1-OLx
290          iMax = sNx+OLx          iMax = sNx+OLx
291          jMin = 1-OLy+1          jMin = 1-OLy
292          jMax = sNy+OLy          jMax = sNy+OLy
293    
   
294  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
295  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
296  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
297    CADJ STORE totphihyd
298    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
299  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
300  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
301  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
302  #endif  #endif
303  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
304    
305    #ifndef DISABLE_DEBUGMODE
306            IF ( debugLevel .GE. debLevB )
307         &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
308    #endif
309    
310  C--     Start of diagnostic loop  C--     Start of diagnostic loop
311          DO k=Nr,1,-1          DO k=Nr,1,-1
312    
# Line 267  C? Patrick, is this formula correct now Line 315  C? Patrick, is this formula correct now
315  C? Do we still need this?  C? Do we still need this?
316  cph kkey formula corrected.  cph kkey formula corrected.
317  cph Needed for rhok, rhokm1, in the case useGMREDI.  cph Needed for rhok, rhokm1, in the case useGMREDI.
318           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  
319  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
320    
321  C--       Integrate continuity vertically for vertical velocity  C--       Integrate continuity vertically for vertical velocity
322            CALL INTEGRATE_FOR_W(  c         CALL INTEGRATE_FOR_W(
323       I                         bi, bj, k, uVel, vVel,  c    I                         bi, bj, k, uVel, vVel,
324       O                         wVel,  c    O                         wVel,
325       I                         myThid )  c    I                         myThid )
326    
327  #ifdef    ALLOW_OBCS  #ifdef    ALLOW_OBCS
328  #ifdef    ALLOW_NONHYDROSTATIC  #ifdef    ALLOW_NONHYDROSTATIC
329  C--       Apply OBC to W if in N-H mode  C--       Apply OBC to W if in N-H mode
330            IF (useOBCS.AND.nonHydrostatic) THEN  c         IF (useOBCS.AND.nonHydrostatic) THEN
331              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  c           CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
332            ENDIF  c         ENDIF
333  #endif    /* ALLOW_NONHYDROSTATIC */  #endif    /* ALLOW_NONHYDROSTATIC */
334  #endif    /* ALLOW_OBCS */  #endif    /* ALLOW_OBCS */
335    
336    C--       Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
337    C--       MOST of THERMODYNAMICS will be disabled
338    #ifndef SINGLE_LAYER_MODE
339    
340  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
341  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
342  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
343            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
344    #ifndef DISABLE_DEBUGMODE
345                IF ( debugLevel .GE. debLevB )
346         &       CALL DEBUG_CALL('FIND_RHO',myThid)
347    #endif
348  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
349  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
350  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
351  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
352              CALL FIND_RHO(              CALL FIND_RHO(
353       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k, k,
354       I        theta, salt,       I        theta, salt,
355       O        rhoK,       O        rhoK,
356       I        myThid )       I        myThid )
357    
358              IF (k.GT.1) THEN              IF (k.GT.1) THEN
359  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
360  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
361  CADJ STORE salt (:,:,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
362  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
363               CALL FIND_RHO(               CALL FIND_RHO(
364       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,
365       I        theta, salt,       I        theta, salt,
366       O        rhoKm1,       O        rhoKm1,
367       I        myThid )       I        myThid )
368              ENDIF              ENDIF
369    #ifndef DISABLE_DEBUGMODE
370                IF ( debugLevel .GE. debLevB )
371         &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)
372    #endif
373              CALL GRAD_SIGMA(              CALL GRAD_SIGMA(
374       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
375       I             rhoK, rhoKm1, rhoK,       I             rhoK, rhoKm1, rhoK,
# Line 318  CADJ STORE salt (:,:,k-1,bi,bj) = comlev Line 377  CADJ STORE salt (:,:,k-1,bi,bj) = comlev
377       I             myThid )       I             myThid )
378            ENDIF            ENDIF
379    
380    #ifdef ALLOW_AUTODIFF_TAMC
381    CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
382    CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
383    #endif /* ALLOW_AUTODIFF_TAMC */
384  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
385  c ==> should use sigmaR !!!  c ==> should use sigmaR !!!
386            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
387    #ifndef DISABLE_DEBUGMODE
388                IF ( debugLevel .GE. debLevB )
389         &       CALL DEBUG_CALL('CALC_IVDC',myThid)
390    #endif
391              CALL CALC_IVDC(              CALL CALC_IVDC(
392       I        bi, bj, iMin, iMax, jMin, jMax, k,       I        bi, bj, iMin, iMax, jMin, jMax, k,
393       I        rhoKm1, rhoK,       I        rhoKm1, rhoK,
# Line 328  c ==> should use sigmaR !!! Line 395  c ==> should use sigmaR !!!
395       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
396            ENDIF            ENDIF
397    
398    #endif /* SINGLE_LAYER_MODE */
399    
400  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
401          ENDDO          ENDDO
402    
403  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
404  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
405  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
406  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
407    
408  #ifdef  ALLOW_OBCS  #ifdef  ALLOW_OBCS
409  C--     Calculate future values on open boundaries  C--     Calculate future values on open boundaries
410          IF (useOBCS) THEN          IF (useOBCS) THEN
411            CALL OBCS_CALC( bi, bj, myTime+deltaT,  #ifndef DISABLE_DEBUGMODE
412              IF ( debugLevel .GE. debLevB )
413         &     CALL DEBUG_CALL('OBCS_CALC',myThid)
414    #endif
415              CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
416       I            uVel, vVel, wVel, theta, salt,       I            uVel, vVel, wVel, theta, salt,
417       I            myThid )       I            myThid )
418          ENDIF          ENDIF
419  #endif  /* ALLOW_OBCS */  #endif  /* ALLOW_OBCS */
420    
421    
422    #ifdef ALLOW_THERM_SEAICE
423           IF (useThermSeaIce) THEN
424    #ifndef DISABLE_DEBUGMODE
425            IF ( debugLevel .GE. debLevB )
426         &    CALL DEBUG_CALL('ICE_FORCING',myThid)
427    #endif
428  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
429  C       relaxation terms, etc.  C       including effects from ice
430          CALL EXTERNAL_FORCING_SURF(          CALL ICE_FORCING(
431       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
432       I             myThid )       I             myThid )
433           ELSE
434    #endif /* ALLOW_THERM_SEAICE */
435    
436    C--     Determines forcing terms based on external fields
437    C       relaxation terms, etc.
438    #ifndef DISABLE_DEBUGMODE
439            IF ( debugLevel .GE. debLevB )
440         &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
441    #endif
442            CALL EXTERNAL_FORCING_SURF(
443         I             bi, bj, iMin, iMax, jMin, jMax,
444         I             myTime, myIter, myThid )
445    
446    #ifdef ALLOW_THERM_SEAICE
447    C--    end of if/else block useThermSeaIce --
448           ENDIF
449    #endif /* ALLOW_THERM_SEAICE */
450    
451  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
452  cph needed for KPP  cph needed for KPP
453  CADJ STORE surfacetendencyU(:,:,bi,bj)  CADJ STORE surfacetendencyU(:,:,bi,bj)
454  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
455  CADJ STORE surfacetendencyV(:,:,bi,bj)  CADJ STORE surfacetendencyV(:,:,bi,bj)
456  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
457  CADJ STORE surfacetendencyS(:,:,bi,bj)  CADJ STORE surfacetendencyS(:,:,bi,bj)
458  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
459  CADJ STORE surfacetendencyT(:,:,bi,bj)  CADJ STORE surfacetendencyT(:,:,bi,bj)
460  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
461  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
462    
463    C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
464    C--     MOST of THERMODYNAMICS will be disabled
465    #ifndef SINGLE_LAYER_MODE
466    
467  #ifdef  ALLOW_GMREDI  #ifdef  ALLOW_GMREDI
468    
469  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
470  CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  cph storing here is needed only for one GMREDI_OPTIONS:
471  CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  cph define GM_BOLUS_ADVEC
472  CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  cph but I've avoided the #ifdef for now, in case more things change
473    CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
474    CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
475    CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
476  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
477    
478  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
479          IF (useGMRedi) THEN          IF (useGMRedi) THEN
480            DO k=1,Nr  #ifndef DISABLE_DEBUGMODE
481              CALL GMREDI_CALC_TENSOR(            IF ( debugLevel .GE. debLevB )
482       I             bi, bj, iMin, iMax, jMin, jMax, k,       &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
483    #endif
484              CALL GMREDI_CALC_TENSOR(
485         I             bi, bj, iMin, iMax, jMin, jMax,
486       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
487       I             myThid )       I             myThid )
           ENDDO  
488  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
489          ELSE          ELSE
490            DO k=1, Nr            CALL GMREDI_CALC_TENSOR_DUMMY(
491              CALL GMREDI_CALC_TENSOR_DUMMY(       I             bi, bj, iMin, iMax, jMin, jMax,
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
492       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
493       I             myThid )       I             myThid )
           ENDDO  
494  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
495          ENDIF          ENDIF
496    
497  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
498  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
499  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
500  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
501  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
502    
503  #endif  /* ALLOW_GMREDI */  #endif  /* ALLOW_GMREDI */
# Line 399  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_ Line 505  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_
505  #ifdef  ALLOW_KPP  #ifdef  ALLOW_KPP
506  C--     Compute KPP mixing coefficients  C--     Compute KPP mixing coefficients
507          IF (useKPP) THEN          IF (useKPP) THEN
508    #ifndef DISABLE_DEBUGMODE
509              IF ( debugLevel .GE. debLevB )
510         &     CALL DEBUG_CALL('KPP_CALC',myThid)
511    #endif
512            CALL KPP_CALC(            CALL KPP_CALC(
513       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myThid )
514  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 413  CADJ STORE KPPghat   (:,:,:,bi,bj) Line 523  CADJ STORE KPPghat   (:,:,:,bi,bj)
523  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
524  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
525  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &   , KPPfrac   (:,:  ,bi,bj)
526  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte
527  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
528    
529  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
530    
531  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
532  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
533  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
534  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
535  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
536  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
537  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
538  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
539  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
540    #endif
541    #ifdef ALLOW_PTRACERS
542    cph-- moved to forward_step to avoid key computation
543    cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
544    cphCADJ &                              key=itdkey, byte=isbyte
545  #endif  #endif
546  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
547    
548  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
549  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  
550          IF ( useAIM ) THEN          IF ( useAIM ) THEN
551           CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  #ifndef DISABLE_DEBUGMODE
552           CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )            IF ( debugLevel .GE. debLevB )
553           CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)       &     CALL DEBUG_CALL('AIM_DO_PHYSICS',myThid)
554    #endif
555             CALL TIMER_START('AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
556             CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
557             CALL TIMER_STOP( 'AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
558          ENDIF          ENDIF
559  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
560    
 #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 */  
   
561  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
562  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
563  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 458  C to be able to exclude this scheme to a Line 569  C to be able to exclude this scheme to a
569  C recomputation. It *is* differentiable, if you need it.  C recomputation. It *is* differentiable, if you need it.
570  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
571  C disable this section of code.  C disable this section of code.
572          IF (multiDimAdvection) THEN          IF (tempMultiDimAdvec) THEN
573           IF (tempStepping .AND.  #ifndef DISABLE_DEBUGMODE
574       &       tempAdvScheme.NE.ENUM_CENTERED_2ND .AND.            IF ( debugLevel .GE. debLevB )
575       &       tempAdvScheme.NE.ENUM_UPWIND_3RD .AND.       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
576       &       tempAdvScheme.NE.ENUM_CENTERED_4TH ) THEN  #endif
577            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
578       U                      theta,gT,       U                      theta,gT,
579       I                      myTime,myIter,myThid)       I                      myTime,myIter,myThid)
580           ENDIF          ENDIF
581           IF (saltStepping .AND.          IF (saltMultiDimAdvec) THEN
582       &       saltAdvScheme.NE.ENUM_CENTERED_2ND .AND.  #ifndef DISABLE_DEBUGMODE
583       &       saltAdvScheme.NE.ENUM_UPWIND_3RD .AND.            IF ( debugLevel .GE. debLevB )
584       &       saltAdvScheme.NE.ENUM_CENTERED_4TH ) THEN       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
585    #endif
586            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
587       U                      salt,gS,       U                      salt,gS,
588       I                      myTime,myIter,myThid)       I                      myTime,myIter,myThid)
          ENDIF  
589          ENDIF          ENDIF
590    C Since passive tracers are configurable separately from T,S we
591    C call the multi-dimensional method for PTRACERS regardless
592    C of whether multiDimAdvection is set or not.
593    #ifdef ALLOW_PTRACERS
594            IF ( usePTRACERS ) THEN
595    #ifndef DISABLE_DEBUGMODE
596              IF ( debugLevel .GE. debLevB )
597         &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
598    #endif
599             CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
600            ENDIF
601    #endif /* ALLOW_PTRACERS */
602  #endif /* DISABLE_MULTIDIM_ADVECTION */  #endif /* DISABLE_MULTIDIM_ADVECTION */
603    
604    #ifndef DISABLE_DEBUGMODE
605           IF ( debugLevel .GE. debLevB )
606         &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
607    #endif
608    
609  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
610          DO k=Nr,1,-1          DO k=Nr,1,-1
611  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
612  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
613  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
614  cph Also, the KappaR? need the index and subscript k!  cph Also, the KappaR? need the index and subscript k!
615           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
616  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
617    
618  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
# Line 506  C--       Get temporary terms used by te Line 634  C--       Get temporary terms used by te
634       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
635       I         myThid)       I         myThid)
636    
637    #ifdef ALLOW_GMREDI
638    
639    C--   Residual transp = Bolus transp + Eulerian transp
640              IF (useGMRedi) THEN
641                CALL GMREDI_CALC_UVFLOW(
642         &            uTrans, vTrans, bi, bj, k, myThid)
643                IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
644         &                    rTrans, bi, bj, k, myThid)
645              ENDIF
646    
647    #ifdef ALLOW_AUTODIFF_TAMC
648    #ifdef GM_BOLUS_ADVEC
649    CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
650    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
651    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
652    #endif
653    #endif /* ALLOW_AUTODIFF_TAMC */
654    
655    #endif /* ALLOW_GMREDI */
656    
657  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
658  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
659  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
# Line 539  C        and step forward storing result Line 687  C        and step forward storing result
687       I         theta, gT,       I         theta, gT,
688       I         myIter, myThid)       I         myIter, myThid)
689           ENDIF           ENDIF
690    
691    #ifdef ALLOW_THERM_SEAICE
692             IF (useThermSeaIce .AND. k.EQ.1) THEN
693               CALL ICE_FREEZE( bi,bj, iMin,iMax,jMin,jMax, myThid )
694             ENDIF
695    #endif
696    
697           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
698             CALL CALC_GS(             CALL CALC_GS(
699       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
# Line 552  C        and step forward storing result Line 707  C        and step forward storing result
707       I         myIter, myThid)       I         myIter, myThid)
708           ENDIF           ENDIF
709  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
710    ceh3 needs an IF ( usePASSIVE_TRACER ) THEN
711           IF ( tr1Stepping ) THEN           IF ( tr1Stepping ) THEN
712             CALL CALC_GTR1(             CALL CALC_GTR1(
713       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
# Line 565  C        and step forward storing result Line 721  C        and step forward storing result
721       I         myIter,myThid)       I         myIter,myThid)
722           ENDIF           ENDIF
723  #endif  #endif
724    #ifdef ALLOW_PTRACERS
725             IF ( usePTRACERS ) THEN
726               CALL PTRACERS_INTEGRATE(
727         I         bi,bj,k,
728         I         xA,yA,uTrans,vTrans,rTrans,maskUp,
729         X         KappaRS,
730         I         myIter,myTime,myThid)
731             ENDIF
732    #endif /* ALLOW_PTRACERS */
733    
734  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
735  C--      Apply open boundary conditions  C--      Apply open boundary conditions
# Line 574  C--      Apply open boundary conditions Line 739  C--      Apply open boundary conditions
739  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
740    
741  C--      Freeze water  C--      Freeze water
742           IF (allowFreezing) THEN           IF ( allowFreezing .AND. .NOT. useSEAICE
743         &       .AND. .NOT.(useThermSeaIce.AND.k.EQ.1) ) THEN
744  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
745  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
746  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
747  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
748              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
749           END IF           ENDIF
750    
751  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
752          ENDDO          ENDDO
753    
754    cswdice -- add ---
755    #ifdef ALLOW_THERM_SEAICE
756    c timeaveraging for ice model values
757    ceh3 This should be wrapped in an IF ( useThermSeaIce ) THEN
758               CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )
759    #endif
760    cswdice --- end add ---
761    
762    
763    
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick? What about this one?  
 cph Keys iikey and idkey dont seem to be needed  
 cph since storing occurs on different tape for each  
 cph impldiff call anyways.  
 cph Thus, common block comlev1_impl isnt needed either.  
 cph Storing below needed in the case useGMREDI.  
         iikey = (ikey-1)*maximpl  
 #endif /* ALLOW_AUTODIFF_TAMC */  
764    
765  C--     Implicit diffusion  C--     Implicit diffusion
766          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
767    
768           IF (tempStepping) THEN           IF (tempStepping) THEN
769  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
770              idkey = iikey + 1  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
771  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
772              CALL IMPLDIFF(              CALL IMPLDIFF(
773       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 613  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bib Line 778  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bib
778    
779           IF (saltStepping) THEN           IF (saltStepping) THEN
780  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
781           idkey = iikey + 2  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
782  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
783              CALL IMPLDIFF(              CALL IMPLDIFF(
784       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 626  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bib Line 790  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bib
790  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
791           IF (tr1Stepping) THEN           IF (tr1Stepping) THEN
792  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
793  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gTr1(:,:,:,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,
# Line 636  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_b Line 800  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_b
800           ENDIF           ENDIF
801  #endif  #endif
802    
803    #ifdef ALLOW_PTRACERS
804    C Vertical diffusion (implicit) for passive tracers
805             IF ( usePTRACERS ) THEN
806               CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
807             ENDIF
808    #endif /* ALLOW_PTRACERS */
809    
810  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
811  C--      Apply open boundary conditions  C--      Apply open boundary conditions
812           IF (useOBCS) THEN           IF (useOBCS) THEN
# Line 648  C--      Apply open boundary conditions Line 819  C--      Apply open boundary conditions
819  C--     End If implicitDiffusion  C--     End If implicitDiffusion
820          ENDIF          ENDIF
821    
822  Ccs-  #ifdef ALLOW_TIMEAVE
823    ceh3 needs an IF ( useTIMEAVE ) THEN
824            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
825              CALL TIMEAVE_CUMUL_1T(ConvectCountTave, ConvectCount,
826         I                           Nr, deltaTclock, bi, bj, myThid)
827            ENDIF
828            useVariableK = useKPP .OR. useGMredi .OR. ivdc_kappa.NE.0.
829            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
830             IF (implicitDiffusion) THEN
831              CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRT,
832         I                        Nr, 3, deltaTclock, bi, bj, myThid)
833             ELSE
834              CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
835         I                        Nr, 3, deltaTclock, bi, bj, myThid)
836             ENDIF
837            ENDIF
838    #endif /* ALLOW_TIMEAVE */
839    
840    #endif /* SINGLE_LAYER_MODE */
841    
842    C--   end bi,bj loops.
843         ENDDO         ENDDO
844        ENDDO        ENDDO
845    
846  #ifdef ALLOW_AIM  #ifndef DISABLE_DEBUGMODE
847        IF ( useAIM ) THEN        If (debugMode) THEN
848         CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
849        ENDIF         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
850         _EXCH_XYZ_R8(gT,myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
851         _EXCH_XYZ_R8(gS,myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
852  #else         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
853        IF (staggerTimeStep.AND.useCubedSphereExchange) THEN         CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
854         _EXCH_XYZ_R8(gT,myThid)         CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
855         _EXCH_XYZ_R8(gS,myThid)         CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
856           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
857    #ifdef ALLOW_PTRACERS
858           IF ( usePTRACERS ) THEN
859             CALL PTRACERS_DEBUG(myThid)
860           ENDIF
861    #endif /* ALLOW_PTRACERS */
862        ENDIF        ENDIF
863  #endif /* ALLOW_AIM */  #endif
864    
865    #ifndef DISABLE_DEBUGMODE
866             IF ( debugLevel .GE. debLevB )
867         &    CALL DEBUG_LEAVE('FORWARD_STEP',myThid)
868    #endif
869    
870        RETURN        RETURN
871        END        END
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.51
  ViewVC Help
Powered by ViewVC 1.1.22