! $Id$
!
! MODULE_DOC: Calculates pressure gradient term for
! MODULE_DOC: polytropic equation of state $p=\text{const}\rho^{\Gamma}$.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lentropy = .false.
! CPARAM logical, parameter :: ltemperature = .false.
! CPARAM logical, parameter :: lthermal_energy = .false.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 0
!
! PENCILS PROVIDED Ma2; fpres(3); tcond; sglnTT(3)
! PENCILS PROVIDED uglnTT; advec_cs2
!
!***************************************************************
module Energy
!
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages
!
  implicit none
!
  include 'energy.h'
!
  logical, pointer :: lpressuregradient_gas
  logical :: lviscosity_heat=.false.
  logical, pointer :: lffree, lrelativistic_eos, lconservative
  real, pointer :: profx_ffree(:),profy_ffree(:),profz_ffree(:)
!
  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_ugradpm=0
  integer :: idiag_thermalpressure=0
  integer :: idiag_ethm=0       ! DIAG_DOC: $\left<\varrho e\right>$
                                ! DIAG_DOC:   \quad(mean thermal
                                ! DIAG_DOC:   [=internal] energy)
  integer :: idiag_pdivum=0     ! DIAG_DOC: $\left<p\nabla\uv\right>$
  integer :: idiag_ufpresm=0
  integer :: idiag_csm=0        ! DIAG_DOC: $\left<c_{\rm s}\right>$
!
  real, dimension(:), pointer :: beta_glnrho_scaled
  real :: gamma
!
  contains
!***********************************************************************
    subroutine register_energy
!
!  No energy equation is being solved; use polytropic equation of state.
!
!  28-mar-02/axel: dummy routine, adapted from entropy.f of 6-nov-01.
!
      use SharedVariables, only: put_shared_variable
!
!  Identify version number.
!
      if (lroot) call svn_id( &
          "$Id$")
!
      call put_shared_variable('lviscosity_heat',lviscosity_heat,caller='register_energy')
!
    endsubroutine register_energy
!***********************************************************************
    subroutine initialize_energy(f)
!
!  Perform any post-parameter-read initialization i.e. calculate derived
!  parameters.
!
!  24-nov-02/tony: coded
!
      use EquationOfState, only: cs0, select_eos_variable, get_gamma_etc
      use SharedVariables, only: get_shared_variable
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      real :: cp
!
!  The following only makes sense if leos=.true.
!
      if (leos) then
        call get_gamma_etc(gamma,cp)
      endif
!
!  Tell the equation of state that we're here and what f variable we use.
!
      if (llocal_iso) then
        if (lroot) call warning('initialize_energy', &
           'llocal_iso=T. Make sure you have the appropriate INITIAL_CONDITION in Makefile.local')
        call select_eos_variable('cs2',-2) !special local isothermal
      elseif (leos) then
        if (gamma == 1.) then
          call select_eos_variable('cs2',-1) !isothermal
        else
          call select_eos_variable('ss',-1) !isentropic => polytropic
        endif
      endif
!
!  Logical variable lpressuregradient_gas shared with hydro modules.
!
      call get_shared_variable('lpressuregradient_gas',lpressuregradient_gas)
!
      call get_shared_variable('beta_glnrho_scaled',beta_glnrho_scaled)
!
      if (ldensity) then
!
!  Check if we are solving the relativistic eos equations.
!  In that case we'd need to get lrelativistic_eos from density.
!
        call get_shared_variable('lrelativistic_eos', lrelativistic_eos)
!
!  Check that cs0 is set correctly when lrelativistic_eos=.true.
!
        if (lrelativistic_eos) then
          if (abs(cs0**2-onethird)>0.01) &
            call warning('initialize_energy','consider putting cs0=1/sqrt(3) for relativistic EoS')
        endif
!
! Check if we are solving the force-free equations in parts of domain
!
        call get_shared_variable('lffree',lffree)
        if (lffree) then
          call get_shared_variable('profx_ffree',profx_ffree,caller='initialize_energy')
          call get_shared_variable('profy_ffree',profy_ffree,caller='initialize_energy')
          call get_shared_variable('profz_ffree',profz_ffree,caller='initialize_energy')
        endif
      else
        allocate(lrelativistic_eos)
        lrelativistic_eos=.false.
      endif
!
!  Check if we are solving for relativistic bulk motions, not just EoS.
!
      if (lhydro.and..not.lhydro_potential.and.iphiuu==0) then
        call get_shared_variable('lconservative', lconservative)
      else
        allocate(lconservative)
        lconservative=.false.
      endif
!
      TT_spec=.false.; ss_spec=.false.

      call keep_compiler_quiet(f)
!
    endsubroutine initialize_energy
!***********************************************************************
    subroutine update_char_vel_energy(f)
!
!  Updates characteristic velocity for slope-limited diffusion.
!
!  25-sep-15/MR+joern: coded
!
      use EquationOfState, only: cs20
!
      real, dimension(mx,my,mz,mfarray), intent(INOUT) :: f
!
      if (lslope_limit_diff) f(2:mx-2,2:my-2,2:mz-2,iFF_char_c) &
                             = max(f(2:mx-2,2:my-2,2:mz-2,iFF_char_c), w_sldchar_ene*sqrt(cs20))
!
    endsubroutine update_char_vel_energy
!***********************************************************************
    subroutine init_energy(f)
!
!  Initialise energy; called from start.f90.
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      call keep_compiler_quiet(f)
!
    endsubroutine init_energy
!***********************************************************************
    subroutine pencil_criteria_energy
!
!  All pencils that the Energy module depends on are specified here.
!
!  20-11-04/anders: coded
!
      if (lhydro.and.lpressuregradient_gas.and..not.lconservative) lpenc_requested(i_fpres)=.true.
      if (leos.and.ldensity.and.lhydro.and.ldt) lpenc_requested(i_cs2)=.true.
      if (any(beta_glnrho_scaled /= 0.)) lpenc_requested(i_cs2)=.true.
!
      if (idiag_ugradpm/=0) then
        lpenc_diagnos(i_rho)=.true.
        lpenc_diagnos(i_uglnrho)=.true.
      endif
!
      if (idiag_thermalpressure/=0) then
        lpenc_diagnos(i_rho)=.true.
        lpenc_diagnos(i_cs2)=.true.
        lpenc_diagnos(i_rcyl_mn)=.true.
      endif
!
      if (idiag_ethm/=0) then
        lpenc_diagnos(i_rho)=.true.
        lpenc_diagnos(i_ee)=.true.
      endif
!
      if (idiag_pdivum/=0) then
        lpenc_diagnos(i_pp)=.true.
        lpenc_diagnos(i_divu)=.true.
      endif
!
      if (idiag_csm/=0) lpenc_diagnos(i_cs2)=.true.
!
    endsubroutine pencil_criteria_energy
!***********************************************************************
    subroutine pencil_interdep_energy(lpencil_in)
!
!  Interdependency among pencils from the Energy module is specified here.
!
!  20-nov-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_fpres)) then
        lpencil_in(i_cs2)=.true.
        if (lstratz) then
          lpencil_in(i_glnrhos)=.true.
        else
          if (.not.lconservative) then
            lpencil_in(i_glnrho)=.true.
          endif
        endif
        if (llocal_iso)  lpencil_in(i_glnTT)=.true.
      endif
!
!  The following only makes sense if leos is true.
!
      if (leos) then
        if (lpencil_in(i_TT1) .and. gamma/=1.) lpencil_in(i_cs2)=.true.
        if (lpencil_in(i_cs2) .and. gamma/=1.) lpencil_in(i_lnrho)=.true.
      endif
!
    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-nov-04/anders: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (pencil_case) :: p
!
      integer :: j
!
      intent(in) :: f
      intent(inout) :: p
! Ma2
      if (lpencil(i_Ma2)) p%Ma2=p%u2/p%cs2
!
!  fpres (=pressure gradient force)
!
      if (lpencil(i_fpres)) then
        if (lstratz) then
          p%fpres = -spread(p%cs2,2,3) * p%glnrhos
        else
          do j=1,3
            if (llocal_iso) then
              p%fpres(:,j)=-p%cs2*(p%glnrho(:,j)+p%glnTT(:,j))
            else
!
!  The relativistic EoS works ok even if cs2 is not 1/3, but
!  it may still be a good idea to put cs0=1/sqrt(3)=0.57735
!
              if (ldensity.and.lrelativistic_eos) then
                !if (.not.lconservative) p%fpres(:,j)=-.75*p%cs2*p%glnrho(:,j)
                if (.not.lconservative) p%fpres(:,j)=-p%cs2/(1 + p%cs2)*p%glnrho(:,j)
              else
                p%fpres(:,j)=-p%cs2*p%glnrho(:,j)
              endif
            endif
!
!  multiply previous p%fpres pencil with profiles
!
            if (ldensity) then
              if (lffree) p%fpres(:,j)=p%fpres(:,j)*profx_ffree*profy_ffree(m)*profz_ffree(n)
            endif
          enddo
        endif
      endif
!
! tcond (dummy)
!
      if (lpencil(i_tcond)) p%tcond=0.
! sglnTT (dummy)
      if (lpencil(i_sglnTT)) p%sglnTT=0.
!
!  ``cs2/dx^2'' for timestep - but only if we are evolving hydrodynamics.
!
      if (lupdate_courant_dt) then
        if (leos.and.ldensity.and.lhydro) then
          p%advec_cs2=p%cs2*dxyz_2
          if (headtt.or.ldebug) print*, 'calc_pencils_energy: max(advec_cs2) =', maxval(p%advec_cs2)
        endif
      endif
!
      call keep_compiler_quiet(f)
!
    endsubroutine calc_pencils_energy
!***********************************************************************
    subroutine energy_before_boundary(f)
!
!  03-apr-20/joern: restructured and fixed slope-limited diffusion
!
      use EquationOfState, only: cs0
!
      real, dimension (mx,my,mz,mfarray), intent(INOUT) :: f
!
!    Slope limited diffusion: update characteristic speed
!    Not staggered yet
!
      if (lslope_limit_diff .and. llast .and. ldensity) then
!      if (lslope_limit_diff) then
        do m=1,my
        do n=1,mz
          f(:,m,n,isld_char)=f(:,m,n,isld_char)+w_sldchar_ene*cs0
        enddo
        enddo
      endif
!
    endsubroutine energy_before_boundary
!***********************************************************************
    subroutine energy_after_boundary(f)
!
!  Dummy routine.
!
      real, dimension (mx,my,mz,mfarray), intent(IN) :: f
!
      call keep_compiler_quiet(f)

    endsubroutine energy_after_boundary
!***********************************************************************
    subroutine denergy_dt(f,df,p)
!
!  Calculate pressure gradient term for isothermal/polytropic equation
!  of state.
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
!
      integer :: j
!
      intent(in) :: f,p
      intent(inout) :: df
!
!  Add isothermal/polytropic pressure term in momentum equation.
!
      if (lhydro.and.lpressuregradient_gas.and..not.lconservative &
          .and.(.not.lhydro_potential)) then

        df(l1:l2,m,n,iux:iuz)=df(l1:l2,m,n,iux:iuz)+p%fpres
!
!  Add pressure force from global density gradient.
!
        if (any(beta_glnrho_scaled /= 0.)) then
          if (headtt) print*, 'denergy_dt: adding global pressure gradient force'
          do j=1,3
            df(l1:l2,m,n,iux-1+j) = df(l1:l2,m,n,iux-1+j) - p%cs2*beta_glnrho_scaled(j)
          enddo
        endif
      endif

      if (lupdate_courant_dt.and.leos.and.ldensity.and.lhydro) advec_cs2 = p%advec_cs2

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

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

      real, dimension(nx) :: ufpres
      integer :: i
!
      call keep_compiler_quiet(f)
!
!  Calculate energy related diagnostics.
!
      if (ldiagnos) then
        if (ldt) then
          if (idiag_dtc/=0) call max_mn_name(sqrt(p%advec_cs2)/cdt,idiag_dtc,l_dt=.true.)
        endif
        if (idiag_ugradpm/=0) call sum_mn_name(p%rho*p%cs2*p%uglnrho,idiag_ugradpm)
        if (idiag_thermalpressure/=0) call sum_lim_mn_name(p%rho*p%cs2,idiag_thermalpressure,p)
        if (idiag_ethm/=0) call sum_mn_name(p%rho*p%ee,idiag_ethm)
        if (idiag_pdivum/=0) call sum_mn_name(p%pp*p%divu,idiag_pdivum)
        call sum_mn_name(p%cs2,idiag_csm,lsqrt=.true.)
        if (idiag_ufpresm/=0) then
          ufpres=0
          do i = 1, 3
            ufpres=ufpres+p%uu(:,i)*p%fpres(:,i)
          enddo
          call sum_mn_name(p%rho*ufpres,idiag_ufpresm)
        endif

      endif

    endsubroutine calc_diagnostics_energy
!***********************************************************************
    subroutine get_slices_energy(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_energy
!***********************************************************************
    subroutine fill_farray_pressure(f)
!
!  18-feb-10/anders: dummy
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      call keep_compiler_quiet(f)
!
    endsubroutine fill_farray_pressure
!***********************************************************************
    subroutine impose_energy_floor(f)
!
!  Dummy subroutine
!
      real, dimension(mx,my,mz,mfarray) :: f
!
      call keep_compiler_quiet(f)
!
    endsubroutine impose_energy_floor
!***********************************************************************
    subroutine dynamical_thermal_diffusion(uc)
!
!  Dummy subroutine
!
      real, intent(in) :: uc
!
      call keep_compiler_quiet(uc)
!
    endsubroutine dynamical_thermal_diffusion
!***********************************************************************
    subroutine read_energy_init_pars(iostat)
!
      integer, intent(out) :: iostat
!
      iostat = 0
!
    endsubroutine read_energy_init_pars
!***********************************************************************
    subroutine write_energy_init_pars(unit)
!
      integer, intent(in) :: unit
!
      call keep_compiler_quiet(unit)
!
    endsubroutine write_energy_init_pars
!***********************************************************************
    subroutine read_energy_run_pars(iostat)
!
      integer, intent(out) :: iostat
!
      iostat = 0
!
    endsubroutine read_energy_run_pars
!***********************************************************************
    subroutine write_energy_run_pars(unit)
!
      integer, intent(in) :: unit
!
      call keep_compiler_quiet(unit)
!
    endsubroutine write_energy_run_pars
!***********************************************************************
    subroutine rprint_energy(lreset,lwrite)
!
!  Reads and registers print parameters relevant to energy.
!
      use Diagnostics, only: parse_name
!
      integer :: iname
      logical :: lreset
      logical, optional :: lwrite
!
!  Reset everything in case of reset
!  (this needs to be consistent with what is defined above!)
!
      if (lreset) then
        idiag_dtc=0; idiag_ugradpm=0; idiag_thermalpressure=0; idiag_ethm=0;
        idiag_pdivum=0; idiag_ufpresm=0; idiag_csm=0
      endif
!
      do iname=1,nname
        call parse_name(iname,cname(iname),cform(iname),'dtc',idiag_dtc)
        call parse_name(iname,cname(iname),cform(iname),'ugradpm',idiag_ugradpm)
        call parse_name(iname,cname(iname),cform(iname),'TTp',idiag_thermalpressure)
        call parse_name(iname,cname(iname),cform(iname),'ethm',idiag_ethm)
        call parse_name(iname,cname(iname),cform(iname),'pdivum',idiag_pdivum)
        call parse_name(iname,cname(iname),cform(iname),'ufpresm',idiag_ufpresm)
        call parse_name(iname,cname(iname),cform(iname),'csm',idiag_csm)
      enddo

      call keep_compiler_quiet(present(lwrite))
!
    endsubroutine rprint_energy
!***********************************************************************
    subroutine energy_after_timestep(f,df,dtsub)
!
      real, dimension(mx,my,mz,mfarray) :: f
      real, dimension(mx,my,mz,mvar) :: df
      real :: dtsub
!
      call keep_compiler_quiet(f,df)
      call keep_compiler_quiet(dtsub)
!
    endsubroutine energy_after_timestep
!***********************************************************************
    subroutine split_update_energy(f)
!
!  Dummy subroutine
!
      real, dimension(mx,my,mz,mfarray), intent(inout) :: f
!
      call keep_compiler_quiet(f)
!
    endsubroutine split_update_energy
!***********************************************************************
    subroutine expand_shands_energy
!
!  Dummy
!
    endsubroutine expand_shands_energy
!***********************************************************************
    subroutine heatcond_TT_2d(TT, hcond, dhcond)
!
! dummy
!
      implicit none
!
      real, dimension(:,:), intent(in) :: TT
      real, dimension(:,:), intent(out) :: hcond
      real, dimension(:,:), optional :: dhcond

      call keep_compiler_quiet(TT,hcond,dhcond)

    endsubroutine heatcond_TT_2d
!***********************************************************************
    subroutine heatcond_TT_1d(TT, hcond, dhcond)
!
! dummy
!
      implicit none
!
      real, dimension(:), intent(in) :: TT
      real, dimension(:), intent(out) :: hcond
      real, dimension(:), optional :: dhcond

      call keep_compiler_quiet(TT,hcond,dhcond)

    endsubroutine heatcond_TT_1d
!***********************************************************************
    subroutine heatcond_TT_0d(TT, hcond, dhcond)
!
! dummy
!
      implicit none
!
      real, intent(in) :: TT
      real, intent(out) :: hcond
      real, optional :: dhcond

      call keep_compiler_quiet(TT,hcond,dhcond)

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

    use Syscalls, only: copy_addr

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

      call keep_compiler_quiet(p_par)
      call copy_addr(lviscosity_heat,p_par(1)) ! bool
      call copy_addr(w_sldchar_ene,p_par(2))

    endsubroutine pushpars2c
!***********************************************************************
endmodule Energy