! $Id$
!
!  This module can replace the energy module by using lnT or T (with
!  ltemperature_nolog=.true.) as dependent variable. For a perfect gas
!  with constant coefficients (no ionization) we have:
!  (1-1/gamma) * cp*T = cs20 * exp( (gamma-1)*ln(rho/rho0)-gamma*s/cp )
!
!  Note that to use lnTT as thermal variable, you may rather want to use
!  energy.f90 with pretend_lnTT=.true. As of March 2007, entropy.f90
!  has way more options and features than temperature_idealgas.f90.
!
!** 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 = .true.
! CPARAM logical, parameter :: lthermal_energy = .false.
!
! MVAR CONTRIBUTION 1
! MAUX CONTRIBUTION 0
!
! PENCILS PROVIDED Ma2; uglnTT; ugTT; fpres(3); tcond; sglnTT(3); d2xlnTT; d2zlnTT, advec_cs2
!
!***************************************************************
module Energy
!
! 12-may-12/MR: made ampl_lnTT a vector; added parameters for initialization
!               by mode to input pars
!
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages
!
  implicit none
!
  include 'energy.h'
!
  real :: radius_lnTT=0.1, widthlnTT=2*epsi
  real, dimension (ninit) :: ampl_lnTT=0.0, ampl_TT=0.0
  real :: lnTT_const=0.0, TT_const=1.0
  real :: Kgperp=0.0, Kgpara=0.0
  real :: chi=impossible, chi_jump=1., chi_z0=0.0, chi_zwidth=0.0, chi_r_reduce=0.0
  real :: zbot=0.0, ztop=0.0
  real :: center1_x=0.0, center1_y=0.0, center1_z=0.0
  real :: r_bcz=0.0, chi_shock=0.0, chi_hyper3=0.0, chi_hyper3_mesh=5.0
  real :: Tbump=0.0, Kmin=0.0, Kmax=0.0, hole_slope=0.0, hole_width=0.0
  real, dimension(nz) :: zmask_temp, zmask_emiss
  real, dimension(nzgrid) :: zmask_temp_global
  real, dimension(2) :: temp_zaver_range=(/-max_real,max_real/)
  real :: mu=1.0
  real :: hcond0=impossible, hcond1=1.0, hcond2=1.0, Fbot=impossible,Ftop=impossible
  real :: luminosity=0.0, wheat=0.1, rcool=0.0, wcool=0.1, cool=0.0
  real :: beta_bouss=-1.0, h_sld_ene=2.0, nlf_sld_ene=1.0
  integer :: temp_zmask_count=1
  integer, parameter :: nheatc_max=3
  logical, pointer :: lpressuregradient_gas
  logical :: ladvection_temperature=.true.
  logical :: lupw_lnTT=.false., lcalc_heat_cool=.false., lheatc_hyper3=.false.
  logical :: lheatc_Kconst=.false., lheatc_Kprof=.false., lheatc_Karctan=.false.
  logical :: lheatc_tensordiffusion=.false., lheatc_hyper3_mesh=.false.
  logical :: lheatc_chiconst=.false., lheatc_chiconst_accurate=.false., lheatc_chi_reduce_ddr=.false.
  logical :: lfreeze_lnTTint=.false., lfreeze_lnTText=.false.
  logical :: lhcond_global=.false., lheatc_chicubicstep=.false.
  logical :: lheatc_shock=.false., lheatc_hyper3_polar=.false.
  logical :: lheatc_Ktherm=.false.
  logical, target :: lheatc_kramers=.false.
  real, target :: hcond0_kramers=0.0, nkramers=0.0
  logical :: lviscosity_heat=.true.
  logical :: lcalc_TTmean=.false.
!  integer :: iglobal_hcond=0
!  integer :: iglobal_glhc=0
  logical :: lenergy_slope_limited=.false.
  logical :: linitial_log=.false.
  logical, pointer :: lreduced_sound_speed   !=>.false.
!  logical, pointer :: lscale_to_cs2top
  character (len=labellen), dimension(nheatc_max) :: iheatcond='nothing'
  character (len=labellen) :: borderss='nothing'
  character (len=labellen) :: div_sld_ene='2nd'
  character (len=labellen), dimension(ninit) :: initlnTT='nothing'
  complex :: coef_lnTT=0.
  character (len=intlen) :: iinit_str
  real    :: kx_lnTT=1.,ky_lnTT=1.,kz_lnTT=1.
  logical :: lADI_mixed=.false., lmultilayer=.false.
  real, pointer :: PrRa   ! preliminary
  real, pointer :: reduce_cs2
  real, target :: mpoly0=1.5, mpoly1=1.5, mpoly2=1.5
!
  real, dimension(nz) :: TTmz, gTTmz 
!
  real, dimension(3) :: gradTT0=(/0.0,0.0,0.0/)
!
!  Init parameters.
!
  namelist /entropy_init_pars/ &
      initlnTT, radius_lnTT, ampl_lnTT, widthlnTT, lnTT_const, TT_const, &
      center1_x, center1_y, center1_z, mpoly0, mpoly1, mpoly2, r_bcz, Fbot, &
      Tbump, Kmin, Kmax, hole_slope, hole_width, ltemperature_nolog, &
      linitial_log, hcond0, luminosity, wheat, coef_lnTT, kx_lnTT, ky_lnTT, kz_lnTT, &
      temp_zaver_range, ampl_TT
!
!  Run parameters.
!
  namelist /entropy_run_pars/ &
      lupw_lnTT, ladvection_temperature, chi, iheatcond, chi_hyper3_mesh, &
      lheatc_chiconst_accurate, hcond0, lcalc_heat_cool, lfreeze_lnTTint, &
      lfreeze_lnTText, widthlnTT, mpoly0, mpoly1, mpoly2, lhcond_global, &
      lviscosity_heat, chi_hyper3, chi_shock, Fbot, Tbump, Kmin, Kmax, &
      hole_slope, hole_width, Kgpara, Kgperp, lADI_mixed, rcool, wcool, &
      cool, beta_bouss, borderss, lmultilayer, lcalc_TTmean, &
      temp_zaver_range, h_sld_ene, nlf_sld_ene, div_sld_ene, &
      gradTT0, w_sldchar_ene, chi_z0, chi_jump, chi_zwidth, &
      hcond0_kramers, nkramers,chi_r_reduce
!
!  Diagnostic variables for print.in
! (needs to be consistent with reset list below)
!
  integer :: idiag_TTmax=0        ! DIAG_DOC: $\max (T)$
  integer :: idiag_gTmax=0        ! DIAG_DOC: $\max (|\nabla T|)$
  integer :: idiag_TTmin=0        ! DIAG_DOC: $\min (T)$
  integer :: idiag_TTm=0          ! DIAG_DOC: $\left< T \right>$
  integer :: idiag_TTzmask=0      ! DIAG_DOC: $\left< T \right>$ for
                                  ! DIAG_DOC: the temp_zaver_range
  integer :: idiag_TT2m=0         ! DIAG_DOC: $\left< T^2 \right>$
  integer :: idiag_TugTm=0        ! DIAG_DOC: $\left< T\uv\cdot\nabla T \right>$
  integer :: idiag_Trms=0         ! DIAG_DOC: $\sqrt{\left< T^2 \right>}$
  integer :: idiag_uxTm=0         ! DIAG_DOC: $\left< u_x T \right>$
  integer :: idiag_uyTm=0         ! DIAG_DOC: $\left< u_y T \right>$
  integer :: idiag_uzTm=0         ! DIAG_DOC: $\left< u_z T \right>$
  integer :: idiag_gT2m=0         ! DIAG_DOC: $\left< (\nabla T)^2 \right>$
  integer :: idiag_guxgTm=0       ! DIAG_DOC: $\left< \nabla u_x \cdot \nabla T \right>$
  integer :: idiag_guygTm=0       ! DIAG_DOC: $\left< \nabla u_y \cdot \nabla T \right>$
  integer :: idiag_guzgTm=0       ! DIAG_DOC: $\left< \nabla u_z \cdot \nabla T \right>$
  integer :: idiag_Tugux_uxugTm=0 ! DIAG_DOC: $\left< T \uv\cdot\nabla u_x + u_x \uv\cdot\nabla T \right>
                                  ! DIAG_DOC:  =\left< \uv\cdot\nabla(u_x T) \right>$
  integer :: idiag_Tuguy_uyugTm=0 ! DIAG_DOC: $\left< T \uv\cdot\nabla u_y + u_y \uv\cdot\nabla T \right>
                                  ! DIAG_DOC:  =\left< \uv\cdot\nabla(u_y T) \right>$
  integer :: idiag_Tuguz_uzugTm=0 ! DIAG_DOC: $\left< T \uv\cdot\nabla u_z + u_z \uv\cdot\nabla T \right>
                                  ! DIAG_DOC:  =\left< \uv\cdot\nabla(u_z T) \right>$
  integer :: idiag_Tdxpm=0        ! DIAG_DOC: $\left< T dp/dx \right>$
  integer :: idiag_Tdypm=0        ! DIAG_DOC: $\left< T dp/dy \right>$
  integer :: idiag_Tdzpm=0        ! DIAG_DOC: $\left< T dp/dz \right>$
!
  integer :: idiag_fradtop=0  ! DIAG_DOC: $<-K{dT\over dz}>_{\text{top}}$
                              ! DIAG_DOC: \quad(top radiative flux)
  integer :: idiag_fradbot=0  ! DIAG_DOC: $<-K{dT\over dz}>_{\text{bot}}$
                              ! DIAG_DOC: \quad(bottom radiative flux)
  integer :: idiag_yHmax=0    ! DIAG_DOC: DOCUMENT ME
  integer :: idiag_yHmin=0    ! DIAG_DOC: DOCUMENT ME
  integer :: idiag_yHm=0      ! DIAG_DOC: DOCUMENT ME
  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_eem=0      ! DIAG_DOC: $\left< e \right> =
                              ! DIAG_DOC:  \left< c_v T \right>$
                              ! DIAG_DOC: \quad(mean internal energy)
  integer :: idiag_ethtot=0   ! DIAG_DOC: $\int_V\varrho e\,dV$
                              ! DIAG_DOC:   \quad(total thermal energy)
  integer :: idiag_ssm=0      ! DIAG_DOC: $\overline{S}$
  integer :: idiag_thcool=0   ! DIAG_DOC: $\tau_{\rm cool}$
  integer :: idiag_ppm=0      ! DIAG_DOC: $\overline{P}$
  integer :: idiag_csm=0      ! DIAG_DOC: $\overline{c}_{\rm s}$
  integer :: idiag_csmax=0    ! DIAG_DOC: $\max (c_{\rm s})$
  integer :: idiag_dtc=0        ! DIAG_DOC: $\delta t/[c_{\delta t}\,\delta_x
                                ! DIAG_DOC:   /\max c_{\rm s}]$
                                ! DIAG_DOC:   \quad(time step relative to
                                ! DIAG_DOC:   acoustic time step;
                                ! DIAG_DOC:   see \S~\ref{time-step})
  integer :: idiag_dtchi=0      ! DIAG_DOC: $\delta t / [c_{\delta t,{\rm v}}\,
                                ! DIAG_DOC:   \delta x^2/\chi_{\rm max}]$
                                ! DIAG_DOC:   \quad(time step relative to time
                                ! DIAG_DOC:   step based on heat conductivity;
                                ! DIAG_DOC:   see \S~\ref{time-step})
!
!
! xy averaged diagnostics given in xyaver.in written every it1d timestep
!
  integer :: idiag_ppmz=0       ! XYAVG_DOC: $\left<p\right>_{xy}$
  integer :: idiag_ppuzmz=0     ! XYAVG_DOC:
  integer :: idiag_TTmz=0       ! XYAVG_DOC: $\left<T\right>_{xy}$
  integer :: idiag_ethmz=0      ! XYAVG_DOC: $\left< e_{\text{th}}
                                ! XYAVG_DOC: \right>_{xy}$
  integer :: idiag_ethuxmz=0    ! XYAVG_DOC:
  integer :: idiag_ethuymz=0    ! XYAVG_DOC:
  integer :: idiag_ethuzmz=0    ! XYAVG_DOC:
  integer :: idiag_fpresxmz=0   ! XYAVG_DOC: $\left<(\nabla p)_x\right>_{xy}$
  integer :: idiag_fpresymz=0   ! XYAVG_DOC: $\left<(\nabla p)_y\right>_{xy}$
  integer :: idiag_fpreszmz=0   ! XYAVG_DOC: $\left<(\nabla p)_z\right>_{xy}$
!
  integer :: idiag_TT2mz=0      ! XYAVG_DOC: $\left<T^2\right>_{xy}$
  integer :: idiag_uxTmz=0      ! XYAVG_DOC: $\left<u_x T\right>_{xy}$
  integer :: idiag_uyTmz=0      ! XYAVG_DOC: $\left<u_y T\right>_{xy}$
  integer :: idiag_uzTmz=0      ! XYAVG_DOC: $\left<u_z T\right>_{xy}$
  integer :: idiag_fradmz=0     ! XYAVG_DOC: $\left<F_{\rm rad}\right>_{xy}$
  integer :: idiag_fconvmz=0    ! XYAVG_DOC: $\left<c_p \varrho u_z 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}$
  integer :: idiag_ethuxmx=0    ! YZAVG_DOC:
!
! variables for slices given in video.in
!
  real, dimension(:,:), allocatable :: pp_xz,pp_yz,pp_xy,pp_xy2,pp_xy3,pp_xy4,pp_xz2
  real, dimension(:,:,:,:,:), allocatable :: pp_r
!
! y averaged diagnostics given in yaver.in
!
  integer :: idiag_TTmxz=0      ! YAVG_DOC: $\left<T\right>_{y}$
  integer :: idiag_EmAIA94mxz=0 ! YAVG_DOC: Emission off AIA 94 channel
                                ! YAVG_DOC: integrated over y direction
  integer :: idiag_EmAIA131mxz=0! YAVG_DOC: Emission off AIA 131 channel
                                ! YAVG_DOC: integrated over y direction
  integer :: idiag_EmAIA171mxz=0! YAVG_DOC: Emission off AIA 171 channel
                                ! YAVG_DOC: integrated over y direction
  integer :: idiag_EmAIA193mxz=0! YAVG_DOC: Emission off AIA 193 channel
                                ! YAVG_DOC: integrated over y direction
  integer :: idiag_EmAIA211mxz=0! YAVG_DOC: Emission off AIA 211 channel
                                ! YAVG_DOC: integrated over y direction
  integer :: idiag_EmAIA304mxz=0! YAVG_DOC: Emission off AIA 304 channel
                                ! YAVG_DOC: integrated over y direction
  integer :: idiag_EmAIA335mxz=0! YAVG_DOC: Emission off AIA 335 channel
                                ! YAVG_DOC: integrated over y direction
  integer :: idiag_EmXRTmxz=0   ! YAVG_DOC: Emission off XRT Al-poly channel
                                ! YAVG_DOC: integrated over y direction
!
! z averaged diagnostics given in zaver.in
!
  integer :: idiag_TTmxy=0      ! ZAVG_DOC: $\left<T\right>_{z}$
  integer :: idiag_EmAIA94mxy=0 ! ZAVG_DOC: Emission off AIA 94 channel
                                ! ZAVG_DOC: integrated over z direction
  integer :: idiag_EmAIA131mxy=0! ZAVG_DOC: Emission off AIA 131 channel
                                ! ZAVG_DOC: integrated over z direction
  integer :: idiag_EmAIA171mxy=0! ZAVG_DOC: Emission off AIA 171 channel
                                ! ZAVG_DOC: integrated over z direction
  integer :: idiag_EmAIA193mxy=0! ZAVG_DOC: Emission off AIA 193 channel
                                ! ZAVG_DOC: integrated over z direction
  integer :: idiag_EmAIA211mxy=0! ZAVG_DOC: Emission off AIA 211 channel
                                ! ZAVG_DOC: integrated over z direction
  integer :: idiag_EmAIA304mxy=0! ZAVG_DOC: Emission off AIA 304 channel
                                ! ZAVG_DOC: integrated over z direction
  integer :: idiag_EmAIA335mxy=0! ZAVG_DOC: Emission off AIA 335 channel
                                ! ZAVG_DOC: integrated over z direction
  integer :: idiag_EmXRTmxy=0   ! ZAVG_DOC: Emission off XRT Al-poly channel
                                ! ZAVG_DOC: integrated over z direction
!
  integer :: ivid_pp=0
!
  real, dimension(nx) :: diffus_chi,diffus_chi3,hcond,cooling
  real :: gamma, gamma1, gamma_m1, cp1
!
  integer :: enum_borderss = 0
!
  contains
!***********************************************************************
    subroutine register_energy
!
!  Initialise variables which should know that we solve an energy
!  equation: ilnTT, etc; increase nvar accordingly.
!
!  06-nov-01/wolf: coded
!  18-may-12/MR: shared variable PrRa fetched from hydro
!  03-apr-20/joern: restructured and fixed slope-limited diffusion
!
      use FArrayManager, only: farray_register_pde, farray_index_append, &
                               farray_register_auxiliary
      use SharedVariables, only: put_shared_variable
!
!  Register TT or lnTT, depending on whether or not ltemperature_nolog
!
      if (ltemperature_nolog) then
        call farray_register_pde('TT',iTT)
        ilnTT=iTT
        !call farray_index_append('ilnTT',ilnTT)
      else
        call farray_register_pde('lnTT',ilnTT)
      endif
!
!  Identify version number.
!
      if (lroot) call svn_id( &
           "$Id$")
!
      if (any(iheatcond=='temperature-slope-limited')) then
        if (dimensionality<3)lisotropic_advection=.true.
        lslope_limit_diff = .true.
        if (isld_char == 0) then
          call farray_register_auxiliary('sld_char',isld_char,communicated=.true.,on_gpu=lgpu)
          if (lroot) write(15,*) 'sld_char = fltarr(mx,my,mz)*one'
          aux_var(aux_count)=',sld_char'
          if (naux+naux_com <  maux+maux_com) aux_var(aux_count)=trim(aux_var(aux_count))//' $'
          aux_count=aux_count+1
        endif
      endif
!
!  Writing files for use with IDL.
!
      if (lroot) then
        if (maux == 0) then
           if (nvar < mvar) write(4,*) ',lnTT $'
           if (nvar == mvar) write(4,*) ',lnTT'
        else
           write(4,*) ',lnTT $'
        endif
        write(15,*) 'lnTT = fltarr(mx,my,mz)*one'
      endif
!
!  Now we share several variables.
!
      call put_shared_variable('hcond0', hcond0, caller='register_energy')
      call put_shared_variable('hcond1', hcond1)
      call put_shared_variable('hcond2', hcond2)
      call put_shared_variable('lmultilayer', lmultilayer)
      call put_shared_variable('widthlnTT', widthlnTT)
      call put_shared_variable('Fbot', Fbot)
      call put_shared_variable('lADI_mixed', lADI_mixed)
      call put_shared_variable('lviscosity_heat',lviscosity_heat)
      call put_shared_variable('mpoly0', mpoly0)
      call put_shared_variable('mpoly1', mpoly1)
      call put_shared_variable('mpoly2', mpoly2)
      call put_shared_variable('hcond0_kramers', hcond0_kramers)
      call put_shared_variable('nkramers', nkramers)
      call put_shared_variable('lheatc_kramers', lheatc_kramers)
      if (lsolid_cells) then
        if (.not. lchemistry) call put_shared_variable('chi', chi)
        call put_shared_variable('ladvection_temperature',ladvection_temperature)
        call put_shared_variable('lheatc_chiconst',lheatc_chiconst)
        call put_shared_variable('lupw_lnTT',lupw_lnTT)
      endif
!
    endsubroutine register_energy
!***********************************************************************
    subroutine initialize_energy(f)
!
!  Called by run.f90 after reading parameters, but before the time loop.
!
!  21-jul-2002/wolf: coded
!
      use FArrayManager, only: farray_register_global
      use Gravity, only: gravz, compute_gravity_star
      use EquationOfState, only : cs2bot, cs2top, select_eos_variable, get_gamma_etc
      use Sub, only: step,der_step
      use SharedVariables, only: get_shared_variable
      use Slices_methods, only: alloc_slice_buffers
      use BorderProfiles, only: request_border_driving
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: dhcond
      logical :: lnothing
      integer :: i
      real :: star_cte, cp

      if (lstart .and. .not.leos) call warning('initialize_energy', &
            'EOS=noeos, but possibly an EQUATION OF STATE is needed')
!
!  Set iTT equal to ilnTT if we are considering non-logarithmic temperature.
!
      if (ltemperature_nolog) then
        call select_eos_variable('TT',iTT)
      else
        call select_eos_variable('lnTT',ilnTT)
      endif
!
      call get_gamma_etc(gamma,cp)
      gamma1=1./gamma; gamma_m1=gamma-1.
      cp1=1./cp
!
!  Freeze temperature.
!
      if (lfreeze_lnTTint) lfreeze_varint(ilnTT)=.true.
      if (lfreeze_lnTText) lfreeze_varext(ilnTT)=.true.
!
!  Check whether we want heat conduction.
!
      lheatc_Kconst= .false.
      lheatc_Kprof= .false.
      lheatc_Karctan= .false.
      lheatc_tensordiffusion=.false.
      lheatc_chiconst = .false.
      lheatc_chicubicstep = .false.
      lheatc_kramers= .false.
      lheatc_chi_reduce_ddr=.false.
!
!  Initialize thermal diffusion.
!
      lheatc_shock=.false.
      lheatc_hyper3=.false.
      lheatc_hyper3_mesh=.false.
      lheatc_hyper3_polar=.false.
      lenergy_slope_limited=.false.
!
!  initialize lnothing. It is needed to prevent multiple output.
!
      lnothing = .false.
!
!  Different choices of heat conduction (if any).
!
      do i=1,nheatc_max
      select case (iheatcond(i))
        case ('K-const')
          lheatc_Kconst=.true.
          call information('initialize_energy', &
          ' heat conduction: K=cst --> gamma*K/rho/TT/cp*div(T*grad lnTT)')
        case ('K-profile')
          lheatc_Kprof=.true.
          lmultilayer=.true.
          call information('initialize_energy','heat conduction: K=K(z) or K=K(r)')
        case ('K-arctan')
          lheatc_Karctan=.true.
          if (.not. ltemperature_nolog) &
            call fatal_error('initialize_energy','K-arctan only valid for TT')
          if (lADI_mixed .and. .not. lADI) &
            call fatal_error('initialize_energy','K-arctan with lADI_mixed=T while lADI=F?')
          call information('initialize_energy','heat conduction: arctan profile')
        case ('chi-const')
          lheatc_chiconst=.true.
          call information('initialize_energy','heat conduction: constant chi')
        case ('chi-no-ddrsph')
          lheatc_chi_reduce_ddr=.true.
          call information('initialize_energy','heat conduction: reduce chi in r direction')
        case ('chi-cubicstep')
          lheatc_chicubicstep=.true.
          call information('initialize_energy','heat conduction: cubic step profile of chi')
        case('K-therm')
          lheatc_Ktherm=.true.
          call information('initialize_energy','heat conduction: temperature dependent K')
        case ('chi-hyper3')
          lheatc_hyper3=.true.
          call information('initialize_energy','hyper conductivity')
        case ('hyper3_mesh','hyper3-mesh')
          lheatc_hyper3_mesh=.true.
          call information('initialize_energy','hyper mesh conductivity')
        case ('hyper3_cyl','hyper3-cyl','hyper3_sph','hyper3-sph')
          lheatc_hyper3_polar=.true.
          call information('initialize_energy', 'hyper conductivity: polar coords')
        case ('shock','chi-shock')
          lheatc_shock=.true.
          call information('initialize_energy','shock conductivity')
        case ('tensor-diffusion')
          lheatc_tensordiffusion=.true.
          if (lroot) print*, 'heat conduction: tensor diffusion'
        case ('kramers')
          lheatc_kramers=.true.
          if (lroot) print*, 'heat conduction: kramers'
        case ('temperature-slope-limited')
          lenergy_slope_limited=.true.
          if (lroot) print*, 'heat conduction: slope limited diffusion'
          if (lroot) print*, 'heat conduction: using ',trim(div_sld_ene),' order'
        case ('nothing')
          if (lroot .and. (.not. lnothing)) print*,'heat conduction: nothing'
        case default
          call fatal_error('initialize_energy','no such iheatcond: '//trim(iheatcond(i)))
        endselect
        lnothing=.true.
      enddo
!
!  Compute and store hcond and dhcond if hcond_global=.true.
!
      if (lhcond_global) then
        call farray_register_global("hcond",iglobal_hcond)
        call farray_register_global("glhc",iglobal_glhc)
        do n=n1,n2
        do m=m1,m2
          hcond  = hcond0*(1. + (hcond1-1.)*    step(x(l1:l2),r_bcz,-widthlnTT))
          dhcond = hcond0*      (hcond1-1.)*der_step(x(l1:l2),r_bcz,-widthlnTT)
          f(l1:l2,m,n,iglobal_hcond)=hcond
          f(l1:l2,m,n,iglobal_glhc)=dhcond
        enddo
        enddo
      endif
!
      if (initlnTT(1)=='gaussian') then
!
!  Needed when one only works with temperature_idealgas to check the radiative
!  diffusion term, i.e. one solves d(TT)/dt=gamma*chi*del2(TT) with bcz='cT'
!  (all other modules are down).
!
        cs2bot=gamma_m1*f(l1,4,n1,ilnTT)
        cs2top=gamma_m1*f(l1,4,n2,ilnTT)
      endif
!
!  Some tricks regarding Fbot and hcond0 when bcz1='c1' (constant flux).
!
      if (bcz12(ilnTT,1)=='c1' .and. lrun) then

        if (.not. leos) &
          call fatal_error('initialize_energy','EOS=noeos, but BC "c1" requires an EQUATION OF STATE')
!
        if (Fbot==impossible .and. hcond0 /= impossible) then
          Fbot=-gamma/gamma_m1*hcond0*gravz/(mpoly0+1.0)
          if (lroot) print*, 'initialize_energy: Calculated Fbot = ', Fbot
        endif
        if (hcond0==impossible .and. Fbot /= impossible) then
          hcond0=-Fbot*gamma_m1/gamma*(mpoly0+1.0)/gravz
          if (lroot) print*, 'initialize_energy: Calculated hcond0 = ', hcond0
        endif
        if (Fbot==impossible .and. hcond0==impossible) &
          call fatal_error('initialize_energy','Both Fbot and hcond0 are unknown')
      endif
!
      if (initlnTT(1)=='star_heat') then

        if (.not. leos) call fatal_error('initialize_energy', &
                        'EOS=noeos, but star heating requires an EQUATION OF STATE')
!
        if (lroot) print*,'star_heat: compute the gravity profile'
        ! compute the gravity profile
        star_cte=(mpoly0+1.)/hcond0*gamma_m1/gamma
        call compute_gravity_star(f, wheat, luminosity, star_cte)
        if (rcool==0.) rcool=r_ext
      endif
!
      if (lmultilayer) then

        if (.not. leos) call fatal_error('initialize_energy', &
                        'EOS=noeos, but multilayer requires an EQUATION OF STATE')
!
        hcond1=(mpoly1+1.)/(mpoly0+1.)
        hcond2=(mpoly2+1.)/(mpoly0+1.)
      endif
!
!  A word of warning...
!
      if (lheatc_Kconst .and. hcond0==0.0) call warning('initialize_energy', 'hcond0 is zero')
      if (lheatc_Ktherm .and. hcond0==0.0) call warning('initialize_energy','hcond0 is zero')
      if (lheatc_Kprof .and. hcond0==0.0) call warning('initialize_energy', 'hcond0 is zero')
      if (lheatc_chiconst .and. chi==0.0) call warning('initialize_energy','chi is zero')
      if (lheatc_chicubicstep .and. chi==0.0) call warning('initialize_energy','chi is zero')
      if (lheatc_kramers .and. hcond0==0.0) call warning('initialize_energy','hcond0 is zero')
      if (lrun) then
        if (lheatc_hyper3 .and. chi_hyper3==0.0) &
            call fatal_error('initialize_energy','Conductivity coefficient chi_hyper3 is zero')
        if (lheatc_shock .and. chi_shock==0.0) &
            call fatal_error('initialize_energy','Conductivity coefficient chi_shock is zero')
      endif
      if (iheatcond(1)=='nothing') then
        if (hcond0/=impossible) call warning('initialize_energy','No heat conduction, but hcond0/=0')
        if (chi/=impossible) call warning('initialize_energy','No heat conduction, but chi/=0')
      endif
      if (lADI_mixed .and. iheatcond(1) /= 'K-arctan') &
        call fatal_error("initialize_energy","lADI_mixed=T while iheatcond /= K-arctan?")
!
      if (llocal_iso) call fatal_error('initialize_energy', &
           'llocal_iso switches on the local isothermal approximation. ' // &
           'Use ENERGY=noenergy in src/Makefile.local')
!
!  Compute mask for z-averaging where z is in temp_zaver_range.
!  Normalize such that the average over the full domain
!  gives still unity.
!
      if (n1 == n2) then
        zmask_temp = 1.
        zmask_temp_global=1.
      else
        where (z(n1:n2) >= temp_zaver_range(1) .and. z(n1:n2)  <= temp_zaver_range(2))
          zmask_temp = 1.
        elsewhere
          zmask_temp = 0.
        endwhere
        where (zglobal(nghost+1:mzgrid-nghost) >= temp_zaver_range(1) .and. &
               zglobal(nghost+1:mzgrid-nghost) <= temp_zaver_range(2))
          zmask_temp_global = 1.
        elsewhere
          zmask_temp_global = 0.
        endwhere
        temp_zmask_count =  max(count(zmask_temp_global ==  1.0),1)
        zmask_temp = zmask_temp*float(nzgrid)/float(temp_zmask_count)
      endif
!
!  debug output
!
      if (lroot.and.ip<14) print*,'zmask_temp=' ,zmask_temp
!
!  Tell the BorderProfiles module if we intend to use border driving, so
!  that the module can request the right pencils.
!
      select case (borderss)
      case ('zero','0','constant','initial-condition')
        call request_border_driving(borderss)
      case ('nothing')
        if (lroot.and.ip<=5) print*, "set_border_entropy: borderss='nothing'"
      case default
        call fatal_error('set_border_entropy','no such borderss: '//trim(borderss))
      endselect
!
!  logical variable lpressuregradient_gas shared with hydro modules
!
      call get_shared_variable('lpressuregradient_gas',lpressuregradient_gas,caller='initialize_energy')
!
!  real variable PrRa shared with hydro modules, used for Boussinesq
!
      if (lboussinesq.and.lviscosity_heat) call get_shared_variable('PrRa',PrRa)
!
!  Check if reduced sound speed is used
!
      if (ldensity) then
        call get_shared_variable('lreduced_sound_speed', lreduced_sound_speed)
        if (lreduced_sound_speed) then
          call get_shared_variable('reduce_cs2',reduce_cs2)
!          call get_shared_variable('lscale_to_cs2top',lscale_to_cs2top)
        endif
      endif
!
      if (ivid_pp/=0) call alloc_slice_buffers(pp_xy,pp_xz,pp_yz,pp_xy2,pp_xy3,pp_xy4,pp_xz2,pp_r)
!
    endsubroutine initialize_energy
!***********************************************************************
    subroutine init_energy(f)
!
!  Initialise lnTT or TT; called from start.f90.
!
!  13-dec-2002/axel+tobi: adapted from init_energy
!
!  initialise energy; called from start.f90
!  07-nov-01/wolf: coded
!  24-nov-02/tony: renamed for consistency (i.e. init_[variable name])
!  12-may-12/MR: initialization with mode added
!
      use General,  only: itoa
      use Sub,      only: blob
      use InitialCondition, only: initial_condition_ss
      use EquationOfState, only: cs2bot, cs2top, cs20, lnrho0, rho0
      use Gravity, only: gravz
      !use Initcond, only: modes
      use Initcond
!
      real, dimension (mx,my,mz,mfarray), intent (inout) :: f
      real, dimension (mz) :: TTz
!
      integer :: j
      logical :: lnothing=.true.
      real :: haut, Rgas, Ttop, alpha, beta, expo, ztop
!
      do j=1,ninit
!
        if (initlnTT(j)/='nothing') then
!
          lnothing=.false.
!
          iinit_str=itoa(j)
!
!  Select between various initial conditions.
!
          select case (initlnTT(j))
          case ('zero', '0'); f(:,:,:,ilnTT) = 0.
!
          case ('const_lnTT'); f(:,:,:,ilnTT)=f(:,:,:,ilnTT)+lnTT_const
!
          case ('const_TT')
            if (ltemperature_nolog) then
              f(:,:,:,iTT)=f(:,:,:,iTT)+TT_const
            else
              f(:,:,:,ilnTT)=f(:,:,:,ilnTT)+log(TT_const)
            endif
            cs2bot=gamma_m1*TT_const
            cs2top=gamma_m1*TT_const
!
! NS: added option for gaussian noise
          case ('gaussian-noise')
                  if (ltemperature_nolog) then
                          call gaunoise(ampl_TT(j),f,iTT,iTT)
                  else
                          call gaunoise(ampl_TT(j),f,ilnTT,ilnTT)
                  endif
          case ('const_dTTdz')
            TTz=TT_const+z
            if (ltemperature_nolog) then
              f(:,:,:,iTT)=f(:,:,:,iTT)+spread(spread(TTz,1,mx),2,my)
            else
              f(:,:,:,ilnTT)=f(:,:,:,ilnTT)+spread(spread(log(TTz),1,mx),2,my)
            endif
            cs2bot=gamma_m1*TTz(n1)
            cs2top=gamma_m1*TTz(n2)
!
          case ('linz_fromBC')
!
!  Linear profile between temperature values at bottom and top supposed to be given
!  in fbcz.
!
            if (ltemperature_nolog) then
              TTz=fbcz(iTT,1)+(z-xyz0(3))/Lxyz(3)*(fbcz(iTT,2)-fbcz(iTT,1))
              f(:,:,:,iTT)=f(:,:,:,iTT)+spread(spread(TTz,1,mx),2,my)
              cs2bot=gamma_m1*fbcz(iTT,1)
              cs2top=gamma_m1*fbcz(iTT,2)
            else
              call not_implemented('init_ss','IC "linz_fromBC" for logarithmic temperature')
            endif
!
          case ('mode')
!
            if (ltemperature_nolog) then
              call modes(ampl_lnTT(j),coef_lnTT,f,iTT,kx_lnTT,ky_lnTT,kz_lnTT)
            else
              call modes(ampl_lnTT(j),coef_lnTT,f,ilnTT,kx_lnTT,ky_lnTT,kz_lnTT)
            endif
!
          case ('single_polytrope'); call single_polytrope(f)
!
          case ('piecew-poly'); call piecew_poly(f)
!
          case ('gaussian')
            do n=n1,n2
              f(l1:l2,4,n,ilnTT)=exp(-(x(l1:l2)/radius_lnTT)**2)*exp(-((z(n)-0.5)/radius_lnTT)**2)
            enddo
            cs2bot=gamma_m1*f(l1,4,n1,ilnTT)
            cs2top=gamma_m1*f(l1,4,n2,ilnTT)
!
          case ('rad_equil')
            call rad_equil(f)
!
          case ('blob_hs')
            if (lroot) print*, 'init_lnTT: hydrostatic blob with ', &
                radius_lnTT, ampl_lnTT(j), center1_x, center1_y, center1_z
            call blob(ampl_lnTT(j),f,ilnTT,radius_lnTT,center1_x, center1_y,center1_z)
            call blob(-ampl_lnTT(j),f,ilnrho,radius_lnTT,center1_x,center1_y,center1_z)
!
          case ('blob')
            if (lroot) print*, 'init_lnTT: blob ', &
                radius_lnTT, ampl_lnTT(j), center1_x, center1_y, center1_z
            call blob(ampl_lnTT(j),f,ilnTT,radius_lnTT,center1_x,center1_y,center1_z)
!
          case ('isothermal')
            if (lroot) print*, 'init_lnTT: isothermal atmosphere'
            if (ltemperature_nolog) then
              f(:,:,:,iTT)  =cs20/gamma_m1
            else
              f(:,:,:,ilnTT)=log(cs20/gamma_m1)
            endif
            haut=-cs20/gamma/gravz
            ztop=xyz1(3)
            if (ldensity_nolog) then
              do n=n1,n2
                f(:,:,n,irho)=rho0*exp((ztop-z(n))/haut)
              enddo
            else
              do n=n1,n2
                f(:,:,n,ilnrho)=lnrho0+(ztop-z(n))/haut
              enddo
            endif
!
          case ('hydro_rad')
            if (lroot) print*, 'init_lnTT: hydrostatic+radiative equilibria'
            if (Fbot==impossible .or. hcond0==impossible) &
                call fatal_error("initialize_lnTT","Fbot or hcond0 not initialized")
            Rgas=(1.-1./gamma)/cp1
            Ttop=cs20/gamma_m1
            beta=-Fbot/hcond0
            alpha=Ttop-beta
            expo=-gravz/beta/Rgas
            do n=n1,n2
              if (ltemperature_nolog) then
                f(:,:,n,iTT)  =beta*z(n)+alpha
              else
                f(:,:,n,ilnTT)=log(beta*z(n)+alpha)
              endif
              f(:,:,n,ilnrho)=lnrho0+(1.+expo)*log((1.+alpha/beta)/(z(n)+alpha/beta))
            enddo
!
          case ('star_heat')
            call star_heat(f)
!
          case default
!
!  Catch unknown values.
!
            call fatal_error('init_energy','no such initlnTT('//trim(iinit_str)//'): '//trim(initlnTT(j)))
!
          endselect
!
          if (lroot) print*,'init_energy: initlnTT('//trim(iinit_str)//') = ',trim(initlnTT(j))
        endif
      enddo
!
!  Interface for user's own initial condition.
!
      if (linitial_condition) call initial_condition_ss(f)
!
      if (lnothing.and.lroot) print*,'init_energy: nothing'
!
      if (ltemperature_nolog.and.linitial_log) f(:,:,:,iTT)=exp(f(:,:,:,ilnTT))
!
    endsubroutine init_energy
!***********************************************************************
    subroutine pencil_criteria_energy
!
!  All pencils that the Energy module depends on are specified here.
!
!  20-11-04/anders: coded
!
      if (lwrite_slices.and.ivid_pp/=0) lpenc_video(i_pp)=.true.
!
!  cs2 affects time step only if the continuity equation is being solved,
!  i.e. not for boussinesq or anelastic.
!
      if (ldt.and.ldensity) lpenc_requested(i_cs2)=.true.
!
      if (lpressuregradient_gas) lpenc_requested(i_fpres)=.true.
      if (lparticles_temperature) lpenc_requested(i_tcond)=.true.
!
      if (lviscosity.and.lviscosity_heat) then
        lpenc_requested(i_cv1)=.true.
        lpenc_requested(i_visc_heat)=.true.
        if (.not.ltemperature_nolog) lpenc_requested(i_TT1)=.true.
      endif
!
      if (ldensity.or.lanelastic) then
        lpenc_requested(i_divu)=.true.
        if (ltemperature_nolog) lpenc_requested(i_TT)=.true.
      endif
!
      if (lcalc_heat_cool) then
        if (ldensity.or.lboussinesq.or.lanelastic) then
          lpenc_requested(i_rho1)=.true.
          lpenc_requested(i_TT)=.true.
          lpenc_requested(i_TT1)=.true.
          lpenc_requested(i_cv1)=.true.
        endif
        if (lgravr) then
          lpenc_requested(i_r_mn)=.true.
          if (lboussinesq) lpenc_requested(i_evr)=.true.
        endif
      endif
!
      if (lheatc_chiconst.or.lheatc_chicubicstep.or.lheatc_chi_reduce_ddr) then
        if (ltemperature_nolog) then
          lpenc_requested(i_del2TT)=.true.
          lpenc_requested(i_gTT)=.true.
        else
          lpenc_requested(i_del2lnTT)=.true.
          lpenc_requested(i_glnTT)=.true.
        endif
        lpenc_requested(i_glnrho)=.true.
        lpenc_requested(i_cp1)=.true.
      endif
      if (lheatc_chicubicstep) lpenc_requested(i_z_mn)=.true.
!
      if (lheatc_Kconst) then
        if (ldensity.or.lboussinesq.or.lanelastic) lpenc_requested(i_rho1)=.true.
        lpenc_requested(i_cp1)=.true.
        if (ltemperature_nolog) then
          lpenc_requested(i_del2TT)=.true.
        else
          lpenc_requested(i_glnTT)=.true.
          lpenc_requested(i_del2lnTT)=.true.
        endif
      endif
!
      if (lheatc_Ktherm) then
        if (ldensity) lpenc_requested(i_rho1)=.true.
        lpenc_requested(i_cp1)=.true.
        if (ltemperature_nolog) then
          lpenc_requested(i_del2TT)=.true.
        else
          lpenc_requested(i_glnTT)=.true.
          lpenc_requested(i_del2lnTT)=.true.
        endif
      endif
!
      if (lheatc_Kprof) then
        lpenc_requested(i_rho1)=.true.
        lpenc_requested(i_cp1)=.true.
        if (ltemperature_nolog) then
          lpenc_requested(i_gTT)=.true.
          lpenc_requested(i_del2TT)=.true.
        else
          lpenc_requested(i_glnTT)=.true.
          lpenc_requested(i_del2lnTT)=.true.
        endif
        if (lgravz) lpenc_requested(i_z_mn)=.true.
        if (lgravr) lpenc_requested(i_r_mn)=.true.
      endif
!
      if (lheatc_Karctan) then
        lpenc_requested(i_rho1)=.true.
        lpenc_requested(i_cp1)=.true.
        lpenc_requested(i_TT)=.true.
        lpenc_requested(i_gTT)=.true.
        if (.not. lADI_mixed) lpenc_requested(i_del2TT)=.true.
      endif
!
      if (lheatc_shock) then
        lpenc_requested(i_glnrho)=.true.
        lpenc_requested(i_shock)=.true.
        lpenc_requested(i_del2lnrho)=.true.
        lpenc_requested(i_gshock)=.true.
        if (ltemperature_nolog) then
          lpenc_requested(i_gTT)=.true.
          lpenc_requested(i_del2TT)=.true.
        else
          lpenc_requested(i_glnTT)=.true.
          lpenc_requested(i_del2lnTT)=.true.
        endif
      endif
!
      if (lheatc_tensordiffusion) then
        lpenc_requested(i_bb)=.true.
        lpenc_requested(i_bij)=.true.
        lpenc_requested(i_rho1)=.true.
        lpenc_requested(i_glnTT)=.true.
        lpenc_requested(i_hlnTT)=.true.
        lpenc_requested(i_cp1)=.true.
      endif
!
      if (lheatc_hyper3) then
        if (ltemperature_nolog) then
          lpenc_requested(i_del6TT)=.true.
        else
          lpenc_requested(i_del6lnTT)=.true.
        endif
      endif
!
      if (lheatc_hyper3_mesh) lpenc_requested(i_TT1)=.true.
!
      if (lheatc_kramers) then

        !lpenc_requested(i_rho)=.true.
        !lpenc_requested(i_lnrho)=.true.
        !lpenc_requested(i_lnTT)=.true.
        lpenc_requested(i_cp1)=.true.
        lpenc_requested(i_rho1)=.true.
        lpenc_requested(i_TT)=.true.
        lpenc_requested(i_glnrho)=.true.
        lpenc_requested(i_glnTT)=.true.
        !lpenc_requested(i_cv1)=.true.
        lpenc_requested(i_del2lnTT)=.true.
      endif
!
      if (ladvection_temperature) then
        if (ltemperature_nolog) then
          lpenc_requested(i_ugTT)=.true.
          lpenc_requested(i_del2TT)=.true.
        else
          lpenc_requested(i_uglnTT)=.true.
        endif
      endif
!
!  Diagnostic pencils.
!
      if ( idiag_TTmax/=0.or.idiag_TTmin/=0 .or.idiag_TTm/=0  .or.idiag_TugTm/=0 .or. &
           idiag_Trms/=0 .or.idiag_uxTm/=0  .or.idiag_uyTm/=0 .or.idiag_uzTm/=0  .or. &
           idiag_TT2mz/=0 .or.idiag_uxTmz/=0.or.idiag_uyTmz/=0.or.idiag_uzTmz/=0 .or. &
           idiag_TT2m/=0 .or.idiag_TTzmask/=0 ) lpenc_diagnos(i_TT)=.true.
!
      if ( idiag_TugTm/=0 .or. idiag_gT2m/=0 .or. &
           idiag_guxgTm/=0 .or. idiag_guygTm/=0 .or. idiag_guzgTm/=0 ) lpenc_diagnos(i_gTT)=.true.
!
      if ( lhydro .or. lhydro_kinematic ) then
        if ( idiag_guxgTm/=0 .or. idiag_guygTm/=0 .or. idiag_guzgTm/=0 ) lpenc_diagnos(i_uij)=.true.
        if ( idiag_Tugux_uxugTm/=0 .or. idiag_Tuguy_uyugTm/=0 .or. idiag_Tuguz_uzugTm/=0 ) then
          lpenc_diagnos(i_ugu)=.true.; lpenc_diagnos(i_ugTT)=.true.
        endif
      endif
!
      if ( idiag_Tdxpm/=0 .or. idiag_Tdypm/=0 .or. idiag_Tdzpm/=0 ) then
        lpenc_diagnos(i_fpres)=.true.
        lpenc_diagnos(i_TT)=.true.
      endif
      if (idiag_gTmax/=0) then
         lpenc_diagnos(i_glnTT) =.true.
         lpenc_diagnos(i_TT) =.true.
      endif
      if (idiag_fradtop/=0.or.idiag_fradbot/=0.or.idiag_fradmz/=0) then
        lpenc_diagnos(i_TT) =.true.
        lpenc_diagnos(i_glnTT) =.true.
      endif
      if (idiag_fconvmz/=0) then
        lpenc_diagnos(i_TT) =.true.
        lpenc_diagnos(i_cp)=.true.
        lpenc_diagnos(i_uu)=.true.
        lpenc_diagnos(i_rho)=.true.
      endif
      if (idiag_yHmax/=0) lpenc_diagnos(i_yH)  =.true.
      if (idiag_yHmin/=0) lpenc_diagnos(i_yH)  =.true.
      if (idiag_yHm/=0)   lpenc_diagnos(i_yH)  =.true.
      if (idiag_ethm/=0.or.idiag_ethmz/=0.or.idiag_ethtot/=0) then
        lpenc_diagnos(i_rho)=.true.
        lpenc_diagnos(i_ee)  =.true.
      endif
      if (idiag_ethuxmz/=0.or.idiag_ethuymz/=0.or.idiag_ethuzmz/=0.or.&
          idiag_ethuxmx/=0) then
        lpenc_diagnos(i_rho)=.true.
        lpenc_diagnos(i_ee) =.true.
        lpenc_diagnos(i_uu) =.true.
      endif
      if (idiag_uxTmz/=0.or.idiag_uyTmz/=0.or.idiag_uzTmz/=0) &
        lpenc_diagnos(i_uu) =.true.
      if (idiag_ssm/=0)    lpenc_diagnos(i_ss)  =.true.
      if (idiag_dtchi/=0)  lpenc_diagnos(i_cs2)=.true.
      if (idiag_csm/=0 .or. idiag_csmax/=0)   lpenc_diagnos(i_cs2)=.true.
      if (idiag_eem/=0)    lpenc_diagnos(i_ee) =.true.
      if (idiag_ppm/=0 .or. idiag_ppmx/=0 .or. idiag_ppmy/=0 .or. &
          idiag_ppmz/=0 .or. idiag_ppuzmz/=0) lpenc_diagnos(i_pp) =.true.
      if (idiag_thcool/=0) lpenc_diagnos(i_rho)=.true.
      if (idiag_TTmx/=0 .or. idiag_TTmy/=0 .or. idiag_TTmz/=0) &
          lpenc_diagnos(i_TT)=.true.
      if (idiag_dtchi/=0) then
        lpenc_diagnos(i_rho1)=.true.
        lpenc_diagnos(i_cv1) =.true.
      endif
      if (idiag_fpresxmz/=0 .or. idiag_fpresymz/=0 .or. &
          idiag_fpreszmz/=0) lpenc_diagnos(i_fpres)=.true.
!
      if (idiag_TTmxy/=0 .or. idiag_TTmxz/=0) lpenc_diagnos2d(i_TT)=.true.
!
      if (idiag_EmAIA94mxy/=0 .or. idiag_EmAIA94mxz/=0 .or.  &
          idiag_EmAIA131mxy/=0 .or. idiag_EmAIA131mxz/=0 .or. &
          idiag_EmAIA171mxy/=0 .or. idiag_EmAIA171mxz/=0 .or. &
          idiag_EmAIA193mxy/=0 .or. idiag_EmAIA193mxz/=0 .or. &
          idiag_EmAIA211mxy/=0 .or. idiag_EmAIA211mxz/=0 .or. &
          idiag_EmAIA304mxy/=0 .or. idiag_EmAIA304mxz/=0 .or. &
          idiag_EmAIA335mxy/=0 .or. idiag_EmAIA335mxz/=0 .or. &
          idiag_EmXRTmxy/=0 .or. idiag_EmXRTmxz/=0 ) then
        lpenc_diagnos2d(i_TT)=.true.
        lpenc_diagnos2d(i_rho)=.true.
      endif
!
    endsubroutine pencil_criteria_energy
!***********************************************************************
    subroutine pencil_interdep_energy(lpencil_in)
!
!  Interdependency among pencils from the Energy module is specified here.
!
!  20-11-04/anders: coded
!
      logical, dimension(npencils) :: lpencil_in
!
      if (lpencil_in(i_Ma2)) then
        lpencil_in(i_u2)=.true.
        lpencil_in(i_cs2)=.true.
      endif
      if (lpencil_in(i_uglnTT)) then
        lpencil_in(i_glnTT)=.true.
        lpencil_in(i_uu)=.true.
      endif
      if (lpencil_in(i_ugTT)) then
        lpencil_in(i_gTT)  =.true.
        lpencil_in(i_uu)=.true.
      endif
!
      if (lpencil_in(i_fpres).and. .not.lboussinesq) then
        lpencil_in(i_cs2)=.true.
        lpencil_in(i_glnrho)=.true.
        lpencil_in(i_glnTT)=.true.
        lpencil_in(i_glnmumol)=.true.
      endif

      if (lpencil_in(i_glnTT).and.ltemperature_nolog.and..not.leos) &
        lpencil_in(i_TT1) =.true.
!
    endsubroutine pencil_interdep_energy
!***********************************************************************
    subroutine calc_pencils_energy(f,p)
!
!  Calculate Energy pencils.
!  Most basic pencils should come first, as others may depend on them.
!
!  20-11-04/anders: coded
!  31-01-18/MR: made calculation of p%gTT corrrect also for log temperature
!
      use Sub, only: u_dot_grad,grad,multmv,del2
      use Deriv, only: der2
!
      real, dimension (mx,my,mz,mfarray), intent (in) :: f
      type (pencil_case), intent (inout) :: p
      integer :: j
      real, dimension(nx,3) :: gpp
      real, dimension(nx) :: temp
!
      if (.not.leos) then
        if (lpencil(i_lnTT)) then
          if (ltemperature_nolog) then
            p%lnTT=alog(f(l1:l2,m,n,iTT))
          else
            p%lnTT=f(l1:l2,m,n,ilnTT)
          endif
        endif
        if (lpencil(i_TT)) then
          if (ltemperature_nolog) then
            p%TT=f(l1:l2,m,n,iTT)
          else
            p%TT=exp(f(l1:l2,m,n,ilnTT))
          endif
        endif
        if (lpencil(i_TT1)) then
          if (ltemperature_nolog) then
            p%TT1=1./f(l1:l2,m,n,iTT)
          else
            p%TT1=exp(-f(l1:l2,m,n,ilnTT))
          endif
        endif
        if (lpencil(i_glnTT)) then
          call grad(f,ilnTT,p%glnTT)
          if (ltemperature_nolog) then
            do j=1,3
              p%glnTT(:,j) = p%glnTT(:,j)*p%TT1
            enddo
          endif
        endif

        if (ltemperature_nolog) then
          if (lpencil(i_del2TT)) call del2(f,iTT,p%del2TT)
        else
          call not_implemented('calc_pencils_energy','del2TT for logarithmic temperature and noeos')
        endif

        if (lpencil(i_hlnTT).or.lpencil(i_del2lnTT).or.lpencil(i_del6TT).or.lpencil(i_del6lnTT)) &
          call not_implemented('calc_pencils_energy','hlnTT, del[26]TT, del[26]lnTT for noeos')

      endif
! Ma2
      if (lpencil(i_Ma2)) p%Ma2=p%u2/p%cs2
! uglnTT
      if (lpencil(i_uglnTT)) call u_dot_grad(f,ilnTT,p%glnTT,p%uu,p%uglnTT,UPWIND=lupw_lnTT)
! ugTT
      if (lpencil(i_ugTT)) call u_dot_grad(f,ilnTT,p%gTT,p%uu,p%ugTT,UPWIND=lupw_lnTT)
! glnTT
      if (lpencil(i_gTT)) then
        call grad(f,ilnTT,p%gTT)
        if (.not.ltemperature_nolog) then
          temp=exp(f(l1:l2,m,n,ilnTT))
          do j=1,3
            p%gTT(:,j) = p%gTT(:,j)*temp
          enddo
        endif
        do j=1,3
          if (gradTT0(j)/=0.) p%gTT(:,j)=p%gTT(:,j)+gradTT0(j)
        enddo
      endif

      if (lpencil(i_d2xlnTT)) call der2(f,ilnTT,p%d2xlnTT,1)
      if (lpencil(i_d2zlnTT)) call der2(f,ilnTT,p%d2zlnTT,3)
!
! fpres
      if (lpencil(i_fpres)) then
        if (lboussinesq) then
!
! subroutine boussinesq has already been called at this instant,
! so f(:,:,:,ipp) contains div(uu), which is in fact the change of div(uu) due to the last
! timestep, hence p.dt
!
            call grad(f,ipp,gpp)
            if (dt==0.) then
              p%fpres=gpp
            else
              p%fpres=gpp/dt
            endif
        else
          do j=1,3
            p%fpres(:,j)=-gamma1*p%cs2*(p%glnrho(:,j)+p%glnTT(:,j)-p%glnmumol(:,j))
          enddo
        endif
      endif
! tcond
      if (lpencil(i_tcond)) then
        if (lheatc_chiconst) then
          p%tcond=chi*p%rho/p%cp1
        elseif (lheatc_Kconst) then
          p%tcond=hcond0
        else
          call not_implemented('calc_pencils_energy','this heatcond with lpencil(i_tcond)')
        endif
      endif
! sglnTT 
      if (lpencil(i_sglnTT)) call multmv(p%sij,p%glnTT,p%sglnTT)
!
!  ``cs2/dx^2'' for timestep
!
      if (lupdate_courant_dt.and.ldensity.and.lhydro) then
        if (lreduced_sound_speed) then
!          if (lscale_to_cs2top) then
!            p%advec_cs2=reduce_cs2*cs2top*dxyz_2
!          else
            p%advec_cs2=reduce_cs2*p%cs2*dxyz_2
!          endif
        else
          p%advec_cs2=p%cs2*dxyz_2
          advec_cs2=p%advec_cs2
        endif
        if (headtt.or.ldebug) print*,'calc_pencils_energy: max(advec_cs2) =',maxval(p%advec_cs2)
      endif
!
    endsubroutine calc_pencils_energy
!***********************************************************************
    subroutine denergy_dt(f,df,p)
!
!  Calculate right hand side of temperature equation.
!
!  lnTT version: DlnTT/Dt = -gamma_m1*divu + gamma*cp1*rho1*TT1*RHS
!    TT version:   DTT/Dt = -gamma_m1*TT*divu + gamma*cp1*rho1*RHS
!
!  13-dec-02/axel+tobi: adapted from energy
!  18-may-12/MR: compression work as heat sink added for boussinesq
!
      use Deriv, only: der6
      use EquationOfState, only: lpres_grad
      use Special, only: special_calc_energy
      use Sub, only: identify_bcs, calc_slope_diff_flux
      use Viscosity, only: calc_viscous_heat
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
!
      real, dimension (nx) :: Hmax=0.0, thdiff, tmp, advec_hypermesh_ss
      integer :: j
!
      intent(inout) :: f,p,df
!
! When initializating  thdiff in the declaration the
! variable gets the SAVE attribute,
! so in the next call thdiff is not initialized anymore.
!
      thdiff = 0.0
!
!  Identify module and boundary conditions.
!
      if (headtt.or.ldebug) print*, 'SOLVE dlnTT_dt'
      if (headtt) then
        if (ltemperature_nolog) then
          print*, 'denergy_dt: TT,cs2=', p%TT(1), p%cs2(1)
          call identify_bcs('TT',iTT)
        else
          print*, 'denergy_dt: lnTT,cs2=', p%lnTT(1), p%cs2(1)
          call identify_bcs('lnTT',ilnTT)
        endif
      endif
!
!  Sound speed squared.
!
      if (headtt) print*, 'denergy_dt: cs20=', p%cs2(1)
!
!  Add pressure gradient term in momentum equation.
!
      if (lpressuregradient_gas) df(l1:l2,m,n,iux:iuz) = df(l1:l2,m,n,iux:iuz) + p%fpres
!
!  Advection term and PdV-work.
!
      if (ladvection_temperature) then
        if (ltemperature_nolog) then
          df(l1:l2,m,n,iTT)   = df(l1:l2,m,n,iTT)   - p%ugTT
        else
          df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) - p%uglnTT
        endif
      endif
!
!  Add divu term.
!
      if (ldensity.or.lanelastic) then
        if (ltemperature_nolog) then
          df(l1:l2,m,n,iTT)   = df(l1:l2,m,n,iTT)   - gamma_m1*p%TT*p%divu
        else
          df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) - gamma_m1*p%divu
        endif
!
! if viscous heating is allowed, compression work must be taken into account, too
!
      elseif (lboussinesq.and.lviscosity.and.lviscosity_heat) then
        if (ltemperature_nolog) then
          df(l1:l2,m,n,iTT)   = df(l1:l2,m,n,iTT)   - p%cv1*PrRa*f(l1:l2,m,n,iTT)*p%uu(:,3)
        else
          df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) - p%cv1*PrRa*p%uu(:,3)
        endif
      endif
!
!  Calculate viscous contribution to temperature.
!
      if (lviscosity.and.lviscosity_heat) call calc_viscous_heat(df,p,Hmax)
!
!  Various heating conduction contributions.
!
      if (lcalc_heat_cool) call calc_heat_cool(f,df,p)
!
!  Thermal conduction
!
      diffus_chi=0.; diffus_chi3=0.
      if (lheatc_chiconst) call calc_heatcond_constchi(df,p)
      if (lheatc_chi_reduce_ddr) call calc_heatcond_chi_no_ddrsph(f,df,p)
      if (lheatc_chicubicstep) call calc_heatcond_cubicstepchi(df,p)
      if (lheatc_Kconst)   call calc_heatcond_constK(df,p)
      if (lheatc_Kprof)    call calc_heatcond(f,df,p)
      if (lheatc_Karctan)  call calc_heatcond_arctan(df,p)
      if (lheatc_tensordiffusion) call calc_heatcond_tensor(df,p)
      if (lheatc_Ktherm) call calc_heatcond_Ktherm (df,p)
      if (lheatc_kramers) call calc_heatcond_kramers (df,p)
!
!     Slope-limited diffusion
!
      if (lenergy_slope_limited.and.llast) then
        if (ltemperature_nolog) then
          call calc_slope_diff_flux(f,iTT,p,h_sld_ene,nlf_sld_ene,tmp,div_sld_ene)
          thdiff=thdiff+tmp
        else
          call calc_slope_diff_flux(f,ilnTT,p,h_sld_ene,nlf_sld_ene,tmp,div_sld_ene)
          thdiff=thdiff+tmp*p%TT1
       endif
     endif
!
!  Hyper diffusion.
!
      if (lheatc_hyper3) then
        if (ltemperature_nolog) then
          thdiff=thdiff+chi_hyper3*p%del6TT
        else
          thdiff=thdiff+chi_hyper3*p%del6lnTT
        endif
        if (lupdate_courant_dt) diffus_chi3=diffus_chi3+chi_hyper3*dxyz_6
        if (headtt) print*,'denergy_dt: chi_hyper3=', chi_hyper3
      endif
!
      if (lheatc_hyper3_mesh) then
        do j=1,3
          call der6(f,ilnTT,tmp,j,IGNOREDX=.true.)
          if (.not.ltemperature_nolog) tmp=tmp*p%TT1
          thdiff = thdiff + chi_hyper3_mesh*pi5_1/60.*tmp*dline_1(:,j)
        enddo
        if (lupdate_courant_dt) then
          advec_hypermesh_ss=chi_hyper3_mesh*pi5_1*sqrt(dxyz_2)
          advec2_hypermesh=advec2_hypermesh+advec_hypermesh_ss**2
        endif
        if (headtt) print*,'denergy_dt: chi_hyper3_mesh=', chi_hyper3_mesh
      endif
!
      if (lheatc_hyper3_polar) then
        do j=1,3
          call der6(f,ilnTT,tmp,j,IGNOREDX=.true.)
          if (.not.ltemperature_nolog) tmp=tmp*p%TT1
          thdiff = thdiff + chi_hyper3*pi4_1*tmp*dline_1(:,j)**2
        enddo
        if (lupdate_courant_dt) diffus_chi3=diffus_chi3+chi_hyper3*pi4_1*dxmin_pencil**4
        if (headtt) print*,'denergy_dt: chi_hyper3=', chi_hyper3
      endif
!
!  Shock diffusion.
!
      if (lheatc_shock) call calc_heatcond_shock(df,p)
!
!  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)
!
!  Add thermal diffusion to temperature equation.
!
      if (ltemperature_nolog) then
        df(l1:l2,m,n,iTT)   = df(l1:l2,m,n,iTT)   + thdiff
      else
        df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + thdiff
      endif
!
!  Boussinesq approximation: - u.grad T_0 added
!
      if (lboussinesq) then
        if (lsphere_in_a_box) then
          df(l1:l2,m,n,iTT) = df(l1:l2,m,n,iTT) - p%r_mn*beta_bouss*( &
                               f(l1:l2,m,n,iux)*p%evr(:,1)+f(l1:l2,m,n,iuy)*p%evr(:,2)+ &
                               f(l1:l2,m,n,iuz)*p%evr(:,3))
        else
!
! background temperature gradient in z direction
!
          df(l1:l2,m,n,iTT) = df(l1:l2,m,n,iTT) - beta_bouss*f(l1:l2,m,n,iuz)
        endif
      endif
!
! Store  pressure gradient in f-array if requested
!
      if (lpres_grad) then
        f(l1:l2,m,n,igpx) = -p%fpres(:,1)*f(l1:l2,m,n,irho)
        f(l1:l2,m,n,igpy) = -p%fpres(:,2)*f(l1:l2,m,n,irho)
      endif
!
!  Information on the timescales.
!
      if (lupdate_courant_dt) then

        maxdiffus=max(maxdiffus,diffus_chi)
        maxdiffus3=max(maxdiffus3,diffus_chi3)
 
        if (headtt.or.ldebug) then
          print*, 'denergy_dt: max(diffus_chi ) =', maxval(diffus_chi)
          print*, 'denergy_dt: max(diffus_chi3) =', maxval(diffus_chi3)
        endif
      endif
!
!  Apply border profile
!
      if (lborder_profiles) call set_border_entropy(f,df,p)

      call calc_diagnostics_energy(f,p)
!
    endsubroutine denergy_dt
!***********************************************************************
    subroutine calc_diagnostics_energy(f,p)

      use Slices_methods, only: store_slices

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

      call keep_compiler_quiet(f)

      call calc_2d_diagnostics_energy(p)
      call calc_1d_diagnostics_energy(p)
      call calc_0d_diagnostics_energy(p)

      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 calc_2d_diagnostics_energy(p)
!
!  2-D averages.
!
      use Diagnostics

      type(pencil_case) :: p      
!
      real, dimension (nx) :: NN2_resp, logT_resp, resp, resp2

      if (l2davgfirst) then

        call zsum_mn_name_xy(p%TT,idiag_TTmxy)
        call ysum_mn_name_xz(p%TT,idiag_TTmxz)
!
!   calculating emission with tabulated response files
!
        if (idiag_EmAIA94mxy/=0  .or. idiag_EmAIA94mxz/=0  .or. &
            idiag_EmAIA131mxy/=0 .or. idiag_EmAIA131mxz/=0 .or. &
            idiag_EmAIA171mxy/=0 .or. idiag_EmAIA171mxz/=0 .or. &
            idiag_EmAIA193mxy/=0 .or. idiag_EmAIA193mxz/=0 .or. &
            idiag_EmAIA211mxy/=0 .or. idiag_EmAIA211mxz/=0 .or. &
            idiag_EmAIA304mxy/=0 .or. idiag_EmAIA304mxz/=0 .or. &
            idiag_EmAIA335mxy/=0 .or. idiag_EmAIA335mxz/=0 .or. &
            idiag_EmXRTmxy/=0    .or. idiag_EmXRTmxz/=0        ) then
!
! particle density squared in SI
!
          NN2_resp = (6./7.*p%rho/m_p*unit_density/unit_mass)**2.
!
! log10 temperature in SI (K)
!
          logT_resp =log10(p%TT*unit_temperature)
!
          if (idiag_EmAIA94mxy/=0 .or. idiag_EmAIA94mxz/=0) then
            call get_AIA_tab_resp('094',logT_resp,NN2_resp,resp)
            resp = resp*unit_length/1e6
            if (idiag_EmAIA94mxy/=0) then
              resp2 = resp*Lz
              call zsum_mn_name_xy(resp2,idiag_EmAIA94mxy)
            endif
            if (idiag_EmAIA94mxz/=0) then
              resp2 = resp*Ly
              call ysum_mn_name_xz(resp2,idiag_EmAIA94mxz)
            endif
          endif
!
          if (idiag_EmAIA131mxy/=0 .or. idiag_EmAIA131mxz/=0) then
            call get_AIA_tab_resp('131',logT_resp,NN2_resp,resp)
            resp = resp*unit_length/1e6
            if (idiag_EmAIA131mxy/=0) then
              resp2 = resp*Lz
              call zsum_mn_name_xy(resp2,idiag_EmAIA131mxy)
            endif
            if (idiag_EmAIA131mxz/=0) then
              resp2 = resp*Ly
              call ysum_mn_name_xz(resp2,idiag_EmAIA131mxz)
            endif
          endif
!
          if (idiag_EmAIA171mxy/=0 .or. idiag_EmAIA171mxz/=0) then
            call get_AIA_tab_resp('171',logT_resp,NN2_resp,resp)
            resp = resp*unit_length/1e6
            if (idiag_EmAIA171mxy/=0) then
              resp2 = resp*Lz
              call zsum_mn_name_xy(resp2,idiag_EmAIA171mxy)
            endif
            if (idiag_EmAIA171mxz/=0) then
              resp2 = resp*Ly
              call ysum_mn_name_xz(resp2,idiag_EmAIA171mxz)
            endif
          endif
!
          if (idiag_EmAIA193mxy/=0 .or. idiag_EmAIA193mxz/=0) then
            call get_AIA_tab_resp('193',logT_resp,NN2_resp,resp)
            resp = resp*unit_length/1e6
            if (idiag_EmAIA193mxy/=0) then
              resp2 = resp*Lz
              call zsum_mn_name_xy(resp2,idiag_EmAIA193mxy)
            endif
            if (idiag_EmAIA193mxz/=0) then
              resp2 = resp*Ly
              call ysum_mn_name_xz(resp2,idiag_EmAIA193mxz)
            endif
          endif
!
          if (idiag_EmAIA211mxy/=0 .or. idiag_EmAIA211mxz/=0) then
            call get_AIA_tab_resp('211',logT_resp,NN2_resp,resp)
              resp = resp*unit_length/1e6
            if (idiag_EmAIA211mxy/=0) then
              resp2 = resp*Lz
              call zsum_mn_name_xy(resp2,idiag_EmAIA211mxy)
            endif
            if (idiag_EmAIA211mxz/=0) then
              resp2 = resp*Ly
              call ysum_mn_name_xz(resp2,idiag_EmAIA211mxz)
            endif
          endif
!
          if (idiag_EmAIA304mxy/=0 .or. idiag_EmAIA304mxz/=0) then
            call get_AIA_tab_resp('304',logT_resp,NN2_resp,resp)
            resp = resp*unit_length/1e6
            if (idiag_EmAIA304mxy/=0) then
              resp2 = resp*Lz
              call zsum_mn_name_xy(resp2,idiag_EmAIA304mxy)
            endif
            if (idiag_EmAIA304mxz/=0) then
              resp2 = resp*Ly
              call ysum_mn_name_xz(resp2,idiag_EmAIA304mxz)
            endif
          endif
!
          if (idiag_EmAIA335mxy/=0 .or. idiag_EmAIA335mxz/=0) then
            call get_AIA_tab_resp('335',logT_resp,NN2_resp,resp)
            resp = resp*unit_length/1e6
            if (idiag_EmAIA335mxy/=0) then
              resp2 = resp*Lz
              call zsum_mn_name_xy(resp2,idiag_EmAIA335mxy)
            endif
            if (idiag_EmAIA335mxz/=0) then
              resp2 = resp*Ly
              call ysum_mn_name_xz(resp2,idiag_EmAIA335mxz)
            endif
          endif
!
          if (idiag_EmXRTmxy/=0 .or. idiag_EmXRTmxz/=0) then
            call get_AIA_tab_resp('XRT',logT_resp,NN2_resp,resp)
            resp = resp*unit_length/1e6
            if (idiag_EmXRTmxy/=0) then
              resp2 = resp*Lz
              call zsum_mn_name_xy(resp2,idiag_EmXRTmxy)
            endif
            if (idiag_EmXRTmxz/=0) then
              resp2 = resp*Ly
              call ysum_mn_name_xz(resp2,idiag_EmXRTmxz)
            endif
          endif
!
        endif

      endif
    endsubroutine calc_2d_diagnostics_energy
!***********************************************************************
    subroutine calc_0d_diagnostics_energy(p)
!
!  Calculate temperature related 0D-diagnostics.
!
      use Diagnostics
      use Sub, only: dot2, dot, dot_mn

      type(pencil_case) :: p

      real, dimension(nx) :: tmp
      real :: fradtop, fradbot

      if (ldiagnos) then
!
        call sum_mn_name(p%TT,idiag_TTm)
        if (idiag_TTzmask/=0) call sum_mn_name(p%TT*zmask_temp(n-n1+1),idiag_TTzmask)
        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%ss,idiag_ssm)
        call sum_mn_name(p%ee,idiag_eem)
        call sum_mn_name(p%pp,idiag_ppm)
        if (idiag_ethm/=0)  call sum_mn_name(p%rho*p%ee,idiag_ethm)
        if (idiag_ethtot/=0) call integrate_mn_name(p%rho*p%ee,idiag_ethtot)
        call sum_mn_name(p%cs2,idiag_csm,lsqrt=.true.)
        call max_mn_name(p%cs2,idiag_csmax,lsqrt=.true.)
        if (idiag_TugTm/=0) call sum_mn_name(p%TT*p%ugTT,idiag_TugTm)
        if (idiag_Trms/=0)  call sum_mn_name(p%TT**2,idiag_Trms,lsqrt=.true.)
        if (idiag_TT2m/=0)  call sum_mn_name(p%TT**2,idiag_TT2m)
        call sum_mn_name(p%uu(:,1)*p%TT,idiag_uxTm)
        call sum_mn_name(p%uu(:,2)*p%TT,idiag_uyTm)
        call sum_mn_name(p%uu(:,3)*p%TT,idiag_uzTm)
        if (idiag_Tugux_uxugTm/=0) call sum_mn_name(p%TT*p%ugu(:,1)+p%uu(:,1)*p%ugTT,idiag_Tugux_uxugTm)
        if (idiag_Tuguy_uyugTm/=0) call sum_mn_name(p%TT*p%ugu(:,2)+p%uu(:,2)*p%ugTT,idiag_Tuguy_uyugTm)
        if (idiag_Tuguz_uzugTm/=0) call sum_mn_name(p%TT*p%ugu(:,3)+p%uu(:,3)*p%ugTT,idiag_Tuguz_uzugTm)
        if (idiag_Tdxpm/=0)  call sum_mn_name(p%TT*p%fpres(:,1),idiag_Tdxpm)
        if (idiag_Tdypm/=0)  call sum_mn_name(p%TT*p%fpres(:,2),idiag_Tdypm)
        if (idiag_Tdzpm/=0)  call sum_mn_name(p%TT*p%fpres(:,3),idiag_Tdzpm)
!
        if (idiag_gT2m/=0) then
          call dot(p%gTT,p%gTT,tmp)
          call sum_mn_name(tmp,idiag_gT2m)
        endif
!
        if (idiag_guxgTm/=0) then
          call dot_mn(p%gTT(:,:),p%uij(:,1,:),tmp)
          call sum_mn_name(tmp,idiag_guxgTm)
        endif
!
        if (idiag_guygTm/=0) then
          call dot_mn(p%gTT(:,:),p%uij(:,2,:),tmp)
          call sum_mn_name(tmp,idiag_guygTm)
        endif
!
        if (idiag_guzgTm/=0) then
          call dot_mn(p%gTT(:,:),p%uij(:,3,:),tmp)
          call sum_mn_name(tmp,idiag_guzgTm)
        endif
!
        if (ldt.and.idiag_dtc/=0) call max_mn_name(sqrt(p%advec_cs2)/cdt,idiag_dtc,l_dt=.true.)

        if (idiag_gTmax/=0) then
          call dot2(p%glnTT,tmp)
          call max_mn_name(p%TT*sqrt(tmp),idiag_gTmax)
        endif
!
        if (idiag_fradtop/=0.and.llast_proc_z) then
          if (n==n2) then
            if (lADI) then
              call heatcond_TT(p%TT,hcond)
            else
              hcond=hcond0
            endif
            fradtop=sum(-hcond*p%TT*p%glnTT(:,3)*dsurfxy)
          endif
          call surf_mn_name(fradtop,idiag_fradtop,n2)
        endif
!
        if (idiag_fradbot/=0.and.lfirst_proc_z) then
          if (n==n1) then
            if (lADI) then
              call heatcond_TT(p%TT,hcond)
            else
              hcond=hcond0
            endif
            fradbot=sum(-hcond*p%TT*p%glnTT(:,3)*dsurfxy)
          endif
          call surf_mn_name(fradbot,idiag_fradbot,n1)
        endif
!
! Cooling power - energy radiated away (luminosity)
!
        if (lcalc_heat_cool .and. idiag_thcool/=0) call sum_lim_mn_name(cooling*p%rho,idiag_thcool,p)
!
        if (ldt) then
          if (idiag_dtchi/=0) call max_mn_name(diffus_chi/cdtv,idiag_dtchi,l_dt=.true.)
        endif
      endif

    endsubroutine calc_0d_diagnostics_energy
!***********************************************************************
    subroutine calc_1d_diagnostics_energy(p)
!
!  1-D averages.
!
      use Diagnostics

      type(pencil_case) :: p

      if (l1davgfirst) then

        if (idiag_fradmz/=0) call xysum_mn_name_z(-hcond0*p%TT*p%glnTT(:,3),idiag_fradmz)
        if (idiag_fconvmz/=0) call xysum_mn_name_z(p%cp*p%rho*p%uu(:,3)*p%TT,idiag_fconvmz)
        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)
        if (idiag_ppuzmz/=0) call xysum_mn_name_z(p%pp*p%uu(:,3),idiag_ppuzmz)
        if (idiag_ethmz/=0) call xysum_mn_name_z(p%rho*p%ee,idiag_ethmz)
        if (idiag_ethuxmx/=0) call yzsum_mn_name_x(p%rho*p%ee*p%uu(:,1),idiag_ethuxmx)
        if (idiag_ethuxmz/=0) call xysum_mn_name_z(p%rho*p%ee*p%uu(:,1),idiag_ethuxmz)
        if (idiag_ethuymz/=0) call xysum_mn_name_z(p%rho*p%ee*p%uu(:,2),idiag_ethuymz)
        if (idiag_ethuzmz/=0) call xysum_mn_name_z(p%rho*p%ee*p%uu(:,3),idiag_ethuzmz)
        call xysum_mn_name_z(p%fpres(:,1),idiag_fpresxmz)
        call xysum_mn_name_z(p%fpres(:,2),idiag_fpresymz)
        call xysum_mn_name_z(p%fpres(:,3),idiag_fpreszmz)
!
        if (idiag_TT2mz/=0) call xysum_mn_name_z(p%TT**2,idiag_TT2mz)
        if (idiag_uxTmz/=0) call xysum_mn_name_z(p%uu(:,1)*p%TT,idiag_uxTmz)
        if (idiag_uyTmz/=0) call xysum_mn_name_z(p%uu(:,2)*p%TT,idiag_uyTmz)
        if (idiag_uzTmz/=0) call xysum_mn_name_z(p%uu(:,3)*p%TT,idiag_uzTmz)
!
      endif

    endsubroutine calc_1d_diagnostics_energy
!***********************************************************************
    subroutine set_border_entropy(f,df,p)
!
!  Calculates the driving term for the border profile
!  of the ss variable.
!
!  28-jul-06/wlad: coded
!
      use BorderProfiles, only: border_driving, set_border_initcond
!
      real, dimension(mx,my,mz,mfarray) :: f
      type (pencil_case) :: p
      real, dimension(mx,my,mz,mvar) :: df
      real, dimension(nx) :: f_target
!
      select case (borderss)
!
      case ('zero','0')
        if (ltemperature_nolog) then
          f_target=0.0
        else
          f_target=1.0
        endif

        if (ltemperature_nolog) then
          call border_driving(f,df,p,f_target,iTT)
        else
          call border_driving(f,df,p,f_target,ilnTT)
        endif

      case ('constant')
        if (ltemperature_nolog) then
          f_target=TT_const
        else
          f_target=lnTT_const
        endif

        if (ltemperature_nolog) then
          call border_driving(f,df,p,f_target,iTT)
        else
          call border_driving(f,df,p,f_target,ilnTT)
        endif

      case ('initial-condition')
        if (ltemperature_nolog) then
          call set_border_initcond(f,iTT,f_target)
        else
          call set_border_initcond(f,ilnTT,f_target)
        endif

        if (ltemperature_nolog) then
          call border_driving(f,df,p,f_target,iTT)
        else
          call border_driving(f,df,p,f_target,ilnTT)
        endif

      case ('nothing')
      endselect
!
    endsubroutine set_border_entropy
!***********************************************************************
    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
!     if (lslope_limit_diff) then
       cs2=0.
       do m=1,my
       do n=1,mz
         if (ltemperature_nolog) then
           cs2 = gamma_m1/cp1*f(:,m,n,iTT)
         else
           cs2 = gamma_m1/cp1*exp(f(:,m,n,ilnTT))
         endif
         f(:,m,n,isld_char)=f(:,m,n,isld_char)+w_sldchar_ene*sqrt(cs2)
       enddo
       enddo
     endif
!
    endsubroutine energy_before_boundary
!***********************************************************************
    subroutine energy_after_boundary(f)
!
!  Calculation of mean quantities.
!
!  17-apr-15/MR: coded
!
      use Sub, only: grad, finalize_aver
      use EquationOfState, only : eoscalc
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      real :: fact
      real, dimension (nx,3):: gradTT
      real, dimension (nx)  :: temp
!
      integer :: nl
!
!  Calculate mean of temperature and its gradient.
!
      if (lcalc_TTmean) then
!
        fact=1./nxygrid
        do n=n1,n2
!
          nl = n-n1+1
          TTmz(nl)=0.; gTTmz(nl)=0.
!
          if (ltemperature_nolog) then
            do m=m1,m2
              TTmz(nl)=TTmz(nl)+sum(f(l1:l2,m,n,iTT))
              call grad(f,iTT,gradTT)
              gTTmz(nl)=gTTmz(nl)+sum(gradTT(:,3))
            enddo
          else
            do m=m1,m2
              temp = exp(f(l1:l2,m,n,ilnTT))
              TTmz(nl)=TTmz(nl)+sum(temp)
              call grad(f,ilnTT,gradTT)
              gTTmz(nl)=gTTmz(nl)+sum(gradTT(:,3)*temp)
            enddo
          endif
!
        enddo
!
        call finalize_aver(nprocxy,12,TTmz)
        call finalize_aver(nprocxy,12,gTTmz)
!
        TTmz  = fact*TTmz
        gTTmz = fact*gTTmz    ! simpler by deriving TTmz!!
!
      endif
!
    endsubroutine energy_after_boundary
!***********************************************************************
    subroutine calc_heatcond_shock(df,p)
!
!  Add shock diffusion to the energy equation,
!
!    De/Dt = ... + div(K*grad(T)) = ... + div(cv*T*Xi*grad(T)) ,
!
!  where e=rho*cv*T and Xi is a regular shock diffusion coefficient.
!
!  01-aug-08/wlad: adapted from entropy
!
      use Diagnostics, only: max_mn_name
      use Sub, only: dot
!
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
      real, dimension (nx) :: thdiff, g2, gshockgTT, gshockglnTT
!
      intent(in) :: p
      intent(inout) :: df
!
      if (headtt) print*, 'calc_heatcond_shock: chi_shock=', chi_shock
!
!  Shock energy diffusivity.
!
      if (ltemperature_nolog) then
        call dot(p%gshock,p%gTT,gshockgTT)
        call dot(p%glnrho,p%gTT,g2)
        thdiff=chi_shock*(p%shock*(p%del2TT+g2)+gshockgTT)
      else
        call dot(p%gshock,p%glnTT,gshockglnTT)
        call dot(p%glnrho+p%glnTT,p%glnTT,g2)
        thdiff=chi_shock*(p%shock*(p%del2lnTT+g2)+gshockglnTT)
      endif
!
      df(l1:l2,m,n,ilntt) = df(l1:l2,m,n,ilntt) + thdiff
!
      if (headtt) print*,'calc_heatcond_shock: added thdiff'
!
      if (lupdate_courant_dt) then
        if (leos_idealgas) then
          diffus_chi=diffus_chi+(gamma*chi_shock*p%shock)*dxyz_2
        else
          diffus_chi=diffus_chi+(chi_shock*p%shock)*dxyz_2
        endif
      endif
!
    endsubroutine calc_heatcond_shock
!***********************************************************************
    subroutine rad_equil(f)
!
!  Compute the radiative and hydrostatic equilibria for a given radiative
!  profile defined in heatcond_TT.
!
!  16-may-07/gastine+dintrans: coded
!
      use Gravity, only: gravz
      use EquationOfState, only: lnrho0,cs20,cs2top,cs2bot,eoscalc,ilnrho_TT
!
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f
      real, dimension (nzgrid) :: temp,lnrho
      real :: hcond, dtemp, dlnrho, ss
      integer :: i,n,iz
!
      if (.not. leos) call fatal_error('rad_equil', &
            'EOS=noeos, but radiative equilibrium requires an EQUATION OF STATE')

      if (.not. ltemperature_nolog) call not_implemented('temperature_idealgas','rad_equil for lnTT')
      if (lroot) print*,'init_energy: rad_equil for kappa-mechanism pb'
!
      if (nzgrid == 1) call not_implemented('rad_equil', "for nzgrid=1")
!
!  Integrate from top to bottom: z(n2) --> z(n1).
!
      temp(nzgrid)=cs20/gamma_m1
      lnrho(nzgrid)=lnrho0
!
!  Calculate the n2-1 gridpoint thanks to a 1st order forward Euler scheme.
!
      call heatcond_TT(temp(nzgrid), hcond)
      dtemp=Fbot/hcond
      temp(nzgrid-1)=temp(nzgrid)+dz*dtemp
      dlnrho=(-gamma/gamma_m1*gravz-dtemp)/temp(nzgrid)
      lnrho(nzgrid-1)=lnrho(nzgrid)+dz*dlnrho
!
!  Now we use a 2nd order centered scheme for the other gridpoints.
!
      do i=nzgrid-1,2,-1
        call heatcond_TT(temp(i), hcond)
        dtemp=Fbot/hcond
        temp(i-1)=temp(i+1)+2.*dz*dtemp
        dlnrho=(-gamma/gamma_m1*gravz-dtemp)/temp(i)
        lnrho(i-1)=lnrho(i+1)+2.*dz*dlnrho
      enddo
!
!  Fill in the density and temperature f-arrays and each z-processor
!  writes its own setup in the file data/proc#/setup.dat
!
      open(unit=11,file=trim(directory)//'/setup.dat')
      write(11,'(5a14)') 'z','rho','temp','ss','hcond'
      do n=1,nz
        iz=ipz*nz+n
        f(:,:,nghost+n,ilnrho)=lnrho(iz)
        f(:,:,nghost+n,ilnTT)=temp(iz)
        call eoscalc(ilnrho_TT,lnrho(iz),temp(iz),ss=ss)
        call heatcond_TT(temp(iz), hcond)
        write(11,'(5e14.5)') z(nghost+n),exp(lnrho(iz)),temp(iz),ss,hcond
      enddo
      close(11)
!
!  Initialize cs2bot by taking into account the new bottom value of temperature
!  Note: cs2top=cs20 already defined in eos_idealgas.
!
      cs2bot=gamma_m1*temp(1)
      print*,'cs2top, cs2bot=', cs2top, cs2bot
!
    endsubroutine rad_equil
!***********************************************************************
    subroutine calc_heat_cool(f,df,p)
!
      use Diagnostics, only: sum_lim_mn_name
      use EquationOfState, only: cs20
      use Sub, only: step
!
      real, dimension(mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
      real, dimension (nx) :: tau, kappa, a1, a3, prof, heat
      real :: a2, kappa0, kappa0_cgs
!
!  Initialize
!
      intent(in) :: p
      intent(inout) :: df
!
      if (.not. leos) call fatal_error('calc_heat_cool', &
            'EOS=noeos, but heatng/cooling requires an EQUATION OF STATE')

      if (headtt) print*,'enter calc_heat_cool', rcool, wcool, cool, cs20
!
      if (lgravr) then
        if (lboussinesq) then
          prof = step(p%r_mn,r_ext,wcool)
          heat = -cool*prof
          prof = 1.-step(p%r_mn,r_int,wcool)
          heat = heat-cool*prof
          df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + heat*f(l1:l2,m,n,ilnTT)
        else
          ! 2-D heating/cooling profiles
          prof = exp(-0.5*(p%r_mn/wheat)**2) * (2*pi*wheat**2)**(-1.)
          heat = luminosity*prof
          prof = step(p%r_mn,rcool,wcool)
          heat = heat - cool*prof*(p%cs2-cs20)/cs20
          df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + p%cv1*p%TT1*heat
        endif
      else
        kappa0_cgs=2e-6  !cm2/g
        kappa0=kappa0_cgs*unit_density*unit_length
        kappa=kappa0*p%TT**2
!
!  Optical Depth tau=kappa*rho*H.
!  If we are using 2D, the pencil value p%rho is actually sigma, the column
!  density, sigma=rho*2*H
!
        if (nzgrid==1) then
          tau = .5*kappa*p%rho
        else
          call not_implemented("calc_heat_cool","opacity for 3D")
          tau = 0.  ! to avoid compiler warnings
        endif
!
!  Analytical gray description of Hubeny (1990)
!  a1 is the optically thick contribution,
!  a3 the optically thin one.
!
        a1=0.375*tau ; a2=0.433013 ; a3=0.25/tau
!
!  Cooling for energy: 2*sigmaSB*p%TT**4/(a1+a2+a3)
!
        cooling = 2*sigmaSB*p%rho1*p%TT**4/(a1+a2+a3)
!
!  This cooling has dimension of energy over time.
!
        df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) - p%cv1*p%TT1*cooling
      endif
!
    endsubroutine calc_heat_cool
!***********************************************************************
    subroutine calc_heatcond_constchi(df,p)
!
!  Calculate the radiative diffusion term for constant chi:
!  lnTT version: cp*chi*Div(rho*TT*glnTT)/(rho*cv*TT)
!           = gamma*chi*(g2.glnTT+g2lnTT) where g2=glnrho+glnTT
!    TT version: cp*chi*Div(rho*gTT)/(rho*cv)
!           = gamma*chi*(g2.gTT+g2TT) where g2=glnrho
!
!  01-mar-07/dintrans: adapted from temperature_ionization
!
      use Diagnostics, only: max_mn_name
      use Sub, only: dot
!
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
      real, dimension (nx) :: g2
!
      intent(in) :: p
      intent(inout) :: df
!
      if (ltemperature_nolog) then
        call dot(p%glnrho,p%gTT,g2)
        g2=g2+p%del2TT
      else
        call dot(p%glnTT+p%glnrho,p%glnTT,g2)
        g2=g2+p%del2lnTT
      endif
!
!  Add heat conduction to RHS of temperature equation.
!
      ! Note that ilnTT=iTT for nonlog temperature.
      df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + gamma*chi*g2
!
!  Check maximum diffusion from thermal diffusion.
!
      if (lupdate_courant_dt) diffus_chi=diffus_chi+gamma*chi*dxyz_2
!
    endsubroutine calc_heatcond_constchi
!***********************************************************************
    subroutine calc_heatcond_chi_no_ddrsph(f,df,p)
!
!  For spherical coordinates only. Same as chi-const, but the
!  r-gradient will be reduced by a factor of 1-chi_r_reduce,
!  i.e. heat conduction in the radial direction is reduced.
!
!  27-mar-25/hongzhe and kuan: coded
!
      use Sub, only: dot
      use Deriv, only: der_x,der2_x
!
      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) :: rr, g2,g2_r, dlnrhodx, dTTdx, d2TTdx2
!
!  For now only implemented for ltemperature_nolog=T
!
      if ((.not.lspherical_coords) .or. (.not.ltemperature_nolog)) &
        call not_implemented('calc_heatcond_chi_no_ddrsph', 'for non-spherical coords or lnTT')
!
      call dot(p%glnrho,p%gTT,g2)
      g2=g2+p%del2TT
!
      rr=x(l1:l2)
      call der_x(f(:,m,n,ilnrho),dlnrhodx)
      call der_x(f(:,m,n,iTT),dTTdx)
      call der2_x(f(:,m,n,iTT),d2TTdx2)
      g2_r = dlnrhodx*dTTdx + 2./rr*dTTdx+d2TTdx2
      g2 = g2 - chi_r_reduce*g2_r
!
      df(l1:l2,m,n,iTT) = df(l1:l2,m,n,iTT) + gamma*chi*g2
!
    endsubroutine calc_heatcond_chi_no_ddrsph
!***********************************************************************
    subroutine calc_heatcond_cubicstepchi(df,p)
!
!  Calculate the radiative diffusion term for chi with cubic step profile:
!
!  lnTT version: cp*chi*Div(rho*TT*glnTT)/(rho*cv*TT)
!           = gamma*[chi*((glnrho+glnTT).glnTT+g2lnTT) + glnTT.gradchi]
!
!    TT version: cp*chi*Div(rho*gTT)/(rho*cv)
!           = gamma*[chi*(glnrho.gTT+g2TT) + gT.gradchi]
!
!  21-oct-18/joern: adapted from calc_heatcond_constchi
!
      use Diagnostics, only: max_mn_name
      use Sub, only: dot, cubic_step, cubic_der_step
!
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
      real, dimension (nx) :: g2, chi_z, gradchi_z
!
      intent(in) :: p
      intent(inout) :: df
!
      if (chi_zwidth == 0.) chi_zwidth = 5.*dz
      chi_z= chi + chi*(chi_jump-1.)*cubic_step(p%z_mn,chi_z0,-chi_zwidth)
!
      gradchi_z=chi*(chi_jump-1.)*cubic_der_step(p%z_mn,chi_z0,-chi_zwidth)
!
      if (ltemperature_nolog) then
        call dot(p%glnrho,p%gTT,g2)
        g2=g2+p%del2TT
      else
        call dot(p%glnTT+p%glnrho,p%glnTT,g2)
        g2=g2+p%del2lnTT
      endif
!
!  Add heat conduction to RHS of temperature equation.
!
      if (ltemperature_nolog) then
        df(l1:l2,m,n,iTT) = df(l1:l2,m,n,iTT) + gamma*(chi_z*g2 + gradchi_z*p%gTT(:,3))
      else
        df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + gamma*(chi_z*g2 + gradchi_z*p%glnTT(:,3))
      endif
!
!  Check maximum diffusion from thermal diffusion.
!
      if (lupdate_courant_dt) diffus_chi=diffus_chi+gamma*chi_z*dxyz_2
!
    endsubroutine calc_heatcond_cubicstepchi
!***********************************************************************
    subroutine calc_heatcond_constK(df,p)
!
!  Calculate the radiative diffusion term for constant K:
!
!  lnTT version: gamma*K/rho/TT/cp*div(TT*grad lnTT)
!                =gamma*K/rho/cp*(gradlnTT.gradlnTT + del2ln TT)
!    TT version: gamma*K/rho/cp*del2(TT)=gamma*chi*del2(TT)
!
!  Note: if ldensity=.false. then rho=1 and chi=K/cp
!
      use Diagnostics, only: max_mn_name
      use Sub, only: dot
!
      real, dimension(mx,my,mz,mvar) :: df
      type (pencil_case)  :: p
      real, dimension(nx) :: g2, chix
!
      intent(in) :: p
      intent(inout) :: df
!
!  Add heat conduction to RHS of temperature equation.
!  Note that rho does not in general need to be unity, even with Boussinesq.
!
      if (ldensity.or.lboussinesq.or.lanelastic) then
        chix=p%rho1*hcond0*p%cp1
      else
        chix=hcond0*p%cp1
      endif
      chix = gamma*chix
!
      if (ltemperature_nolog) then
        df(l1:l2,m,n,iTT)   = df(l1:l2,m,n,iTT)   + chix*p%del2TT
      else
        call dot(p%glnTT,p%glnTT,g2)
        df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + chix*(g2 + p%del2lnTT)
      endif
!
!  Check maximum diffusion from thermal diffusion.
!
      if (lupdate_courant_dt) diffus_chi=diffus_chi+chix*dxyz_2
!
    endsubroutine calc_heatcond_constK
!***********************************************************************
    subroutine calc_heatcond_Ktherm(df,p)
!
!  Calculate the radiative diffusion term for K=cte:
!
!  lnTT version: gamma*K/rho/TT/cp*div(T*grad lnTT)
!                =gamma*K/rho/cp*(gradlnTT.gradlnTT + del2ln TT)
!    TT version: gamma*K/rho/cp*del2(TT)=gamma*chi*del2(TT)
!
!  Note: if ldensity=.false. then rho=1 and chi=K/cp
!
      use EquationOfState, only: rho0
      use Diagnostics, only: max_mn_name
      use Sub, only: dot
!
      real, dimension(mx,my,mz,mvar) :: df
      type (pencil_case)  :: p
      real, dimension(nx) :: g2, chix,hcondTT
!
      intent(in) :: p
      intent(inout) :: df
!
      hcondTT=hcond0*sqrt(exp(p%lnTT))
!
!  Add heat conduction to RHS of temperature equation.
!  If ldensity=F, we need to divide by rho0, which can be /= 1.
!
      if (ldensity) then
        chix=p%rho1*hcondTT*p%cp1
      else
        chix=hcondTT*p%cp1/rho0
      endif
      chix = gamma*chix
!
      if (ltemperature_nolog) then
        df(l1:l2,m,n,iTT)   = df(l1:l2,m,n,iTT)   + chix*p%del2TT
      else
        call dot(p%glnTT,p%glnTT,g2)
        df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + chix*(g2 + p%del2lnTT)
      endif
!
!  Check maximum diffusion from thermal diffusion.
!
      if (lupdate_courant_dt) diffus_chi=diffus_chi+chix*dxyz_2
!
    endsubroutine calc_heatcond_Ktherm
!***********************************************************************
    subroutine calc_heatcond_kramers(df,p)
!
!  Heat conduction using Kramers' opacity law
!
!  23-feb-11/pete: coded
!  24-aug-15/MR: bounds for chi introduced
!
      use Diagnostics
      use Debug_IO, only: output_pencil
      use Sub, only: dot
      use General, only: notanumber
!
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
!
      intent(in) :: p
      intent(inout) :: df
!
      real, dimension(nx) :: thdiff,Krho1,g2
!
!  Diffusion of the form
!      cv*rho*D log(T)/Dt = ... + nab.(K*gradT)/T ,
!  where
!      K = K_0*(T**6.5/rho**2)**n.
!  In reality n=1, but we may need to use n\=1 for numerical reasons.
!
      Krho1 = hcond0_kramers*p%rho1**(2.*nkramers+1.)*p%TT**(6.5*nkramers)   ! = K/rho
      !Krho1 = hcond0_kramers*exp(-p%lnrho*(2.*nkramers+1.)+p%lnTT*(6.5*nkramers))   ! = K/rho
      !if (chimax_kramers>0.) &
      !  Krho1 = max(min(Krho1,chimax_kramers/p%cp1),chimin_kramers/p%cp1)
      call dot(-2.*nkramers*p%glnrho+(6.5*nkramers+1)*p%glnTT,p%glnTT,g2)
      thdiff = Krho1*(p%del2lnTT+g2)

      df(l1:l2,m,n,ilntt) = df(l1:l2,m,n,ilntt) + thdiff

      if (lupdate_courant_dt) diffus_chi=diffus_chi+(gamma*p%cp1*Krho1)*dxyz_2

    endsubroutine calc_heatcond_kramers
!***********************************************************************
    subroutine calc_heatcond_arctan(df,p)
!
!  Radiative diffusion with an arctan profile for the conductivity
!
!  Calculate gamma/(rho*cp)*div(K * grad TT)=
!      gamma*K/(rho*cp)*(grad LnK.grad TT + del2 TT)
!
!  16-may-07/gastine+dintrans: coded
!  01-mar-10/dintrans: introduced a mixed version with the ADI scheme that only
!  computes *during the explicit step* the term
!  gamma/(rho*cp)*grad(K).grad(TT) with grad(K)=dK/dT*grad(TT),
!  this term being less restrictive for the explicit timestep
!
      use Diagnostics, only: max_mn_name
      use Sub, only: dot, multsv
!
      real, dimension(mx,my,mz,mvar) :: df
      real, dimension (nx)   :: dhcond, g1, chix
      real, dimension (nx,3) :: gLnhcond=0.
      type (pencil_case)     :: p
!
      intent(in) :: p
      intent(inout) :: df
!
      call heatcond_TT(p%TT, hcond, dhcond)
!
!  must specify the new bottom value of hcond for the 'c1' BC
!     if (n == n1) hcond0=hcond(1)
      if (lADI_mixed) then
        call dot(p%gTT, p%gTT, g1)
        df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + gamma*p%rho1*p%cp1*dhcond*g1
        chix=0.
      else
!
!  grad LnK=grad_T Ln K.grad(TT)
!
        dhcond=dhcond/hcond
        call multsv(dhcond, p%gTT, gLnhcond)
        call dot(gLnhcond, p%gTT, g1)
        chix=p%rho1*p%cp1*hcond
        df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + gamma*chix*(g1+p%del2TT)
      endif
!
!  Check maximum diffusion from thermal diffusion.
!
      if (lupdate_courant_dt) diffus_chi=diffus_chi+gamma*chix*dxyz_2
!
    endsubroutine calc_heatcond_arctan
!***********************************************************************
    subroutine calc_heatcond(f,df,p)
!
!  Calculate the radiative diffusion term for a variable K:
!    ivar=lnTT --> 1/(rho*cv*T)*div(K*grad TT)
!    ivar=TT   --> 1/(rho*cv)*div(K*grad TT)
!
!  12-Mar-07/dintrans: coded
!
      use Diagnostics, only: max_mn_name
      use Sub, only: dot, step, der_step
      use Gravity, only: z1, z2
!
      real, dimension(mx,my,mz,mfarray) :: f
      real, dimension(mx,my,mz,mvar) :: df
      type (pencil_case) :: p
!
      real, dimension(nx) :: g2, dhcond, chix
      real, dimension (nx,3) :: glhc=0.,glnThcond
!
      intent(in) :: f,p
      intent(inout) :: df
!
      if (lhcond_global) then
        hcond=f(l1:l2,m,n,iglobal_hcond)
        glhc(:,1)=f(l1:l2,m,n,iglobal_glhc)
      else
        if (lgravz) then
          hcond = 1. + (hcond1-1.)*step(p%z_mn,z1,-widthlnTT) &
                     + (hcond2-1.)*step(p%z_mn,z2,widthlnTT)
          hcond = hcond0*hcond
          glhc(:,3) = (hcond1-1.)*der_step(p%z_mn,z1,-widthlnTT) &
                    + (hcond2-1.)*der_step(p%z_mn,z2,widthlnTT)
          glhc(:,3) = hcond0*glhc(:,3)
        elseif (lcylindrical_coords) then
          hcond = 1. + (hcond1-1.)*step(rcyl_mn,r_bcz,-widthlnTT)
          hcond = hcond0*hcond
          glhc(:,1) = hcond0*(hcond1-1.)*der_step(rcyl_mn,r_bcz,-widthlnTT)
        elseif (lgravr) then
          hcond = 1. + (hcond1-1.)*step(p%r_mn,r_bcz,-widthlnTT) &
                     + (hcond2-1.)*step(p%r_mn,r_ext,widthlnTT)
          hcond = hcond0*hcond
          dhcond=(hcond1-1.)*der_step(p%r_mn,r_bcz,-widthlnTT) &
                +(hcond2-1.)*der_step(p%r_mn,r_ext,widthlnTT)
          dhcond=hcond0*dhcond
          glhc(:,1) = x(l1:l2)/p%r_mn*dhcond
          glhc(:,2) = y(m)/p%r_mn*dhcond
          glhc(:,3) = z(n)/p%r_mn*dhcond
        endif
      endif
!
      if (ltemperature_nolog) then
        glnThcond = glhc/spread(hcond,2,3)              ! grad ln(hcond)
        call dot(p%gTT,glnThcond,g2)
        g2 = g2 + p%del2TT
      else
        glnThcond = p%glnTT + glhc/spread(hcond,2,3)    ! grad ln(T*hcond)
        call dot(p%glnTT,glnThcond,g2)
        g2 = g2 + p%del2lnTT
      endif
!
!  Add heat conduction to RHS of temperature equation.
!
      chix=p%rho1*hcond*p%cp1
      df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + gamma*chix*g2
!
!  Check maximum diffusion from thermal diffusion.
!
      if (lupdate_courant_dt) diffus_chi=diffus_chi+gamma*chix*dxyz_2
!
    endsubroutine calc_heatcond
!***********************************************************************
    subroutine calc_heatcond_tensor(df,p)
!
!  Calculates heat conduction parallel and perpendicular (isotropic)
!  to magnetic field lines.
!
!  25-aug-09/bing: moved from denergy_dt to here
!
      use Diagnostics, only: max_mn_name
      use Sub, only: dot,dot2,tensor_diffusion_coef
!
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
!
      real, dimension (nx) :: cosbgT,gT2,b2
      real, dimension (nx) :: vKpara,vKperp,rhs
!
      vKpara(:) = Kgpara
      vKperp(:) = Kgperp
!
      call tensor_diffusion_coef(p%glnTT,p%hlnTT,p%bij,p%bb,vKperp,vKpara,rhs,llog=.true.)
!
      df(l1:l2,m,n,ilnTT)=df(l1:l2,m,n,ilnTT)+rhs*p%rho1*gamma*p%cp1
!
      call dot(p%bb,p%glnTT,cosbgT)
      call dot2(p%glnTT,gT2)
      call dot2(p%bb,b2)
!
      where ((gT2<=tini).or.(b2<=tini))
        cosbgT=0.
      elsewhere
        cosbgT=cosbgT/sqrt(gT2*b2)
      endwhere
!
      if (lupdate_courant_dt) diffus_chi=diffus_chi+cosbgT*gamma*Kgpara*p%rho1*p%cp1*dxyz_2
!
    endsubroutine calc_heatcond_tensor
!***********************************************************************
    subroutine get_AIA_tab_resp(channel,lgT,part_den2,respon)
!
!  Determine the indices for the temperature used in the tabulated
!  response functions of AIA on SDO.
!
!  8-may-19/joern: coded
!
      real, dimension  (nx) :: respon
      real, dimension  (nx) :: part_den2
      real, dimension  (nx) :: lgT
!
      real                  :: a0, a1, a2
!
      character (len=3) :: channel
!
!  the response functions are approximated by gaussian curves
!  using the amplitude (a0), center (a1), width (the standard deviation) (a2) of gaussian
!  following:
!        func(x) = a0 * exp((x-a1)**2/2*a2**2)
!
   select case (trim(channel))
!
!   AIA channels
!
     case ('094')
       a0=  3.68132d-27
       a1=      6.86406
       a2=     0.171703
!
     case ('131')
       a0=  5.69305d-26
       a1=      5.75804
       a2=     0.225604
!
     case ('171')
       a0=  1.29671d-24
       a1=      5.93474
       a2=     0.173903
!
     case ('193')
       a0=  5.42492d-25
       a1=      6.17694
       a2=     0.120203
!
     case ('211')
       a0=  1.71252d-25
       a1=      6.26945
       a2=     0.118102
!
     case ('304')
       a0=  4.26987d-25
       a1=      4.92302
       a2=     0.128603
!
     case ('335')
       a0=  6.38811d-27
       a1=      5.33373
       a2=     0.471810
!
!  XRT channel'
!  we use Al-poly
     case ('XRT')
       a0=  2.72928d-25
       a1=      6.93243
       a2=     0.317607
!
     case default
       a0= 1.
       a1= 1.
       a2=0.1
!
  endselect
!
!  Construct the reponse function
!
      respon=a0*(exp(-1.*((lgT-a1)**2.)/(2.*a2**2.)))
!
!  two seconds exposure time
!  mulitplying with the particle density squared
!  getting the units right
!
      respon=2*respon*part_den2*1e-4
!
    endsubroutine get_AIA_tab_resp
!***********************************************************************
    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 energy.
!
!   1-jun-02/axel: adapted from magnetic fields
!
      use Diagnostics, only: parse_name, set_type
      use FArrayManager, only: farray_index_append
!
      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_TTmax=0; idiag_TTzmask=0; idiag_TTmin=0; idiag_TTm=0; idiag_fradtop=0
        idiag_TugTm=0; idiag_Trms=0; idiag_fradbot=0
        idiag_uxTm=0; idiag_uyTm=0; idiag_uzTm=0; idiag_gT2m=0
        idiag_guxgTm=0; idiag_guygTm=0; idiag_guzgTm=0
        idiag_Tugux_uxugTm=0; idiag_Tuguy_uyugTm=0; idiag_Tuguz_uzugTm=0
        idiag_Tdxpm=0; idiag_Tdypm=0; idiag_Tdzpm=0
        idiag_yHmax=0; idiag_yHmin=0; idiag_yHm=0; idiag_gTmax=0
        idiag_ethm=0; idiag_ssm=0; idiag_thcool=0
        idiag_dtchi=0; idiag_dtc=0
        idiag_eem=0; idiag_ppm=0; idiag_csm=0; idiag_csmax=0
        idiag_ppmx=0; idiag_ppmy=0; idiag_ppmz=0; idiag_ppuzmz=0
        idiag_TTmx=0; idiag_TTmy=0; idiag_TTmz=0; idiag_ethuxmx=0
        idiag_TT2mz=0; idiag_uxTmz=0; idiag_uyTmz=0; idiag_uzTmz=0
        idiag_ethmz=0; idiag_ethuxmz=0; idiag_ethuymz=0; idiag_ethuzmz=0
        idiag_TTmxy=0; idiag_TTmxz=0
        idiag_fpresxmz=0; idiag_fpresymz=0; idiag_fpreszmz=0; idiag_fradmz=0
        idiag_ethtot=0; idiag_fconvmz=0; idiag_TT2m=0
        ivid_pp=0
        idiag_EmAIA94mxy=0; idiag_EmAIA94mxz=0
        idiag_EmAIA131mxy=0; idiag_EmAIA131mxz=0
        idiag_EmAIA171mxy=0; idiag_EmAIA171mxz=0
        idiag_EmAIA193mxy=0; idiag_EmAIA193mxz=0
        idiag_EmAIA211mxy=0; idiag_EmAIA211mxz=0
        idiag_EmAIA304mxy=0; idiag_EmAIA304mxz=0
        idiag_EmAIA335mxy=0; idiag_EmAIA335mxz=0
        idiag_EmXRTmxy=0; idiag_EmXRTmxz=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),'TTmax',idiag_TTmax)
        call parse_name(iname,cname(iname),cform(iname),'TTzmask',idiag_TTzmask)
        call parse_name(iname,cname(iname),cform(iname),'gTmax',idiag_gTmax)
        call parse_name(iname,cname(iname),cform(iname),'TTmin',idiag_TTmin)
        call parse_name(iname,cname(iname),cform(iname),'TTm',idiag_TTm)
        call parse_name(iname,cname(iname),cform(iname),'TugTm',idiag_TugTm)
        call parse_name(iname,cname(iname),cform(iname),'Trms',idiag_Trms)
        call parse_name(iname,cname(iname),cform(iname),'TT2m',idiag_TT2m)
        call parse_name(iname,cname(iname),cform(iname),'uxTm',idiag_uxTm)
        call parse_name(iname,cname(iname),cform(iname),'uyTm',idiag_uyTm)
        call parse_name(iname,cname(iname),cform(iname),'uzTm',idiag_uzTm)
        call parse_name(iname,cname(iname),cform(iname),'gT2m',idiag_gT2m)
        call parse_name(iname,cname(iname),cform(iname),'guxgTm',idiag_guxgTm)
        call parse_name(iname,cname(iname),cform(iname),'guygTm',idiag_guygTm)
        call parse_name(iname,cname(iname),cform(iname),'guzgTm',idiag_guzgTm)
        call parse_name(iname,cname(iname),cform(iname),'Tugux_uxugT',idiag_Tugux_uxugTm)
        call parse_name(iname,cname(iname),cform(iname),'Tuguy_uyugT',idiag_Tuguy_uyugTm)
        call parse_name(iname,cname(iname),cform(iname),'Tuguz_uzugT',idiag_Tuguz_uzugTm)
        call parse_name(iname,cname(iname),cform(iname),'Tdxpm',idiag_Tdxpm)
        call parse_name(iname,cname(iname),cform(iname),'Tdypm',idiag_Tdypm)
        call parse_name(iname,cname(iname),cform(iname),'Tdzpm',idiag_Tdzpm)
        call parse_name(iname,cname(iname),cform(iname),'fradtop',idiag_fradtop)
        call parse_name(iname,cname(iname),cform(iname),'fradbot',idiag_fradbot)
        call parse_name(iname,cname(iname),cform(iname),'ethm',idiag_ethm)
        call parse_name(iname,cname(iname),cform(iname),'ethtot',idiag_ethtot)
        call parse_name(iname,cname(iname),cform(iname),'ssm',idiag_ssm)
        call parse_name(iname,cname(iname),cform(iname),'dtchi',idiag_dtchi)
        call parse_name(iname,cname(iname),cform(iname),'dtc',idiag_dtc)
        call parse_name(iname,cname(iname),cform(iname),'eem',idiag_eem)
        call parse_name(iname,cname(iname),cform(iname),'ppm',idiag_ppm)
        call parse_name(iname,cname(iname),cform(iname),'csmax',idiag_csmax)
        call parse_name(iname,cname(iname),cform(iname),'thcool',idiag_thcool)
      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)
        call parse_name(inamex,cnamex(inamex),cformx(inamex),'ethuxmx',idiag_ethuxmx)
      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)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'ppuzmz',idiag_ppuzmz)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'ethmz',idiag_ethmz)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'ethuxmz',idiag_ethuxmz)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'ethuymz',idiag_ethuymz)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'ethuzmz',idiag_ethuzmz)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'fpresxmz',idiag_fpresxmz)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'fpresymz',idiag_fpresymz)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'fpreszmz',idiag_fpreszmz)
!
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'TT2mz',idiag_TT2mz)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'uxTmz',idiag_uxTmz)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'uyTmz',idiag_uyTmz)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'uzTmz',idiag_uzTmz)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'fradmz',idiag_fradmz)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'fconvmz',idiag_fconvmz)
!
      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)
        call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy),'EmAIA94mxy',idiag_EmAIA94mxy)
        call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy),'EmAIA131mxy',idiag_EmAIA131mxy)
        call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy),'EmAIA171mxy',idiag_EmAIA171mxy)
        call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy),'EmAIA193mxy',idiag_EmAIA193mxy)
        call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy),'EmAIA211mxy',idiag_EmAIA211mxy)
        call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy),'EmAIA304mxy',idiag_EmAIA304mxy)
        call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy),'EmAIA335mxy',idiag_EmAIA335mxy)
        call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy),'EmXRTmxy',idiag_EmXRTmxy)
      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)
        call parse_name(inamexz,cnamexz(inamexz),cformxz(inamexz),'EmAIA94mxz',idiag_EmAIA94mxz)
        call parse_name(inamexz,cnamexz(inamexz),cformxz(inamexz),'EmAIA131mxz',idiag_EmAIA131mxz)
        call parse_name(inamexz,cnamexz(inamexz),cformxz(inamexz),'EmAIA171mxz',idiag_EmAIA171mxz)
        call parse_name(inamexz,cnamexz(inamexz),cformxz(inamexz),'EmAIA193mxz',idiag_EmAIA193mxz)
        call parse_name(inamexz,cnamexz(inamexz),cformxz(inamexz),'EmAIA211mxz',idiag_EmAIA211mxz)
        call parse_name(inamexz,cnamexz(inamexz),cformxz(inamexz),'EmAIA304mxz',idiag_EmAIA304mxz)
        call parse_name(inamexz,cnamexz(inamexz),cformxz(inamexz),'EmAIA335mxz',idiag_EmAIA335mxz)
        call parse_name(inamexz,cnamexz(inamexz),cformxz(inamexz),'EmXRTmxz',idiag_EmXRTmxz)
      enddo
!
!  check for those quantities for which we want video slices
!     
      if (lwrite_slices) then 
        where(cnamev=='TT'.or.cnamev=='lnTT') cformv='DEFINED'
        do iname=1,nnamev
          call parse_name(iname,cnamev(iname),cformv(iname),'pp',ivid_pp)
        enddo
      endif
!
!  Write column where which variable is stored.
!
      if (lwr) then
        if (ltemperature_nolog) then
          call farray_index_append('ilnTT', 0)
        else
          call farray_index_append('iTT', 0)
        endif
      endif
!
    endsubroutine rprint_energy
!***********************************************************************
    subroutine get_slices_energy(f,slices)
!
      use Slices_methods, only: assign_slices_scal, process_slices, log2d, exp2d
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (slice_data) :: slices
!
!  Loop over slices
!
      select case (trim(slices%name))
!
!  Temperature.
!
        case ('TT')
          if (iTT>0) then
            call assign_slices_scal(slices,f,iTT)
          else
            call assign_slices_scal(slices,f,ilnTT)
            call process_slices(slices,exp2d)
          endif
!  lnTT
        case ('lnTT')
          call assign_slices_scal(slices,f,ilnTT)
          if (iTT>0) call process_slices(slices,log2d)
!  Pressure
        case ('pp')
          call assign_slices_scal(slices,pp_xy,pp_xz,pp_yz,pp_xy2,pp_xy3,pp_xy4,pp_xz2,pp_r)
!
      endselect
!
    endsubroutine get_slices_energy
!***********************************************************************
    subroutine single_polytrope(f)
!
!  04-aug-07/dintrans: a single polytrope with index mpoly0
!
      use Gravity, only: gravz
      use EquationOfState, only: cs20, lnrho0, cs2bot, cs2top
!
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f
      real :: beta, zbot, ztop, T0, temp
!
!  beta is the (negative) temperature gradient
!  beta = -(g/cp) /[(1-1/gamma)*(m+1)]
!  gamma*(Rgas/mu)T0 = cs2(ad) = cp*T0*gamma_m1,
!  so T0 = cs20*cp1/gamma_m1
!
      if (.not. leos) call fatal_error('single_polytrope', &
                      'EOS=noeos, but polytrope requires an EQUATION OF STATE')

      beta=-cp1*gravz/(mpoly0+1.)*gamma/gamma_m1
      ztop=xyz0(3)+Lxyz(3)
      zbot=xyz0(3)
      T0=cs20*cp1/gamma_m1
      print*, 'polytrope: mpoly0, beta, T0=', mpoly0, beta, T0
!
      do imn=1,ny*nz
        n=nn(imn)
        m=mm(imn)
        temp=T0+beta*(ztop-z(n))
        if (ltemperature_nolog) then
          f(:,m,n,iTT)  =temp
        else
          f(:,m,n,ilnTT)=log(temp)
        endif
        f(:,m,n,ilnrho)=lnrho0+mpoly0*log(temp/T0)
      enddo
      cs2bot=gamma_m1*(T0+beta*(ztop-zbot))
      cs2top=cs20
!
    endsubroutine single_polytrope
!***********************************************************************
    subroutine piecew_poly(f)
!
!  Computes piecewise polytropic and hydrostatic atmosphere.
!  Adapted from single_polytrope.
!
!  19-jan-10/bing: coded
!  The layout is the same than in entropy.f90:
!  ------------ ztop Ttop
!     mpoly2
!  ------------ z2   T2, lnrho2
!     mpoly0
!  ------------ z1   T1, lnrho1
!     mpoly1
!  ------------ zbot
!
      use Gravity, only: gravz, z1, z2
      use EquationOfState, only: cs2top, cs2bot, lnrho0
!
      real, dimension(mx,my,mz,mfarray) :: f
      real :: Ttop, T1, T2, beta0, beta1, beta2, temp
      real :: lnrhotop, lnrho1, lnrho2, ztop
      integer :: i
!
      if (.not. leos) &
        call fatal_error('piecew_poly','EOS=noeos, but polytrope requires an EQUATION OF STATE')

!  Top boundary values.
!
      Ttop=cs2top*cp1/gamma_m1
      lnrhotop = lnrho0
      ztop=xyz0(3)+Lxyz(3)
!
!  Temperature gradients.
!
      beta0 =-cp1*gravz/(mpoly0+1.)*gamma/gamma_m1
      beta1 =-cp1*gravz/(mpoly1+1.)*gamma/gamma_m1
      beta2 =-cp1*gravz/(mpoly2+1.)*gamma/gamma_m1
!
      T2 = Ttop + beta2*(ztop-z2)
      T1 = T2   + beta0*(z2-z1)
!
      lnrho2 = lnrhotop+mpoly2*log(T2/Ttop)
      lnrho1 = lnrho2  +mpoly0*log(T1/T2)
!
      do  i=n2,n1,-1
        if (z(i) >= z2)                 temp = Ttop + beta2*(ztop-z(i))
        if (z(i) < z2 .and. z(i) >= z1) temp = T2   + beta0*(z2-z(i))
        if (z(i) < z1)                  temp = T1   + beta1*(z1-z(i))
!
        if (ltemperature_nolog) then
          f(:,:,i,iTT)  =temp
        else
          f(:,:,i,ilnTT)=log(temp)
        endif
!
        if (z(i) >= z2) f(:,:,i,ilnrho)=lnrhotop+mpoly2*log(temp/Ttop)
        if (z(i) < z2 .and. z(i) >= z1 ) f(:,:,i,ilnrho)=lnrho2+mpoly0*log(temp/T2)
        if (z(i) < z1) f(:,:,i,ilnrho)=lnrho1+mpoly1*log(temp/T1)
      enddo
!
! one also needs to refresh cs2bot in case of a 'cT' BC for the temperature
!
      cs2bot=gamma_m1*(T1 + beta1*(z1-xyz0(3)))
!
    endsubroutine piecew_poly
!***********************************************************************
    subroutine star_heat(f)
!
!  Initialize energy for two superposed polytropes with a central heating
!
!  04-fev-2011/dintrans: coded
!
      use EquationOfState, only: rho0, lnrho0, get_soundspeed, eoscalc, ilnrho_TT
      use Sub, only: step, interp1, erfunc
!
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f
      integer, parameter   :: nr=100
      integer              :: i,l,iter
      real, dimension (nr) :: r, lnrho, temp, lumi, g, hcond
      real                 :: u,r_mn,lnrho_r,temp_r,cs2,ss
      real                 :: rhotop, rbot,rt_old,rt_new,rhobot
      real                 :: rb_old,rb_new,crit,r_max

      if (.not. leos) &
        call fatal_error('star_heat','EOS=noeos, but star heating requires an EQUATION OF STATE')
!
!  Define the radial grid r=[0,r_max], luminosity and gravity
!
      r_max=sqrt(xyz1(1)**2+xyz1(2)**2+xyz1(3)**2)
      r(1)=0. ; lumi(1)=0. ; g(1)=0.
      do i=2,nr
        r(i)=r_max*real(i-1)/(nr-1)
        u=r(i)/sqrt(2.)/wheat
        lumi(i)=luminosity*(1.-exp(-u**2))
        g(i)=-lumi(i)/(2.*pi*r(i))*(mpoly0+1.)/hcond0*gamma_m1/gamma
      enddo
!
      hcond1=(mpoly1+1.)/(mpoly0+1.)
      hcond2=(mpoly2+1.)/(mpoly0+1.)
      hcond = 1. + (hcond1-1.)*step(r,r_bcz,-widthlnTT) &
                 + (hcond2-1.)*step(r,r_ext,widthlnTT)
      hcond = hcond0*hcond
!
      rbot=rho0
      rt_old=0.01*rbot
      rt_new=0.012*rbot
      rhotop=rt_old
      call strat_heat(nr, r, lumi, g, hcond, temp, lnrho, rhotop, rhobot)
      print*, 'find rhobot=', rhobot
      rb_old=rhobot
!
      rhotop=rt_new
      call strat_heat(nr, r, lumi, g, hcond, temp, lnrho, rhotop, rhobot)
      if (lroot) print*, 'find rhobot=', rhobot
      rb_new=rhobot
!
      do iter=1,10
        rhotop=rt_old+(rt_new-rt_old)/(rb_new-rb_old)*(rbot-rb_old)
!
        crit=abs(rhotop/rt_new-1.)
        if (crit<=1e-4) exit
        call strat_heat(nr, r, lumi, g, hcond, temp, lnrho, rhotop, rhobot)
!
!  Update new estimates.
!
        rt_old=rt_new
        rb_old=rb_new
        rt_new=rhotop
        rb_new=rhobot
      enddo
      if (lroot) print*,'- iteration completed: rhotop,crit=',rhotop,crit
!
!  One needs to refresh rho0 and lnrho0 because the density top value
!  has changed --> important for the future EOS calculations (ss, ...)
!
      lnrho0=lnrho(nr)
      rho0=exp(lnrho0)
      print*,'new values for lnrho0 and rho0:', lnrho0, rho0
!
      do imn=1,ny*nz
        n=nn(imn)
        m=mm(imn)
        do l=l1,l2
          r_mn=sqrt(x(l)**2+y(m)**2+z(n)**2)
          lnrho_r=interp1(r,lnrho,nr,r_mn)
          temp_r=interp1(r,temp,nr,r_mn)
          f(l,m,n,ilnrho)=lnrho_r
          f(l,m,n,ilnTT)=temp_r
        enddo
      enddo
!
      if (lroot) then
        print*,'--> writing initial setup to data/proc0/setup.dat'
        open(unit=11,file=trim(directory)//'/setup.dat')
        write(11,'(a1,a5,6a12)') '#','r','rho','ss','cs2','grav','lumi','hcond'
        do i=nr,1,-1
          u=r(i)/sqrt(2.)/wheat
          call get_soundspeed(temp(i),cs2)
          call eoscalc(ilnrho_TT,lnrho(i),temp(i),ss=ss)
          write(11,'(f6.3,6e12.3)') r(i),exp(lnrho(i)),ss,cs2,g(i),lumi(i), hcond(i)
        enddo
        close(11)
      endif
!
    endsubroutine star_heat
!***********************************************************************
    subroutine strat_heat(nr,r,lumi,g,hcond,temp,lnrho,rhotop,rhobot)
!
      use EquationOfState, only: cs20
      use Sub, only: interp1
!
      integer              :: nr, i
      real, dimension (nr) :: r, lnrho, temp, lumi, g, hcond
      real                 :: dr,dtemp,dlnrho
      real                 :: rhotop,rhobot,lnrhobot
!
      if (.not. leos) call fatal_error('strat_heat', &
            'EOS=noeos, but stratified heating requires an EQUATION OF STATE')

      temp(nr)=cs20/gamma_m1 ; lnrho(nr)=alog(rhotop)
      dr=r(2)
      do i=nr-1,1,-1
        if (r(i+1) > r_ext) then
          ! Isothermal exterior: mpoly2 but force T=cte
          dtemp=0.
          dlnrho=-gamma*g(i+1)/cs20
        elseif (r(i+1) > r_bcz) then
          ! Convection zone: mpoly0
! adiabatic stratification
!          dtemp=-g(i+1)
!          dlnrho=3./2.*dtemp/temp(i+1)
          dtemp=lumi(i+1)/(2.*pi*r(i+1))/hcond(i+1)
          dlnrho=mpoly0*dtemp/temp(i+1)
        else
          ! Radiative zone: mpoly1
          dtemp=lumi(i+1)/(2.*pi*r(i+1))/hcond(i+1)
          dlnrho=mpoly1*dtemp/temp(i+1)
        endif
        temp(i)=temp(i+1)+dtemp*dr
        lnrho(i)=lnrho(i+1)+dlnrho*dr
      enddo
!
      lnrhobot=interp1(r,lnrho,nr,r_ext)
      rhobot=exp(lnrhobot)
!
    endsubroutine strat_heat
!***********************************************************************
    subroutine dynamical_thermal_diffusion(uc)
!
!  Dummy subroutine
!
      real, intent(in) :: uc
!
      call keep_compiler_quiet(uc)
      call not_implemented('dynamical_thermal_diffusion','')
!
    endsubroutine dynamical_thermal_diffusion
!***********************************************************************
    subroutine expand_shands_energy
!
!  Expands shorthand labels of temperature diagnostics.
!
!  16-may-12/MR: coded
!
      use Diagnostics, only: expand_cname
      use General, only: lextend_vector
!
      integer :: nname_prev

      if (nname>0) then
        nname_prev=nname
        call expand_cname(cname,nname,'uuTm', 'u', cform)
        call expand_cname(cname,nname,'gugTm','gu',cform)
        call expand_cname(cname,nname,'Tdpm', 'Td',cform)
        if (nname>nname_prev) then
          if (.not.lextend_vector(fname,nname)) & ! sanity check
            call fatal_error('expand_shands_energy','could not extend fname')
        endif
      endif
!
    endsubroutine expand_shands_energy
!***********************************************************************
    subroutine heatcond_TT_2d(TT, hcond, dhcond)
!
! 07-Sep-07/gastine: computes 2-D radiative conductivity hcond(T) with
! its derivative dhcond=dhcond(T)/dT for radiative conductivity hole (kappa-mechanism problem).
!
      real, dimension(:,:), intent(in) :: TT
      real, dimension(:,:), intent(out) :: hcond
      real, dimension(:,:), optional :: dhcond
!
      real :: hole_alpha

      hole_alpha=(Kmax-Kmin)/(pi/2.+atan(hole_slope*hole_width**2))
      hcond=hole_slope*(TT-Tbump-hole_width)*(TT-Tbump+hole_width)
      if (present(dhcond)) dhcond=2.*hole_alpha/(1.+hcond**2)*hole_slope*(TT-Tbump)
      hcond=Kmax+hole_alpha*(-pi/2.+atan(hcond))
!
    endsubroutine heatcond_TT_2d
!***********************************************************************
    subroutine heatcond_TT_1d(TT, hcond, dhcond)
!
! 18-Sep-07/dintrans: computes 1-D radiative conductivity with
! its derivative dhcond=dhcond(T)/dT for radiative conductivity hole (kappa-mechanism problem).
!
      real, dimension(:), intent(in) :: TT
      real, dimension(:), intent(out) :: hcond
      real, dimension(:), optional :: dhcond
!
      real :: hole_alpha

      hole_alpha=(Kmax-Kmin)/(pi/2.+atan(hole_slope*hole_width**2))
      hcond=hole_slope*(TT-Tbump-hole_width)*(TT-Tbump+hole_width)
      if (present(dhcond)) dhcond=2.*hole_alpha/(1.+hcond**2)*hole_slope*(TT-Tbump)
      hcond=Kmax+hole_alpha*(-pi/2.+atan(hcond))
!
    endsubroutine heatcond_TT_1d
!***********************************************************************
    subroutine heatcond_TT_0d(TT, hcond, dhcond)
!
! 07-Sep-07/gastine: computes the radiative conductivity hcond(T) with
! its derivative dhcond=dhcond(T)/dT for radiative conductivity hole (kappa-mechanism problem).
!
      real, intent(in) :: TT
      real, intent(out) :: hcond
      real, optional :: dhcond
!
      real :: hole_alpha

      hole_alpha=(Kmax-Kmin)/(pi/2.+atan(hole_slope*hole_width**2))
      hcond=hole_slope*(TT-Tbump-hole_width)*(TT-Tbump+hole_width)
      if (present(dhcond)) dhcond=2.*hole_alpha/(1.+hcond**2)*hole_slope*(TT-Tbump)
      hcond=Kmax+hole_alpha*(-pi/2.+atan(hcond))
!
    endsubroutine heatcond_TT_0d
!***********************************************************************
    subroutine update_char_vel_energy(f)
!
!  Updates characteristic veelocity for slope-limited diffusion.
!
!  9-jul-18/joern: adapted from update_char_vel_energy in entropy.f90
!
      use EquationOfState, only: eoscalc
!      use General, only: staggered_mean_scal
      use General, only: staggered_max_scal
!
      real, dimension(mx,my,mz,mfarray), intent(INOUT) :: f
!
      real, dimension(mx) :: cs2

      if (.not. leos) call fatal_error('update_char_vel_energy', &
                      'EOS=noeos, but sound speed requires an EQUATION OF STATE')
!
!  Calculate sound speed and store temporarily in first slot of diffusive fluxes.
!
        do n=1,mz; do m=1,my
          call eoscalc(f,mx,cs2=cs2)
          f(:,m,n,iFF_diff) = sqrt(cs2)   ! sqrt needed as we need the speed.
        enddo; enddo
!
!        call staggered_mean_scal(f,iFF_diff,iFF_char_c,w_sldchar_ene)
        call staggered_max_scal(f,iFF_diff,iFF_char_c,w_sldchar_ene)
!
    endsubroutine update_char_vel_energy
!***********************************************************************
    subroutine pushpars2c(p_par)

    use Syscalls, only: copy_addr
    use General, only: string_to_enum

    integer, parameter :: n_pars=100
    integer(KIND=ikind8), dimension(n_pars) :: p_par

    call copy_addr(chi,p_par(1))
    call copy_addr(widthlntt,p_par(2))
    call copy_addr(lntt_const,p_par(3))
    call copy_addr(tt_const,p_par(4))
    call copy_addr(kgperp,p_par(5))
    call copy_addr(kgpara,p_par(6))
    call copy_addr(chi_jump,p_par(7))
    call copy_addr(chi_z0,p_par(8))
    call copy_addr(chi_r_reduce,p_par(9))
    call copy_addr(r_bcz,p_par(10))
    call copy_addr(chi_shock,p_par(11))
    call copy_addr(chi_hyper3,p_par(12))
    call copy_addr(chi_hyper3_mesh,p_par(13))
    call copy_addr(tbump,p_par(14))
    call copy_addr(kmin,p_par(15))
    call copy_addr(kmax,p_par(16))
    call copy_addr(hole_slope,p_par(17))
    call copy_addr(hole_width,p_par(18))
    call copy_addr(hcond0,p_par(19))
    call copy_addr(hcond1,p_par(20))
    call copy_addr(hcond2,p_par(21))
    call copy_addr(luminosity,p_par(22))
    call copy_addr(wheat,p_par(23))
    call copy_addr(rcool,p_par(24))
    call copy_addr(wcool,p_par(25))
    call copy_addr(cool,p_par(26))
    call copy_addr(beta_bouss,p_par(27))
    call copy_addr(h_sld_ene,p_par(28))
    call copy_addr(nlf_sld_ene,p_par(29))
    call copy_addr(ladvection_temperature,p_par(30)) ! bool
    call copy_addr(lupw_lntt,p_par(31)) ! bool
    call copy_addr(lcalc_heat_cool,p_par(32)) ! bool
    call copy_addr(lheatc_hyper3,p_par(33)) ! bool
    call copy_addr(lheatc_kconst,p_par(34)) ! bool
    call copy_addr(lheatc_kprof,p_par(35)) ! bool
    call copy_addr(lheatc_karctan,p_par(36)) ! bool
    call copy_addr(lheatc_tensordiffusion,p_par(37)) ! bool
    call copy_addr(lheatc_hyper3_mesh,p_par(38)) ! bool
    call copy_addr(lheatc_chiconst,p_par(39)) ! bool
    call copy_addr(lheatc_chi_reduce_ddr,p_par(40)) ! bool
    call copy_addr(lhcond_global,p_par(41)) ! bool
    call copy_addr(lheatc_chicubicstep,p_par(42)) ! bool
    call copy_addr(lheatc_shock,p_par(43)) ! bool
    call copy_addr(lheatc_hyper3_polar,p_par(44)) ! bool
    call copy_addr(lheatc_ktherm,p_par(45)) ! bool
    call copy_addr(lheatc_kramers,p_par(46)) ! bool
    call copy_addr(hcond0_kramers,p_par(47))
    call copy_addr(nkramers,p_par(48))
    call copy_addr(lviscosity_heat,p_par(49)) ! bool
    call copy_addr(lenergy_slope_limited,p_par(52)) ! bool
    call copy_addr(ladi_mixed,p_par(53)) ! bool
    call copy_addr(gamma,p_par(54))
    call copy_addr(gamma1,p_par(55))
    call copy_addr(gamma_m1,p_par(56))
    call copy_addr(gradtt0,p_par(57)) ! real3
    call string_to_enum(enum_borderss,borderss)
    call copy_addr(enum_borderss,p_par(58)) ! int
 
    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
