model/src/temp_integrate.F

4794324d0b Jean*0001 #include "PACKAGES_CONFIG.h"
                0002 #include "CPP_OPTIONS.h"
517dbdc414 Jean*0003 #ifdef ALLOW_AUTODIFF
                0004 # include "AUTODIFF_OPTIONS.h"
                0005 #endif
fb8688a588 Jean*0006 #ifdef ALLOW_GENERIC_ADVDIFF
                0007 # include "GAD_OPTIONS.h"
                0008 #endif
4794324d0b Jean*0009
                0010 CBOP
                0011 C     !ROUTINE: TEMP_INTEGRATE
                0012 C     !INTERFACE:
                0013       SUBROUTINE TEMP_INTEGRATE(
0b191c5f5a Jean*0014      I           bi, bj, recip_hFac,
fb8688a588 Jean*0015      I           uFld, vFld, wFld,
                0016      U           KappaRk,
4794324d0b Jean*0017      I           myTime, myIter, myThid )
                0018 C     !DESCRIPTION: \bv
                0019 C     *==========================================================*
                0020 C     | SUBROUTINE TEMP_INTEGRATE
bd27360393 Jean*0021 C     | o Calculate tendency for temperature and integrates
                0022 C     |   forward in time. The temperature array is updated here
                0023 C     |   while adjustments (filters, conv.adjustment) are applied
                0024 C     |   later, in S/R TRACERS_CORRECTION_STEP.
4794324d0b Jean*0025 C     *==========================================================*
e874fa47e5 Jean*0026 C     | A procedure called APPLY_FORCING_T is called from
4794324d0b Jean*0027 C     | here. These procedures can be used to add per problem
                0028 C     | heat flux source terms.
                0029 C     | Note: Although it is slightly counter-intuitive the
                0030 C     |       EXTERNAL_FORCING routine is not the place to put
                0031 C     |       file I/O. Instead files that are required to
                0032 C     |       calculate the external source terms are generally
                0033 C     |       read during the model main loop. This makes the
                0034 C     |       logistics of multi-processing simpler and also
                0035 C     |       makes the adjoint generation simpler. It also
                0036 C     |       allows for I/O to overlap computation where that
                0037 C     |       is supported by hardware.
                0038 C     | Aside from the problem specific term the code here
                0039 C     | forms the tendency terms due to advection and mixing
                0040 C     | The baseline implementation here uses a centered
                0041 C     | difference form for the advection term and a tensorial
                0042 C     | divergence of a flux form for the diffusive term. The
fb8688a588 Jean*0043 C     | diffusive term is formulated so that isopycnal mixing
                0044 C     | and GM-style subgrid-scale terms can be incorporated by
                0045 C     | simply setting the diffusion tensor terms appropriately.
4794324d0b Jean*0046 C     *==========================================================*
                0047 C     \ev
                0048
                0049 C     !USES:
                0050       IMPLICIT NONE
                0051 C     == GLobal variables ==
                0052 #include "SIZE.h"
                0053 #include "EEPARAMS.h"
                0054 #include "PARAMS.h"
0b191c5f5a Jean*0055 #include "GRID.h"
                0056 #include "DYNVARS.h"
4794324d0b Jean*0057 #include "RESTART.h"
                0058 #ifdef ALLOW_GENERIC_ADVDIFF
fb8688a588 Jean*0059 # include "GAD.h"
                0060 # include "GAD_SOM_VARS.h"
4794324d0b Jean*0061 #endif
0b191c5f5a Jean*0062 #ifdef ALLOW_TIMEAVE
                0063 # include "TIMEAVE_STATV.h"
                0064 #endif
7c50f07931 Mart*0065 #ifdef ALLOW_AUTODIFF_TAMC
4794324d0b Jean*0066 # include "tamc.h"
                0067 #endif
                0068
                0069 C     !INPUT/OUTPUT PARAMETERS:
                0070 C     == Routine arguments ==
0b191c5f5a Jean*0071 C     bi, bj,    :: tile indices
                0072 C     recip_hFac :: reciprocal of cell open-depth factor (@ next iter)
                0073 C     uFld,vFld  :: Local copy of horizontal velocity field
                0074 C     wFld       :: Local copy of vertical velocity field
                0075 C     KappaRk    :: Vertical diffusion for Tempertature
                0076 C     myTime     :: current time
                0077 C     myIter     :: current iteration number
                0078 C     myThid     :: my Thread Id. number
fb8688a588 Jean*0079       INTEGER bi, bj
0b191c5f5a Jean*0080       _RS recip_hFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0081       _RL uFld      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0082       _RL vFld      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0083       _RL wFld      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0084       _RL KappaRk   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
4794324d0b Jean*0085       _RL     myTime
                0086       INTEGER myIter
                0087       INTEGER myThid
                0088 CEOP
                0089
                0090 #ifdef ALLOW_GENERIC_ADVDIFF
31067d68ad Jean*0091 #ifdef ALLOW_DIAGNOSTICS
                0092 C     !FUNCTIONS:
                0093       LOGICAL  DIAGNOSTICS_IS_ON
                0094       EXTERNAL DIAGNOSTICS_IS_ON
                0095 #endif /* ALLOW_DIAGNOSTICS */
                0096
4794324d0b Jean*0097 C     !LOCAL VARIABLES:
0b191c5f5a Jean*0098 C     iMin, iMax :: 1rst index loop range
                0099 C     jMin, jMax :: 2nd  index loop range
                0100 C     k          :: vertical index
                0101 C     kM1        :: =k-1 for k>1, =1 for k=1
                0102 C     kUp        :: index into 2 1/2D array, toggles between 1|2
                0103 C     kDown      :: index into 2 1/2D array, toggles between 2|1
                0104 C     xA         :: Tracer cell face area normal to X
                0105 C     yA         :: Tracer cell face area normal to X
                0106 C     maskUp     :: Land/water mask for Wvel points (interface k)
                0107 C     uTrans     :: Zonal volume transport through cell face
                0108 C     vTrans     :: Meridional volume transport through cell face
                0109 C     rTrans     ::   Vertical volume transport at interface k
                0110 C     rTransKp   :: Vertical volume transport at inteface k+1
                0111 C     fZon       :: Flux of temperature (T) in the zonal direction
                0112 C     fMer       :: Flux of temperature (T) in the meridional direction
                0113 C     fVer       :: Flux of temperature (T) in the vertical direction
                0114 C                   at the upper(U) and lower(D) faces of a cell.
1ece1facde Jean*0115 C     gT_loc     :: Temperature tendency (local to this S/R)
31067d68ad Jean*0116 C     gtForc     :: Temperature forcing tendency
e874fa47e5 Jean*0117 C     gt_AB      :: Adams-Bashforth temperature tendency increment
0b191c5f5a Jean*0118 C   useVariableK :: T when vertical diffusion is not constant
fb8688a588 Jean*0119       INTEGER iMin, iMax, jMin, jMax
4794324d0b Jean*0120       INTEGER i, j, k
                0121       INTEGER kUp, kDown, kM1
                0122       _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0123       _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0124       _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0125       _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0126       _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0127       _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0128       _RL rTransKp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0b191c5f5a Jean*0129       _RL fZon    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0130       _RL fMer    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0131       _RL fVer    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
1ece1facde Jean*0132       _RL gT_loc  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
31067d68ad Jean*0133       _RL gtForc  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
4794324d0b Jean*0134       _RL gt_AB   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31067d68ad Jean*0135 #ifdef ALLOW_DIAGNOSTICS
                0136       LOGICAL diagForcing, diagAB_tend
                0137 #endif
4794324d0b Jean*0138       LOGICAL calcAdvection
                0139       INTEGER iterNb
                0140 #ifdef ALLOW_ADAMSBASHFORTH_3
972c0130ec Jean*0141       INTEGER m2
4794324d0b Jean*0142 #endif
0b191c5f5a Jean*0143 #ifdef ALLOW_TIMEAVE
                0144       LOGICAL useVariableK
                0145 #endif
7c50f07931 Mart*0146 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0147 C     tkey :: tape key (tile dependent)
                0148 C     kkey :: tape key (level and tile dependent)
                0149       INTEGER tkey, kkey
7c50f07931 Mart*0150 #endif
0b191c5f5a Jean*0151
fb8688a588 Jean*0152 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0153
                0154       iterNb = myIter
0b191c5f5a Jean*0155       IF (staggerTimeStep) iterNb = myIter - 1
4794324d0b Jean*0156
d3548146a8 Jean*0157 C-    Loop ranges for daughter routines:
                0158 c     iMin = 1
                0159 c     iMax = sNx
                0160 c     jMin = 1
                0161 c     jMax = sNy
                0162 C     Regarding model dynamics, only needs to get correct tracer tendency
                0163 C     (gT_loc) in tile interior (1:sNx,1:sNy);
                0164 C     However, for some diagnostics, we may want to get valid tendency
                0165 C     extended over 1 point in tile halo region (0:sNx+1,0:sNy=1).
                0166       iMin = 0
                0167       iMax = sNx+1
                0168       jMin = 0
                0169       jMax = sNy+1
                0170
31067d68ad Jean*0171 #ifdef ALLOW_DIAGNOSTICS
                0172       diagForcing = .FALSE.
                0173       diagAB_tend = .FALSE.
                0174       IF ( useDiagnostics .AND. tempForcing )
                0175      &     diagForcing = DIAGNOSTICS_IS_ON( 'gT_Forc ', myThid )
                0176       IF ( useDiagnostics .AND. AdamsBashforthGt )
                0177      &     diagAB_tend = DIAGNOSTICS_IS_ON( 'AB_gT   ', myThid )
                0178 #endif
                0179
4794324d0b Jean*0180 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0181       tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
4794324d0b Jean*0182 #endif /* ALLOW_AUTODIFF_TAMC */
                0183
ff1bd1eafd Jean*0184 C-    Apply AB on T :
                0185       IF ( AdamsBashforth_T ) THEN
                0186 C     compute T^n+1/2 (stored in gtNm) extrapolating T forward in time
                0187 #ifdef ALLOW_ADAMSBASHFORTH_3
                0188 c         m1 = 1 + MOD(iterNb+1,2)
                0189 c         m2 = 1 + MOD( iterNb ,2)
                0190           CALL ADAMS_BASHFORTH3(
                0191      I                           bi, bj, 0, Nr,
                0192      I                           theta(1-OLx,1-OLy,1,bi,bj),
                0193      U                           gtNm, gt_AB,
                0194      I                           tempStartAB, iterNb, myThid )
                0195 #else /* ALLOW_ADAMSBASHFORTH_3 */
                0196           CALL ADAMS_BASHFORTH2(
                0197      I                           bi, bj, 0, Nr,
                0198      I                           theta(1-OLx,1-OLy,1,bi,bj),
                0199      U                           gtNm1(1-OLx,1-OLy,1,bi,bj), gt_AB,
                0200      I                           tempStartAB, iterNb, myThid )
                0201 #endif /* ALLOW_ADAMSBASHFORTH_3 */
                0202       ENDIF
                0203
fb8688a588 Jean*0204 C-    Tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
                0205       DO k=1,Nr
                0206        DO j=1-OLy,sNy+OLy
                0207         DO i=1-OLx,sNx+OLx
1ece1facde Jean*0208          gT_loc(i,j,k) = 0. _d 0
fb8688a588 Jean*0209         ENDDO
                0210        ENDDO
                0211       ENDDO
4794324d0b Jean*0212       DO j=1-OLy,sNy+OLy
                0213        DO i=1-OLx,sNx+OLx
0b191c5f5a Jean*0214          fVer(i,j,1) = 0. _d 0
                0215          fVer(i,j,2) = 0. _d 0
4794324d0b Jean*0216        ENDDO
                0217       ENDDO
fb8688a588 Jean*0218 #ifdef ALLOW_AUTODIFF
                0219       DO k=1,Nr
                0220        DO j=1-OLy,sNy+OLy
                0221         DO i=1-OLx,sNx+OLx
                0222          kappaRk(i,j,k) = 0. _d 0
                0223         ENDDO
                0224        ENDDO
                0225       ENDDO
4e4ad91a39 Jean*0226 #endif /* ALLOW_AUTODIFF */
                0227 #ifdef ALLOW_AUTODIFF_TAMC
598aebfcee Mart*0228 CADJ STORE wFld(:,:,:)         = comlev1_bibj, key=tkey, kind=isbyte
                0229 CADJ STORE theta(:,:,:,bi,bj)  = comlev1_bibj, key=tkey, kind=isbyte
ff1bd1eafd Jean*0230 # ifdef ALLOW_ADAMSBASHFORTH_3
598aebfcee Mart*0231 CADJ STORE gtNm(:,:,:,bi,bj,1) = comlev1_bibj, key=tkey, kind=isbyte
                0232 CADJ STORE gtNm(:,:,:,bi,bj,2) = comlev1_bibj, key=tkey, kind=isbyte
ff1bd1eafd Jean*0233 # else
598aebfcee Mart*0234 CADJ STORE gtNm1(:,:,:,bi,bj)  = comlev1_bibj, key=tkey, kind=isbyte
ff1bd1eafd Jean*0235 # endif
4e4ad91a39 Jean*0236 #endif /* ALLOW_AUTODIFF_TAMC */
4794324d0b Jean*0237
fb8688a588 Jean*0238 #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
0b191c5f5a Jean*0239       CALL CALC_3D_DIFFUSIVITY(
fb8688a588 Jean*0240      I         bi, bj, iMin, iMax, jMin, jMax,
                0241      I         GAD_TEMPERATURE, useGMredi, useKPP,
                0242      O         kappaRk,
                0243      I         myThid )
d9efc637ba Mart*0244 # ifdef ALLOW_AUTODIFF_TAMC
                0245 CADJ STORE kappaRk = comlev1_bibj, key = tkey, kind = isbyte
                0246 # endif
fb8688a588 Jean*0247 #endif /* INCLUDE_CALC_DIFFUSIVITY_CALL */
                0248
                0249 #ifndef DISABLE_MULTIDIM_ADVECTION
                0250 C--     Some advection schemes are better calculated using a multi-dimensional
                0251 C       method in the absence of any other terms and, if used, is done here.
                0252 C
                0253 C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
                0254 C The default is to use multi-dimensinal advection for non-linear advection
                0255 C schemes. However, for the sake of efficiency of the adjoint it is necessary
                0256 C to be able to exclude this scheme to avoid excessive storage and
                0257 C recomputation. It *is* differentiable, if you need it.
                0258 C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
                0259 C disable this section of code.
                0260 #ifdef GAD_ALLOW_TS_SOM_ADV
                0261 # ifdef ALLOW_AUTODIFF_TAMC
598aebfcee Mart*0262 C This looks awkward from an F77 point of view, but TAF does the right thing.
                0263 CADJ STORE som_T(:,:,:,bi,bj,:) = comlev1_bibj, key=tkey, kind=isbyte
fb8688a588 Jean*0264 # endif
                0265       IF ( tempSOM_Advection ) THEN
                0266 # ifdef ALLOW_DEBUG
                0267         IF (debugMode) CALL DEBUG_CALL('GAD_SOM_ADVECT',myThid)
                0268 # endif
                0269         CALL GAD_SOM_ADVECT(
e874fa47e5 Jean*0270      I             tempImplVertAdv,
                0271      I             tempAdvScheme, tempVertAdvScheme, GAD_TEMPERATURE,
                0272      I             dTtracerLev, uFld, vFld, wFld, theta,
fb8688a588 Jean*0273      U             som_T,
1ece1facde Jean*0274      O             gT_loc,
fb8688a588 Jean*0275      I             bi, bj, myTime, myIter, myThid )
                0276       ELSEIF (tempMultiDimAdvec) THEN
                0277 #else /* GAD_ALLOW_TS_SOM_ADV */
                0278       IF (tempMultiDimAdvec) THEN
                0279 #endif /* GAD_ALLOW_TS_SOM_ADV */
                0280 # ifdef ALLOW_DEBUG
                0281         IF (debugMode) CALL DEBUG_CALL('GAD_ADVECTION',myThid)
                0282 # endif
                0283         CALL GAD_ADVECTION(
e874fa47e5 Jean*0284      I             tempImplVertAdv,
                0285      I             tempAdvScheme, tempVertAdvScheme, GAD_TEMPERATURE,
                0286      I             dTtracerLev, uFld, vFld, wFld, theta,
1ece1facde Jean*0287      O             gT_loc,
fb8688a588 Jean*0288      I             bi, bj, myTime, myIter, myThid )
                0289       ENDIF
                0290 #endif /* DISABLE_MULTIDIM_ADVECTION */
                0291
                0292 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0293
                0294 C-    Start vertical index (k) loop (Nr:1)
                0295       calcAdvection = tempAdvection .AND. .NOT.tempMultiDimAdvec
4794324d0b Jean*0296       DO k=Nr,1,-1
                0297 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0298         kkey = k + (tkey-1)*Nr
4794324d0b Jean*0299 #endif /* ALLOW_AUTODIFF_TAMC */
                0300         kM1  = MAX(1,k-1)
                0301         kUp  = 1+MOD(k+1,2)
                0302         kDown= 1+MOD(k,2)
                0303
                0304 #ifdef ALLOW_AUTODIFF_TAMC
d9efc637ba Mart*0305 CADJ STORE fVer(:,:,:)         = comlev1_bibj_k, key=kkey, kind = isbyte
                0306 CADJ STORE gT_loc(:,:,k)       = comlev1_bibj_k, key=kkey, kind = isbyte
4794324d0b Jean*0307 # ifdef ALLOW_ADAMSBASHFORTH_3
d9efc637ba Mart*0308 CADJ STORE gtNm(:,:,k,bi,bj,1) = comlev1_bibj_k, key=kkey, kind = isbyte
                0309 CADJ STORE gtNm(:,:,k,bi,bj,2) = comlev1_bibj_k, key=kkey, kind = isbyte
4794324d0b Jean*0310 # else
d9efc637ba Mart*0311 CADJ STORE gtNm1(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, kind = isbyte
4794324d0b Jean*0312 # endif
                0313 #endif /* ALLOW_AUTODIFF_TAMC */
                0314         CALL CALC_ADV_FLOW(
                0315      I                uFld, vFld, wFld,
                0316      U                rTrans,
                0317      O                uTrans, vTrans, rTransKp,
                0318      O                maskUp, xA, yA,
                0319      I                k, bi, bj, myThid )
                0320
31067d68ad Jean*0321 C--   Collect forcing term in local array gtForc:
                0322         DO j=1-OLy,sNy+OLy
                0323          DO i=1-OLx,sNx+OLx
                0324           gtForc(i,j) = 0. _d 0
                0325          ENDDO
                0326         ENDDO
                0327         IF ( tempForcing ) THEN
                0328           CALL APPLY_FORCING_T(
                0329      U                        gtForc,
                0330      I                        iMin,iMax,jMin,jMax, k, bi,bj,
                0331      I                        myTime, myIter, myThid )
                0332 #ifdef ALLOW_DIAGNOSTICS
                0333           IF ( diagForcing ) THEN
                0334             CALL DIAGNOSTICS_FILL(gtForc,'gT_Forc ',k,1,2,bi,bj,myThid)
                0335           ENDIF
                0336 #endif /* ALLOW_DIAGNOSTICS */
                0337         ENDIF
                0338
4794324d0b Jean*0339 #ifdef ALLOW_ADAMSBASHFORTH_3
972c0130ec Jean*0340 c       m1 = 1 + MOD(iterNb+1,2)
4794324d0b Jean*0341         m2 = 1 + MOD( iterNb ,2)
                0342         CALL GAD_CALC_RHS(
                0343      I           bi, bj, iMin,iMax,jMin,jMax, k, kM1, kUp, kDown,
                0344      I           xA, yA, maskUp, uFld(1-OLx,1-OLy,k),
                0345      I           vFld(1-OLx,1-OLy,k), wFld(1-OLx,1-OLy,k),
                0346      I           uTrans, vTrans, rTrans, rTransKp,
                0347      I           diffKhT, diffK4T, KappaRk(1-OLx,1-OLy,k), diffKr4T,
cb30c1718d Jean*0348      I           theta(1-OLx,1-OLy,1,bi,bj),
                0349      I           gtNm(1-OLx,1-OLy,1,bi,bj,m2), dTtracerLev,
4794324d0b Jean*0350      I           GAD_TEMPERATURE, tempAdvScheme, tempVertAdvScheme,
                0351      I           calcAdvection, tempImplVertAdv, AdamsBashforth_T,
46918f1b26 Jean*0352      I           tempVertDiff4, useGMRedi, useKPP, temp_stayPositive,
0b191c5f5a Jean*0353      O           fZon, fMer,
1ece1facde Jean*0354      U           fVer, gT_loc,
4794324d0b Jean*0355      I           myTime, myIter, myThid )
                0356 #else /* ALLOW_ADAMSBASHFORTH_3 */
                0357         CALL GAD_CALC_RHS(
                0358      I           bi, bj, iMin,iMax,jMin,jMax, k, kM1, kUp, kDown,
                0359      I           xA, yA, maskUp, uFld(1-OLx,1-OLy,k),
                0360      I           vFld(1-OLx,1-OLy,k), wFld(1-OLx,1-OLy,k),
                0361      I           uTrans, vTrans, rTrans, rTransKp,
                0362      I           diffKhT, diffK4T, KappaRk(1-OLx,1-OLy,k), diffKr4T,
cb30c1718d Jean*0363      I           theta(1-OLx,1-OLy,1,bi,bj),
                0364      I           gtNm1(1-OLx,1-OLy,1,bi,bj), dTtracerLev,
4794324d0b Jean*0365      I           GAD_TEMPERATURE, tempAdvScheme, tempVertAdvScheme,
                0366      I           calcAdvection, tempImplVertAdv, AdamsBashforth_T,
46918f1b26 Jean*0367      I           tempVertDiff4, useGMRedi, useKPP, temp_stayPositive,
0b191c5f5a Jean*0368      O           fZon, fMer,
1ece1facde Jean*0369      U           fVer, gT_loc,
4794324d0b Jean*0370      I           myTime, myIter, myThid )
46918f1b26 Jean*0371 #endif /* ALLOW_ADAMSBASHFORTH_3 */
4794324d0b Jean*0372
                0373 C--   External thermal forcing term(s) inside Adams-Bashforth:
31067d68ad Jean*0374         IF ( tempForcing .AND. tracForcingOutAB.NE.1 ) THEN
                0375           DO j=1-OLy,sNy+OLy
                0376            DO i=1-OLx,sNx+OLx
1ece1facde Jean*0377             gT_loc(i,j,k) = gT_loc(i,j,k) + gtForc(i,j)
31067d68ad Jean*0378            ENDDO
                0379           ENDDO
                0380         ENDIF
4794324d0b Jean*0381
                0382         IF ( AdamsBashforthGt ) THEN
                0383 #ifdef ALLOW_ADAMSBASHFORTH_3
                0384           CALL ADAMS_BASHFORTH3(
                0385      I                          bi, bj, k, Nr,
1ece1facde Jean*0386      U                          gT_loc, gtNm,
                0387      O                          gt_AB,
4794324d0b Jean*0388      I                          tempStartAB, iterNb, myThid )
                0389 #else
                0390           CALL ADAMS_BASHFORTH2(
                0391      I                          bi, bj, k, Nr,
1ece1facde Jean*0392      U                          gT_loc, gtNm1(1-OLx,1-OLy,1,bi,bj),
                0393      O                          gt_AB,
4794324d0b Jean*0394      I                          tempStartAB, iterNb, myThid )
                0395 #endif
                0396 #ifdef ALLOW_DIAGNOSTICS
31067d68ad Jean*0397           IF ( diagAB_tend ) THEN
4794324d0b Jean*0398             CALL DIAGNOSTICS_FILL(gt_AB,'AB_gT   ',k,1,2,bi,bj,myThid)
                0399           ENDIF
                0400 #endif /* ALLOW_DIAGNOSTICS */
                0401         ENDIF
                0402
                0403 C--   External thermal forcing term(s) outside Adams-Bashforth:
31067d68ad Jean*0404         IF ( tempForcing .AND. tracForcingOutAB.EQ.1 ) THEN
                0405           DO j=1-OLy,sNy+OLy
                0406            DO i=1-OLx,sNx+OLx
1ece1facde Jean*0407             gT_loc(i,j,k) = gT_loc(i,j,k) + gtForc(i,j)
31067d68ad Jean*0408            ENDDO
                0409           ENDDO
                0410         ENDIF
4794324d0b Jean*0411
                0412 #ifdef NONLIN_FRSURF
d9efc637ba Mart*0413 # ifdef ALLOW_AUTODIFF_TAMC
                0414 CADJ STORE gT_loc(:,:,k) = comlev1_bibj_k, key = kkey, kind = isbyte
                0415 # endif
4794324d0b Jean*0416         IF (nonlinFreeSurf.GT.0) THEN
                0417           CALL FREESURF_RESCALE_G(
                0418      I                            bi, bj, k,
1ece1facde Jean*0419      U                            gT_loc,
4794324d0b Jean*0420      I                            myThid )
                0421          IF ( AdamsBashforthGt ) THEN
                0422 #ifdef ALLOW_ADAMSBASHFORTH_3
                0423 # ifdef ALLOW_AUTODIFF_TAMC
d9efc637ba Mart*0424 CADJ STORE gtNm(:,:,k,bi,bj,1) = comlev1_bibj_k, key=kkey, kind=isbyte
                0425 CADJ STORE gtNm(:,:,k,bi,bj,2) = comlev1_bibj_k, key=kkey, kind=isbyte
4794324d0b Jean*0426 # endif
                0427           CALL FREESURF_RESCALE_G(
                0428      I                            bi, bj, k,
1ece1facde Jean*0429      U                            gtNm(1-OLx,1-OLy,1,bi,bj,1),
4794324d0b Jean*0430      I                            myThid )
                0431           CALL FREESURF_RESCALE_G(
                0432      I                            bi, bj, k,
1ece1facde Jean*0433      U                            gtNm(1-OLx,1-OLy,1,bi,bj,2),
4794324d0b Jean*0434      I                            myThid )
                0435 #else
d9efc637ba Mart*0436 # ifdef ALLOW_AUTODIFF_TAMC
                0437 CADJ STORE gtNm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, kind = isbyte
                0438 # endif
4794324d0b Jean*0439           CALL FREESURF_RESCALE_G(
                0440      I                            bi, bj, k,
1ece1facde Jean*0441      U                            gtNm1(1-OLx,1-OLy,1,bi,bj),
4794324d0b Jean*0442      I                            myThid )
                0443 #endif
                0444          ENDIF
                0445         ENDIF
                0446 #endif /* NONLIN_FRSURF */
                0447
                0448 C-    end of vertical index (k) loop (Nr:1)
                0449       ENDDO
                0450
0b191c5f5a Jean*0451 #ifdef ALLOW_DOWN_SLOPE
                0452       IF ( useDOWN_SLOPE ) THEN
d9efc637ba Mart*0453 # ifdef ALLOW_AUTODIFF_TAMC
                0454 CADJ STORE recip_hFac = comlev1_bibj, key = tkey, kind = isbyte
                0455 # endif
0b191c5f5a Jean*0456         IF ( usingPCoords ) THEN
                0457           CALL DWNSLP_APPLY(
                0458      I                  GAD_TEMPERATURE, bi, bj, kSurfC,
6fca64f237 Jean*0459      I                  theta(1-OLx,1-OLy,1,bi,bj),
1ece1facde Jean*0460      U                  gT_loc,
6fca64f237 Jean*0461      I                  recip_hFac, recip_rA, recip_drF,
                0462      I                  dTtracerLev, myTime, myIter, myThid )
0b191c5f5a Jean*0463         ELSE
                0464           CALL DWNSLP_APPLY(
                0465      I                  GAD_TEMPERATURE, bi, bj, kLowC,
6fca64f237 Jean*0466      I                  theta(1-OLx,1-OLy,1,bi,bj),
1ece1facde Jean*0467      U                  gT_loc,
6fca64f237 Jean*0468      I                  recip_hFac, recip_rA, recip_drF,
                0469      I                  dTtracerLev, myTime, myIter, myThid )
0b191c5f5a Jean*0470         ENDIF
                0471       ENDIF
                0472 #endif /* ALLOW_DOWN_SLOPE */
                0473
1ece1facde Jean*0474 C-    Integrate forward in time, storing in gT_loc:  gT <= T + dt*gT
6fca64f237 Jean*0475       CALL TIMESTEP_TRACER(
                0476      I                  bi, bj, dTtracerLev,
                0477      I                  theta(1-OLx,1-OLy,1,bi,bj),
1ece1facde Jean*0478      U                  gT_loc,
6fca64f237 Jean*0479      I                  myTime, myIter, myThid )
                0480
0b191c5f5a Jean*0481 C--   Implicit vertical advection & diffusion
                0482
d9efc637ba Mart*0483 #ifdef ALLOW_AUTODIFF_TAMC
598aebfcee Mart*0484 CADJ STORE gT_loc (:,:,:)     = comlev1_bibj, key = tkey, kind = isbyte
d9efc637ba Mart*0485 #endif /* ALLOW_AUTODIFF_TAMC */
0b191c5f5a Jean*0486 #ifdef INCLUDE_IMPLVERTADV_CODE
5685a5c914 Jean*0487       IF ( tempImplVertAdv .OR. implicitDiffusion ) THEN
                0488 C     to recover older (prior to 2016-10-05) results:
                0489 c     IF ( tempImplVertAdv ) THEN
0b191c5f5a Jean*0490 #ifdef ALLOW_AUTODIFF_TAMC
598aebfcee Mart*0491 CADJ STORE wFld(:,:,:)        = comlev1_bibj, key = tkey, kind = isbyte
                0492 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = tkey, kind = isbyte
                0493 CADJ STORE recip_hFac(:,:,:)  = comlev1_bibj, key = tkey, kind = isbyte
0b191c5f5a Jean*0494 #endif /* ALLOW_AUTODIFF_TAMC */
                0495         CALL GAD_IMPLICIT_R(
                0496      I         tempImplVertAdv, tempVertAdvScheme, GAD_TEMPERATURE,
1ece1facde Jean*0497      I         dTtracerLev, kappaRk, recip_hFac, wFld,
                0498      I         theta(1-OLx,1-OLy,1,bi,bj),
                0499      U         gT_loc,
0b191c5f5a Jean*0500      I         bi, bj, myTime, myIter, myThid )
                0501       ELSEIF ( implicitDiffusion ) THEN
                0502 #else /* INCLUDE_IMPLVERTADV_CODE */
                0503       IF     ( implicitDiffusion ) THEN
                0504 #endif /* INCLUDE_IMPLVERTADV_CODE */
                0505         CALL IMPLDIFF(
                0506      I         bi, bj, iMin, iMax, jMin, jMax,
                0507      I         GAD_TEMPERATURE, kappaRk, recip_hFac,
1ece1facde Jean*0508      U         gT_loc,
0b191c5f5a Jean*0509      I         myThid )
                0510       ENDIF
                0511
                0512 #ifdef ALLOW_TIMEAVE
d8d1486ca1 Jean*0513       useVariableK = useKPP .OR. usePP81 .OR. useKL10 .OR. useMY82
                0514      &       .OR. useGGL90 .OR. useGMredi .OR. ivdc_kappa.NE.0.
0b191c5f5a Jean*0515       IF ( taveFreq.GT.0. .AND. useVariableK
                0516      &                    .AND.implicitDiffusion ) THEN
1ece1facde Jean*0517         CALL TIMEAVE_CUMUL_DIF_1T( TdiffRtave,
                0518      I                        gT_loc, kappaRk,
                0519      I                        Nr, 3, deltaTClock, bi, bj, myThid )
0b191c5f5a Jean*0520       ENDIF
                0521 #endif /* ALLOW_TIMEAVE */
                0522
bd27360393 Jean*0523       IF ( AdamsBashforth_T ) THEN
                0524 C-    Save current tracer field (for AB on tracer) and then update tracer
ff1bd1eafd Jean*0525 #ifdef ALLOW_ADAMSBASHFORTH_3
bd27360393 Jean*0526         CALL CYCLE_AB_TRACER(
cb30c1718d Jean*0527      I             bi, bj, gT_loc,
                0528      U             theta(1-OLx,1-OLy,1,bi,bj),
                0529      O             gtNm(1-OLx,1-OLy,1,bi,bj,m2),
bd27360393 Jean*0530      I             myTime, myIter, myThid )
                0531 #else /* ALLOW_ADAMSBASHFORTH_3 */
ff1bd1eafd Jean*0532         CALL CYCLE_AB_TRACER(
                0533      I             bi, bj, gT_loc,
                0534      U             theta(1-OLx,1-OLy,1,bi,bj),
                0535      O             gtNm1(1-OLx,1-OLy,1,bi,bj),
                0536      I             myTime, myIter, myThid )
bd27360393 Jean*0537 #endif /* ALLOW_ADAMSBASHFORTH_3 */
ff1bd1eafd Jean*0538       ELSE
bd27360393 Jean*0539 C-    Update tracer fields:  T(n) = T**
                0540         CALL CYCLE_TRACER(
                0541      I             bi, bj,
cb30c1718d Jean*0542      O             theta(1-OLx,1-OLy,1,bi,bj),
                0543      I             gT_loc, myTime, myIter, myThid )
bd27360393 Jean*0544       ENDIF
                0545
4794324d0b Jean*0546 #endif /* ALLOW_GENERIC_ADVDIFF */
                0547
                0548       RETURN
                0549       END