#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif
#if (defined ALLOW_PTRACERS) && (!defined ALLOW_LONGSTEP)
# define DO_PTRACERS_HERE
#endif

#ifdef ALLOW_AUTODIFF
# ifdef ALLOW_GMREDI
#  include "GMREDI_OPTIONS.h"
# endif
# ifdef ALLOW_KPP
#  include "KPP_OPTIONS.h"
# endif
#endif /* ALLOW_AUTODIFF */

CBOP
C     !ROUTINE: THERMODYNAMICS
C     !INTERFACE:
      SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE THERMODYNAMICS
C     | o Controlling routine for the prognostic part of the
C     |   thermo-dynamics.
C     *===========================================================
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "RESTART.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "FFIELDS.h"
#ifdef ALLOW_GENERIC_ADVDIFF
# include "GAD.h"
#endif
#ifdef DO_PTRACERS_HERE
# include "PTRACERS_SIZE.h"
# include "PTRACERS_PARAMS.h"
# include "PTRACERS_FIELDS.h"
#endif

#ifdef ALLOW_AUTODIFF
# include "tamc.h"
# include "tamc_keys.h"
# include "EOS.h"
# ifdef ALLOW_KPP
#  include "KPP.h"
# endif
# ifdef ALLOW_GMREDI
#  include "GMREDI.h"
# endif
# ifdef ALLOW_EBM
#  include "EBM.h"
# endif
# ifdef ALLOW_SALT_PLUME
#  include "SALT_PLUME.h"
# endif
#endif /* ALLOW_AUTODIFF */

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myTime :: Current time in simulation
C     myIter :: Current iteration number in simulation
C     myThid :: Thread number for this instance of the routine.
      _RL myTime
      INTEGER myIter
      INTEGER myThid

#ifdef ALLOW_GENERIC_ADVDIFF
# ifdef ALLOW_MONITOR
C     !FUNCTIONS:
      LOGICAL  DIFFERENT_MULTIPLE
      EXTERNAL DIFFERENT_MULTIPLE
# endif /* ALLOW_MONITOR */

C     !LOCAL VARIABLES:
C     == Local variables
C     uFld,vFld,wFld :: Local copy of velocity field (3 components)
C     kappaRk        :: Total diffusion in vertical, all levels, 1 tracer
C     recip_hFacNew  :: reciprocal of futur time-step hFacC
C     bi, bj         :: Tile indices
C     i, j, k        :: loop indices
      _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL wFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL kappaRk (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RS recip_hFacNew(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      INTEGER bi, bj
      INTEGER i, j, k
# ifdef ALLOW_MONITOR
      LOGICAL monOutputCFL
      _RL wrTime
      _RL trAdvCFL(3,nSx,nSy)
# endif /* ALLOW_MONITOR */
CEOP

#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
#endif

# ifdef ALLOW_MONITOR
      monOutputCFL = .FALSE.
      IF ( monitorSelect.GE.2 ) THEN
        wrTime = myTime
        IF ( .NOT.staggerTimeStep ) wrTime = myTime + deltaTClock
        monOutputCFL =
     &       DIFFERENT_MULTIPLE( monitorFreq, wrTime, deltaTClock )
      ENDIF
# endif /* ALLOW_MONITOR */

#ifdef ALLOW_AUTODIFF_TAMC
C--   dummy statement to end declaration part
      ikey = 1
      itdkey = 1
#endif /* ALLOW_AUTODIFF_TAMC */

#ifdef ALLOW_AUTODIFF_TAMC
C--   HPF directive to help TAMC
CHPF$ INDEPENDENT
#endif /* ALLOW_AUTODIFF_TAMC */

C-- Compute correction at the surface for Lin Free Surf.
#ifdef ALLOW_AUTODIFF
      TsurfCor = 0. _d 0
      SsurfCor = 0. _d 0
#endif
      IF ( linFSConserveTr ) THEN
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE theta,salt,wvel = comlev1, key = ikey_dynamics, byte=isbyte
#endif
       CALL CALC_WSURF_TR( theta, salt, wVel,
     &                     myTime, myIter, myThid )
      ENDIF
#ifdef ALLOW_LAYERS
      IF ( useLayers ) THEN
        CALL LAYERS_WSURF_TR( theta, salt, wVel,
     &                     myTime, myIter, myThid )
      ENDIF
#endif /* ALLOW_LAYERS */

#ifdef DO_PTRACERS_HERE
#ifdef ALLOW_AUTODIFF
      DO k=1,PTRACERS_numInUse
        meanSurfCorPTr(k) = 0.0 _d 0
      ENDDO
#endif
      IF ( PTRACERS_calcSurfCor ) THEN
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE ptracer = comlev1, key = ikey_dynamics, byte=isbyte
#endif
       CALL PTRACERS_CALC_WSURF_TR(wVel,myTime,myIter,myThid)
      ENDIF
#endif /* DO_PTRACERS_HERE */

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)

#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
     &                      + act3*max1*max2
     &                      + act4*max1*max2*max3
#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 k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           recip_hFacNew(i,j,k) = 0. _d 0
C This is currently also used by IVDC and Diagnostics
           kappaRk(i,j,k)    = 0. _d 0
          ENDDO
         ENDDO
        ENDDO

C--     Compute new reciprocal hFac for implicit calculation
#ifdef NONLIN_FRSURF
        IF ( nonlinFreeSurf.GT.0 ) THEN
         IF ( select_rStar.GT.0 ) THEN
# ifndef DISABLE_RSTAR_CODE
          DO k=1,Nr
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             recip_hFacNew(i,j,k) = recip_hFacC(i,j,k,bi,bj)
     &                            / rStarExpC(i,j,bi,bj)
            ENDDO
           ENDDO
          ENDDO
# endif /* DISABLE_RSTAR_CODE */
         ELSEIF ( selectSigmaCoord.NE.0 ) THEN
# ifndef DISABLE_SIGMA_CODE
          DO k=1,Nr
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             recip_hFacNew(i,j,k) = recip_hFacC(i,j,k,bi,bj)
     &        /( 1. _d 0 + dEtaHdt(i,j,bi,bj)*deltaTFreeSurf
     &                    *dBHybSigF(k)*recip_drF(k)
     &                    *recip_hFacC(i,j,k,bi,bj)
     &         )
            ENDDO
           ENDDO
          ENDDO
# endif /* DISABLE_RSTAR_CODE */
         ELSE
          DO k=1,Nr
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             IF ( k.EQ.kSurfC(i,j,bi,bj) ) THEN
              recip_hFacNew(i,j,k) = 1. _d 0 / hFac_surfC(i,j,bi,bj)
             ELSE
              recip_hFacNew(i,j,k) = recip_hFacC(i,j,k,bi,bj)
             ENDIF
            ENDDO
           ENDDO
          ENDDO
         ENDIF
        ELSE
#endif /* NONLIN_FRSURF */
          DO k=1,Nr
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             recip_hFacNew(i,j,k) = _recip_hFacC(i,j,k,bi,bj)
            ENDDO
           ENDDO
          ENDDO
#ifdef NONLIN_FRSURF
        ENDIF
#endif /* NONLIN_FRSURF */

C--   Set up 3-D velocity field that we use to advect tracers:
C-    just do a local copy:
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           uFld(i,j,k) = uVel(i,j,k,bi,bj)
           vFld(i,j,k) = vVel(i,j,k,bi,bj)
           wFld(i,j,k) = wVel(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDDO
#ifdef ALLOW_GMREDI
C-    add Bolus velocity to Eulerian-mean velocity:
        IF (useGMRedi) THEN
          CALL GMREDI_RESIDUAL_FLOW(
     U                  uFld, vFld, wFld,
     I                  bi, bj, myIter, myThid )
        ENDIF
#endif /* ALLOW_GMREDI */
#ifdef ALLOW_MONITOR
        IF ( monOutputCFL  ) THEN
          CALL MON_CALC_ADVCFL_TILE( Nr, bi, bj,
     I                         uFld, vFld, wFld, dTtracerLev,
     O                         trAdvCFL(1,bi,bj),
     I                         myIter, myThid )
c        WRITE(standardMessageUnit,'(A,I8,2I3,A,1P3E14.6)')
c    &     ' trAdv_CFL: it,bi,bj=', myIter,bi,bj,
c    &     ' , CFL =', (trAdvCFL(i,bi,bj),i=1,3)
        ENDIF
#endif /* ALLOW_MONITOR */

C-    Apply AB on T,S : moved inside TEMP/SALT_INTEGRATE

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE recip_hFacNew(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
CADJ STORE uFld (:,:,:)       = comlev1_bibj, key=itdkey, byte=isbyte
CADJ STORE vFld (:,:,:)       = comlev1_bibj, key=itdkey, byte=isbyte
CADJ STORE wFld (:,:,:)       = comlev1_bibj, key=itdkey, byte=isbyte
CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
# ifdef ALLOW_SALT_PLUME
CADJ STORE saltPlumeFlux(:,:,bi,bj)   =
CADJ &     comlev1_bibj, key=itdkey,kind = isbyte
CADJ STORE saltPlumeDepth(:,:,bi,bj)   =
CADJ &     comlev1_bibj, key=itdkey,kind = isbyte
# endif
# if ((defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)) && (defined ALLOW_GMREDI)
#  ifdef GM_NON_UNITY_DIAGONAL
CADJ STORE kux(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
CADJ STORE kvy(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
#  endif
#  ifdef GM_EXTRA_DIAGONAL
CADJ STORE kuz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
CADJ STORE kvz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
#  endif
# endif
#endif /* ALLOW_AUTODIFF_TAMC */

C--     Calculate active tracer tendencies and step forward in time.
C       Active tracer arrays are updated while adjustments (filters,
C       conv.adjustment) are applied later in TRACERS_CORRECTION_STEP.

        IF ( tempStepping ) THEN
#ifdef ALLOW_DEBUG
          IF (debugMode) CALL DEBUG_CALL('TEMP_INTEGRATE',myThid)
#endif
          CALL TEMP_INTEGRATE(
     I         bi, bj, recip_hFacNew,
     I         uFld, vFld, wFld,
     U         kappaRk,
     I         myTime, myIter, myThid )
        ENDIF

        IF ( saltStepping ) THEN
#ifdef ALLOW_DEBUG
          IF (debugMode) CALL DEBUG_CALL('SALT_INTEGRATE',myThid)
#endif
          CALL SALT_INTEGRATE(
     I         bi, bj, recip_hFacNew,
     I         uFld, vFld, wFld,
     U         kappaRk,
     I         myTime, myIter, myThid )
        ENDIF

#ifdef DO_PTRACERS_HERE
C--     Calculate passive tracer tendencies and step forward in time.
C       Passive tracer arrays are updated while adjustments (filters,
C       conv.adjustment) are applied later in TRACERS_CORRECTION_STEP.
C       Also apply open boundary conditions for each passive tracer
        IF ( usePTRACERS ) THEN
#ifdef ALLOW_DEBUG
          IF (debugMode) CALL DEBUG_CALL('PTRACERS_INTEGRATE',myThid)
#endif
          CALL PTRACERS_INTEGRATE(
     I          bi, bj, recip_hFacNew,
     I          uFld, vFld, wFld,
     U          kappaRk,
     I          myTime, myIter, myThid )
        ENDIF
#endif /* DO_PTRACERS_HERE */

#ifdef   ALLOW_OBCS
C--   Apply open boundary conditions
        IF ( useOBCS ) THEN
          CALL OBCS_APPLY_TS( bi, bj, 0, theta, salt, myThid )
        ENDIF
#endif   /* ALLOW_OBCS */

#ifdef ALLOW_FRICTION_HEATING
#ifdef ALLOW_DIAGNOSTICS
        IF ( addFrictionHeating .AND. useDiagnostics ) THEN
          CALL DIAGNOSTICS_FILL_RS( frictionHeating, 'HeatDiss',
     &                              0, Nr, 1, bi, bj, myThid )
        ENDIF
#endif /* ALLOW_DIAGNOSTICS */
C-    Reset frictionHeating to zero
        IF ( addFrictionHeating .AND. .NOT.staggerTimeStep ) THEN
          DO k=1,Nr
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
              frictionHeating(i,j,k,bi,bj) = 0. _d 0
            ENDDO
           ENDDO
          ENDDO
        ENDIF
#endif /* ALLOW_FRICTION_HEATING */

C--   end bi,bj loops.
       ENDDO
      ENDDO

#ifdef ALLOW_MONITOR
      IF ( monOutputCFL ) THEN
        CALL MON_CALC_ADVCFL_GLOB(
     I                       trAdvCFL, myIter, myThid )
      ENDIF
#endif /* ALLOW_MONITOR */

#ifdef ALLOW_DEBUG
      IF ( debugLevel.GE.debLevD ) THEN
       CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
       CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
       CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
       CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
       CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
#ifndef ALLOW_ADAMSBASHFORTH_3
       CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
       CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
#endif
#ifdef DO_PTRACERS_HERE
       IF ( usePTRACERS ) THEN
         CALL PTRACERS_DEBUG(myThid)
       ENDIF
#endif /* DO_PTRACERS_HERE */
      ENDIF
#endif /* ALLOW_DEBUG */

#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
#endif

#endif /* ALLOW_GENERIC_ADVDIFF */

      RETURN
      END
