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

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

Parent Directory | 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.54 by edhill,
Fri Oct 24 05:29:35 2003 UTC
 Line 1
  C $Header$
  C $Name$
+ #include "PACKAGES_CONFIG.h"
  #include "CPP_OPTIONS.h"
+ #ifdef ALLOW_AUTODIFF_TAMC
+ # ifdef ALLOW_GMREDI
+ #  include "GMREDI_OPTIONS.h"
+ # endif
+ # ifdef ALLOW_KPP
+ #  include "KPP_OPTIONS.h"
+ # endif
+ #ifdef ALLOW_PTRACERS
+ # include "PTRACERS_OPTIONS.h"
+ #endif
+ #endif /* ALLOW_AUTODIFF_TAMC */
+ CBOP
+ C     !ROUTINE: THERMODYNAMICS
+ C     !INTERFACE:
        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
- C     /==========================================================\
+ C     !DESCRIPTION: \bv
- C     | SUBROUTINE THERMODYNAMICS                                |
+ C     *==========================================================*
- C     | o Controlling routine for the prognostic part of the     |
+ C     | SUBROUTINE THERMODYNAMICS
- C     |   thermo-dynamics.                                       |
+ C     | o Controlling routine for the prognostic part of the
- C     |==========================================================|
+ C     |   thermo-dynamics.
- C     \==========================================================/
+ C     *===========================================================
-       IMPLICIT NONE
+ 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     |
+ C     *==========================================================*
+ C     \ev
+ C     !USES:
+       IMPLICIT NONE
  C     == Global variables ===
  #include "SIZE.h"
  #include "EEPARAMS.h"
  #include "PARAMS.h"
  #include "DYNVARS.h"
  #include "GRID.h"
+ #include "GAD.h"
  #ifdef ALLOW_PASSIVE_TRACER
  #include "TR1.h"
  #endif
+ #ifdef ALLOW_PTRACERS
+ #include "PTRACERS.h"
+ #endif
+ #ifdef ALLOW_TIMEAVE
+ #include "TIMEAVE_STATV.h"
+ #endif
  #ifdef ALLOW_AUTODIFF_TAMC
  # include "tamc.h"
  # include "tamc_keys.h"
  # include "FFIELDS.h"
+ # include "EOS.h"
  # ifdef ALLOW_KPP
  #  include "KPP.h"
  # endif
-Line 34 
 C     == Global variables ===
+Line 104 
 C     == Global variables ===
  # endif
  #endif /* ALLOW_AUTODIFF_TAMC */
- #ifdef ALLOW_TIMEAVE
+ C     !INPUT/OUTPUT PARAMETERS:
- #include "TIMEAVE_STATV.h"
- #endif
  C     == Routine arguments ==
  C     myTime - Current time in simulation
  C     myIter - Current iteration number in simulation
-Line 46 
 C     myThid - Thread number for this in
+Line 113 
 C     myThid - Thread number for this in
        INTEGER myIter
        INTEGER myThid
+ C     !LOCAL VARIABLES:
  C     == Local variables
  C     xA, yA                 - Per block temporaries holding face areas
  C     uTrans, vTrans, rTrans - Per block temporaries holding flow
-Line 59 
 C                                      i
+Line 127 
 C                                      i
  C                                      so we need an fVer for each
  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
  C     KappaRT,       - Total diffusion in vertical for T and S.
  C     KappaRS          (background + spatially varying, isopycnal term).
+ C     useVariableK   = T when vertical diffusion is not constant
  C     iMin, iMax     - Ranges and sub-block indices on which calculations
  C     jMin, jMax       are applied.
  C     bi, bj
  C     k, kup,        - Index for layer above and below. kup and kDown
  C     kDown, km1       are switched with layer to be the appropriate
  C                      index into fVerTerm.
- C     tauAB - Adams-Bashforth timestepping weight: 0=forward ; 1/2=Adams-Bashf.
        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
-Line 84 
 C     tauAB - Adams-Bashforth timesteppi
+Line 147 
 C     tauAB - Adams-Bashforth timesteppi
        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
        _RL fVerS   (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)
        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
-Line 94 
 C     tauAB - Adams-Bashforth timesteppi
+Line 156 
 C     tauAB - Adams-Bashforth timesteppi
        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
-       _RL tauAB
+ C     This is currently used by IVDC and Diagnostics
- C This is currently used by IVDC and Diagnostics
        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
+       LOGICAL useVariableK
        INTEGER iMin, iMax
        INTEGER jMin, jMax
        INTEGER bi, bj
        INTEGER i, j
        INTEGER k, km1, kup, kDown
+       INTEGER iTracer
- 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---
+ #ifndef DISABLE_DEBUGMODE
+          IF ( debugLevel .GE. debLevB )
+      &    CALL DEBUG_ENTER('FORWARD_STEP',myThid)
+ #endif
  #ifdef ALLOW_AUTODIFF_TAMC
  C--   dummy statement to end declaration part
        ikey = 1
+       itdkey = 1
  #endif /* ALLOW_AUTODIFF_TAMC */
- 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
  #ifdef ALLOW_AUTODIFF_TAMC
  C--   HPF directive to help TAMC
  CHPF$ INDEPENDENT
-Line 195 
 CHPF$ INDEPENDENT
+Line 189 
 CHPF$ INDEPENDENT
  #ifdef ALLOW_AUTODIFF_TAMC
  C--    HPF directive to help TAMC
  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
- CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA
+ CHPF$&                  ,utrans,vtrans,xA,yA
  CHPF$&                  ,KappaRT,KappaRS
  CHPF$&                  )
  #endif /* ALLOW_AUTODIFF_TAMC */
-Line 205 
 CHPF$&                  )
+Line 199 
 CHPF$&                  )
  #ifdef ALLOW_AUTODIFF_TAMC
            act1 = bi - myBxLo(myThid)
            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
            act2 = bj - myByLo(myThid)
            max2 = myByHi(myThid) - myByLo(myThid) + 1
            act3 = myThid - 1
            max3 = nTx*nTy
            act4 = ikey_dynamics - 1
+           itdkey = (act1 + 1) + act2*max1
-           ikey = (act1 + 1) + act2*max1
       &                      + act3*max1*max2
       &                      + act4*max1*max2*max3
  #endif /* ALLOW_AUTODIFF_TAMC */
- C--     Set up work arrays that need valid initial values
+ 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
+           rhoKM1 (i,j)   = 0. _d 0
+           phiSurfX(i,j)  = 0. _d 0
+           phiSurfY(i,j)  = 0. _d 0
            rTrans (i,j)   = 0. _d 0
            fVerT  (i,j,1) = 0. _d 0
            fVerT  (i,j,2) = 0. _d 0
-Line 236 
 C--     Set up work arrays that need val
+Line 239 
 C--     Set up work arrays that need val
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
  C This is currently also used by IVDC and Diagnostics
+            sigmaX(i,j,k) = 0. _d 0
+            sigmaY(i,j,k) = 0. _d 0
+            sigmaR(i,j,k) = 0. _d 0
             ConvectCount(i,j,k) = 0.
-            KappaRT(i,j,k) = 0. _d 0
+            KappaRT(i,j,k)    = 0. _d 0
-            KappaRS(i,j,k) = 0. _d 0
+            KappaRS(i,j,k)    = 0. _d 0
+ C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
+            gT(i,j,k,bi,bj)   = 0. _d 0
+            gS(i,j,k,bi,bj)   = 0. _d 0
+ # ifdef ALLOW_PASSIVE_TRACER
+ ceh3 needs an IF ( use PASSIVE_TRACER) THEN
+            gTr1(i,j,k,bi,bj) = 0. _d 0
+ # endif
+ # ifdef ALLOW_PTRACERS
+ ceh3 this should have an   IF ( usePTRACERS ) THEN
+            DO iTracer=1,PTRACERS_numInUse
+             gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
+            ENDDO
+ # endif
+ #ifdef ALLOW_AUTODIFF_TAMC
+ cph all the following init. are necessary for TAF
+ cph although some of these are re-initialised later.
+ # ifdef ALLOW_GMREDI
+            Kwx(i,j,k,bi,bj)  = 0. _d 0
+            Kwy(i,j,k,bi,bj)  = 0. _d 0
+            Kwz(i,j,k,bi,bj)  = 0. _d 0
+ #  ifdef GM_NON_UNITY_DIAGONAL
+            Kux(i,j,k,bi,bj)  = 0. _d 0
+            Kvy(i,j,k,bi,bj)  = 0. _d 0
+ #  endif
+ #  ifdef GM_EXTRA_DIAGONAL
+            Kuz(i,j,k,bi,bj)  = 0. _d 0
+            Kvz(i,j,k,bi,bj)  = 0. _d 0
+ #  endif
+ #  ifdef GM_BOLUS_ADVEC
+            GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
+            GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
+ #  endif
+ #  ifdef GM_VISBECK_VARIABLE_K
+            VisbeckK(i,j,bi,bj)   = 0. _d 0
+ #  endif
+ # endif /* ALLOW_GMREDI */
+ #endif /* ALLOW_AUTODIFF_TAMC */
            ENDDO
           ENDDO
          ENDDO
-         iMin = 1-OLx+1
+         iMin = 1-OLx
          iMax = sNx+OLx
-         jMin = 1-OLy+1
+         jMin = 1-OLy
          jMax = sNy+OLy
  #ifdef ALLOW_AUTODIFF_TAMC
- CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
+ CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
- CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
+ CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
+ CADJ STORE totphihyd
+ CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
+ #ifdef ALLOW_KPP
+ CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
+ CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
+ #endif
  #endif /* ALLOW_AUTODIFF_TAMC */
+ #ifndef DISABLE_DEBUGMODE
+         IF ( debugLevel .GE. debLevB )
+      &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
+ #endif
  C--     Start of diagnostic loop
          DO k=Nr,1,-1
-Line 262 
 C? Patrick, is this formula correct now
+Line 315 
 C? Patrick, is this formula correct now
  C? Do we still need this?
  cph kkey formula corrected.
  cph Needed for rhok, rhokm1, in the case useGMREDI.
-          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
  #endif /* ALLOW_AUTODIFF_TAMC */
  C--       Integrate continuity vertically for vertical velocity
-           CALL INTEGRATE_FOR_W(
+ c         CALL INTEGRATE_FOR_W(
-      I                         bi, bj, k, uVel, vVel,
+ c    I                         bi, bj, k, uVel, vVel,
-      O                         wVel,
+ c    O                         wVel,
-      I                         myThid )
+ c    I                         myThid )
  #ifdef    ALLOW_OBCS
  #ifdef    ALLOW_NONHYDROSTATIC
  C--       Apply OBC to W if in N-H mode
-           IF (useOBCS.AND.nonHydrostatic) THEN
+ c         IF (useOBCS.AND.nonHydrostatic) THEN
-             CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
+ c           CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
-           ENDIF
+ c         ENDIF
  #endif    /* ALLOW_NONHYDROSTATIC */
  #endif    /* ALLOW_OBCS */
+ C--       Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
+ C--       MOST of THERMODYNAMICS will be disabled
+ #ifndef SINGLE_LAYER_MODE
  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
+ #ifndef DISABLE_DEBUGMODE
+             IF ( debugLevel .GE. debLevB )
+      &       CALL DEBUG_CALL('FIND_RHO',myThid)
+ #endif
  #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        bi, bj, iMin, iMax, jMin, jMax, k, k,
       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        bi, bj, iMin, iMax, jMin, jMax, k-1, k,
       I        theta, salt,
       O        rhoKm1,
       I        myThid )
              ENDIF
+ #ifndef DISABLE_DEBUGMODE
+             IF ( debugLevel .GE. debLevB )
+      &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)
+ #endif
              CALL GRAD_SIGMA(
       I             bi, bj, iMin, iMax, jMin, jMax, k,
       I             rhoK, rhoKm1, rhoK,
-Line 313 
 CADJ STORE salt (:,:,k-1,bi,bj) = comlev
+Line 377 
 CADJ STORE salt (:,:,k-1,bi,bj) = comlev
       I             myThid )
            ENDIF
+ #ifdef ALLOW_AUTODIFF_TAMC
+ CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
+ CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
+ #endif /* ALLOW_AUTODIFF_TAMC */
  C--       Implicit Vertical Diffusion for Convection
  c ==> should use sigmaR !!!
            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
+ #ifndef DISABLE_DEBUGMODE
+             IF ( debugLevel .GE. debLevB )
+      &       CALL DEBUG_CALL('CALC_IVDC',myThid)
+ #endif
              CALL CALC_IVDC(
       I        bi, bj, iMin, iMax, jMin, jMax, k,
       I        rhoKm1, rhoK,
-Line 323 
 c ==> should use sigmaR !!!
+Line 395 
 c ==> should use sigmaR !!!
       I        myTime, myIter, myThid)
            ENDIF
+ #endif /* SINGLE_LAYER_MODE */
  C--     end of diagnostic k loop (Nr:1)
          ENDDO
  #ifdef ALLOW_AUTODIFF_TAMC
  cph avoids recomputation of integrate_for_w
- CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
+ CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
  #ifdef  ALLOW_OBCS
  C--     Calculate future values on open boundaries
          IF (useOBCS) THEN
-           CALL OBCS_CALC( bi, bj, myTime+deltaT,
+ #ifndef DISABLE_DEBUGMODE
+           IF ( debugLevel .GE. debLevB )
+      &     CALL DEBUG_CALL('OBCS_CALC',myThid)
+ #endif
+           CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
       I            uVel, vVel, wVel, theta, salt,
       I            myThid )
          ENDIF
  #endif  /* ALLOW_OBCS */
+ #ifdef ALLOW_THERM_SEAICE
+        IF (useThermSeaIce) THEN
+ #ifndef DISABLE_DEBUGMODE
+         IF ( debugLevel .GE. debLevB )
+      &    CALL DEBUG_CALL('ICE_FORCING',myThid)
+ #endif
  C--     Determines forcing terms based on external fields
- C       relaxation terms, etc.
+ C       including effects from ice
-         CALL EXTERNAL_FORCING_SURF(
+         CALL ICE_FORCING(
       I             bi, bj, iMin, iMax, jMin, jMax,
       I             myThid )
+        ELSE
+ #endif /* ALLOW_THERM_SEAICE */
+ C--     Determines forcing terms based on external fields
+ C       relaxation terms, etc.
+ #ifndef DISABLE_DEBUGMODE
+         IF ( debugLevel .GE. debLevB )
+      &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
+ #endif
+         CALL EXTERNAL_FORCING_SURF(
+      I             bi, bj, iMin, iMax, jMin, jMax,
+      I             myTime, myIter, myThid )
+ #ifdef ALLOW_THERM_SEAICE
+ C--    end of if/else block useThermSeaIce --
+        ENDIF
+ #endif /* ALLOW_THERM_SEAICE */
  #ifdef ALLOW_AUTODIFF_TAMC
  cph needed for KPP
  CADJ STORE surfacetendencyU(:,:,bi,bj)
- CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
+ CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
  CADJ STORE surfacetendencyV(:,:,bi,bj)
- CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
+ CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
  CADJ STORE surfacetendencyS(:,:,bi,bj)
- CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
+ CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
  CADJ STORE surfacetendencyT(:,:,bi,bj)
- CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
+ CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
+ C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
+ C--     MOST of THERMODYNAMICS will be disabled
+ #ifndef SINGLE_LAYER_MODE
  #ifdef  ALLOW_GMREDI
  #ifdef ALLOW_AUTODIFF_TAMC
- CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte
+ cph storing here is needed only for one GMREDI_OPTIONS:
- CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte
+ cph define GM_BOLUS_ADVEC
- CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte
+ cph but I've avoided the #ifdef for now, in case more things change
+ CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
+ CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
+ CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
          IF (useGMRedi) THEN
-           DO k=1,Nr
+ #ifndef DISABLE_DEBUGMODE
-             CALL GMREDI_CALC_TENSOR(
+           IF ( debugLevel .GE. debLevB )
-      I             bi, bj, iMin, iMax, jMin, jMax, k,
+      &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
+ #endif
+           CALL GMREDI_CALC_TENSOR(
+      I             bi, bj, iMin, iMax, jMin, jMax,
       I             sigmaX, sigmaY, sigmaR,
       I             myThid )
-           ENDDO
  #ifdef ALLOW_AUTODIFF_TAMC
          ELSE
-           DO k=1, Nr
+           CALL GMREDI_CALC_TENSOR_DUMMY(
-             CALL GMREDI_CALC_TENSOR_DUMMY(
+      I             bi, bj, iMin, iMax, jMin, jMax,
-      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 Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
- CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte
+ CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
- CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte
+ CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
  #endif  /* ALLOW_GMREDI */
-Line 394 
 CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_
+Line 505 
 CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_
  #ifdef  ALLOW_KPP
  C--     Compute KPP mixing coefficients
          IF (useKPP) THEN
+ #ifndef DISABLE_DEBUGMODE
+           IF ( debugLevel .GE. debLevB )
+      &     CALL DEBUG_CALL('KPP_CALC',myThid)
+ #endif
            CALL KPP_CALC(
       I                  bi, bj, myTime, myThid )
  #ifdef ALLOW_AUTODIFF_TAMC
-Line 405 
 C--     Compute KPP mixing coefficients
+Line 520 
 C--     Compute KPP mixing coefficients
  #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
+ CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
  #endif  /* ALLOW_KPP */
  #ifdef ALLOW_AUTODIFF_TAMC
- CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
+ CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
- CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
+ CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
- CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
+ CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
- 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
  #ifdef ALLOW_PASSIVE_TRACER
- CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
+ CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
+ #endif
+ #ifdef ALLOW_PTRACERS
+ cph-- moved to forward_step to avoid key computation
+ cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
+ cphCADJ &                              key=itdkey, byte=isbyte
  #endif
  #endif /* ALLOW_AUTODIFF_TAMC */
  #ifdef ALLOW_AIM
  C       AIM - atmospheric intermediate model, physics package code.
- C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics
          IF ( useAIM ) THEN
-          CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
+ #ifndef DISABLE_DEBUGMODE
-          CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )
+           IF ( debugLevel .GE. debLevB )
-          CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
+      &     CALL DEBUG_CALL('AIM_DO_PHYSICS',myThid)
+ #endif
+          CALL TIMER_START('AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
+          CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
+          CALL TIMER_STOP( 'AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
          ENDIF
  #endif /* ALLOW_AIM */
+ #ifndef DISABLE_MULTIDIM_ADVECTION
+ C--     Some advection schemes are better calculated using a multi-dimensional
+ C       method in the absence of any other terms and, if used, is done here.
+ C
+ C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
+ C The default is to use multi-dimensinal advection for non-linear advection
+ C schemes. However, for the sake of efficiency of the adjoint it is necessary
+ C to be able to exclude this scheme to avoid excessive storage and
+ C recomputation. It *is* differentiable, if you need it.
+ C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
+ C disable this section of code.
+         IF (tempMultiDimAdvec) THEN
+ #ifndef DISABLE_DEBUGMODE
+           IF ( debugLevel .GE. debLevB )
+      &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
+ #endif
+           CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
+      U                      theta,gT,
+      I                      myTime,myIter,myThid)
+         ENDIF
+         IF (saltMultiDimAdvec) THEN
+ #ifndef DISABLE_DEBUGMODE
+           IF ( debugLevel .GE. debLevB )
+      &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
+ #endif
+           CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
+      U                      salt,gS,
+      I                      myTime,myIter,myThid)
+         ENDIF
+ C Since passive tracers are configurable separately from T,S we
+ C call the multi-dimensional method for PTRACERS regardless
+ C of whether multiDimAdvection is set or not.
+ #ifdef ALLOW_PTRACERS
+         IF ( usePTRACERS ) THEN
+ #ifndef DISABLE_DEBUGMODE
+           IF ( debugLevel .GE. debLevB )
+      &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
+ #endif
+          CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
+         ENDIF
+ #endif /* ALLOW_PTRACERS */
+ #endif /* DISABLE_MULTIDIM_ADVECTION */
+ #ifndef DISABLE_DEBUGMODE
+        IF ( debugLevel .GE. debLevB )
+      &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
+ #endif
  C--     Start of thermodynamics loop
          DO k=Nr,1,-1
-Line 443 
 C--     Start of thermodynamics loop
+Line 610 
 C--     Start of thermodynamics loop
  C? Patrick Is this formula correct?
  cph Yes, but I rewrote it.
  cph Also, the KappaR? need the index and subscript k!
-          kkey = (ikey-1)*Nr + k
+          kkey = (itdkey-1)*Nr + k
  #endif /* ALLOW_AUTODIFF_TAMC */
  C--       km1    Points to level above k (=k-1)
-Line 465 
 C--       Get temporary terms used by te
+Line 632 
 C--       Get temporary terms used by te
       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
       I         myThid)
+ #ifdef ALLOW_GMREDI
+ C--   Residual transp = Bolus transp + Eulerian transp
+           IF (useGMRedi) THEN
+             CALL GMREDI_CALC_UVFLOW(
+      &            uTrans, vTrans, bi, bj, k, myThid)
+             IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
+      &                    rTrans, bi, bj, k, myThid)
+           ENDIF
  #ifdef ALLOW_AUTODIFF_TAMC
- CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
+ #ifdef GM_BOLUS_ADVEC
- CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
+ CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
+ CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
+ CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
+ #endif
  #endif /* ALLOW_AUTODIFF_TAMC */
+ #endif /* ALLOW_GMREDI */
  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
  C--      Calculate the total vertical diffusivity
           CALL CALC_DIFFUSIVITY(
-Line 477 
 C--      Calculate the total vertical di
+Line 659 
 C--      Calculate the total vertical di
       I        maskUp,
       O        KappaRT,KappaRS,
       I        myThid)
+ # ifdef ALLOW_AUTODIFF_TAMC
+ CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
+ CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
+ # endif /* ALLOW_AUTODIFF_TAMC */
  #endif
            iMin = 1-OLx+2
-Line 492 
 C        and step forward storing result
+Line 678 
 C        and step forward storing result
       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
       I         KappaRT,
       U         fVerT,
-      I         myTime, myThid)
+      I         myTime,myIter,myThid)
-            tauAB = 0.5d0 + abEps
             CALL TIMESTEP_TRACER(
-      I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
+      I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
       I         theta, gT,
-      U         gTnm1,
       I         myIter, myThid)
           ENDIF
+ #ifdef ALLOW_THERM_SEAICE
+          IF (useThermSeaIce .AND. k.EQ.1) THEN
+            CALL ICE_FREEZE( bi,bj, iMin,iMax,jMin,jMax, myThid )
+          ENDIF
+ #endif
           IF ( saltStepping ) THEN
             CALL CALC_GS(
       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
       I         KappaRS,
       U         fVerS,
-      I         myTime, myThid)
+      I         myTime,myIter,myThid)
-            tauAB = 0.5d0 + abEps
             CALL TIMESTEP_TRACER(
-      I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
+      I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
       I         salt, gS,
-      U         gSnm1,
       I         myIter, myThid)
           ENDIF
  #ifdef ALLOW_PASSIVE_TRACER
+ ceh3 needs an IF ( usePASSIVE_TRACER ) THEN
           IF ( tr1Stepping ) THEN
             CALL CALC_GTR1(
       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
       I         KappaRT,
       U         fVerTr1,
-      I         myTime, myThid)
+      I         myTime,myIter,myThid)
-            tauAB = 0.5d0 + abEps
             CALL TIMESTEP_TRACER(
-      I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
+      I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
       I         Tr1, gTr1,
-      U         gTr1NM1,
+      I         myIter,myThid)
-      I         myIter, myThid)
           ENDIF
  #endif
+ #ifdef ALLOW_PTRACERS
+          IF ( usePTRACERS ) THEN
+            CALL PTRACERS_INTEGRATE(
+      I         bi,bj,k,
+      I         xA,yA,uTrans,vTrans,rTrans,maskUp,
+      X         KappaRS,
+      I         myIter,myTime,myThid)
+          ENDIF
+ #endif /* ALLOW_PTRACERS */
  #ifdef   ALLOW_OBCS
  C--      Apply open boundary conditions
           IF (useOBCS) THEN
-            CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
+            CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
           END IF
  #endif   /* ALLOW_OBCS */
  C--      Freeze water
-          IF (allowFreezing) THEN
+          IF ( allowFreezing .AND. .NOT. useSEAICE
+      &       .AND. .NOT.(useThermSeaIce.AND.k.EQ.1) ) THEN
  #ifdef ALLOW_AUTODIFF_TAMC
- CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
+ CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
  CADJ &   , key = kkey, byte = isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
-          END IF
+          ENDIF
  C--     end of thermodynamic k loop (Nr:1)
          ENDDO
+ cswdice -- add ---
+ #ifdef ALLOW_THERM_SEAICE
+ c timeaveraging for ice model values
+ ceh3 This should be wrapped in an IF ( useThermSeaIce ) THEN
+            CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )
+ #endif
+ cswdice --- end add ---
- #ifdef ALLOW_AUTODIFF_TAMC
- C? Patrick? What about this one?
- cph Keys iikey and idkey don't seem to be needed
- 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 KappaRT(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
- CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
+ CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
              CALL IMPLDIFF(
       I         bi, bj, iMin, iMax, jMin, jMax,
       I         deltaTtracer, KappaRT, recip_HFacC,
-      U         gTNm1,
+      U         gT,
       I         myThid )
           ENDIF
           IF (saltStepping) THEN
  #ifdef ALLOW_AUTODIFF_TAMC
-          idkey = iikey + 2
+ CADJ STORE KappaRS(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
- CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
+ CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
              CALL IMPLDIFF(
       I         bi, bj, iMin, iMax, jMin, jMax,
       I         deltaTtracer, KappaRS, recip_HFacC,
-      U         gSNm1,
+      U         gS,
       I         myThid )
           ENDIF
  #ifdef ALLOW_PASSIVE_TRACER
           IF (tr1Stepping) THEN
  #ifdef ALLOW_AUTODIFF_TAMC
- CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
+ CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
            CALL IMPLDIFF(
       I      bi, bj, iMin, iMax, jMin, jMax,
       I      deltaTtracer, KappaRT, recip_HFacC,
-      U      gTr1Nm1,
+      U      gTr1,
       I      myThid )
           ENDIF
  #endif
+ #ifdef ALLOW_PTRACERS
+ C Vertical diffusion (implicit) for passive tracers
+          IF ( usePTRACERS ) THEN
+            CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
+          ENDIF
+ #endif /* ALLOW_PTRACERS */
  #ifdef   ALLOW_OBCS
  C--      Apply open boundary conditions
           IF (useOBCS) THEN
             DO K=1,Nr
-              CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
+              CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
             ENDDO
           END IF
  #endif   /* ALLOW_OBCS */
-Line 613 
 C--      Apply open boundary conditions
+Line 818 
 C--      Apply open boundary conditions
  C--     End If implicitDiffusion
          ENDIF
- Ccs-
+ #ifdef ALLOW_TIMEAVE
+ ceh3 needs an IF ( useTIMEAVE ) THEN
+         IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
+           CALL TIMEAVE_CUMUL_1T(ConvectCountTave, ConvectCount,
+      I                           Nr, deltaTclock, bi, bj, myThid)
+         ENDIF
+         useVariableK = useKPP .OR. useGMredi .OR. ivdc_kappa.NE.0.
+         IF (taveFreq.GT.0. .AND. useVariableK ) THEN
+          IF (implicitDiffusion) THEN
+           CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRT,
+      I                        Nr, 3, deltaTclock, bi, bj, myThid)
+          ELSE
+           CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
+      I                        Nr, 3, deltaTclock, bi, bj, myThid)
+          ENDIF
+         ENDIF
+ #endif /* ALLOW_TIMEAVE */
+ #endif /* SINGLE_LAYER_MODE */
+ C--   end bi,bj loops.
         ENDDO
        ENDDO
- #ifdef ALLOW_AIM
+ #ifndef DISABLE_DEBUGMODE
-       IF ( useAIM ) THEN
+       If (debugMode) THEN
-        CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )
+        CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
-       ENDIF
+        CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
-        _EXCH_XYZ_R8(gTnm1,myThid)
+        CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
-        _EXCH_XYZ_R8(gSnm1,myThid)
+        CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
- #else
+        CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
-       IF (staggerTimeStep.AND.useCubedSphereExchange) THEN
+        CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
-        _EXCH_XYZ_R8(gTnm1,myThid)
+        CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
-        _EXCH_XYZ_R8(gSnm1,myThid)
+        CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
+        CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
+ #ifdef ALLOW_PTRACERS
+        IF ( usePTRACERS ) THEN
+          CALL PTRACERS_DEBUG(myThid)
+        ENDIF
+ #endif /* ALLOW_PTRACERS */
        ENDIF
- #endif /* ALLOW_AIM */
+ #endif
+ #ifndef DISABLE_DEBUGMODE
+          IF ( debugLevel .GE. debLevB )
+      &    CALL DEBUG_LEAVE('FORWARD_STEP',myThid)
+ #endif
        RETURN
        END

 Legend:



Removed from v.1.2
 


changed lines


 
Added in v.1.54
 Legend:



Removed from v.1.2
 


changed lines


 
Added in v.1.54
-Removed from v.1.2
+Added in v.1.54

	ViewVC Help
Powered by ViewVC 1.1.22