C $Header: /u/gcmpack/MITgcm/pkg/longstep/longstep_thermodynamics.F,v 1.9 2013/12/27 16:01:16 jmc Exp $
C $Name:  $

#include "LONGSTEP_OPTIONS.h"

#ifdef ALLOW_AUTODIFF_TAMC
# ifdef ALLOW_GMREDI
#  include "GMREDI_OPTIONS.h"
# endif
#endif /* ALLOW_AUTODIFF_TAMC */

CBOP
C     !ROUTINE: LONGSTEP_THERMODYNAMICS
C     !INTERFACE:
      SUBROUTINE LONGSTEP_THERMODYNAMICS( myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE LONGSTEP_THERMODYNAMICS
C     | o Controlling routine for the prognostics of passive tracers
C     |   with longer time step.
C     *===========================================================
C     | This is a copy of THERMODYNAMICS, but only with the
C     | parts relevant to ptracers, and dynamics fields replaced
C     | by their longstep averages.
C     | When THERMODYNAMICS is changed, this routine probably has
C     | to be changed too :(
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"
#ifdef ALLOW_GENERIC_ADVDIFF
# include "GAD.h"
#endif
#include "LONGSTEP_PARAMS.h"
#include "LONGSTEP.h"
#ifdef ALLOW_PTRACERS
# include "PTRACERS_SIZE.h"
# include "PTRACERS_PARAMS.h"
# include "PTRACERS_FIELDS.h"
#endif

#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
# include "FFIELDS.h"
# include "EOS.h"
# ifdef ALLOW_GMREDI
#  include "GMREDI.h"
# endif
# ifdef ALLOW_EBM
#  include "EBM.h"
# endif
#endif /* ALLOW_AUTODIFF_TAMC */

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_LONGSTEP
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     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         :: 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 iMin, iMax
      INTEGER jMin, jMax
      INTEGER bi, bj
      INTEGER i, j, k
CEOP

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

C     time for a ptracer time step?
      IF ( LS_doTimeStep ) THEN

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

      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
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 */

        iMin = 1-OLx
        iMax = sNx+OLx
        jMin = 1-OLy
        jMax = sNy+OLy

C     need to recompute surfaceForcingPtr using LS_fwFlux
        CALL LONGSTEP_FORCING_SURF(
     I        bi, bj, iMin, iMax, jMin, jMax,
     I        myTime,myIter,myThid )

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) = LS_uVel(i,j,k,bi,bj)
           vFld(i,j,k) = LS_vVel(i,j,k,bi,bj)
           wFld(i,j,k) = LS_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_PTRACERS
C--     Calculate passive tracer tendencies
C       and step forward, storing result in gPtr
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 /* ALLOW_PTRACERS */

C--   end bi,bj loops.
       ENDDO
      ENDDO

#ifdef ALLOW_DEBUG
      IF ( debugLevel.GE.debLevD ) THEN
       CALL DEBUG_STATS_RL(Nr,LS_uVel,'LS_Uvel (THERMODYNAMICS)',myThid)
       CALL DEBUG_STATS_RL(Nr,LS_vVel,'LS_Vvel (THERMODYNAMICS)',myThid)
       CALL DEBUG_STATS_RL(Nr,LS_wVel,'LS_Wvel (THERMODYNAMICS)',myThid)
       CALL DEBUG_STATS_RL(Nr,LS_theta,'LS_Theta (THERMODYNAMICS)',
     &                     myThid)
       CALL DEBUG_STATS_RL(Nr,LS_salt,'LS_Salt (THERMODYNAMICS)',myThid)
#ifdef ALLOW_PTRACERS
       IF ( usePTRACERS ) THEN
         CALL PTRACERS_DEBUG(myThid)
       ENDIF
#endif /* ALLOW_PTRACERS */
      ENDIF
#endif /* ALLOW_DEBUG */

C     LS_doTimeStep
      ENDIF

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

#endif /* ALLOW_LONGSTEP */

      RETURN
      END