! $Id$
!
!  This module can replace the energy module by using the thermal energy
!  eth as dependent variable. For a perfect gas we have
!
!    deth/dt + div(u*eth) = -P*div(u)
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lentropy = .false.
! CPARAM logical, parameter :: ltemperature = .false.
! CPARAM logical, parameter :: lthermal_energy = .true.
!
! MVAR CONTRIBUTION 1
!
! PENCILS PROVIDED Ma2; fpres(3); ugeths; transpeth; sglnTT(3)
!
!***************************************************************
module Energy
!
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages
!
  implicit none
!
  include 'energy.h'
!
  real :: eth_left, eth_right, widtheth, eth_const=1.0
  real :: chi=0.0, chi_shock=0.0, chi_hyper3_mesh=0.
  real :: kappa_rosseland=1.0
  real :: energy_floor = 0.0, temperature_floor = 0.0
  real :: rk_eps = 1.e-3
  real :: nu_z=0., cs_z=0.
  real :: TTref = 0., tau_cool = 1.
  real :: xjump_mid=0.,yjump_mid=0.,zjump_mid=0.
  integer :: rk_nmax = 100
  integer :: njeans = 4
  logical :: lsplit_update=.false.
  logical :: lviscosity_heat=.true.
  logical :: lupw_eth=.false.
  logical :: lcheck_negative_energy=.false.
  logical :: lKI02 = .false.
  logical :: lSD93 = .false.
  logical :: lconst_cooling_time = .false.
  logical :: ljeans_floor=.false.
  logical :: lchi_rosseland=.false.
  logical, pointer :: lpressuregradient_gas
  logical :: lenergy_slope_limited=.false.
  character (len=labellen), dimension(ninit) :: initeth='nothing'
  character(len=labellen) :: feedback = 'linear'
!
!  Input parameters.
!
  namelist /entropy_init_pars/ &
      initeth, eth_left, eth_right, widtheth, eth_const,&
      xjump_mid, yjump_mid, zjump_mid
!
!  Run parameters.
!
  namelist /entropy_run_pars/ &
      lviscosity_heat, lupw_eth, &
      chi, chi_shock, chi_hyper3_mesh, &
      lchi_rosseland, kappa_rosseland, &
      energy_floor, temperature_floor, lcheck_negative_energy, &
      rk_eps, rk_nmax, lKI02, lSD93, nu_z, cs_z, &
      lconst_cooling_time, TTref, tau_cool, &
      ljeans_floor, njeans, w_sldchar_ene
!
!  Diagnostic variables (needs to be consistent with reset list below).
!
  integer :: idiag_TTmax=0    ! DIAG_DOC: $\max (T)$
  integer :: idiag_TTmin=0    ! DIAG_DOC: $\min (T)$
  integer :: idiag_ppm=0      ! DIAG_DOC: $\left< p \right>$
  integer :: idiag_TTm=0      ! DIAG_DOC: $\left<T\right>$
  integer :: idiag_pdivum=0
  integer :: idiag_ethm=0     ! DIAG_DOC: $\left< e_{\text{th}}\right> =
                              ! DIAG_DOC:  \left< c_v \rho T \right> $
                              ! DIAG_DOC: \quad(mean thermal energy)
  integer :: idiag_ethtot=0   ! DIAG_DOC: $\int_V e_{\text{th}}\,dV$
                              ! DIAG_DOC:   \quad(total thermal energy)
  integer :: idiag_ethmin=0   ! DIAG_DOC: $\mathrm{min} e_\text{th}$
  integer :: idiag_ethmax=0   ! DIAG_DOC: $\mathrm{max} e_\text{th}$
  integer :: idiag_eem=0      ! DIAG_DOC: $\left< e \right> =
                              ! DIAG_DOC:  \left< c_v T \right>$
                              ! DIAG_DOC: \quad(mean internal energy)
  integer :: idiag_etot=0     ! DIAG_DOC: $\langle e_\textrm{th} + \rho u^2 / 2\rangle$
!
! xy averaged diagnostics given in xyaver.in
!
  integer :: idiag_ppmz=0     ! XYAVG_DOC: $\left<p\right>_{xy}$
  integer :: idiag_TTmz=0     ! XYAVG_DOC: $\left<T\right>_{xy}$
!
! xz averaged diagnostics given in xzaver.in
!
  integer :: idiag_ppmy=0     ! XZAVG_DOC: $\left<p\right>_{xz}$
  integer :: idiag_TTmy=0     ! XZAVG_DOC: $\left<T\right>_{xz}$
!
! yz averaged diagnostics given in yzaver.in
!
  integer :: idiag_ppmx=0     ! YZAVG_DOC: $\left<p\right>_{yz}$
  integer :: idiag_TTmx=0     ! YZAVG_DOC: $\left<T\right>_{yz}$
!
! y averaged diagnostics given in yaver.in
!
  integer :: idiag_TTmxz=0    ! YAVG_DOC: $\left<T\right>_{y}$
!
! z averaged diagnostics given in zaver.in
!
  integer :: idiag_TTmxy=0    ! ZAVG_DOC: $\left<T\right>_{z}$
!
! variables for slices given in video.in
!
  integer :: ivid_pp=0
  real, dimension(:,:), allocatable :: pp_xz,pp_yz,pp_xy,pp_xy2,pp_xy3,pp_xy4,pp_xz2
  real, dimension(:,:,:,:,:), allocatable :: pp_r
!
! General variables for operator split terms.
!
  real :: cv1 = 0.0, cv1_temp = 0.0
!
! Coefficients for the KI02 terms
!
  real, parameter :: KI_heat = 2.0e-26                             ! heating rate [erg/s]
  real, parameter :: KI_v1 = 1.e7, KI_v2 = 1.4e-2                  ! coefficients of the two terms in Lambda / Gamma [cm^3]
  real, parameter :: KI_T1 = 1.184e5, KI_T2 = 1.e3, KI_T3 = 92.    ! coefficients for the temperatures [K]
  real :: KI_a0 = 0., KI_a1 = 0., KI_a2 = 0.
!
! Sutherland & Dopita (1993) cooling function.
!
  real, dimension(:), allocatable :: SD_logTT, SD_logLambda
  integer :: SD_nt = 0
  real :: SD_a0 = 0.
!
! Coefficient for Jeans energy floor
!
  real :: Jeans_c0 = 0.
!
! Background stratification
!
  real, dimension(mz) :: eth0z = 0.0, dlneth0dz = 0.0
  real, dimension(mz) :: rho0z = 0.0
  real, dimension(nx) :: diffus_chi, diffus_chi3
!
  real :: gamma, gamma_m1

  contains
!***********************************************************************
    subroutine register_energy
!
!  Initialise variables which should know that we solve an energy equation.
!
!  04-nov-10/anders+evghenii: adapted
!
      use FArrayManager, only: farray_register_pde
      use SharedVariables, only: put_shared_variable
!
      call farray_register_pde('eth',ieth)
!
!  Identify version number.
!
      if (lroot) call svn_id( &
           "$Id$")
!
      call put_shared_variable('lviscosity_heat',lviscosity_heat, caller='register_energy')
!
    endsubroutine register_energy
!***********************************************************************
    subroutine initialize_energy(f)
!
!  Called by run.f90 after reading parameters, but before the time loop.
!
!  04-nov-10/anders+evghenii: adapted
!  03-oct-11/ccyang: add initialization for KI02
!
      use EquationOfState, only: select_eos_variable, get_stratz, getmu, cs0, cs20, get_gamma_etc
      use SharedVariables, only: get_shared_variable
      use Slices_methods, only: alloc_slice_buffers
!
      real, dimension (mx,my,mz,mfarray), intent (inout) :: f
!
      real, pointer :: tsg
!
      integer :: istat
      integer :: i, j, k
      real :: mu
      real :: c0, c1, cp, cv
!
      call select_eos_variable('eth',ieth)
!
!  logical variable lpressuregradient_gas shared with hydro modules
!
      call get_shared_variable('lpressuregradient_gas',lpressuregradient_gas,caller='initialize_energy')

      if (.not.leos_idealgas) call fatal_error('initialize_energy','currently assumes EOS=eos_idealgas')
!
      call get_gamma_etc(gamma,cp,cv)
      gamma_m1=gamma-1.
      cv1=1./cv
!
!  Decide if operator splitting is required.
!
      lsplit_update = lconst_cooling_time .or. lKI02 .or. lSD93
      if (lsplit_update .and. .not. ldensity) &
        call fatal_error('initialize_energy','density is required for split_update_energy')
!
!  General variables required by split_update_energy.
!
      if (lsplit_update) then
        cv1_temp = cv1 * unit_temperature
        if (lstratz) cv1_temp = cv1_temp * cs20 / (gamma * gamma_m1)
      endif
!
!  Initialize the KI02 terms.
!
      KI02: if (lKI02) then
        call getmu(mu_tmp=mu)
        c1 = log10(mu * m_u)
        KI_a0 = log10(KI_heat) + log10(unit_time) - log10(unit_energy)
        c0 = KI_a0 - 3. * log10(unit_length) - 2. * c1
        KI_a0 = KI_a0 - c1
        SI: if (unit_system == 'SI') then
          KI_a0 = KI_a0 - 7.
          c0 = c0 - 13.
        endif SI
        KI_a0 = 10.**KI_a0
        c0 = 10.**c0
        disk: if (nzgrid == 1) then
          c0 = c0 * nu_z / (2. * sqrtpi)
          if (cs_z > 0.) c0 = c0 / cs_z
        endif disk
        KI_a1 = c0 * KI_v1
        KI_a2 = c0 * KI_v2
      endif KI02
!
!  Initialize the SD93 terms.
!
      if (lSD93) call init_cooling_SD93
!
!  Initialize the coefficient of the Jeans energy floor.
!
      Jeans: if (ljeans_floor) then
        if (nzgrid /= 1) then
          call not_implemented('initialize_energy', '3D Jeans floor')
        else
          Jeans_c0 = real(njeans) * G_Newton * dxmax / (gamma * gamma_m1)
        endif
      endif Jeans
!
!  Get background stratification, if any.
!
      if (lstratz) call get_stratz(z, rho0z, dlneth0dz, eth0z)  ! dlnrho0dz = dlneth0dz
!
      if (llocal_iso) call fatal_error('initialize_energy', &
          'llocal_iso switches on the local isothermal approximation. ' // &
          'Use ENERGY=noenergy in src/Makefile.local')
!
      if (ivid_pp/=0) call alloc_slice_buffers(pp_xy,pp_xz,pp_yz,pp_xy2,pp_xy3,pp_xy4,pp_xz2,pp_r)

      call keep_compiler_quiet(f)
!
    endsubroutine initialize_energy
!***********************************************************************
    subroutine init_energy(f)
!
!  Initialise thermal energy.
!
!  04-nov-10/anders+evghenii: adapted
!
      use General, only: itoa
      use Initcond, only: jump
      use EquationOfState, only: rho0, cs20
      use InitialCondition, only: initial_condition_ss
!
      real, dimension (mx,my,mz,mfarray), intent (inout) :: f
!
      integer :: j
      logical :: lnothing=.true.
      character (len=intlen) :: iinit_str
!
      do j=1,ninit
!
        if (initeth(j)/='nothing') then
!
          lnothing=.false.
!
          iinit_str=itoa(j)
!
!  Select between various initial conditions.
!
          select case (initeth(j))
!
          case ('zero', '0'); f(:,:,:,ieth) = 0.0
!
          case ('xjump'); call jump(f,ieth,eth_left,eth_right,widtheth,xjump_mid,yjump_mid,zjump_mid,'x')
!
          case ('const_eth'); f(:,:,:,ieth)=f(:,:,:,ieth)+eth_const
!
          case ('basic_state')
            if (lstratz) then
              f(:,:,:,ieth) = 0.0
            else
              f(:,:,:,ieth) = rho0 * cs20 / (gamma * gamma_m1)
            endif
!
          case default
!
!  Catch unknown values.
!
            call fatal_error('init_energy','no such initeth: '//trim(initeth(j)))
!
          endselect
        endif
      enddo
!
!  Call for the InitialCondition facility.
!
      if (linitial_condition) call initial_condition_ss(f)
!
    endsubroutine init_energy
!***********************************************************************
    subroutine pencil_criteria_energy
!
!  All pencils that the Energy module depends on are specified here.
!
!  04-nov-10/anders+evghenii: adapted
!
      lpenc_requested(i_divu) = .true.
      lpenc_requested(i_pp) = .true.
      stratz: if (lstratz) then
        lpenc_requested(i_eths) = .true.
        lpenc_requested(i_ugeths) = .true.
        lpenc_requested(i_uu) = .true.
      else stratz
        weno: if (lweno_transport) then
          lpenc_requested(i_transpeth) = .true.
        else weno
          lpenc_requested(i_geth) = .true.
          lpenc_requested(i_eth) = .true.
          lpenc_requested(i_uu) = .true.
        endif weno
        lpenc_requested(i_pp) = .true.
      endif stratz
      if (lhydro) lpenc_requested(i_fpres)=.true.
      if (ldt) lpenc_requested(i_cs2)=.true.
      if (lviscosity.and.lviscosity_heat) lpenc_requested(i_visc_heat)=.true.
      if (chi/=0.0) then
        lpenc_requested(i_rho)=.true.
        lpenc_requested(i_cp)=.true.
        lpenc_requested(i_del2TT)=.true.
        lpenc_requested(i_grho)=.true.
        lpenc_requested(i_gTT)=.true.
      endif
      if (lchi_rosseland) then
        lpenc_requested(i_rho)=.true.
        lpenc_requested(i_rho1)=.true.
        lpenc_requested(i_del2TT)=.true.
        lpenc_requested(i_grho)=.true.
        lpenc_requested(i_gTT)=.true.
      endif
!
      shock: if (chi_shock /= 0.0) then
        lpenc_requested(i_shock) = .true.
        lpenc_requested(i_gshock) = .true.
        stratz1: if (lstratz) then
          lpenc_requested(i_geths) = .true.
        else stratz1
          lpenc_requested(i_geth) = .true.
          lpenc_requested(i_del2eth) = .true.
        endif stratz1
      endif shock
!
!  Diagnostic pencils.
!
      if (idiag_ethm/=0 .or. idiag_ethmin/=0 .or. idiag_ethmax/=0 .or. idiag_ethtot/=0) &
         lpenc_diagnos(i_eth)=.true.
      if (idiag_eem/=0) lpenc_diagnos(i_ee)=.true.
      if (idiag_etot /= 0) then
        lpenc_diagnos(i_eth) = .true.
        lpenc_diagnos(i_ekin) = .true.
      endif
      if (idiag_ppm/=0) lpenc_diagnos(i_pp)=.true.
      if (idiag_pdivum/=0) then
        lpenc_diagnos(i_pp)=.true.
        lpenc_diagnos(i_divu)=.true.
      endif
!
      if (idiag_TTm/=0   .or. idiag_TTmin/=0 .or. idiag_TTmax/=0 .or. &
          idiag_TTmxy/=0 .or. idiag_TTmxz/=0 .or. idiag_TTmx/=0  .or. &
          idiag_TTmy/=0  .or. idiag_TTmz/=0 ) lpenc_diagnos(i_TT)=.true.
!
    endsubroutine pencil_criteria_energy
!***********************************************************************
    subroutine pencil_interdep_energy(lpencil_in)
!
!  Interdependency among pencils from the Energy module is specified here.
!
!  04-nov-10/anders+evghenii: adapted
!
      logical, dimension(npencils) :: lpencil_in
!
      fpres: if (lpencil_in(i_fpres)) then
        lpencil_in(i_rho1) = .true.
        stratz: if (lstratz) then
          lpencil_in(i_geths) = .true.
        else stratz
          lpencil_in(i_geth) = .true.
        endif stratz
      endif fpres
!
      ugeths: if (lpencil_in(i_ugeths)) then
        lpencil_in(i_geths) = .true.
        lpencil_in(i_uu) = .true.
      endif ugeths
!
    endsubroutine pencil_interdep_energy
!***********************************************************************
    subroutine calc_pencils_energy(f,p)
!
!  Calculate Energy pencils.
!  This routine is called after  calc_pencils_eos
!  Most basic pencils should come first, as others may depend on them.
!
!  04-nov-10/anders+evghenii: adapted
!  14-feb-11/bing: moved eth dependend pecnils to eos_idealgas
!
      use EquationOfState, only: cs20
      use Sub, only: u_dot_grad
      use WENO_transport, only: weno_transp
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (pencil_case) :: p
!
      real, dimension(nx) :: penc
      integer :: j
!
! fpres
!
      fpres: if (lpencil(i_fpres)) then
        stratz: if (lstratz) then
          penc = gamma_m1 * eth0z(n) * p%rho1
          p%fpres = -p%geths
          p%fpres(:,3) = p%fpres(:,3) + (f(l1:l2,m,n,irho) - f(l1:l2,m,n,ieth)) * dlneth0dz(n)
          forall(j = 1:3) p%fpres(:,j) = penc * p%fpres(:,j)
        else stratz
          forall(j = 1:3) p%fpres(:,j) = -gamma_m1 * p%rho1 * p%geth(:,j)
        endif stratz
      endif fpres
!
! ugeths
!
      if (lpencil(i_ugeths)) call u_dot_grad(f, ieth, p%geths, p%uu, p%ugeths)
!
! transpeth
!
      if (lpencil(i_transpeth)) call weno_transp(f,m,n,ieth,-1,iux,iuy,iuz,p%transpeth,dx_1,dy_1,dz_1)
! sglnTT 
      if (lpencil(i_sglnTT)) &
        call not_implemented('calc_pencils_energy','pencil sglnTT for thermal_energy')
!
    endsubroutine calc_pencils_energy
!***********************************************************************
    subroutine denergy_dt(f,df,p)
!
!  Calculate right hand side of energy equation.
!
!    deth/dt + div(u*eth) = -P*div(u)
!
!  04-nov-10/anders+evghenii: coded
!  02-aug-11/ccyang: add mesh hyper-diffusion
!
      use Special, only: special_calc_energy
      use Sub, only: identify_bcs, u_dot_grad, del2
      use Viscosity, only: calc_viscous_heat
      use Deriv, only: der6
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
      real, dimension(mx,my,mz,mvar), intent(inout) :: df
      type(pencil_case), intent(in) :: p
!
      real, dimension(nx) :: Hmax=0.0, ugeth, d6eth, penc
      integer :: j
!
!  Identify module and boundary conditions.
!
      if (headtt.or.ldebug) print*, 'denergy_dt: solve deth_dt'
      if (headtt) call identify_bcs('eth',ieth)
!
!  Sound speed squared.
!
      if (headtt) print*, 'denergy_dt: cs20=', p%cs2(1)
!
!  ``cs2/dx^2'' for timestep
!
      if (ldensity.and.lhydro.and.lfirst.and.ldt) then
        advec_cs2=p%cs2*dxyz_2
        if (headtt.or.ldebug) print*,'denergy_dt: max(advec_cs2) =',maxval(advec_cs2)
      endif
!
!  Add pressure gradient term in momentum equation.
!
      if (lhydro) df(l1:l2,m,n,iux:iuz) = df(l1:l2,m,n,iux:iuz) + p%fpres
!
!  Entry possibility for "personal" entries.
!  In that case you'd need to provide your own "special" routine.
!
      if (lspecial) call special_calc_energy(f,df,p)
!
      if (lstratz) then

        df(l1:l2,m,n,ieth) = df(l1:l2,m,n,ieth) - p%ugeths - p%eths*(gamma*p%divu + p%uu(:,3)*dlneth0dz(n))

      else
!
!  Add energy transport term.
!
        if (lweno_transport) then
          df(l1:l2,m,n,ieth) = df(l1:l2,m,n,ieth) - p%transpeth
        else
          call u_dot_grad(f, ieth, p%geth, p%uu, ugeth, UPWIND=lupw_eth)
          df(l1:l2,m,n,ieth) = df(l1:l2,m,n,ieth) - p%eth*p%divu - ugeth
        endif
!
!  Add P*dV work.
!
        df(l1:l2,m,n,ieth) = df(l1:l2,m,n,ieth) - p%pp*p%divu
!
      endif
!
!  Calculate viscous contribution to temperature.
!
      if (lviscosity.and.lviscosity_heat) call calc_viscous_heat(df,p,Hmax)
!
!  Thermal energy diffusion.
!
      diffus_chi=0; diffus_chi3=0.
      if (chi/=0.0) then
        df(l1:l2,m,n,ieth) = df(l1:l2,m,n,ieth) + chi * p%cp * (p%rho*p%del2TT + sum(p%grho*p%gTT, 2))
        if (lfirst .and. ldt) diffus_chi = diffus_chi + gamma*chi*dxyz_2
      endif
!
!  Shock diffusion
!
      shock: if (chi_shock /= 0.0) then
        stratz1: if (lstratz) then
          call del2(f, ieth, penc)
          df(l1:l2,m,n,ieth) = df(l1:l2,m,n,ieth) + chi_shock*(p%shock*penc + sum(p%gshock*p%geths, 2))
        else stratz1
          df(l1:l2,m,n,ieth) = df(l1:l2,m,n,ieth) + chi_shock*(p%shock*p%del2eth + sum(p%gshock*p%geth,2))
        endif stratz1
        if (lfirst .and. ldt) diffus_chi = diffus_chi + chi_shock * p%shock * dxyz_2
      endif shock
!
!  Mesh hyper-diffusion
!
      if (chi_hyper3_mesh/=0.0) then
        do j=1,3
          call der6(f, ieth, d6eth, j, IGNOREDX=.true.)
          df(l1:l2,m,n,ieth) = df(l1:l2,m,n,ieth) + chi_hyper3_mesh * d6eth * dline_1(:,j)
        enddo
        if (lfirst .and. ldt) diffus_chi3 = diffus_chi3 + chi_hyper3_mesh*sum(dline_1,2)                  
      endif
!
!  Radiative diffusion (Rosseland approximation) through thermal energy diffusion.
!
      if (lchi_rosseland) then
        df(l1:l2,m,n,ieth) = df(l1:l2,m,n,ieth) + ( 16.0*sigmaSB/(3.0*kappa_rosseland))*p%TT*p%TT*p%rho1*( &
                             3.0*sum(p%gTT*p%gTT, 2) - p%TT*p%rho1*sum(p%grho*p%gTT, 2) + p%TT*p%del2TT )
!       This expression has an extra TT/eth = cv1*rho1  to convert to the right units - needed because of the sigmaSB
!       The timestep gets the factor TT/eth = cv1*rho1
        if (lfirst .and. ldt) diffus_chi = diffus_chi &
            +(16.0*sigmaSB/(3.0*kappa_rosseland))*p%TT*p%TT*p%TT/p%rho *cv1*p%rho1 * dxyz_2
      endif
      if (lfirst .and. ldt) then
        maxdiffus=max(maxdiffus,diffus_chi)
        maxdiffus3=max(maxdiffus3,diffus_chi3)
      endif
!
!  Diagnostics.
!
      call calc_diagnostics_energy(f,p)
!
    endsubroutine denergy_dt
!***********************************************************************
    subroutine calc_diagnostics_energy(f,p)

      use Diagnostics
      use Slices_methods, only: store_slices

      real, dimension (mx,my,mz,mfarray) :: f
      type(pencil_case) :: p

      call keep_compiler_quiet(f)

      if (ldiagnos) then
        call sum_mn_name(p%TT,idiag_TTm)
        call max_mn_name(p%TT,idiag_TTmax)
        if (idiag_TTmin/=0)  call max_mn_name(-p%TT,idiag_TTmin,lneg=.true.)
        call sum_mn_name(p%eth,idiag_ethm)
        call integrate_mn_name(p%eth,idiag_ethtot)
        if (idiag_ethmin/=0) call max_mn_name(-p%eth,idiag_ethmin,lneg=.true.)
        call max_mn_name(p%eth,idiag_ethmax)
        call sum_mn_name(p%ee,idiag_eem)
        call sum_mn_name(p%pp,idiag_ppm)
        if (idiag_etot /= 0) call sum_mn_name(p%eth + p%ekin, idiag_etot)
        if (idiag_pdivum/=0) call sum_mn_name(p%pp*p%divu,idiag_pdivum)
      endif
!
      if (l1davgfirst) then
        call yzsum_mn_name_x(p%pp,idiag_ppmx)
        call xzsum_mn_name_y(p%pp,idiag_ppmy)
        call xysum_mn_name_z(p%pp,idiag_ppmz)
        call yzsum_mn_name_x(p%TT,idiag_TTmx)
        call xzsum_mn_name_y(p%TT,idiag_TTmy)
        call xysum_mn_name_z(p%TT,idiag_TTmz)
      endif
!
!  2-D averages.
!
      if (l2davgfirst) then
        call zsum_mn_name_xy(p%TT,idiag_TTmxy)
        call ysum_mn_name_xz(p%TT,idiag_TTmxz)
      endif
!
      if (lvideo.and.lfirst) then
        if (ivid_pp/=0) call store_slices(p%pp,pp_xy,pp_xz,pp_yz,pp_xy2,pp_xy3,pp_xy4,pp_xz2,pp_r)
      endif

    endsubroutine calc_diagnostics_energy
!***********************************************************************
    subroutine energy_before_boundary(f)
!
!  Actions to take before boundary conditions are set.
!
!   1-apr-20/joern: coded
!
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f
      real, dimension (mx) :: cs2
!
!    Slope limited diffusion: update characteristic speed
!    Not staggered yet
!
     if (lslope_limit_diff .and. llast) then
       cs2=0.
       do m=1,my
       do n=1,mz
         if (ldensity_nolog) then
           cs2 = gamma*gamma_m1*f(:,m,n,ieth)/f(:,m,n,irho)
         else
           cs2 = gamma*gamma_m1*f(:,m,n,ieth)*exp(-f(:,m,n,ilnrho))
         endif
         f(:,m,n,isld_char)=f(:,m,n,isld_char)+w_sldchar_ene*sqrt(cs2)
       enddo
       enddo
     endif
!
    endsubroutine energy_before_boundary
!***********************************************************************
    subroutine read_energy_init_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=entropy_init_pars, IOSTAT=iostat)
!
    endsubroutine read_energy_init_pars
!***********************************************************************
    subroutine write_energy_init_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=entropy_init_pars)
!
    endsubroutine write_energy_init_pars
!***********************************************************************
    subroutine read_energy_run_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=entropy_run_pars, IOSTAT=iostat)
!
    endsubroutine read_energy_run_pars
!***********************************************************************
    subroutine write_energy_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=entropy_run_pars)
!
    endsubroutine write_energy_run_pars
!***********************************************************************
    subroutine rprint_energy(lreset,lwrite)
!
!  Reads and registers print parameters relevant to entropy.
!
!  04-nov-10/anders+evghenii: adapted from temperature_idealgas.f90
!
      use Diagnostics, only: parse_name
!
      logical :: lreset
      logical, optional :: lwrite
!
      integer :: iname, inamex, inamey, inamez, inamexy, inamexz
      logical :: lwr
!
      lwr = .false.
      if (present(lwrite)) lwr=lwrite
!
!  Reset everything in case of reset.
!  (this needs to be consistent with what is defined above!)
!
      if (lreset) then
        idiag_TTm=0; idiag_TTmax=0; idiag_TTmin=0
        idiag_ethm=0; idiag_ethmin=0; idiag_ethmax=0; idiag_eem=0
        idiag_etot = 0; idiag_ethtot=0
        idiag_pdivum=0; idiag_ppm=0
        ivid_pp=0
      endif
!
!  iname runs through all possible names that may be listed in print.in
!
      do iname=1,nname
        call parse_name(iname,cname(iname),cform(iname),'TTm',idiag_TTm)
        call parse_name(iname,cname(iname),cform(iname),'TTmax',idiag_TTmax)
        call parse_name(iname,cname(iname),cform(iname),'TTmin',idiag_TTmin)
        call parse_name(iname,cname(iname),cform(iname),'ethm',idiag_ethm)
        call parse_name(iname,cname(iname),cform(iname),'ethmin',idiag_ethmin)
        call parse_name(iname,cname(iname),cform(iname),'ethmax',idiag_ethmax)
        call parse_name(iname,cname(iname),cform(iname),'eem',idiag_eem)
        call parse_name(iname,cname(iname),cform(iname),'etot',idiag_etot)
        call parse_name(iname,cname(iname),cform(iname),'ppm',idiag_ppm)
        call parse_name(iname,cname(iname),cform(iname),'pdivum',idiag_pdivum)
        call parse_name(iname,cname(iname),cform(iname),'ethtot',idiag_ethtot)
      enddo
!
!  Check for those quantities for which we want yz-averages.
!
      do inamex=1,nnamex
        call parse_name(inamex,cnamex(inamex),cformx(inamex),'ppmx',idiag_ppmx)
        call parse_name(inamex,cnamex(inamex),cformx(inamex),'TTmx',idiag_TTmx)
      enddo
!
!  Check for those quantities for which we want xz-averages.
!
      do inamey=1,nnamey
        call parse_name(inamey,cnamey(inamey),cformy(inamey),'ppmy',idiag_ppmy)
        call parse_name(inamey,cnamey(inamey),cformy(inamey),'TTmy',idiag_TTmy)
      enddo
!
!  Check for those quantities for which we want xy-averages.
!
      do inamez=1,nnamez
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'ppmz',idiag_ppmz)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'TTmz',idiag_TTmz)
      enddo
!
!  Check for those quantities for which we want z-averages.
!
      do inamexy=1,nnamexy
        call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy),'TTmxy',idiag_TTmxy)
      enddo
!
!  Check for those quantities for which we want y-averages.
!
      do inamexz=1,nnamexz
        call parse_name(inamexz,cnamexz(inamexz),cformxz(inamexz),'TTmxz',idiag_TTmxz)
      enddo
!
!  check for those quantities for which we want video slices
!
      if (lwrite_slices) then 
        where(cnamev=='TT'.or.cnamev=='lnTT') cformv='DEFINED'
      endif
      do iname=1,nnamev
        call parse_name(iname,cnamev(iname),cformv(iname),'pp',ivid_pp)
      enddo
!
    endsubroutine rprint_energy
!***********************************************************************
    subroutine get_slices_energy(f,slices)
!
!  04-nov-10/anders+evghenii: adapted
!
      use Slices_methods, only: assign_slices_scal, process_slices, exp2d
      use EquationOfState, only: eoscalc, irho_eth, ilnrho_eth
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (slice_data) :: slices
!
      real, dimension(nx) :: penc
      integer :: ieosvars
      integer :: m, n
!
!  Loop over slices
!
      slice_name: select case (trim(slices%name))
!
!  Pressure
!
        case ('pp') slice_name
          call assign_slices_scal(slices,pp_xy,pp_xz,pp_yz,pp_xy2,pp_xy3,pp_xy4,pp_xz2,pp_r)
!
!  Temperature
!
        case ('TT', 'lnTT') slice_name
          density: if (ldensity_nolog) then
            ieosvars = irho_eth
          else density
            ieosvars = ilnrho_eth
          endif density
!
          do m = m1, m2
            if (lwrite_slice_xy) &
              call eoscalc(ieosvars, f(l1:l2,m,iz_loc,ilnrho), f(l1:l2,m,iz_loc,ieth), iz=iz_loc, lnTT=slices%xy(:,m-nghost))
            if (lwrite_slice_xy2) &
              call eoscalc(ieosvars, f(l1:l2,m,iz2_loc,ilnrho), f(l1:l2,m,iz2_loc,ieth), iz=iz2_loc, lnTT=slices%xy(:,m-nghost))
            if (lwrite_slice_xy3) &
              call eoscalc(ieosvars, f(l1:l2,m,iz3_loc,ilnrho), f(l1:l2,m,iz3_loc,ieth), iz=iz3_loc, lnTT=slices%xy(:,m-nghost))
            if (lwrite_slice_xy4) &
              call eoscalc(ieosvars, f(l1:l2,m,iz4_loc,ilnrho), f(l1:l2,m,iz4_loc,ieth), iz=iz4_loc, lnTT=slices%xy(:,m-nghost))
          enddo
!
          if (lwrite_slice_xz.or.lwrite_slice_xz2) then
            do n = n1, n2
              if (lwrite_slice_xz) &
                call eoscalc(ieosvars, f(l1:l2,iy_loc,n,ilnrho), f(l1:l2,iy_loc,n,ieth), iz=n, lnTT=slices%xz(:,n-nghost))
              if (lwrite_slice_xz2) &
                call eoscalc(ieosvars, f(l1:l2,iy2_loc,n,ilnrho), f(l1:l2,iy2_loc,n,ieth), iz=n, lnTT=slices%xz2(:,n-nghost))
            enddo
          endif
!
          if (lwrite_slice_yz) then
            do n = n1, n2
              do m = m1, m2
                call eoscalc(ieosvars, f(ix_loc,m,n,ilnrho), f(ix_loc,m,n,ieth), iz=n, lnTT=slices%yz(m-nghost,n-nghost))
              enddo
            enddo
          endif
!
          if (trim(slices%name) == 'TT') call process_slices(slices,exp2d)
!
          slices%ready = .true.
!
      endselect slice_name
!
    endsubroutine get_slices_energy
!***********************************************************************
    subroutine impose_energy_floor(f)
!
!  Trap any negative energy or impose a floor in minimum thermal energy.
!
!  29-aug-11/ccyang: coded
!
      real, dimension(mx,my,mz,mfarray) :: f
!
      real, dimension(nx) :: eth, rho, eth1
      integer :: i, j, k
      real :: deth
!
!  Impose the energy floor.
!
      efloor: if (energy_floor > 0.0) then
        stratz: if (lstratz) then
          zscan: do k = 1, mz
            deth = energy_floor / eth0z(k) - 1.0
            where(f(:,:,k,ieth) < deth) f(:,:,k,ieth) = deth
          enddo zscan
        else stratz
          where(f(:,:,:,ieth) < energy_floor) f(:,:,:,ieth) = energy_floor
        endif stratz
      endif efloor
!
!  Apply Jeans energy floor.
!
      jfloor: if (ljeans_floor) then
        zscan1: do k = n1, n2
          yscan1: do j = m1, m2
            stratz1: if (lstratz) then
              eth = eth0z(k) * (1.0 + f(l1:l2,j,k,ieth))
              rho = rho0z(k) * (1.0 + f(l1:l2,j,k,irho))
              eth1 = eth
              call jeans_floor(eth, rho)
              where(eth1 /= eth) f(l1:l2,j,k,ieth) = eth / eth0z(k) - 1.0
            else stratz1
              eth = f(l1:l2,j,k,ieth)
              rho = f(l1:l2,j,k,irho)
              if (.not. ldensity_nolog) rho = exp(rho)
              call jeans_floor(eth, rho)
              f(l1:l2,j,k,ieth) = eth
            endif stratz1
          enddo yscan1
        enddo zscan1
      endif jfloor
!
!  Stop the code if negative energy exists.
!
      chkneg: if (lcheck_negative_energy) then
        negeth: if (any(f(l1:l2,m1:m2,n1:n2,ieth) <= merge(-1.0, 0.0, lstratz))) then
          zscan2: do k = n1, n2
            yscan2: do j = m1, m2
              xscan2: do i = l1, l2
                if (lstratz) then
                  if (f(i,j,k,ieth) <= -1.0) print 10, f(i,j,k,ieth), x(i), y(j), z(k)
                else
                  if (f(i,j,k,ieth) <= 0.0) print 20, f(i,j,k,ieth), x(i), y(j), z(k)
                endif
                10 format (1x, 'deth = ',es13.6, ' at x = ',es13.6, ', y = ',es13.6, ', z = ',es13.6)
                20 format (1x, 'eth = ', es13.6, ' at x = ', es13.6, ', y = ', es13.6, ', z = ', es13.6)
              enddo xscan2
            enddo yscan2
          enddo zscan2
          call fatal_error('impose_energy_floor', 'negative energy detected')
        endif negeth
      endif chkneg
!
    endsubroutine impose_energy_floor
!***********************************************************************
    subroutine dynamical_thermal_diffusion(uc)
!
!  Dynamically set thermal diffusion coefficient given fixed mesh Reynolds number.
!
!  02-aug-11/ccyang: coded
!
!  Input Argument
!      uc
!          Characteristic velocity of the system.
!
      real, intent(in) :: uc
!
!  Hyper-diffusion coefficient
!
      if (chi_hyper3_mesh /= 0.0) chi_hyper3_mesh = pi5_1 * uc / re_mesh / sqrt(real(dimensionality))
!
    endsubroutine dynamical_thermal_diffusion
!***********************************************************************
    subroutine split_update_energy(f)
!
!  Update the thermal energy by integrating the operator split energy
!  terms.
!
!  19-jan-13/ccyang: coded
!
      use Boundcond, only: zero_ghosts, update_ghosts
      use Density, only: impose_density_floor
!
      real, dimension(mx,my,mz,mfarray), intent(inout) :: f
!
      integer, dimension(mx,my,mz) :: status
      real, dimension(mx,my,mz) :: delta_eth
      character(len=256) :: message
      integer :: i, j, k
!
      split_update: if (lsplit_update) then
!
!  Impose density and energy floors.
!
        call impose_density_floor(f)
        call impose_energy_floor(f)
!
!  Update the ghost cells.
!
        call zero_ghosts(f, ieth)
        call update_ghosts(f, ieth)
        gcrho: if (ldensity_nolog) then
          call zero_ghosts(f, irho)
          call update_ghosts(f, irho)
        else gcrho
          call zero_ghosts(f, ilnrho)
          call update_ghosts(f, ilnrho)
        endif gcrho
!
!  Update the energy.
!
        deth: if (lstratz) then
          zscan: do k = 1, mz
            call get_delta_eth(real(t), dt, eth0z(k) * (1.0 + f(:,:,k,ieth)), rho0z(k) * (1.0 + f(:,:,k,irho)), &
                               delta_eth(:,:,k), status(:,:,k))
            delta_eth(:,:,k) = delta_eth(:,:,k) / eth0z(k)
          enddo zscan
        elseif (ldensity_nolog) then deth
          call get_delta_eth(real(t), dt, f(:,:,:,ieth), f(:,:,:,irho), delta_eth, status)
        else deth
          call get_delta_eth(real(t), dt, f(:,:,:,ieth), exp(f(:,:,:,ilnrho)), delta_eth, status)
        endif deth
!
        if (any(status < 0)) then
          write(message,10) count(status == -1), ' cells had underflows and ', &
                            count(status == -2), ' cells reached maximum iterations'
          10 format(1x, i3, a, i3, a)
          call warning('split_update_energy', trim(message))
          do k = n1, n2; do j = m1, m2; do i = l1, l2
            if (status(i,j,k) == -2) print 20, f(i,j,k,irho), get_temperature(f(i,j,k,ieth), f(i,j,k,irho))
            20 format (1x, 'rho = ', es13.6, ' & TT = ', es13.6)
          enddo; enddo; enddo
        endif
!
        if (ldebug) then
          where(status < 0) status = rk_nmax
          print *, 'Minimum, maximum, and average numbers of iterations = ', &
                   minval(status), maxval(status), real(sum(status)) / real(nw)
        endif
!
        f(l1:l2,m1:m2,n1:n2,ieth) = f(l1:l2,m1:m2,n1:n2,ieth) + delta_eth(l1:l2,m1:m2,n1:n2)
      endif split_update
!
    endsubroutine
!***********************************************************************
    pure subroutine const_cooling_time(eth, rho, dt)
!
!  Exactly integrate the cooling term dTT/dt = (TTref - TT) / tau_cool.
!  Ideal equation of state is assumed.
!
!  31-jan-13/ccyang: coded.
!
      real, intent(inout) :: eth
      real, intent(in) :: rho, dt
!
      real :: temp, a, b, c, x
!
      c = cv1 / rho
      temp = c * eth
      x = dt / tau_cool
      a = exp(-x)
      b = 1.0 - a
      if (b == 0.0) b = x * (1.0 - 0.5 * x)
      temp = temp * a + TTref * b
      eth = temp / c
!
    endsubroutine const_cooling_time
!***********************************************************************
    elemental subroutine get_delta_eth(t, dt, eth, rho, delta_eth, status)
!
!  Integrate the operator split energy terms that are local:
!
!    de/dt = H(e(t),rho)
!
!  08-aug-11/ccyang: coded
!
      real, intent(in) :: t, dt, eth, rho
      real, intent(out) :: delta_eth
      integer, intent(out), optional :: status
!
!  Cash-Karp parameters for embedded Runga-Kutta method
!
      real, parameter :: safety = 0.9, errcon = (5. / safety)**5
      real, parameter :: a2 = .2, a3 = .3, a4 = .6, a5 = 1., a6 = .875
      real, parameter :: b2 = .2
      real, dimension(2), parameter :: b3 = (/ .075, .225 /)
      real, dimension(3), parameter :: b4 = (/ .3, -.9, 1.2 /)
      real, dimension(4), parameter :: b5 = (/ -11./54., 2.5, -70./27., 35./27. /)
      real, dimension(5), parameter :: b6 = (/ 1631./55296., 175./512., 575./13824., 44275./110592., 253./4096. /)
      real, dimension(6), parameter :: c = (/ 37./378., 0., 250./621., 125./594., 0., 512./1771. /)
      real, dimension(6), parameter :: d = c - (/ 2825./27648., 0., 18575./48384., 13525./55296., 277./14336., .25 /)
!
      real, dimension(6) :: de
      logical :: last, lovershoot
      integer :: nok, nbad, status1
      real :: tf, t1, dt1, eth1, ethj, deth, error
!
!  Initialization
!
      delta_eth = 0.
!
      tf = abs(t + dt)
      t1 = t
      dt1 = dt
      eth1 = eth
      nok = 0
      nbad = 0
      last = .true.
!
      rkck: do
!
!  Embedded Runga-Kutta step
!
        de(1) = dt1 * calc_heat_split(t1, eth1, rho)
        de(2) = dt1 * calc_heat_split(t1 + a2 * dt1, eth1 + b2 * de(1), rho)
        de(3) = dt1 * calc_heat_split(t1 + a3 * dt1, eth1 + sum(b3 * de(:2)), rho)
        de(4) = dt1 * calc_heat_split(t1 + a4 * dt1, eth1 + sum(b4 * de(:3)), rho)
        de(5) = dt1 * calc_heat_split(t1 + a5 * dt1, eth1 + sum(b5 * de(:4)), rho)
        de(6) = dt1 * calc_heat_split(t1 + a6 * dt1, eth1 + sum(b6 * de(:5)), rho)
!
        deth = sum(c * de)
!
!  Check overshoot towards negative energy.
!
        lovershoot = .false.
        floor: if (ljeans_floor) then
          ethj = eth1 + deth
          call jeans_floor(ethj, rho)
          if (ethj > eth1 + deth) then
            lovershoot = .true.
            deth = ethj - eth1
            error = 0.03125
          endif
        elseif (eth1 + deth <= 0.) then floor
          lovershoot = .true.
          error = 32.
        endif floor
!
!  Estimate the error.
!
        if (.not. lovershoot) error = abs(sum(d * de)) / (rk_eps * (eth1 + abs(de(1))))
!
!  Time step control
!
        passed: if (error <= 1.) then
          nok = nok + 1
          delta_eth = delta_eth + deth
!
          done: if (last) then
            if (lovershoot) then
              status1 = 1
            else
              status1 = 0
            endif
            exit rkck
          else if (t1 + dt1 == t1) then done
            status1 = -1
            exit rkck
          endif done
!
          t1 = t1 + dt1
          eth1 = eth + delta_eth
!
          if (error > errcon) then
            dt1 = safety * dt1 * error**(-0.2)
          else
            dt1 = 5. * dt1
          endif
          if (abs(t1 + dt1) >= tf) then
            dt1 = t + dt - t1
            last = .true.
          endif
        else passed
          nbad = nbad + 1
          dt1 = sign(max(abs(safety*dt1*error**(-0.25)), 0.1*abs(dt1)), dt1)
          last = .false.
        endif passed
!
        if (nok + nbad > rk_nmax) then
          status1 = -2
          exit rkck
        endif
!
      enddo rkck
!
      if (present(status)) then
        if (ldebug) then
          status = nok + nbad
        else
          status = status1
        endif
      endif
!
!  Cases for which there exists exact solutions.
!
      const: if (lconst_cooling_time) then
        eth1 = eth + delta_eth
        call const_cooling_time(eth1, rho, dt)
        delta_eth = eth1 - eth
      endif const
!
    endsubroutine get_delta_eth
!***********************************************************************
    elemental function calc_heat_split(t, eth, rho) result(heat)
!
!  Calculate the total heat generated from the operator split energy
!  terms.
!
!  06-aug-11/ccyang: coded
!
      real, intent(in) :: t, eth, rho
      real :: heat,eth1
!
      if (t==t) eth1 = eth     ! to quiet the compiler, can be later removed if t is ever used.
      if (ljeans_floor) call jeans_floor(eth1, rho)
!
      heat = 0.
      if (lKI02 .and. lSD93) then
        heat = heat - cool_KI02_SD93(eth1, rho)
      elseif (lKI02) then
        heat = heat + heat_KI02(eth1, rho)
      elseif (lSD93) then
        heat = heat - cool_SD93(eth1, rho)
      endif
!
    endfunction calc_heat_split
!***********************************************************************
    elemental real function heat_KI02(eth, rho)
!
!  Add Koyama & Inutsuka (2002) heating and cooling terms.
!
!  03-oct-11/ccyang: coded
!
      use EquationOfState, only: get_soundspeed
!
      real, intent(in) :: eth, rho
!
      real :: temp, cs2
      real :: c0, c1, c2
!
!  Find the temperature and its derived values.
!
      heat_KI02 = 0.
      temp = get_temperature(eth, rho)
      if (temp <= 0.) return
      c1 = exp(-KI_T1 / (temp + KI_T2))
      c2 = exp(-KI_T3 / temp)
!
!  Calculate the net heat.
!    For 3D run: rho is volume mass density
!    For 2D run: rho is surface mass density
!    cs_z, if > 0, fix the disk vertical scale height
!
      c0 = rho
      if (nzgrid == 1 .and. cs_z <= 0.) then
        call get_soundspeed(real(temp / unit_temperature), cs2)
        c0 = rho / sqrt(cs2)
      endif
      heat_KI02 = rho * (KI_a0 - c0 * (KI_a1 * c1 + KI_a2 * sqrt(temp) * c2))
!
    endfunction
!***********************************************************************
    subroutine init_cooling_SD93
!
!  Initializes the tabulated cooling function of Sutherland and Dopita (1993).
!
!  02-feb-13/ccyang: coded.
!   8-nov-13/axel: changed some calls to f95 standard
!
      use Mpicomm
      use Syscalls, only: get_env_var
      use EquationOfState, only: getmu
!
      integer, parameter :: lun = 1
      integer, parameter :: ncol = 12
      real, dimension(ncol) :: col
      character(len=256) :: msg, src
      integer :: stat
      real :: mu
      integer :: i
!
!  Find the number of table entries.
!
      if (lroot) then
!-      call get_environment_variable('PENCIL_HOME', src)
        call get_env_var('PENCIL_HOME', src)
        src = trim(src) // '/src/cooling_SD93_Z00.dat'
        msg=''
!-      open (unit=lun, file=src, action='read', iostat=stat, iomsg=msg)
        open (unit=lun, file=src, action='read', iostat=stat)
        if (stat /= 0) call fatal_error('init_cooling_SD93','cannot open the cooling table; '// &
                                        trim(msg),force=.true.)
        nline: do
!-        read (lun,*,iostat=stat,iomsg=msg) col
          read (lun,*,iostat=stat) col
          if (stat < 0) exit nline
          if (stat > 0) call fatal_error('init_cooling_SD93','error in reading the cooling table; '// &
                                         trim(msg),force=.true.)
          SD_nt = SD_nt + 1
        enddo nline
!-      close (unit=lun, iostat=stat, iomsg=msg)
        close (unit=lun, iostat=stat)
        if (stat /= 0) call fatal_error('init_cooling_SD93', 'cannot close the cooling table; '// &
                                        trim(msg),force=.true.)
      endif
      call mpibcast_int(SD_nt,comm=MPI_COMM_WORLD)
!
!  Allocate the cooling table.
!
      allocate (SD_logTT(SD_nt), SD_logLambda(SD_nt), stat=stat)
      if (stat /= 0) call fatal_error('init_cooling_SD93', 'cannot allocate the cooling table')
!
!  Read the cooling table.
!
      table: if (lroot) then
        open (unit=lun, file=src, action='read')
        reading: do i = 1, SD_nt
          read (lun,*) col
          SD_logTT(i) = col(1)
          SD_logLambda(i) = col(5)
        enddo reading
        close (unit=lun)
      endif table
      call mpibcast_real(SD_logTT, SD_nt, comm=MPI_COMM_WORLD)
      call mpibcast_real(SD_logLambda, SD_nt,comm=MPI_COMM_WORLD)
!
      if (lroot) print *, 'init_cooling_SD93: read the Sutherland & Dopita (1993) cooling table.'
!
!  Save the conversion factor.
!
      call getmu(mu_tmp=mu)
      SD_a0 = exp(log(unit_time) - log(unit_energy) - 3. * log(unit_length) - 2. * log(mu * m_u))
      if (unit_system == 'SI') SD_a0 = 1.e-13 * SD_a0
      disk: if (nzgrid == 1) then
        SD_a0 = 0.5 * SD_a0 * nu_z / sqrtpi
        if (cs_z > 0.) SD_a0 = SD_a0 / cs_z
      endif disk
!
    endsubroutine init_cooling_SD93
!***********************************************************************
    elemental real function cool_SD93(eth, rho)
!
!  Add Sutherland & Dopita (1993) cooling function.
!
!  02-feb-13/ccyang: coded.
!
      use General, only: spline
      use EquationOfState, only: get_soundspeed
!
      real, intent(in) :: eth, rho
!
      logical :: err
      real, dimension(1) :: logLambda
      real :: temp, cs2
!
!  Find the temperature.
!
      cool_SD93 = 0.
      temp = get_temperature(eth, rho)
      if (temp <= 0.) return
!
!  Interpolate the cooling table.
!
      call spline(SD_logTT, SD_logLambda, (/log10(temp)/), logLambda, SD_nt, 1, err)
      if (err) return
!
!  Calculate the net cooling rate.
!    For 3D run: rho is volume mass density
!    For 2D run: rho is surface mass density
!    cs_z, if > 0, fix the disk vertical scale height
!
      cool_SD93 = SD_a0 * rho**2 * 10.**logLambda(1)
      disk: if (nzgrid == 1 .and. cs_z <= 0.) then
        call get_soundspeed(real(temp / unit_temperature), cs2)
        cool_SD93 = cool_SD93 / sqrt(cs2)
      endif disk
!
    endfunction cool_SD93
!***********************************************************************
    elemental real function cool_KI02_SD93(eth, rho) result(cool)
!
!  Smoothly connects the KI02 and the SD93 heating/cooling functions.
!
!  04-feb-13/ccyang: coded.
!
      real, intent(in) :: eth, rho
!
      real, parameter :: logTTmin = 4., logTTmax = 4.2
      real, parameter :: dlogTT = logTTmax - logTTmin
      real :: temp, logTT, cool1, cool2
!
!  Find the temperature.
!
      temp = get_temperature(eth, rho)
      if (temp <= 0.) then
        cool = 0.
        return
      endif
      logTT = log10(temp)
!
!  Linearly combine the two heating/cooling functions.
!
      step: if (logTT < logTTmin) then
        cool = -heat_KI02(eth, rho)
      elseif (logTT > logTTmax) then step
        cool = cool_SD93(eth, rho)
      else step
        cool1 = -heat_KI02(eth, rho)
        cool2 = cool_SD93(eth, rho)
        temp = cos(pi * (logTT - logTTmin) / dlogTT)
        cool = 0.5 * ((1. + temp) * cool1 + (1. - temp) * cool2)
      endif step
!
    endfunction cool_KI02_SD93
!***********************************************************************
    elemental subroutine jeans_floor(eth, rho)
!
!  Apply a floor to energy where Jeans length is unresolved.
!
!  29-aug-11/ccyang: coded.
!
      real, intent(inout) :: eth
      real, intent(in) :: rho
!
      real :: eth_min
!
!  Find the Jeans energy floor.
!
      if (nzgrid /= 1) then
!       3D run: rho is volume mass density
      else
!       2D run: rho is surface mass density
        eth_min = Jeans_c0 * rho**2
      endif
!
!  Impose the floor where needed.
!
      if (eth < eth_min) eth = eth_min
!
    endsubroutine
!***********************************************************************
    pure real function get_temperature(eth, rho) result(temp)
!
!  Finds the absolute temperature in Kelvins assuming the ideal-gas EOS.
!
!  09-sep-14/ccyang: coded.
!
      real, intent(in) :: eth, rho
!
!  Find the temperature.
!
      if (lstratz) then
        temp = cv1_temp * (1.0 + eth) / (1.0 + rho)
      else
        temp = cv1_temp * eth / rho
      endif
!
    endfunction get_temperature
!***********************************************************************
    subroutine update_char_vel_energy(f)
!
! TB implemented.
!
!   25-sep-15/MR+joern: coded
!
      real, dimension (mx,my,mz,mfarray), intent(in) :: f

      call keep_compiler_quiet(f)

      call warning('update_char_vel_energy', &
                   'characteristic velocity not yet implemented for thermal energy')

    endsubroutine update_char_vel_energy
!***********************************************************************
    subroutine pushpars2c(p_par)

    use Syscalls, only: copy_addr
 
    integer, parameter :: n_pars=1
    integer(KIND=ikind8), dimension(n_pars) :: p_par

    call copy_addr(chi,p_par(1))

    endsubroutine pushpars2c
!***********************************************************************
!********************************************************************
!********************************************************************
!************        DO NOT DELETE THE FOLLOWING        *************
!********************************************************************
!**  This is an automatically generated include file that creates  **
!**  copies dummy routines from nospecial.f90 for any Special      **
!**  routines not implemented in this file                         **
!**                                                                **
    include 'energy_dummies.inc'
!***********************************************************************
endmodule Energy
