!$Id$
!  This module can replace the entropy module by using _T_ as dependent
!  variable. For a perfect gas with constant coefficients (no ionization)
!  we have (1-1/gamma) * cp*T = cs02 * exp( (gamma-1)*ln(rho/rho0)-gamma*s/cp )
!
!  At a later point we may want to rename the module Entropy into Energy

!** 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.
!
! MVAR CONTRIBUTION 1
! MAUX CONTRIBUTION 0
!
! PENCILS PROVIDED Ma2; uglnTT; fpres(3)
!
!***************************************************************
module Entropy

  use Cparam
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages
  use EquationOfState, only: mpoly0,mpoly1
  use Interstellar

  implicit none
  public :: ADI_constK, ADI_Kprof


  include 'entropy.h'

  real :: radius_lnTT=0.1,ampl_lnTT=0.,widthlnTT=2*epsi
  real :: lnTT_left=1.0,lnTT_right=1.0,lnTT_const=0.0,TT_const=1.0
  real :: kx_lnTT=1.0,ky_lnTT=1.0,kz_lnTT=1.0
  real :: chi=impossible,heat_uniform=0.0
  real :: zbot=0.0,ztop=0.0
  real :: tau_heat_cor=-1.0,tau_damp_cor=-1.0,zcor=0.0,TT_cor=0.0
  real :: center1_x=0., center1_y=0., center1_z=0.
  real :: r_bcz=0.
  real :: Tbump=0.,Kmin=0.,Kmax=0.
  integer, parameter :: nheatc_max=1
  logical :: lpressuregradient_gas=.true.,ladvection_temperature=.true.
  logical :: lupw_lnTT=.false.
  logical :: lheatc_Kconst=.false.,lheatc_Kprof=.false.,lheatc_Karctan=.false.
  logical :: lheatc_chiconst=.false.,lheatc_chiconst_accurate=.false.
  character (len=labellen) :: iheatcond='nothing'
  logical :: lhcond_global=.false.
  logical :: lviscosity_heat=.true.
  integer :: iglobal_hcond=0
  integer :: iglobal_glhc=0

  character (len=labellen), dimension(ninit) :: initlnTT='nothing'
  character (len=intlen) :: iinit_str

! Delete (or use) me asap!
  real :: hcond0=impossible, hcond1=1.,Fbot,FbotKbot,Ftop,Kbot,FtopKtop
  logical :: lmultilayer=.false.

! input parameters
  namelist /entropy_init_pars/ &
      initlnTT,radius_lnTT,ampl_lnTT,widthlnTT, &
      lnTT_left,lnTT_right,lnTT_const,TT_const, &
      kx_lnTT,ky_lnTT,kz_lnTT,center1_x,center1_y,center1_z, &
      mpoly0,mpoly1,r_bcz, &
      Fbot,Tbump,Kmin,Kmax

! run parameters
  namelist /entropy_run_pars/ &
      lupw_lnTT,lpressuregradient_gas,ladvection_temperature, &
      heat_uniform,chi,iheatcond,tau_heat_cor,tau_damp_cor,zcor,TT_cor, &
      lheatc_chiconst_accurate,hcond0, &
      Tbump,Kmin,Kmax, &
      widthlnTT,mpoly0,mpoly1, &
      lhcond_global,lviscosity_heat, &
      Fbot,Tbump,Kmin,Kmax
!
! other variables (needs to be consistent with reset list below)
  integer :: idiag_TTmax=0,idiag_TTmin=0,idiag_TTm=0
  integer :: idiag_ethm=0,idiag_ssm=0
  integer :: idiag_dtchi=0,idiag_dtc=0
  integer :: idiag_eem=0,idiag_ppm=0,idiag_csm=0
 
  contains

!***********************************************************************
    subroutine register_entropy()
!
!  initialise variables which should know that we solve an entropy
!  equation: ilnTT, etc; increase nvar accordingly
!
!  6-nov-01/wolf: coded
!
      use FArrayManager
!
      call farray_register_pde('lnTT',ilnTT)
!
!  Identify version number.
!
      if (lroot) call svn_id( &
           "$Id$")
!
!  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
!                                       
    endsubroutine register_entropy
!***********************************************************************
    subroutine initialize_entropy(f)
!
!  called by run.f90 after reading parameters, but before the time loop
!
!  21-jul-2002/wolf: coded
!
      use Cdata
      use FArrayManager
      use Gravity, only: g0
      use EquationOfState
      use Sub, only: step,der_step
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx)   :: hcond,dhcond
      logical :: lnothing
!
      if (.not. leos) then
         call fatal_error('initialize_entropy','EOS=noeos but temperature_TT requires an EQUATION OF STATE for the fluid')
      endif
!
      call select_eos_variable('TT',ilnTT)
!
!  Check whether we want heat conduction
!
      lheatc_Kconst= .false.
      lheatc_Kprof= .false.
      lheatc_Karctan= .false.
      lheatc_chiconst = .false.
      lnothing = .false.
!
      select case (iheatcond)
        case ('K-const')
          lheatc_Kconst=.true.
          call information('initialize_TT', &
          ' heat conduction: K=cst --> gamma*K/(rho*cp)*lapl(T)')
          if (initlnTT(1).ne.'rad_equil') &
            Fbot=gamma/(gamma-1.)*hcond0*g0/(mpoly0+1.)
        case ('K-profile')
          lheatc_Kprof=.true.
! 
!  TODO..... ailleurs !
!
          hcond1=(mpoly1+1.)/(mpoly0+1.)
          Fbot=gamma/(gamma-1.)*hcond0*g0/(mpoly0+1.)
          call information('initialize_entropy',' heat conduction: K=K(r)')
        case ('K-arctan')
          lheatc_Karctan=.true.         
          call information('initialize_entropy',' heat conduction: arctan profile')
        case ('chi-const')
          lheatc_chiconst=.true.
          call information('initialize_entropy',' heat conduction: constant chi')
        case ('nothing')
          if (lroot .and. (.not. lnothing)) print*,'heat conduction: nothing'
        case default
          if (lroot) then
            write(unit=errormsg,fmt=*)  &
                'No such value iheatcond = ', trim(iheatcond)
            call fatal_error('initialize_entropy',errormsg)
          endif
       endselect
       lnothing=.true.
!
!  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 = 1. + (hcond1-1.)*step(x(l1:l2),r_bcz,-widthlnTT)
           hcond = hcond0*hcond
           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
!
!  A word of warning...
!
      if (lheatc_Kconst .and. hcond0==0.0) then
        call warning('initialize_entropy', 'hcond0 is zero!')
      endif
      if (lheatc_Kprof .and. hcond0==0.0) then
        call warning('initialize_entropy', 'hcond0 is zero!')
      endif
      if (lheatc_chiconst .and. chi==0.0) then
        call warning('initialize_entropy','chi is zero!')
      endif
      if (iheatcond=='nothing') then
        if (hcond0 /= impossible) call warning('initialize_entropy', 'No heat conduction, but hcond0 /= 0')
        if (chi /= impossible) call warning('initialize_entropy', 'No heat conduction, but chi /= 0')
      endif
!
    endsubroutine initialize_entropy
!***********************************************************************
    subroutine init_ss(f)
!
!  initialise lnTT; called from start.f90
!
!  13-dec-2002/axel+tobi: adapted from init_ss
!
!  initialise entropy; called from start.f90
!  07-nov-2001/wolf: coded 
!  24-nov-2002/tony: renamed for consistancy (i.e. init_[variable name])
!
      use General,  only: itoa
      use Sub,      only: blob
      use Initcond, only: jump
      use InitialCondition, only: initial_condition_ss
!
      real, dimension (mx,my,mz,mfarray), intent (inout) :: f
      logical :: lnothing=.true.
      integer :: j
!
      do j=1,ninit
!
        if (initlnTT(j)/='nothing') then
!
          lnothing=.false.
!
          iinit_str=itoa(j)
!
!  select different initial conditions
!
          select case (initlnTT(j))
          case ('zero', '0'); f(:,:,:,ilnTT) = 0.
          case ('const_lnTT'); f(:,:,:,ilnTT)=f(:,:,:,ilnTT)+lnTT_const
          case ('const_TT'); f(:,:,:,ilnTT)=f(:,:,:,ilnTT)+log(TT_const)
          case ('single_polytrope'); call single_polytrope(f)
          case ('rad_equil')
            call rad_equil(f)
            if (ampl_lnTT.ne.0.) then
              print*,'add a bubble with:',ampl_lnTT,radius_lnTT,center1_x,center1_y,center1_z
              call blob(ampl_lnTT,f,ilnTT,radius_lnTT,center1_x,center1_y,center1_z)
              call blob(-ampl_lnTT,f,ilnrho,radius_lnTT,center1_x,center1_y,center1_z)
            endif
          case ('bubble_hs')
!         print*,'init_lnTT: put bubble in hydrostatic equilibrium: radius_lnTT,ampl_lnTT=',radius_lnTT,ampl_lnTT,center1_x,center1_y,center1_z
            call blob(ampl_lnTT,f,ilnTT,radius_lnTT,center1_x,center1_y,center1_z)
            call blob(-ampl_lnTT,f,ilnrho,radius_lnTT,center1_x,center1_y,center1_z)
!
          case default
          !
          !  Catch unknown values
          !
            write(unit=errormsg,fmt=*) 'No such value for init_TT(' &
                //trim(iinit_str)//'): ',trim(initlnTT(j))
            call fatal_error('init_TT',errormsg)

          endselect

          if (lroot) print*,'init_TT: init_TT(' &
              //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_ss: nothing'
!
    endsubroutine init_ss
!***********************************************************************
    subroutine pencil_criteria_entropy()
!
      use Cdata
!
      if (ldt) lpenc_requested(i_cs2)=.true.
      if (lpressuregradient_gas) lpenc_requested(i_fpres)=.true.
!
      if (lviscosity.and.lviscosity_heat) then
        lpenc_requested(i_cv1)=.true.
        lpenc_requested(i_visc_heat)=.true.
      endif
!
      if (ldensity) lpenc_requested(i_divu)=.true.
!
      if (ladvection_temperature) lpenc_requested(i_uglnTT)=.true.
!
      if (lheatc_chiconst) lpenc_requested(i_del2lnTT)=.true.
!
      if (lheatc_Kconst) then
        lpenc_requested(i_rho1)=.true.
        lpenc_requested(i_del2lnTT)=.true.
        lpenc_requested(i_cp1)=.true.
      endif
!
      if (lheatc_Kprof) then
        lpenc_requested(i_rho1)=.true.
        lpenc_requested(i_glnTT)=.true.
        lpenc_requested(i_del2lnTT)=.true.
        lpenc_requested(i_cp1)=.true.
      endif
!
      if (lheatc_Karctan) then
        lpenc_requested(i_rho1)=.true.
        lpenc_requested(i_TT)=.true.
        lpenc_requested(i_glnTT)=.true.
        lpenc_requested(i_del2lnTT)=.true.
        lpenc_requested(i_cp1)=.true.
      endif
!
!  Diagnostics
!
      if (idiag_TTmax/=0) lpenc_diagnos(i_TT)  =.true.
      if (idiag_TTmin/=0) lpenc_diagnos(i_TT)  =.true.
      if (idiag_TTm/=0)   lpenc_diagnos(i_TT)  =.true.
      if (idiag_ethm/=0) then
                          lpenc_diagnos(i_rho1)=.true.
                          lpenc_diagnos(i_ee)  =.true.
      endif
      if (idiag_ssm/=0)   lpenc_diagnos(i_ss)  =.true.
      if (idiag_dtchi/=0) then
                          lpenc_diagnos(i_rho1)=.true.
                          lpenc_diagnos(i_cv1) =.true.
      endif
      if (idiag_dtchi/=0)  lpenc_diagnos(i_cs2)=.true.
      if (idiag_csm/=0)    lpenc_diagnos(i_cs2)=.true.
      if (idiag_eem/=0)    lpenc_diagnos(i_ee) =.true.
      if (idiag_ppm/=0)    lpenc_diagnos(i_pp) =.true.
!
    endsubroutine pencil_criteria_entropy
!***********************************************************************
    subroutine pencil_interdep_entropy(lpencil_in)
!
!  Interdependency among pencils from the Entropy 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)) lpencil_in(i_glnTT)=.true.
!
      if (lpencil_in(i_fpres)) then
        lpencil_in(i_glnrho)=.true.
        lpencil_in(i_glnTT)=.true.
        lpencil_in(i_TT)=.true.
      endif
!
    endsubroutine pencil_interdep_entropy
!*********************************************************************** 
    subroutine calc_pencils_entropy(f,p)
!
!  Calculate Entropy pencils.
!  Most basic pencils should come first, as others may depend on them.
!
!  20-11-04/anders: coded
! 
      use EquationOfState
      use Sub

      real, dimension (mx,my,mz,mfarray), intent (in) :: f
      type (pencil_case), intent (inout) :: p
      integer :: j
!
!  Mach Speed
!
      if (lpencil(i_Ma2)) p%Ma2=p%u2/p%cs2
!
!  Temperature advection
!  (Needs to be here because of lupw_lnTT)
!
      if (lpencil(i_uglnTT)) &
        call u_dot_grad(f,ilnTT,p%glnTT,p%uu,p%uglnTT,UPWIND=lupw_lnTT)
!
! fpres
!
      if (lpencil(i_fpres)) then
        do j=1,3
          p%fpres(:,j)=-gamma_m1*gamma1*(p%TT*p%glnrho(:,j) + p%glnTT(:,j))
        enddo
      endif
!
    endsubroutine calc_pencils_entropy
!**********************************************************************
    subroutine dss_dt(f,df,p)
!
!  calculate right hand side of entropy equation
!  heat condution is currently disabled until old stuff,
!  which in now in calc_heatcond, has been reinstalled.
!  DTT/Dt = -gamma_m1*TT*divu + gamma*cp1*rho1*RHS
!
!  13-dec-02/axel+tobi: adapted from entropy
!
      use Cdata
      use Diagnostics
      use Mpicomm
      use Sub
      use Viscosity, only: calc_viscous_heat
      use EquationOfState, only: gamma_m1
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
      real, dimension(nx) :: Hmax=0.
!
      intent(inout) :: f,p
      intent(out) :: df
!
!  identify module and boundary conditions
!
      if (headtt.or.ldebug) print*,'SOLVE dlnTT_dt'
      if (headtt) call identify_bcs('lnTT',ilnTT)
      if (headtt) print*,'dss_dt: lnTT,cs2=', p%lnTT(1), p%cs2(1)
!
!  sound speed squared
!
      if (headtt) print*,'dss_dt: cs20=',p%cs2(1)
!
!  ``cs2/dx^2'' for timestep
!
      if (lfirst.and.ldt) advec_cs2=p%cs2*dxyz_2
      if (headtt.or.ldebug) print*,'dss_dt: max(advec_cs2) =',maxval(advec_cs2)
!
!  subtract pressure gradient term in momentum equation
!
      if (lhydro.and.lpressuregradient_gas) &
         df(l1:l2,m,n,iux:iuz) = df(l1:l2,m,n,iux:iuz) + p%fpres(:,iux:iuz)
!
!  advection term
!
      if (ladvection_temperature) &
         df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) - p%uglnTT
!
!  Calculate viscous contribution to temperature
!
      if (lviscosity.and.lviscosity_heat) call calc_viscous_heat(df,p,Hmax)
!
!  Thermal conduction
!
      if (lheatc_chiconst) call calc_heatcond_constchi(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)
!
!  Need to add left-hand-side of the continuity equation (see manual)
!  Check this
      if (ldensity) &
        df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) - gamma_m1*p%TT*p%divu
!
!  Calculate entropy related diagnostics
!
      if (ldiagnos) then
        if (idiag_TTmax/=0) call max_mn_name(p%TT,idiag_TTmax)
        if (idiag_TTmin/=0) call max_mn_name(-p%TT,idiag_TTmin,lneg=.true.)
        if (idiag_TTm/=0)   call sum_mn_name(p%TT,idiag_TTm)
        if (idiag_ethm/=0)   call sum_mn_name(p%ee/p%rho1,idiag_ethm)
        if (idiag_ssm/=0)   call sum_mn_name(p%ss,idiag_ssm)
        if (idiag_dtc/=0) &
          call max_mn_name(sqrt(advec_cs2)/cdt,idiag_dtc,l_dt=.true.)
        if (idiag_eem/=0) call sum_mn_name(p%ee,idiag_eem)
        if (idiag_ppm/=0) call sum_mn_name(p%pp,idiag_ppm)
        if (idiag_csm/=0) call sum_mn_name(p%cs2,idiag_csm,lsqrt=.true.)
      endif
!
    endsubroutine dss_dt
!***********************************************************************
    subroutine single_polytrope(f)
!
! 04-aug-2007/dintrans: a simple polytrope with index mpoly0
!
      use Cdata
      use Gravity, only: gravz
      use EquationOfState, only: cs20, lnrho0, gamma, gamma_m1, get_cp1, &
                                 cs2bot, cs2top
!
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f
      real :: beta, zbot, ztop, cp1, T0, temp
!
!  beta is the (negative) temperature gradient
!  beta = -(g/cp) /[(1-1/gamma)*(m+1)]
!
      call get_cp1(cp1)
      beta=-cp1*gravz/(mpoly0+1.)*gamma/gamma_m1
      ztop=xyz0(3)+Lxyz(3)
      zbot=xyz0(3)
      T0=cs20/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))
        f(:,m,n,ilnTT)=temp
        f(:,m,n,ilnrho)=lnrho0+mpoly0*log(temp)-mpoly0*log(T0)
      enddo
      cs2bot=gamma_m1*(T0+beta*(ztop-zbot))
      cs2top=cs20
!
    endsubroutine single_polytrope
!***********************************************************************
    subroutine rad_equil(f)
!
! 16-mai-2007/tgastine+dintrans: compute the radiative and hydrostatic 
! equilibria for a given radiative profile (here a hole for the moment)
!
      use Cdata
      use Gravity, only: gravz
      use EquationOfState, only:cs20,lnrho0,cs2top,cs2bot,gamma,gamma_m1
!
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f
      real, dimension (mz) :: temp,lnrho
      real :: arg,hcond,dtemp,dlnrho
      real :: alp,sig,ecart
      integer :: i
!
      sig=7.
      ecart=0.4
      alp=(Kmax-Kmin)/(pi/2.+atan(sig*ecart**2))
      print*,'sig, ecart, alp=',sig, ecart, alp
      print*,'Kmin, Kmax, Fbot, Tbump=', Kmin, Kmax, Fbot, Tbump
!
! integrate from the top to the bottom
!
      temp(n2)=cs20/gamma_m1
      lnrho(n2)=lnrho0
      f(:,:,n2,ilnTT)=cs20/gamma_m1
      f(:,:,n2,ilnrho)=lnrho0
      do i=n2-1,n1,-1
        arg=sig*(temp(i+1)-Tbump-ecart)*(temp(i+1)-Tbump+ecart)
        hcond=Kmax+alp*(-pi/2.+atan(arg))
        dtemp=Fbot/hcond
        temp(i)=temp(i+1)+dz*dtemp
        dlnrho=2d0*(-gamma/gamma_m1*gravz-dtemp)/(7.d0/6.d0*temp(i)+5.d0/6.d0*temp(i+1))
        lnrho(i)=lnrho(i+1)+dz*dlnrho
        f(:,:,i,ilnTT)=temp(i)
        f(:,:,i,ilnrho)=lnrho(i)
      enddo
! initialize cs2bot by taking into account the new bottom value of temperature
! note: cs2top is already defined in eos_init by assuming cs2top=cs20
!      cs2bot=gamma_m1*temp(n1)
      print*,'cs2top, cs2bot=', cs2top, cs2bot
!
      if (lroot) then
        print*,'--> write the initial setup in data/proc0/setup.dat'
        open(unit=11,file=trim(directory)//'/setup.dat')
        write(11,'(4a14)') 'z','rho','temp','hcond'
        do i=n2,n1,-1
          arg=sig*(temp(i)-Tbump-ecart)*(temp(i)-Tbump+ecart)
          hcond=Kmax+alp*(-pi/2.+atan(arg))
          write(11,'(4e14.5)') z(i),exp(lnrho(i)),temp(i),hcond
        enddo
        close(11)
      endif
!
    endsubroutine rad_equil
!***********************************************************************
    subroutine calc_heatcond_constchi(df,p)
!
!  01-mar-07/dintrans: adapted from temperature_ionization
!
!  calculate chi*grad(rho*T*glnTT)/(rho*TT)
!           =chi*(g2.glnTT+g2lnTT) where g2=glnrho+glnTT
!
      use Diagnostics
      use EquationOfState, only: gamma
      use Sub

      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
      real, dimension (nx) :: g2
!
      call dot(p%glnTT+p%glnrho,p%glnTT,g2)
!
!  Add heat conduction to RHS of temperature equation
!
      df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + gamma*chi*(g2 + p%del2lnTT)
!
!  check maximum diffusion from thermal diffusion
!
      if (lfirst.and.ldt) then
        diffus_chi=diffus_chi+gamma*chi*dxyz_2
        if (ldiagnos.and.idiag_dtchi/=0) then
          call max_mn_name(diffus_chi/cdtv,idiag_dtchi,l_dt=.true.)
        endif
      endif

    endsubroutine calc_heatcond_constchi
!***********************************************************************
    subroutine calc_heatcond_constK(df,p)
!
!  calculate gamma*K/rho/cp*div(grad TT)= gamma*K/rho/cp*del2 TT
!  Be careful! here lnTT means TT...
!
      use Diagnostics
      use EquationOfState, only: gamma
      use Sub

      real, dimension(mx,my,mz,mvar) :: df
      type (pencil_case) :: p
      real, dimension(nx) :: hcond,chix
!
      hcond=hcond0
!
!  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*p%del2lnTT
!
!  check maximum diffusion from thermal diffusion
!
      if (lfirst.and.ldt) then
        diffus_chi=diffus_chi+gamma*chix*dxyz_2
        if (ldiagnos.and.idiag_dtchi/=0) then
          call max_mn_name(diffus_chi/cdtv,idiag_dtchi,l_dt=.true.)
        endif
      endif

    endsubroutine calc_heatcond_constK
!***********************************************************************
    subroutine calc_heatcond_arctan(df,p)
!
! 16-mai-2007/tgastine+dintrans: radiative diffusion with an arctan
!  profile for the conductivity
!  calculate gamma/rho*cp*div(K T*grad lnTT)=
!    gamma*K/rho*cp*(gradlnTT.gradlnTT + del2ln TT + gradlnTT.gradlnK)
!
      use Diagnostics
      use EquationOfState, only: gamma
      use HDF5IO, only: output_profile
      use Sub

      real, dimension(mx,my,mz,mvar) :: df
      type (pencil_case) :: p
      real, dimension(nx) :: arg,hcond,chiT,g2,chix
      real, dimension (nx,3) :: glnhcond=0.
      real :: alp,sig,ecart
!
! todo: must be defined in the namelist (used by rad_equil and _arctan...)
      sig=7.
      ecart=0.4
      alp=(Kmax-Kmin)/(pi/2.+atan(sig*ecart**2))
!
      arg=sig*(p%TT-Tbump-ecart)*(p%TT-Tbump+ecart)
      hcond=Kmax+alp*(-pi/2.+atan(arg))
      chiT=2./hcond*sig*alp*(p%TT-Tbump)/(1.+arg**2)  ! d(ln K)/dT
!
      glnhcond(:,3)=chiT*p%glnTT(:,3)
!      call dot(p%glnTT,p%glnTT,g1)
      call dot(p%glnTT,glnhcond,g2)
!
!  Write out hcond z-profile (during first time step only)
!
      if (m==m1) then
        call output_profile('hcond',(/z(n)/),(/hcond(1)/),'z', lsave_name=(n==n1))
        call output_profile('glnhcond',(/z(n)/),(/glnhcond(1,3)/),'z', lsave_name=(n==n1))
        call output_profile('K_T',(/z(n)/),(/chiT/),'z', lsave_name=(n==n1))
      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*(g1+g2+p%del2lnTT)
      df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + gamma*chix*(g2+p%del2lnTT)
!
!  check maximum diffusion from thermal diffusion
!
      if (lfirst.and.ldt) then
        diffus_chi=diffus_chi+gamma*chix*dxyz_2
        if (ldiagnos.and.idiag_dtchi/=0) then
          call max_mn_name(diffus_chi/cdtv,idiag_dtchi,l_dt=.true.)
        endif
      endif

    endsubroutine calc_heatcond_arctan
!***********************************************************************
    subroutine calc_heatcond(f,df,p)
!
!  12-Mar-2007/dintrans: coded
!  calculate gamma*K/rho*cp*div(T*grad lnTT)= 
!              gamma*K/rho*cp*(gradlnTT.gradln(hcond*TT) + del2ln TT)
!
!
      use Diagnostics
      use EquationOfState, only: gamma
      use Sub

      real, dimension(mx,my,mz,mfarray) :: f
      real, dimension(mx,my,mz,mvar) :: df
      type (pencil_case) :: p
      real, dimension(nx) :: g2,hcond,chix
      real, dimension (nx,3) :: glhc=0.,glnThcond
      integer :: i
      logical :: lwrite_hcond=.true.
      save :: lwrite_hcond
!
      if (lhcond_global) then
        hcond=f(l1:l2,m,n,iglobal_hcond)
        glhc(:,1)=f(l1:l2,m,n,iglobal_glhc)
      else
        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)
      endif
      if (lroot .and. lwrite_hcond) then
        open(1,file=trim(directory)//'/hcond.dat',position='append')
        write(1,'(3e14.5)') (rcyl_mn(i),hcond(i),glhc(i,1),i=1,nx)
        close(1)
        lwrite_hcond=.false.
      endif
!
      glnThcond = p%glnTT + glhc/spread(hcond,2,3)    ! grad ln(T*hcond)
      call dot(p%glnTT,glnThcond,g2)
!
!  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 + p%del2lnTT)
!
!  check maximum diffusion from thermal diffusion
!
      if (lfirst.and.ldt) then
        diffus_chi=diffus_chi+gamma*chix*dxyz_2
        if (ldiagnos.and.idiag_dtchi/=0) then
          call max_mn_name(diffus_chi/cdtv,idiag_dtchi,l_dt=.true.)
        endif
      endif

    endsubroutine calc_heatcond
!***********************************************************************
    subroutine read_entropy_init_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=entropy_init_pars, IOSTAT=iostat)
!
    endsubroutine read_entropy_init_pars
!***********************************************************************
    subroutine write_entropy_init_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=entropy_init_pars)
!
    endsubroutine write_entropy_init_pars
!***********************************************************************
    subroutine read_entropy_run_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=entropy_run_pars, IOSTAT=iostat)
!
    endsubroutine read_entropy_run_pars
!***********************************************************************
    subroutine write_entropy_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=entropy_run_pars)
!
    endsubroutine write_entropy_run_pars
!***********************************************************************
    subroutine rprint_entropy(lreset,lwrite)
!
!  reads and registers print parameters relevant to entropy
!
!   1-jun-02/axel: adapted from magnetic fields
!
      use Diagnostics
      use FArrayManager, only: farray_index_append
!
      integer :: iname
      logical :: lreset,lwr
      logical, optional :: lwrite
!
      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_TTmin=0; idiag_TTm=0
        idiag_ethm=0; idiag_ssm=0
        idiag_dtchi=0; idiag_dtc=0
        idiag_eem=0; idiag_ppm=0; idiag_csm=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),'TTmin',idiag_TTmin)
        call parse_name(iname,cname(iname),cform(iname),'TTm',idiag_TTm)
        call parse_name(iname,cname(iname),cform(iname),'ethm',idiag_ethm)
        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),'csm',idiag_csm)
      enddo
!
!  write column where which variable is stored
!
      if (lwr) then
        call farray_index_append('nname',nname)
        call farray_index_append('ilnTT',ilnTT)
        call farray_index_append('iss',iss)
        call farray_index_append('i_TTmax',idiag_TTmax)
        call farray_index_append('i_TTmin',idiag_TTmin)
        call farray_index_append('i_TTm',idiag_TTm)
        call farray_index_append('i_ethm',idiag_ethm)
        call farray_index_append('i_ssm',idiag_ssm)
        call farray_index_append('i_dtchi',idiag_dtchi)
        call farray_index_append('i_dtc',idiag_dtc)
        call farray_index_append('i_eem',idiag_eem)
        call farray_index_append('i_ppm',idiag_ppm)
        call farray_index_append('i_csm',idiag_csm)
      endif
!
    endsubroutine rprint_entropy
!***********************************************************************
    subroutine get_slices_entropy(f,slices)
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (slice_data) :: slices
!
      call keep_compiler_quiet(f)
      call keep_compiler_quiet(slices%ready)
!
    endsubroutine get_slices_entropy
!***********************************************************************
    subroutine ADI_constK(finit,f)
       
      use Cdata
      use Cparam
      use EquationOfState, only: gamma, gamma_m1, cs2bot, cs2top

      implicit none

      integer :: i,j
      real, dimension(mx,my,mz,mfarray) :: finit,f
      real, dimension(mx,mz) :: finter,source,rho
      real, dimension(nx)    :: a,b,c
      real, dimension(nz)    :: rhs,work
      real    :: aalpha, bbeta
!
      source=(f(:,4,:,ilnTT)-finit(:,4,:,ilnTT))/dt
      rho=exp(f(:,4,:,ilnrho))
      rho=rho/gamma
!
!  lignes en implicite
!
      do j=n1,n2
        a=-hcond0*dt/(2.*rho(l1:l2,j)*dx**2)
!
        b=1.+hcond0*dt/(rho(l1:l2,j)*dx**2)
!
        c=-hcond0*dt/(2.*rho(l1:l2,j)*dx**2) 
!
        rhs=finit(l1:l2,4,j,ilnTT)+hcond0*dt/(2.*rho(l1:l2,j)*dz**2) &
            *(finit(l1:l2,4,j+1,ilnTT)-2.*finit(l1:l2,4,j,ilnTT)+ &
            finit(l1:l2,4,j-1,ilnTT))+dt/2.*source(l1:l2,j)
!
        aalpha=c(nx)
        bbeta=a(1)
        c(nx)=0.
        a(1)=0.
        call cyclic(a,b,c,aalpha,bbeta,rhs,work,nx)
        finter(l1:l2,j)=work(1:nx)
      enddo
!
      call BC_CT(finter)
!
!  colonnes en implicite
!
      do i=l1,l2
        a=-hcond0*dt/(2.*rho(i,n1:n2)*dz**2)
!
        b=1.+hcond0*dt/(rho(i,n1:n2)*dz**2)
!
        c=-hcond0*dt/(2.*rho(i,n1:n2)*dz**2)
!
        rhs=finter(i,n1:n2)+hcond0*dt/(2.*rho(i,n1:n2)*dx**2) &
           *(finter(i+1,n1:n2)-2.*finter(i,n1:n2)+finter(i-1,n1:n2)) &
           +dt/2.*source(i,n1:n2)
!
        c(nz)=0.
        a(1)=0.
        b(1)=1.
        b(nz)=1.
        c(1)=0.
        a(nz)=0.
        rhs(1)=cs2bot/gamma_m1
        rhs(nx)=cs2top/gamma_m1
        call tridag(a,b,c,rhs,work,nz)
        f(i,4,n1:n2,ilnTT)=work(1:nz)
      enddo
!
      call BC_CT(f(:,4,:,ilnTT))
!
      endsubroutine ADI_constK
!**************************************************************
      subroutine ADI_Kprof(finit,f)
       
      use Cdata
      use Cparam
      use EquationOfState, only: gamma,cs2bot,cs2top

      implicit none

      integer :: i,j
      real    :: aalpha, bbeta
      real, dimension(mx,my,mz,mfarray) :: finit,f
      real, dimension(mx,mz) :: source,rho,chiprof,dchi,valinter,val
      real, dimension(nx)    :: a,b,c
      real, dimension(nz)    :: rhs,work

      source=(f(:,4,:,ilnTT)-finit(:,4,:,ilnTT))/dt
      rho=exp(f(:,4,:,ilnrho))
      rho=rho/gamma
      call hcond_ADI(f,chiprof,dchi)
!
!  lignes en implicite
!
      do j=n1,n2
       a=-dt/(4d0*rho(l1:l2,j)*dx**2)*(dchi(l1-1:l2-1,j) &
         *(finit(l1-1:l2-1,4,j,ilnTT)-finit(l1:l2,4,j,ilnTT)) &
         +chiprof(l1-1:l2-1,j)+chiprof(l1:l2,j))
!
       b=1d0+dt/(4d0*rho(l1:l2,j)*dx**2)*(dchi(l1:l2,j) &
         *(2d0*finit(l1:l2,4,j,ilnTT)-finit(l1-1:l2-1,4,j,ilnTT) &
         -finit(l1+1:l2+1,4,j,ilnTT))+2d0*chiprof(l1:l2,j) &
       +chiprof(l1+1:l2+1,j)+chiprof(l1-1:l2-1,j))
!
       c=-dt/(4d0*rho(l1:l2,j)*dx**2)*(dchi(l1+1:l2+1,j) &
          *(finit(l1+1:l2+1,4,j,ilnTT)-finit(l1:l2,4,j,ilnTT)) &
       +chiprof(l1:l2,j)+chiprof(l1+1:l2+1,j))
!
       rhs=1d0/(2d0*rho(l1:l2,j)*dz**2)*((chiprof(l1:l2,j+1) &
           +chiprof(l1:l2,j))*(finit(l1:l2,4,j+1,ilnTT)-finit(l1:l2,4,j,ilnTT))&
           -(chiprof(l1:l2,j)+chiprof(l1:l2,j-1)) &
       *(finit(l1:l2,4,j,ilnTT)-finit(l1:l2,4,j-1,ilnTT)))
!
       rhs=rhs+1d0/(2d0*rho(l1:l2,j)*dx**2)*((chiprof(l1+1:l2+1,j) &
         +chiprof(l1:l2,j))*(finit(l1+1:l2+1,4,j,ilnTT)-finit(l1:l2,4,j,ilnTT))&
           -(chiprof(l1:l2,j)+chiprof(l1-1:l2-1,j)) &
           *(finit(l1:l2,4,j,ilnTT)-finit(l1-1:l2-1,4,j,ilnTT)))
!
       aalpha=c(nx)
       bbeta=a(1)
       c(nx)=0d0
       a(1)=0d0
       call cyclic(a,b,c,aalpha,bbeta,rhs,work,nx)
       valinter(l1:l2,j)=work(1:nx)
      enddo
!
!  colonnes en implicite
!
      do i=l1,l2
       a=-dt/(4d0*rho(i,n1:n2)*dz**2)*(dchi(i,n1-1:n2-1) &
         *(finit(i,4,n1-1:n2-1,ilnTT)-finit(i,4,n1:n2,ilnTT))&
         +chiprof(i,n1-1:n2-1)+chiprof(i,n1:n2))
!
       b=1d0+dt/(4d0*rho(i,n1:n2)*dz**2)*(dchi(i,n1:n2)* &
         (2d0*finit(i,4,n1:n2,ilnTT)-finit(i,4,n1-1:n2-1,ilnTT) &
         -finit(i,4,n1+1:n2+1,ilnTT))+2d0*chiprof(i,n1:n2) &
         +chiprof(i,n1+1:n2+1)+chiprof(i,n1-1:n2-1))
!
       c=-dt/(4d0*rho(i,n1:n2)*dz**2)*(dchi(i,n1+1:n2+1) &
         *(finit(i,4,n1+1:n2+1,ilnTT)-finit(i,4,n1:n2,ilnTT))&
         +chiprof(i,n1:n2)+chiprof(i,n1+1:n2+1))
!
       rhs=valinter(i,n1:n2)
!
       c(nz)=0d0
       a(1)=0d0
! Constant flux at the bottom
!       b(1)=-1.
!       c(1)=0.
!       c(1)=1.
!       rhs(1)=0.
! Constant temperature at the top
       b(nx)=1.
       a(nx)=0.
       rhs(nz)=0.
!
       call tridag(a,b,c,rhs,work,nz)
       val(i,n1:n2)=work(1:nz)
      enddo
!
      f(:,4,:,ilnTT)=finit(:,4,:,ilnTT)+dt*val+dt*source
!
!      f(:,:,n1,ilnTT)=cs2bot/(gamma-1d0)
      f(:,:,n2,ilnTT)=cs2top/(gamma-1d0)
!      call BC_CT(f)
      call BC_flux(f,chiprof)
      call hcond_ADI(f,chiprof,dchi)
!
      endsubroutine ADI_Kprof
!**************************************************************
      subroutine BC_CT(f_2d)

      implicit none

      real, dimension(mx,mz) :: f_2d

! z-direction
      f_2d(:,n1-1)=2.*f_2d(:,n1)-f_2d(:,n1+1)
      f_2d(:,n2+1)=2.*f_2d(:,n2)-f_2d(:,n2-1)
! x-direction
      f_2d(1:l1-1,:)=f_2d(l2i:l2,:)
      f_2d(l2+1:mx,:)=f_2d(l1:l1i,:)

      endsubroutine BC_CT
!**************************************************************
      subroutine BC_flux(f,chiprof)

      real, dimension(mx,my,mz,mfarray) :: f
      real, dimension(mx,mz) :: chiprof
      integer :: i

      do i=1,nghost
        f(:,:,n1-i,ilnTT)=f(:,:,n1+i,ilnTT)+2*i*dz*Fbot/chiprof(10,n1+i)
      enddo
      f(:,:,n2+1,ilnTT)=2*f(:,:,n2,ilnTT)-f(:,:,n2-1,ilnTT)
      ! x-direction
      f(1:l1-1,:,:,ilnTT)=f(l2i:l2,:,:,ilnTT)
      f(l2+1:mx,:,:,ilnTT)=f(l1:l1i,:,:,ilnTT)

      endsubroutine BC_flux
!**************************************************************
      subroutine hcond_ADI(f,chiprof,dchi)

      use Sub, only: write_zprof
      real, dimension(mx,my,mz,mfarray) :: f
      real , dimension(mx,mz) :: chiprof, dchi
!      double precision :: chi0
      real, dimension(mx,mz) :: arg
      real :: alp,sig,ecart
!
!      chi0=1d-3
!      chiprof=chi0
!      dchi=0d0
      sig=7.
      ecart=0.4
      alp=(Kmax-Kmin)/(pi/2.+atan(sig*ecart**2))
      arg=sig*(f(:,4,:,ilnTT)-Tbump-ecart)*(f(:,4,:,ilnTT)-Tbump+ecart)
      chiprof=Kmax+alp*(-pi/2.+atan(arg))
      dchi=2d0*alp/(1d0+arg**2)*sig*(f(:,4,:,ilnTT)-Tbump)
!
!      call write_zprof('hcond',chiprof(10,n1:n2))
!      call write_zprof('dhcond/dT',dchi(10,n1:n2))

      endsubroutine hcond_ADI
!**************************************************************
      subroutine tridag(a,b,c,r,u,n)

      integer :: j,n
      integer, parameter :: NMAX=500
      real    :: bet
      real, dimension(n) :: a,b,c,r,u
      real, dimension(NMAX) :: gam
!
      if (b(1).eq.0.) pause 'tridag: rewrite equations'
      bet=b(1)
      u(1)=r(1)/bet
      do 11 j=2,n
        gam(j)=c(j-1)/bet
        bet=b(j)-a(j)*gam(j)
        if (bet.eq.0.) pause 'tridag failed'
        u(j)=(r(j)-a(j)*u(j-1))/bet
11    continue
      do 12 j=n-1,1,-1
        u(j)=u(j)-gam(j+1)*u(j+1)
12    continue
!
      return
      endsubroutine tridag
!**************************************************************
      subroutine cyclic(a,b,c,alpha,beta,r,x,n)
!
      implicit none
!
      integer :: i,n
      integer, parameter    :: NMAX=500
      real    :: alpha, beta,gamma,fact      
      real, dimension(n)    :: a,b,c,r,x
      real, dimension(NMAX) :: bb,u,z
!
      if (n.le.2)pause 'n too small in cyclic'
      if (n.gt.NMAX)pause 'NMAX too small in cyclic'
      gamma=-b(1)
      bb(1)=b(1)-gamma
      bb(n)=b(n)-alpha*beta/gamma
      do 11 i=2,n-1
        bb(i)=b(i)
11    continue
      call tridag(a,bb,c,r,x,n)
      u(1)=gamma
      u(n)=alpha
      do 12 i=2,n-1
        u(i)=0.
12    continue
      call tridag(a,bb,c,u,z,n)
      fact=(x(1)+beta*x(n)/gamma)/(1.+z(1)+beta*z(n)/gamma)
      do 13 i=1,n
        x(i)=x(i)-fact*z(i)
13    continue
!
      return
      endsubroutine cyclic
!***************************************************************
    subroutine calc_heatcond_ADI(finit,f)
!
      implicit none
!
      real, dimension(mx,my,mz,mfarray) :: finit,f
!
      call keep_compiler_quiet(finit)
      call keep_compiler_quiet(f)
!
    endsubroutine calc_heatcond_ADI
!***********************************************************************
endmodule Entropy