mod_hd_phys.t

Go to the documentation of this file.

!> Hydrodynamics physics module
module mod_hd_phys
  use mod_thermal_conduction, only: tc_fluid
  use mod_radiative_cooling, only: rc_fluid
  use mod_thermal_emission, only: te_fluid
  use mod_physics
  use mod_comm_lib, only: mpistop
  implicit none
  private
 
  !> Whether an energy equation is used
  logical, public, protected              :: hd_energy = .true.
 
  !> Whether thermal conduction is added
  logical, public, protected              :: hd_thermal_conduction = .false.
  type(tc_fluid), allocatable, public :: tc_fl
  type(te_fluid), allocatable, public :: te_fl_hd
 
  !> Whether radiative cooling is added
  logical, public, protected              :: hd_radiative_cooling = .false.
  type(rc_fluid), allocatable, public :: rc_fl
 
  !> Whether dust is added
  logical, public, protected              :: hd_dust = .false.
 
  !> Whether dust is added using and implicit update in IMEX
  logical, public, protected              :: hd_dust_implicit = .false.
 
  !> Whether radiation-gas interaction is handled using flux limited diffusion
  logical, public, protected              :: hd_radiation_fld = .false.
 
  !> Formalism to treat radiation: either fld or afld (anisotropic fld)
  character(len=8), public :: hd_radiation_fld_formalism = 'fld'
 
  !> Whether viscosity is added
  logical, public, protected              :: hd_viscosity = .false.
 
  !> Whether gravity is added
  logical, public, protected              :: hd_gravity = .false.
 
  !> Whether particles module is added
  logical, public, protected              :: hd_particles = .false.
 
  !> Whether rotating frame is activated
  logical, public, protected              :: hd_rotating_frame = .false.
 
  !> Whether CAK radiation line force is activated
  logical, public, protected              :: hd_cak_force = .false.
 
  !> Number of tracer species
  integer, public, protected              :: hd_n_tracer = 0
 
  !> Whether plasma is partially ionized
  logical, public, protected              :: hd_partial_ionization = .false.
 
  !> Index of the density (in the w array)
  integer, public, protected              :: rho_
 
  !> Indices of the momentum density
  integer, allocatable, public, protected :: mom(:)
 
  !> Indices of the momentum density for the form of better vectorization
  integer, public, protected              :: ^c&m^C_
 
  !> Indices of the tracers
  integer, allocatable, public, protected :: tracer(:)
 
  !> Index of the energy density (-1 if not present)
  integer, public, protected              :: e_
 
  !> Index of the gas pressure (-1 if not present) should equal e_
  integer, public, protected              :: p_
 
  !> Index of the radiation energy (when fld active)
  integer, public, protected              :: r_e
 
  !> Indices of temperature
  integer, public, protected :: te_
 
  !> Index of the cutoff temperature for the TRAC method
  integer, public, protected              :: tcoff_
 
  !> The adiabatic index
  double precision, public                :: hd_gamma = 5.d0/3.0d0
 
  !> gamma minus one and its inverse
  double precision :: gamma_1, inv_gamma_1
 
  !> The adiabatic constant
  double precision, public                :: hd_adiab = 1.0d0
 
  !> The small_est allowed energy
  double precision, protected             :: small_e
 
  !> The smallest allowed radiation energy (when fld active)
  double precision, public, protected     :: small_r_e
 
  !> Whether TRAC method is used
  logical, public, protected              :: hd_trac = .false.
  integer, public, protected              :: hd_trac_type = 1
 
  !> Helium abundance over Hydrogen
  double precision, public, protected  :: he_abundance=0.1d0
  !> Ionization fraction of H
  !> H_ion_fr = H+/(H+ + H)
  double precision, public, protected  :: h_ion_fr=1d0
  !> Ionization fraction of He
  !> He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
  double precision, public, protected  :: he_ion_fr=1d0
  !> Ratio of number He2+ / number He+ + He2+
  !> He_ion_fr2 = He2+/(He2+ + He+)
  double precision, public, protected  :: he_ion_fr2=1d0
  ! used for eq of state when it is not defined by units,
  ! the units do not contain terms related to ionization fraction
  ! and it is p = RR * rho * T
  double precision, public, protected  :: rr=1d0
  ! remove the below flag  and assume default value = .false.
  ! when eq state properly implemented everywhere
  ! and not anymore through units
  logical, public, protected :: eq_state_units = .true.
 
  procedure(sub_get_pthermal), pointer :: hd_get_rfactor   => null()
  ! Public methods
  public :: hd_phys_init
  public :: hd_kin_en
  public :: hd_get_pthermal
  public :: hd_get_csound2
  public :: hd_to_conserved
  public :: hd_to_primitive
  public :: hd_check_params
  public :: hd_check_w
  public :: hd_e_to_ei
  public :: hd_ei_to_e
  public :: hd_get_rfactor
  ! Begin: following relevant for radiative hydro using FLD
  ! first three are local and of interest for mod_usr applications
  public :: hd_get_pradiation
  public :: hd_get_ptot
  public :: hd_get_trad
  ! the following used in FLD modules
  !    as pointer phys_get_tgas
  public :: hd_get_temperature_from_etot
  !    as pointer phys_set_mg_bounds
  public :: hd_set_mg_bounds
  ! End: following relevant for radiative hydro using FLD
 
contains
 
  !> Read this module's parameters from a file
  subroutine hd_read_params(files)
    use mod_global_parameters
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /hd_list/ hd_energy, hd_n_tracer, hd_gamma, hd_adiab, &
    hd_dust, hd_dust_implicit, hd_thermal_conduction, hd_radiative_cooling, hd_viscosity, &
    hd_gravity, he_abundance,h_ion_fr, he_ion_fr, he_ion_fr2, eq_state_units, &
    si_unit, hd_particles, hd_rotating_frame, hd_trac, &
    hd_trac_type, hd_cak_force, hd_partial_ionization, &
    hd_radiation_fld,hd_radiation_fld_formalism
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status="old")
       read(unitpar, hd_list, end=111)
111    close(unitpar)
    end do
 
  end subroutine hd_read_params
 
  !> Write this module's parameters to a snapsoht
  subroutine hd_write_info(fh)
    use mod_global_parameters
    integer, intent(in)                 :: fh
    integer, parameter                  :: n_par = 1
    double precision                    :: values(n_par)
    character(len=name_len)             :: names(n_par)
    integer, dimension(MPI_STATUS_SIZE) :: st
    integer                             :: er
 
    call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
 
    names(1) = "gamma"
    values(1) = hd_gamma
    call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
    call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
  end subroutine hd_write_info
 
  !> Initialize the module
  subroutine hd_phys_init()
    use mod_global_parameters
    use mod_thermal_conduction
    use mod_radiative_cooling
    use mod_dust, only: dust_init
    use mod_viscosity, only: viscosity_init
    use mod_gravity, only: gravity_init
    use mod_particles, only: particles_init
    use mod_rotating_frame, only:rotating_frame_init
    use mod_cak_force, only: cak_init
    use mod_supertimestepping, only: sts_init, add_sts_method,&
            set_conversion_methods_to_head, set_error_handling_to_head
    use mod_ionization_degree
    use mod_usr_methods, only: usr_rfactor
    use mod_fld
    use mod_afld
 
    integer :: itr, idir
 
    call hd_read_params(par_files)
 
    physics_type = "hd"
    phys_energy  = hd_energy
    phys_total_energy  = hd_energy
    phys_internal_e = .false.
    phys_gamma = hd_gamma
    phys_partial_ionization=hd_partial_ionization
 
    phys_trac=hd_trac
    if(phys_trac) then
      if(ndim .eq. 1) then
        if(hd_trac_type .gt. 2) then
          hd_trac_type=1
          if(mype==0) write(*,*) 'WARNING: set hd_trac_type=1'
        end if
        phys_trac_type=hd_trac_type
      else
        phys_trac=.false.
        if(mype==0) write(*,*) 'WARNING: set hd_trac=F when ndim>=2'
      end if
    end if
 
    ! set default gamma for polytropic/isothermal process
    if(.not.hd_energy) then
      if(hd_thermal_conduction) then
        hd_thermal_conduction=.false.
        if(mype==0) write(*,*) 'WARNING: set hd_thermal_conduction=F when hd_energy=F'
      end if
      if(hd_radiative_cooling) then
        hd_radiative_cooling=.false.
        if(mype==0) write(*,*) 'WARNING: set hd_radiative_cooling=F when hd_energy=F'
      end if
    end if
    if(.not.eq_state_units) then
      if(hd_partial_ionization) then
        hd_partial_ionization=.false.
        if(mype==0) write(*,*) 'WARNING: set hd_partial_ionization=F when eq_state_units=F'
      end if
    end if
    use_particles = hd_particles
 
    allocate(start_indices(number_species),stop_indices(number_species))
 
    ! set the index of the first flux variable for species 1
    start_indices(1)=1
    ! Determine flux variables
    rho_ = var_set_rho()
 
    allocate(mom(ndir))
    mom(:) = var_set_momentum(ndir)
    m^c_=mom(^c);
 
    ! Set index of energy variable
    if (hd_energy) then
       e_ = var_set_energy()
       p_ = e_
    else
       e_ = -1
       p_ = -1
    end if
 
    if(hd_radiation_fld)then
       if(hd_cak_force)then
          if(mype==0) then
            write(*,*)'Warning: CAK force addition together with FLD radiation'
          endif
       endif
       if(hd_radiative_cooling)then
          if(mype==0) then
            write(*,*)'Warning: Optically thin cooling together with FLD radiation'
          endif
       endif
       if(hd_dust.and.hd_dust_implicit)then
          call mpistop('implicit dust addition not compatible with FLD radiation')
       endif
       if(si_unit)then
          call mpistop('using FLD implies the use of cgs units')
       endif
       if(.not.hd_energy)then
          call mpistop('using FLD implies the use of an energy equation, set hd_energy=T')
       else
          !> set added variable and equation for radiation energy
          r_e = var_set_radiation_energy()
          phys_set_mg_bounds       => hd_set_mg_bounds
          phys_get_tgas            => hd_get_temperature_from_etot
          !> Initiate radiation-closure module
          select case (hd_radiation_fld_formalism)
          case('fld')
           call fld_init(he_abundance, hd_gamma)
          case('afld')
           {^ifoned
          call mpistop('using anisotropic FLD implies multidimensional setup')
           }
           call afld_init(he_abundance, hd_gamma)
          case default
           call mpistop('Radiation formalism unknown')
          end select
       endif
    else
      r_e=-1
    endif
 
    phys_get_dt              => hd_get_dt
    phys_get_cmax            => hd_get_cmax
    phys_get_tcutoff         => hd_get_tcutoff
    phys_get_cbounds         => hd_get_cbounds
    phys_get_flux            => hd_get_flux
    phys_add_source_geom     => hd_add_source_geom
    phys_add_source          => hd_add_source
    phys_to_conserved        => hd_to_conserved
    phys_to_primitive        => hd_to_primitive
    phys_check_params        => hd_check_params
    phys_check_w             => hd_check_w
    phys_get_pthermal        => hd_get_pthermal
    phys_get_v               => hd_get_v
    phys_get_rho             => hd_get_rho
    phys_write_info          => hd_write_info
    phys_handle_small_values => hd_handle_small_values
 
    ! derive units from basic units
    call hd_physical_units()
 
    if (hd_dust) then
        call dust_init(rho_, mom(:), e_)
    endif
 
    allocate(tracer(hd_n_tracer))
 
    ! Set starting index of tracers
    do itr = 1, hd_n_tracer
       tracer(itr) = var_set_fluxvar("trc", "trp", itr, need_bc=.false.)
    end do
 
    !  set temperature as an auxiliary variable to get ionization degree
    if(hd_partial_ionization) then
      te_ = var_set_auxvar('Te','Te')
    else
      te_ = -1
    end if
 
    ! set number of variables which need update ghostcells
    nwgc=nwflux+nwaux
 
    ! set the index of the last flux variable for species 1
    stop_indices(1)=nwflux
 
    if(hd_trac) then
      tcoff_ = var_set_wextra()
      iw_tcoff=tcoff_
    else
      tcoff_ = -1
    end if
 
    ! choose Rfactor in ideal gas law
    if(hd_partial_ionization) then
      hd_get_rfactor=>rfactor_from_temperature_ionization
      phys_update_temperature => hd_update_temperature
    else if(associated(usr_rfactor)) then
      hd_get_rfactor=>usr_rfactor
    else
      hd_get_rfactor=>rfactor_from_constant_ionization
    end if
 
    ! initialize thermal conduction module
    if (hd_thermal_conduction) then
      if (.not. hd_energy) &
           call mpistop("thermal conduction needs hd_energy=T")
 
      call sts_init()
      call tc_init_params(hd_gamma)
 
      allocate(tc_fl)
      call tc_get_hd_params(tc_fl,tc_params_read_hd)
      call add_sts_method(hd_get_tc_dt_hd,hd_sts_set_source_tc_hd,e_,1,e_,1,.false.)
      call set_conversion_methods_to_head(hd_e_to_ei, hd_ei_to_e)
      call set_error_handling_to_head(hd_tc_handle_small_e)
      tc_fl%get_temperature_from_conserved => hd_get_temperature_from_etot
      tc_fl%get_temperature_from_eint => hd_get_temperature_from_eint
      tc_fl%get_rho => hd_get_rho
      tc_fl%e_ = e_
    end if
 
    ! Initialize radiative cooling module
    if (hd_radiative_cooling) then
      if (.not. hd_energy) &
           call mpistop("radiative cooling needs hd_energy=T")
      call radiative_cooling_init_params(hd_gamma,he_abundance)
      allocate(rc_fl)
      call radiative_cooling_init(rc_fl,rc_params_read)
      rc_fl%get_rho => hd_get_rho
      rc_fl%get_pthermal => hd_get_pthermal
      rc_fl%get_var_Rfactor => hd_get_rfactor
      rc_fl%e_ = e_
      rc_fl%Tcoff_ = tcoff_
    end if
    allocate(te_fl_hd)
    te_fl_hd%get_rho=> hd_get_rho
    te_fl_hd%get_pthermal=> hd_get_pthermal
    te_fl_hd%get_var_Rfactor => hd_get_rfactor
{^ifthreed
    phys_te_images => hd_te_images
}
    ! Initialize viscosity module
    if (hd_viscosity) call viscosity_init(phys_wider_stencil)
 
    ! Initialize gravity module
    if (hd_gravity) call gravity_init()
 
    ! Initialize rotating_frame module
    if (hd_rotating_frame) call rotating_frame_init()
 
    ! Initialize CAK radiation force module
    if (hd_cak_force) call cak_init(hd_gamma)
 
    ! Initialize particles module
    if (hd_particles) then
       call particles_init()
    end if
 
    ! Check whether custom flux types have been defined
    if (.not. allocated(flux_type)) then
       allocate(flux_type(ndir, nw))
       flux_type = flux_default
    else if (any(shape(flux_type) /= [ndir, nw])) then
       call mpistop("phys_check error: flux_type has wrong shape")
    end if
 
    nvector      = 1 ! No. vector vars
    allocate(iw_vector(nvector))
    iw_vector(1) = mom(1) - 1
    ! initialize ionization degree table
    if(hd_partial_ionization) call ionization_degree_init()
 
  end subroutine hd_phys_init
 
{^ifthreed
  subroutine hd_te_images
    use mod_global_parameters
    use mod_thermal_emission
    select case(convert_type)
      case('EIvtiCCmpi','EIvtuCCmpi')
        call get_euv_image(unitconvert,te_fl_hd)
      case('ESvtiCCmpi','ESvtuCCmpi')
        call get_euv_spectrum(unitconvert,te_fl_hd)
      case('SIvtiCCmpi','SIvtuCCmpi')
        call get_sxr_image(unitconvert,te_fl_hd)
      case('WIvtiCCmpi','WIvtuCCmpi')
        call get_whitelight_image(unitconvert,te_fl_hd)
      case default
        call mpistop("Error in synthesize emission: Unknown convert_type")
      end select
  end subroutine hd_te_images
}
!!start th cond
  ! wrappers for STS functions in thermal_conductivity module
  ! which take as argument the tc_fluid (defined in the physics module)
  subroutine  hd_sts_set_source_tc_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
    use mod_global_parameters
    use mod_fix_conserve
    use mod_thermal_conduction, only: sts_set_source_tc_hd
    integer, intent(in) :: ixi^l, ixo^l, igrid, nflux
    double precision, intent(in) ::  x(ixi^s,1:ndim)
    double precision, intent(inout) ::  wres(ixi^s,1:nw), w(ixi^s,1:nw)
    double precision, intent(in) :: my_dt
    logical, intent(in) :: fix_conserve_at_step
    call sts_set_source_tc_hd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl)
  end subroutine hd_sts_set_source_tc_hd
 
  function hd_get_tc_dt_hd(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
    !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para/rho)
    !and   tc_k_para can depend on T=p/rho
    use mod_global_parameters
    use mod_thermal_conduction, only: get_tc_dt_hd
 
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
    double precision, intent(in) :: w(ixi^s,1:nw)
    double precision :: dtnew
 
    dtnew=get_tc_dt_hd(w,ixi^l,ixo^l,dx^d,x,tc_fl)
  end function hd_get_tc_dt_hd
 
  subroutine hd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
    ! move this in a different  routine as in mhd if needed in more places
    use mod_global_parameters
    use mod_small_values
 
    integer, intent(in)             :: ixi^l,ixo^l
    double precision, intent(inout) :: w(ixi^s,1:nw)
    double precision, intent(in)    :: x(ixi^s,1:ndim)
    integer, intent(in)    :: step
 
    integer :: idir
    logical :: flag(ixi^s,1:nw)
    character(len=140) :: error_msg
 
    flag=.false.
    where(w(ixo^s,e_)<small_e) flag(ixo^s,e_)=.true.
    if(any(flag(ixo^s,e_))) then
      select case (small_values_method)
      case ("replace")
        where(flag(ixo^s,e_)) w(ixo^s,e_)=small_e
      case ("average")
        call small_values_average(ixi^l, ixo^l, w, x, flag, e_)
      case default
        ! small values error shows primitive variables
        w(ixo^s,e_)=w(ixo^s,e_)*(hd_gamma - 1.0d0)
        do idir = 1, ndir
           w(ixo^s, iw_mom(idir)) = w(ixo^s, iw_mom(idir))/w(ixo^s,rho_)
        end do
        write(error_msg,"(a,i3)") "Thermal conduction step ", step
        call small_values_error(w, x, ixi^l, ixo^l, flag, error_msg)
      end select
    end if
  end subroutine hd_tc_handle_small_e
 
    ! fill in tc_fluid fields from namelist
    subroutine tc_params_read_hd(fl)
      use mod_global_parameters, only: unitpar,par_files
      type(tc_fluid), intent(inout) :: fl
      integer                      :: n
      logical :: tc_saturate=.false.
      double precision :: tc_k_para=0d0
 
      namelist /tc_list/ tc_saturate, tc_k_para
 
      do n = 1, size(par_files)
         open(unitpar, file=trim(par_files(n)), status="old")
         read(unitpar, tc_list, end=111)
111      close(unitpar)
      end do
      fl%tc_saturate = tc_saturate
      fl%tc_k_para = tc_k_para
 
    end subroutine tc_params_read_hd
 
  subroutine hd_get_rho(w,x,ixI^L,ixO^L,rho)
    use mod_global_parameters
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s,1:nw),x(ixi^s,1:ndim)
    double precision, intent(out) :: rho(ixi^s)
 
    rho(ixo^s) = w(ixo^s,rho_)
 
  end subroutine hd_get_rho
 
!!end th cond
!!rad cool
    subroutine rc_params_read(fl)
      use mod_global_parameters, only: unitpar,par_files
      use mod_constants, only: bigdouble
      use mod_basic_types, only: std_len
      type(rc_fluid), intent(inout) :: fl
      integer                      :: n
      ! list parameters
      integer :: ncool = 4000
    
      !> Name of cooling curve
      character(len=std_len)  :: coolcurve='JCcorona'
    
      !> Fixed temperature not lower than tlow
      logical    :: tfix=.false.
    
      !> Lower limit of temperature
      double precision   :: tlow=bigdouble
    
      !> Add cooling source in a split way (.true.) or un-split way (.false.)
      logical    :: rc_split=.false.
 
 
      namelist /rc_list/ coolcurve, ncool, tlow, tfix, rc_split
  
      do n = 1, size(par_files)
        open(unitpar, file=trim(par_files(n)), status="old")
        read(unitpar, rc_list, end=111)
111     close(unitpar)
      end do
 
      fl%ncool=ncool
      fl%coolcurve=coolcurve
      fl%tlow=tlow
      fl%Tfix=tfix
      fl%rc_split=rc_split
    end subroutine rc_params_read
!! end rad cool
 
  subroutine hd_check_params
    use mod_global_parameters
    use mod_dust, only: dust_check_params, dust_implicit_update, dust_evaluate_implicit
 
    if (.not. hd_energy) then
       if (hd_gamma <= 0.0d0) call mpistop ("Error: hd_gamma <= 0")
       if (hd_adiab < 0.0d0) call mpistop  ("Error: hd_adiab < 0")
       small_pressure= hd_adiab*small_density**hd_gamma
    else
       if (hd_gamma <= 0.0d0 .or. hd_gamma == 1.0d0) &
            call mpistop ("Error: hd_gamma <= 0 or hd_gamma == 1.0")
       small_e = small_pressure/(hd_gamma - 1.0d0)
       inv_gamma_1=1.d0/(hd_gamma-1.d0)
       small_r_e = small_pressure/(hd_gamma - 1.0d0)
    end if
 
    if (hd_dust) call dust_check_params()
 
    if(hd_dust_implicit) then
        if(.not.use_imex_scheme)then
           call mpistop('select IMEX scheme for implicit dust update')
        endif
        ! implicit dust update
        phys_implicit_update => dust_implicit_update
        phys_evaluate_implicit => dust_evaluate_implicit
    endif  
    if(hd_radiation_fld) then
        if(.not.use_imex_scheme)then
           call mpistop('select IMEX scheme for FLD radiation use')
        endif
        if(use_multigrid)then
           call hd_set_mg_bounds()
        else
           call mpistop('multigrid must have BCs for IMEX and FLD radiation use')
        endif
    endif
 
  end subroutine hd_check_params
 
  !> Set the boundaries for the diffusion of E
  subroutine hd_set_mg_bounds
    use mod_global_parameters
    use mod_multigrid_coupling
    use mod_usr_methods
 
    integer :: ib
 
    ! Set boundary conditions for the multigrid solver
    do ib = 1, 2*ndim
       select case (typeboundary(r_e, ib))
       case (bc_symm)
          ! d/dx u = 0
          mg%bc(ib, mg_iphi)%bc_type = mg_bc_neumann
          mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
       case (bc_asymm)
          ! u = 0
          mg%bc(ib, mg_iphi)%bc_type = mg_bc_dirichlet
          mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
       case (bc_cont)
          ! d/dx u = 0
          ! mg%bc(iB, mg_iphi)%bc_type = mg_bc_continuous
          mg%bc(ib, mg_iphi)%bc_type = mg_bc_neumann
          mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
       case (bc_periodic)
          ! Nothing to do here
       case (bc_noinflow)
          call usr_special_mg_bc(ib)
       case (bc_special)
          call usr_special_mg_bc(ib)
       case default
          call mpistop("divE_multigrid warning: unknown b.c. ")
       end select
    end do
  end subroutine hd_set_mg_bounds
 
  subroutine hd_physical_units
    use mod_global_parameters
    double precision :: mp,kb
    double precision :: a,b
    ! Derive scaling units
    if(si_unit) then
      mp=mp_si
      kb=kb_si
    else
      mp=mp_cgs
      kb=kb_cgs
    end if
    if(eq_state_units) then
      a=1d0+4d0*he_abundance
      if(hd_partial_ionization) then
        b=1d0+h_ion_fr+he_abundance*(he_ion_fr*(he_ion_fr2+1d0)+1d0)
      else
        b=2d0+3d0*he_abundance
      end if
      rr=1d0
    else
      a=1d0
      b=1d0
      rr=(1d0+h_ion_fr+he_abundance*(he_ion_fr*(he_ion_fr2+1d0)+1d0))/(1d0+4d0*he_abundance)
    end if
    if(unit_density/=1.d0 .or. unit_numberdensity/=1.d0) then
      if(unit_density/=1.d0) then
        unit_numberdensity=unit_density/(a*mp)
      else if(unit_numberdensity/=1.d0) then
        unit_density=a*mp*unit_numberdensity
      end if
      if(unit_temperature/=1.d0) then
        unit_pressure=b*unit_numberdensity*kb*unit_temperature
        unit_velocity=sqrt(unit_pressure/unit_density)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_pressure/=1.d0) then
        unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
        unit_velocity=sqrt(unit_pressure/unit_density)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_velocity/=1.d0) then
        unit_pressure=unit_density*unit_velocity**2
        unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_time/=1.d0) then
        unit_velocity=unit_length/unit_time
        unit_pressure=unit_density*unit_velocity**2
        unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
      end if
    else if(unit_temperature/=1.d0) then
      ! units of temperature and velocity are dependent
      if(unit_pressure/=1.d0) then
        unit_numberdensity=unit_pressure/(b*unit_temperature*kb)
        unit_density=a*mp*unit_numberdensity
        unit_velocity=sqrt(unit_pressure/unit_density)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      end if
    else if(unit_pressure/=1.d0) then
      if(unit_velocity/=1.d0) then
        unit_density=unit_pressure/unit_velocity**2
        unit_numberdensity=unit_density/(a*mp)
        unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_time/=0.d0) then
        unit_velocity=unit_length/unit_time
        unit_density=unit_pressure/unit_velocity**2
        unit_numberdensity=unit_density/(a*mp)
        unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
      end if
    end if
    unit_mass = unit_density * unit_length**3
 
    !> Units for radiative flux and opacity, latter is used in FLD
    unit_radflux = unit_velocity*unit_pressure
    unit_opacity = one/(unit_density*unit_length)
 
  end subroutine hd_physical_units
 
  !> Returns logical argument flag where values are ok
  subroutine hd_check_w(primitive, ixI^L, ixO^L, w, flag)
    use mod_global_parameters
    use mod_dust, only: dust_check_w
 
    logical, intent(in)          :: primitive
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, nw)
    logical, intent(inout)       :: flag(ixi^s,1:nw)
    double precision             :: tmp(ixi^s)
 
    flag=.false.
 
    if (hd_energy) then
       if (primitive) then
          where(w(ixo^s, e_) < small_pressure) flag(ixo^s,e_) = .true.
       else
          tmp(ixo^s)=(hd_gamma-1.0d0)*(w(ixo^s,e_)-&
           half*(^c&w(ixo^s,m^c_)**2+)/w(ixo^s,rho_))
          where(tmp(ixo^s) < small_pressure) flag(ixo^s,e_) = .true.
       endif
       if(hd_radiation_fld)then
          where(w(ixo^s, r_e) < small_r_e) flag(ixo^s,r_e) = .true.
       endif
    end if
 
    where(w(ixo^s, rho_) < small_density) flag(ixo^s,rho_) = .true.
 
    if(hd_dust) call dust_check_w(ixi^l,ixo^l,w,flag)
 
  end subroutine hd_check_w
 
  !> Transform primitive variables into conservative ones
  subroutine hd_to_conserved(ixI^L, ixO^L, w, x)
    use mod_global_parameters
    use mod_dust, only: dust_to_conserved
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
    integer :: ix^d
 
    {do ix^db=ixomin^db,ixomax^db\}
      if (hd_energy) then
         ! Calculate total energy from pressure and kinetic energy
         w(ix^d,e_)=w(ix^d, e_)*inv_gamma_1+&
          half*(^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)
      end if
      ! Convert velocity to momentum
      ^c&w(ix^d,m^c_)=w(ix^d,rho_)*w(ix^d,m^c_)\
    {end do\}
 
    if (hd_dust) then
      call dust_to_conserved(ixi^l, ixo^l, w, x)
    end if
 
  end subroutine hd_to_conserved
 
  !> Transform conservative variables into primitive ones
  subroutine hd_to_primitive(ixI^L, ixO^L, w, x)
    use mod_global_parameters
    use mod_dust, only: dust_to_primitive
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
    double precision                :: inv_rho
    integer :: ix^d
 
    if (fix_small_values) then
      call hd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'hd_to_primitive')
    end if
 
   {do ix^db=ixomin^db,ixomax^db\}
      inv_rho = 1.d0/w(ix^d,rho_)
      ! Convert momentum to velocity
      ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
      ! Calculate pressure = (gamma-1) * (e-ek)
      if(hd_energy) then
         ! Compute pressure
        w(ix^d,p_)=(hd_gamma-1.d0)*(w(ix^d,e_)&
                  -half*w(ix^d,rho_)*(^c&w(ix^d,m^c_)**2+))
      end if
   {end do\}
 
    ! Convert dust momentum to dust velocity
    if (hd_dust) then
      call dust_to_primitive(ixi^l, ixo^l, w, x)
    end if
 
  end subroutine hd_to_primitive
 
  !> Transform internal energy to total energy
  subroutine hd_ei_to_e(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
    ! Calculate total energy from internal and kinetic energy
    w(ixo^s,e_)=w(ixo^s,e_)+half*(^c&w(ixo^s,m^c_)**2+)/w(ixo^s,rho_)
 
  end subroutine hd_ei_to_e
 
  !> Transform total energy to internal energy
  subroutine hd_e_to_ei(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
    ! Calculate ei = e - ek
    w(ixo^s,e_)=w(ixo^s,e_)-half*(^c&w(ixo^s,m^c_)**2+)/w(ixo^s,rho_)
 
  end subroutine hd_e_to_ei
 
  !> Calculate v_i = m_i / rho within ixO^L
  subroutine hd_get_v_idim(w, x, ixI^L, ixO^L, idim, v)
    use mod_global_parameters
    integer, intent(in)           :: ixi^l, ixo^l, idim
    double precision, intent(in)  :: w(ixi^s, nw), x(ixi^s, 1:ndim)
    double precision, intent(out) :: v(ixi^s)
 
    v(ixo^s) = w(ixo^s, mom(idim)) / w(ixo^s, rho_)
  end subroutine hd_get_v_idim
 
  !> Calculate velocity vector v_i = m_i / rho within ixO^L
  subroutine hd_get_v(w,x,ixI^L,ixO^L,v)
    use mod_global_parameters
 
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s,nw), x(ixi^s,1:^nd)
    double precision, intent(out) :: v(ixi^s,1:ndir)
 
    integer :: idir
 
    do idir=1,ndir
      v(ixo^s,idir) = w(ixo^s, mom(idir)) / w(ixo^s, rho_)
    end do
 
  end subroutine hd_get_v
 
  !> Calculate cmax_idim = csound + abs(v_idim) within ixO^L
  subroutine hd_get_cmax(w, x, ixI^L, ixO^L, idim, cmax)
    use mod_global_parameters
    use mod_dust, only: dust_get_cmax_prim
    use mod_usr_methods, only: usr_set_pthermal
 
    integer, intent(in)                       :: ixi^l, ixo^l, idim
    ! w in primitive form
    double precision, intent(in)              :: w(ixi^s, nw), x(ixi^s, 1:ndim)
    double precision, intent(inout)           :: cmax(ixi^s)
 
    if(hd_energy) then
      cmax(ixo^s)=dabs(w(ixo^s,mom(idim)))+dsqrt(hd_gamma*w(ixo^s,p_)/w(ixo^s,rho_))
    else
      if (.not. associated(usr_set_pthermal)) then
        cmax(ixo^s) = hd_adiab * w(ixo^s, rho_)**hd_gamma
      else
        call usr_set_pthermal(w,x,ixi^l,ixo^l,cmax)
      end if
      cmax(ixo^s)=dabs(w(ixo^s,mom(idim)))+dsqrt(hd_gamma*cmax(ixo^s)/w(ixo^s,rho_))
    end if
 
    if (hd_dust) then
      call dust_get_cmax_prim(w, x, ixi^l, ixo^l, idim, cmax)
    end if
  end subroutine hd_get_cmax
 
  !> get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
  subroutine hd_get_tcutoff(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
    use mod_global_parameters
    integer, intent(in) :: ixi^l,ixo^l
    double precision, intent(in) :: x(ixi^s,1:ndim)
    ! in primitive form
    double precision, intent(inout) :: w(ixi^s,1:nw)
    double precision, intent(out) :: tco_local, tmax_local
 
    double precision, parameter :: trac_delta=0.25d0
    double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
    double precision :: ltr(ixi^s),ltrc,ltrp,tcoff(ixi^s)
    integer :: jxo^l,hxo^l
    integer :: jxp^l,hxp^l,ixp^l
    logical :: lrlt(ixi^s)
 
    {^ifoned
    call hd_get_rfactor(w,x,ixi^l,ixi^l,te)
    te(ixi^s)=w(ixi^s,p_)/(te(ixi^s)*w(ixi^s,rho_))
 
    tco_local=zero
    tmax_local=maxval(te(ixo^s))
    select case(hd_trac_type)
    case(0)
      w(ixi^s,tcoff_)=3.d5/unit_temperature
    case(1)
      hxo^l=ixo^l-1;
      jxo^l=ixo^l+1;
      lts(ixo^s)=0.5d0*dabs(te(jxo^s)-te(hxo^s))/te(ixo^s)
      lrlt=.false.
      where(lts(ixo^s) > trac_delta)
        lrlt(ixo^s)=.true.
      end where
      if(any(lrlt(ixo^s))) then
        tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
      end if
    case(2)
      !> iijima et al. 2021, LTRAC method
      ltrc=1.5d0
      ltrp=2.5d0
      ixp^l=ixo^l^ladd1;
      hxo^l=ixo^l-1;
      jxo^l=ixo^l+1;
      hxp^l=ixp^l-1;
      jxp^l=ixp^l+1;
      lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
      ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
      w(ixo^s,tcoff_)=te(ixo^s)*&
        (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
    case default
      call mpistop("hd_trac_type not allowed for 1D simulation")
    end select
    }
  end subroutine hd_get_tcutoff
 
  !> Calculate cmax_idim = csound + abs(v_idim) within ixO^L
  subroutine hd_get_cbounds(wLC, wRC, wLp, wRp, x, ixI^L, ixO^L, idim,Hspeed,cmax, cmin)
    use mod_global_parameters
    use mod_dust, only: dust_get_cmax
    use mod_variables
 
    integer, intent(in)             :: ixi^l, ixo^l, idim
    ! conservative left and right status
    double precision, intent(in)    :: wlc(ixi^s, nw), wrc(ixi^s, nw)
    ! primitive left and right status
    double precision, intent(in)    :: wlp(ixi^s, nw), wrp(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
    double precision, intent(inout) :: cmax(ixi^s,1:number_species)
    double precision, intent(inout), optional :: cmin(ixi^s,1:number_species)
    double precision, intent(in)    :: hspeed(ixi^s,1:number_species)
 
    double precision :: wmean(ixi^s,nw)
    double precision, dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
    integer :: ix^d
 
    select case(boundspeed)
    case (1)
      ! This implements formula (10.52) from "Riemann Solvers and Numerical
      ! Methods for Fluid Dynamics" by Toro.
 
      tmp1(ixo^s)=dsqrt(wlp(ixo^s,rho_))
      tmp2(ixo^s)=dsqrt(wrp(ixo^s,rho_))
      tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
      umean(ixo^s)=(wlp(ixo^s,mom(idim))*tmp1(ixo^s)+wrp(ixo^s,mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
 
      if(hd_energy) then
        csoundl(ixo^s)=hd_gamma*wlp(ixo^s,p_)/wlp(ixo^s,rho_)
        csoundr(ixo^s)=hd_gamma*wrp(ixo^s,p_)/wrp(ixo^s,rho_)
      else
        call hd_get_csound2(wlc,x,ixi^l,ixo^l,csoundl)
        call hd_get_csound2(wrc,x,ixi^l,ixo^l,csoundr)
      end if
 
      dmean(ixo^s) = (tmp1(ixo^s)*csoundl(ixo^s)+tmp2(ixo^s)*csoundr(ixo^s)) * &
           tmp3(ixo^s) + 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2 * &
           (wrp(ixo^s,mom(idim))-wlp(ixo^s,mom(idim)))**2
 
      dmean(ixo^s)=dsqrt(dmean(ixo^s))
      if(present(cmin)) then
        cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
        cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
        if(h_correction) then
          {do ix^db=ixomin^db,ixomax^db\}
            cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
            cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
          {end do\}
        end if
      else
        cmax(ixo^s,1)=dabs(umean(ixo^s))+dmean(ixo^s)
      end if
 
      if (hd_dust) then
        wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
        call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
      end if
 
    case (2)
      wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
      tmp1(ixo^s)=wmean(ixo^s,mom(idim))/wmean(ixo^s,rho_)
      call hd_get_csound2(wmean,x,ixi^l,ixo^l,csoundr)
      csoundr(ixo^s) = dsqrt(csoundr(ixo^s))
 
      if(present(cmin)) then
        cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
        cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
        if(h_correction) then
          {do ix^db=ixomin^db,ixomax^db\}
            cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
            cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
          {end do\}
        end if
      else
        cmax(ixo^s,1)=dabs(tmp1(ixo^s))+csoundr(ixo^s)
      end if
 
      if (hd_dust) then
        call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
      end if
    case (3)
      ! Miyoshi 2005 JCP 208, 315 equation (67)
      if(hd_energy) then
        csoundl(ixo^s)=hd_gamma*wlp(ixo^s,p_)/wlp(ixo^s,rho_)
        csoundr(ixo^s)=hd_gamma*wrp(ixo^s,p_)/wrp(ixo^s,rho_)
      else
        call hd_get_csound2(wlc,x,ixi^l,ixo^l,csoundl)
        call hd_get_csound2(wrc,x,ixi^l,ixo^l,csoundr)
      end if
      csoundl(ixo^s)=max(dsqrt(csoundl(ixo^s)),dsqrt(csoundr(ixo^s)))
      if(present(cmin)) then
        cmin(ixo^s,1)=min(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))-csoundl(ixo^s)
        cmax(ixo^s,1)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
        if(h_correction) then
          {do ix^db=ixomin^db,ixomax^db\}
            cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
            cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
          {end do\}
        end if
      else
        cmax(ixo^s,1)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
      end if
      if (hd_dust) then
        wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
        call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
      end if
    end select
 
  end subroutine hd_get_cbounds
 
  !> Calculate the square of the thermal sound speed csound2 within ixO^L.
  !> csound2=gamma*p/rho
  subroutine hd_get_csound2(w,x,ixI^L,ixO^L,csound2)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: w(ixi^s,nw)
    double precision, intent(in)    :: x(ixi^s,1:ndim)
    double precision, intent(out)   :: csound2(ixi^s)
 
    call hd_get_pthermal(w,x,ixi^l,ixo^l,csound2)
    csound2(ixo^s)=hd_gamma*csound2(ixo^s)/w(ixo^s,rho_)
 
  end subroutine hd_get_csound2
 
  !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho) within ixO^L
  subroutine hd_get_pthermal(w, x, ixI^L, ixO^L, pth)
    use mod_global_parameters
    use mod_usr_methods, only: usr_set_pthermal
    use mod_small_values, only: trace_small_values
 
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out):: pth(ixi^s)
    integer                      :: iw, ix^d
 
    if (hd_energy) then
       pth(ixo^s) = (hd_gamma - 1.0d0) * (w(ixo^s, e_) - &
            hd_kin_en(w, ixi^l, ixo^l))
    else
       if (.not. associated(usr_set_pthermal)) then
          pth(ixo^s) = hd_adiab * w(ixo^s, rho_)**hd_gamma
       else
          call usr_set_pthermal(w,x,ixi^l,ixo^l,pth)
       end if
    end if
 
    if (fix_small_values) then
      {do ix^db= ixo^lim^db\}
         if(pth(ix^d)<small_pressure) then
            pth(ix^d)=small_pressure
         endif
      {enddo^d&\}
    else if (check_small_values) then
      {do ix^db= ixo^lim^db\}
         if(pth(ix^d)<small_pressure) then
           write(*,*) "Error: small value of gas pressure",pth(ix^d),&
                " encountered when call hd_get_pthermal"
           write(*,*) "Iteration: ", it, " Time: ", global_time
           write(*,*) "Location: ", x(ix^d,:)
           write(*,*) "Cell number: ", ix^d
           do iw=1,nw
             write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
           end do
           ! use erroneous arithmetic operation to crash the run
           if(trace_small_values) write(*,*) dsqrt(pth(ix^d)-bigdouble)
           write(*,*) "Saving status at the previous time step"
           crash=.true.
         end if
      {enddo^d&\}
    end if
 
  end subroutine hd_get_pthermal
 
  !> Calculate radiation pressure within ixO^L
  subroutine hd_get_pradiation(w, x, ixI^L, ixO^L, prad)
    use mod_global_parameters
    use mod_fld
    use mod_afld
 
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out):: prad(ixo^s, 1:ndim, 1:ndim)
 
    select case (hd_radiation_fld_formalism)
    case('fld')
      call fld_get_radpress(w, x, ixi^l, ixo^l, prad, nghostcells)
    case('afld')
      call afld_get_radpress(w, x, ixi^l, ixo^l, prad, nghostcells)
    case default
      call mpistop('Radiation formalism unknown')
    end select
  end subroutine hd_get_pradiation
 
  !> calculates the sum of the gas pressure and max Prad tensor element
  subroutine hd_get_ptot(w, x, ixI^L, ixO^L, ptot)
    use mod_global_parameters
 
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision             :: pth(ixi^s)
    double precision             :: prad_tensor(ixo^s, 1:ndim, 1:ndim)
    double precision             :: prad_max(ixo^s)
    double precision, intent(out):: ptot(ixi^s)
    integer :: ix^d
 
    call hd_get_pthermal(w, x, ixi^l, ixo^l, pth)
    call hd_get_pradiation(w, x, ixi^l, ixo^l, prad_tensor)
 
    {do ix^d = ixomin^d,ixomax^d\}
      prad_max(ix^d) = maxval(prad_tensor(ix^d,:,:))
    {enddo\}
 
    ptot(ixo^s) = pth(ixo^s) + prad_max(ixo^s)
 
  end subroutine hd_get_ptot
 
  !> Calculates radiation temperature
  ! note: const_rad_a is assuming cgs units
  subroutine hd_get_trad(w, x, ixI^L, ixO^L, trad)
    use mod_global_parameters
    use mod_constants
 
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out):: trad(ixi^s)
 
    trad(ixi^s) = (w(ixi^s,r_e)*unit_pressure&
    /const_rad_a)**(1.d0/4.d0)/unit_temperature
 
  end subroutine hd_get_trad
 
  !> Calculate temperature=p/rho when in e_ the  total energy is stored
  subroutine hd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out):: res(ixi^s)
 
    double precision :: r(ixi^s)
 
    call hd_get_rfactor(w,x,ixi^l,ixo^l,r)
    call hd_get_pthermal(w, x, ixi^l, ixo^l, res)
    res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,rho_))
  end subroutine hd_get_temperature_from_etot
 
  !> Calculate temperature=p/rho when in e_ the  internal energy is stored
  subroutine hd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out):: res(ixi^s)
 
    double precision :: r(ixi^s)
 
    call hd_get_rfactor(w,x,ixi^l,ixo^l,r)
    res(ixo^s) = (hd_gamma - 1.0d0) * w(ixo^s, e_)/(w(ixo^s,rho_)*r(ixo^s))
  end subroutine hd_get_temperature_from_eint
 
  ! Calculate flux f_idim[iw]
  subroutine hd_get_flux(wC, w, x, ixI^L, ixO^L, idim, f)
    use mod_global_parameters
    use mod_dust, only: dust_get_flux_prim
 
    integer, intent(in)             :: ixi^l, ixo^l, idim
    ! conservative w
    double precision, intent(in)    :: wc(ixi^s, 1:nw)
    ! primitive w
    double precision, intent(in)    :: w(ixi^s, 1:nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
    double precision, intent(out)   :: f(ixi^s, nwflux)
 
    double precision                :: pth(ixi^s)
    integer                         :: ix^db
 
    if (hd_energy) then
     {do ix^db=ixomin^db,ixomax^db\}
        f(ix^d,rho_)=w(ix^d,mom(idim))*w(ix^d,rho_)
        ! Momentum flux is v_i*m_i, +p in direction idim
        ^c&f(ix^d,m^c_)=w(ix^d,mom(idim))*wc(ix^d,m^c_)\
        f(ix^d,mom(idim))=f(ix^d,mom(idim))+w(ix^d,p_)
        ! Energy flux is v_i*(e + p)
        f(ix^d,e_)=w(ix^d,mom(idim))*(wc(ix^d,e_)+w(ix^d,p_))
     {end do\}
    else
      call hd_get_pthermal(wc, x, ixi^l, ixo^l, pth)
     {do ix^db=ixomin^db,ixomax^db\}
        f(ix^d,rho_)=w(ix^d,mom(idim))*w(ix^d,rho_)
        ! Momentum flux is v_i*m_i, +p in direction idim
        ^c&f(ix^d,m^c_)=w(ix^d,mom(idim))*wc(ix^d,m^c_)\
        f(ix^d,mom(idim))=f(ix^d,mom(idim))+pth(ix^d)
     {end do\}
    end if
 
    if(hd_radiation_fld)then
     {do ix^db=ixomin^db,ixomax^db\}
        ! advection of radiation enery v_i*r_e
        f(ix^d,r_e)=w(ix^d,mom(idim))*wc(ix^d,r_e)
     {end do\}
    endif
 
    do ix1 = 1, hd_n_tracer
       f(ixo^s, tracer(ix1)) = w(ixo^s,mom(idim)) * w(ixo^s, tracer(ix1))
    end do
 
    ! Dust fluxes
    if (hd_dust) then
      call dust_get_flux_prim(w, x, ixi^l, ixo^l, idim, f)
    end if
 
  end subroutine hd_get_flux
 
  !> Add geometrical source terms to w
  !>
  !> Notice that the expressions of the geometrical terms depend only on ndir,
  !> not ndim. Eg, they are the same in 2.5D and in 3D, for any geometry.
  !>
  subroutine hd_add_source_geom(qdt, dtfactor, ixI^L, ixO^L, wCT, wprim, w, x)
    use mod_global_parameters
    use mod_usr_methods, only: usr_set_surface, usr_set_pthermal
    use mod_rotating_frame, only: rotating_frame_add_source
    use mod_dust, only: dust_n_species, dust_mom, dust_rho
    use mod_geometry
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: qdt, dtfactor, x(ixi^s, 1:ndim)
    double precision, intent(inout) :: wct(ixi^s, 1:nw), wprim(ixi^s,1:nw),w(ixi^s, 1:nw)
    double precision :: pth(ixi^s), source(ixi^s), minrho
    integer                         :: iw,idir, h1x^l{^nooned, h2x^l}
    integer :: mr_,mphi_ ! Polar var. names
    integer :: irho, ifluid, n_fluids
    double precision :: exp_factor(ixi^s), del_exp_factor(ixi^s), exp_factor_primitive(ixi^s)
 
    if (hd_dust) then
       n_fluids = 1 + dust_n_species
    else
       n_fluids = 1
    end if
 
    select case (coordinate)
 
    case(cartesian_expansion)
      !the user provides the functions of exp_factor and del_exp_factor
      if(associated(usr_set_surface)) call usr_set_surface(ixi^l,x,block%dx,exp_factor,del_exp_factor,exp_factor_primitive)
      if(hd_energy) then
        source(ixo^s)=wprim(ixo^s, p_)
      else
        if(.not. associated(usr_set_pthermal)) then
          source(ixo^s)=hd_adiab * wprim(ixo^s, rho_)**hd_gamma
        else
          call usr_set_pthermal(wct,x,ixi^l,ixo^l,source)
        end if
      end if
      source(ixo^s) = source(ixo^s)*del_exp_factor(ixo^s)/exp_factor(ixo^s)
      w(ixo^s,mom(1)) = w(ixo^s,mom(1)) + qdt*source(ixo^s)
 
    case (cylindrical)
       do ifluid = 0, n_fluids-1
          ! s[mr]=(pthermal+mphi**2/rho)/radius
          if (ifluid == 0) then
            ! gas
            irho  = rho_
            mr_   = mom(r_)
            if(phi_>0) mphi_ = mom(phi_)
            if(hd_energy) then
              source(ixo^s)=wprim(ixo^s, p_)
            else
              if(.not. associated(usr_set_pthermal)) then
                source(ixo^s)=hd_adiab * wprim(ixo^s, rho_)**hd_gamma
              else
                call usr_set_pthermal(wct,x,ixi^l,ixo^l,source)
              end if
            end if
            minrho = 0.0d0
          else
            ! dust : no pressure
            irho  = dust_rho(ifluid)
            mr_   = dust_mom(r_, ifluid)
            if(phi_>0) mphi_ = dust_mom(phi_, ifluid)
            source(ixi^s) = zero
            minrho = 0.0d0
          end if
          if(phi_ > 0) then
            where (wct(ixo^s, irho) > minrho)
              source(ixo^s) = source(ixo^s) + wct(ixo^s,mphi_)*wprim(ixo^s,mphi_)
              w(ixo^s, mr_) = w(ixo^s, mr_) + qdt*source(ixo^s)/x(ixo^s,r_)
            end where
            ! s[mphi]=(-mphi*vr)/radius
            where (wct(ixo^s, irho) > minrho)
              source(ixo^s) = -wct(ixo^s, mphi_) * wprim(ixo^s, mr_)
              w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * source(ixo^s) / x(ixo^s, r_)
            end where
          else
            ! s[mr]=2pthermal/radius
            w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * source(ixo^s) / x(ixo^s, r_)
          end if
       end do
    case (spherical)
       if (hd_dust) then
          call mpistop("Dust geom source terms not implemented yet with spherical geometries")
       end if
       mr_   = mom(r_)
       if(phi_>0) mphi_ = mom(phi_)
       h1x^l=ixo^l-kr(1,^d); {^nooned h2x^l=ixo^l-kr(2,^d);}
       if(hd_energy) then
         pth(ixo^s)=wprim(ixo^s, p_)
       else
         if(.not. associated(usr_set_pthermal)) then
           pth(ixo^s)=hd_adiab * wprim(ixo^s, rho_)**hd_gamma
         else
           call usr_set_pthermal(wct,x,ixi^l,ixo^l,pth)
         end if
       end if
       ! s[mr]=((vtheta**2+vphi**2)*rho+2*p)/r
       source(ixo^s) = pth(ixo^s) * x(ixo^s, 1) &
            *(block%surfaceC(ixo^s, 1) - block%surfaceC(h1x^s, 1)) &
            /block%dvolume(ixo^s)
       do idir = 2, ndir
         source(ixo^s) = source(ixo^s) + wprim(ixo^s, mom(idir))**2 * wprim(ixo^s, rho_)
       end do
       w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * source(ixo^s) / x(ixo^s, 1)
 
       {^nooned
       ! s[mtheta]=-(vr*vtheta*rho)/r+cot(theta)*(vphi**2*rho+p)/r
       source(ixo^s) = pth(ixo^s) * x(ixo^s, 1) &
            * (block%surfaceC(ixo^s, 2) - block%surfaceC(h2x^s, 2)) &
            / block%dvolume(ixo^s)
       if (ndir == 3) then
         source(ixo^s) = source(ixo^s) + (wprim(ixo^s, mom(3))**2 * wprim(ixo^s, rho_)) / tan(x(ixo^s, 2))
       end if
       source(ixo^s) = source(ixo^s) - (wprim(ixo^s, mom(2)) * wprim(ixo^s, mr_)) * wprim(ixo^s, rho_)
       w(ixo^s, mom(2)) = w(ixo^s, mom(2)) + qdt * source(ixo^s) / x(ixo^s, 1)
 
       if (ndir == 3) then
         ! s[mphi]=-(vphi*vr/rho)/r-cot(theta)*(vtheta*vphi/rho)/r
         source(ixo^s) = -(wprim(ixo^s, mom(3)) * wprim(ixo^s, mr_)) * wprim(ixo^s, rho_)&
                        - (wprim(ixo^s, mom(2)) * wprim(ixo^s, mom(3))) * wprim(ixo^s, rho_) / tan(x(ixo^s, 2))
         w(ixo^s, mom(3)) = w(ixo^s, mom(3)) + qdt * source(ixo^s) / x(ixo^s, 1)
       end if
       }
    end select
 
    if (hd_rotating_frame) then
       if (hd_dust) then
          call mpistop("Rotating frame not implemented yet with dust")
       else
          call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
       end if
    end if
 
  end subroutine hd_add_source_geom
 
  ! w[iw]= w[iw]+qdt*S[wCT, qtC, x] where S is the source based on wCT within ixO
  subroutine hd_add_source(qdt,dtfactor, ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
    use mod_global_parameters
    use mod_radiative_cooling, only: radiative_cooling_add_source
    use mod_dust, only: dust_add_source, dust_mom, dust_rho, dust_n_species
    use mod_viscosity, only: viscosity_add_source
    use mod_usr_methods, only: usr_gravity
    use mod_gravity, only: gravity_add_source, grav_split
    use mod_cak_force, only: cak_add_source
 
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: qdt, dtfactor
    double precision, intent(in)    :: wct(ixi^s, 1:nw),wctprim(ixi^s,1:nw), x(ixi^s, 1:ndim)
    double precision, intent(inout) :: w(ixi^s, 1:nw)
    logical, intent(in)             :: qsourcesplit
    logical, intent(inout)          :: active
 
    double precision :: gravity_field(ixi^s, 1:ndim)
    integer :: idust, idim
 
    if(hd_dust .and. .not. hd_dust_implicit) then
      call dust_add_source(qdt,ixi^l,ixo^l,wct,w,x,qsourcesplit,active)
    end if
 
    if(hd_radiative_cooling) then
      call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,wctprim,w,x,&
           qsourcesplit,active, rc_fl)
    end if
 
    if(hd_viscosity) then
      call viscosity_add_source(qdt,ixi^l,ixo^l,wct,wctprim,w,x,&
           hd_energy,qsourcesplit,active)
    end if
 
    if (hd_gravity) then
      call gravity_add_source(qdt,ixi^l,ixo^l,wct,wctprim,w,x,&
           hd_energy,qsourcesplit,active)
 
      if (hd_dust .and. qsourcesplit .eqv. grav_split) then
         active = .true.
 
         call usr_gravity(ixi^l, ixo^l, wct, x, gravity_field)
         do idust = 1, dust_n_species
            do idim = 1, ndim
               w(ixo^s, dust_mom(idim, idust)) = w(ixo^s, dust_mom(idim, idust)) &
                    + qdt * gravity_field(ixo^s, idim) * wct(ixo^s, dust_rho(idust))
            end do
         end do
      end if
    end if
 
    if (hd_cak_force) then
      call cak_add_source(qdt,ixi^l,ixo^l,wct,w,x,hd_energy,qsourcesplit,active)
    end if
 
    ! This is where the radiation force and heating/cooling are added
    if (hd_radiation_fld) then
       call hd_add_radiation_source(qdt,ixi^l,ixo^l,wct,w,x,qsourcesplit,active)
    endif
 
    if(hd_partial_ionization) then
      if(.not.qsourcesplit) then
        active = .true.
        call hd_update_temperature(ixi^l,ixo^l,wct,w,x)
      end if
    end if
 
  end subroutine hd_add_source
 
  subroutine hd_add_radiation_source(qdt,ixI^L,ixO^L,wCT,w,x,qsourcesplit,active)
    use mod_constants
    use mod_global_parameters
    use mod_usr_methods
    use mod_fld
    use mod_afld
 
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: qdt, x(ixi^s,1:ndim)
    double precision, intent(in)    :: wct(ixi^s,1:nw)
    double precision, intent(inout) :: w(ixi^s,1:nw)
    logical, intent(in) :: qsourcesplit
    logical, intent(inout) :: active
    double precision :: cmax(ixi^s)
 
    select case(hd_radiation_fld_formalism)
    case('fld')
      call fld_get_diffcoef_central(w, wct, x, ixi^l, ixo^l)
      ! radiation force
      call get_fld_rad_force(qdt,ixi^l,ixo^l,wct,w,x,&
        hd_energy,qsourcesplit,active)
      call hd_handle_small_values(.true., w, x, ixi^l, ixo^l, 'fld_add_radiation')
    case('afld')
      call afld_get_diffcoef_central(w, wct, x, ixi^l, ixo^l)
      ! radiation force
      call get_afld_rad_force(qdt,ixi^l,ixo^l,wct,w,x,&
        hd_energy,qsourcesplit,active)
      call hd_handle_small_values(.true., w, x, ixi^l, ixo^l, 'afld_add_radiation')
      ! photon tiring, heating and cooling
      call get_afld_energy_interact(qdt,ixi^l,ixo^l,wct,w,x,&
        hd_energy,qsourcesplit,active)
    case default
      call mpistop('Radiation formalism unknown')
    end select
 
  end subroutine hd_add_radiation_source
 
  subroutine hd_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
    use mod_global_parameters
    use mod_dust, only: dust_get_dt
    use mod_viscosity, only: viscosity_get_dt
    use mod_gravity, only: gravity_get_dt
    use mod_cak_force, only: cak_get_dt
    use mod_fld, only: fld_radforce_get_dt
    use mod_afld, only: afld_radforce_get_dt
 
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: dx^d, x(ixi^s, 1:^nd)
    double precision, intent(in)    :: w(ixi^s, 1:nw)
    double precision, intent(inout) :: dtnew
 
    dtnew = bigdouble
 
    if(hd_dust) then
      call dust_get_dt(w, ixi^l, ixo^l, dtnew, dx^d, x)
    end if
 
    if(hd_viscosity) then
      call viscosity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
    end if
 
    if(hd_gravity) then
      call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
   end if
 
   if (hd_cak_force) then
     call cak_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
   end if
 
   if(hd_radiation_fld) then
      select case(hd_radiation_fld_formalism)
        case('fld')
          call fld_radforce_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
        case('afld')
          call afld_radforce_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
        case default
          call mpistop('Radiation formalism unknown')
      end select
   endif
 
  end subroutine hd_get_dt
 
  function hd_kin_en(w, ixI^L, ixO^L, inv_rho) result(ke)
    use mod_global_parameters, only: nw, ndim
    integer, intent(in)                    :: ixi^l, ixo^l
    double precision, intent(in)           :: w(ixi^s, nw)
    double precision                       :: ke(ixo^s)
    double precision, intent(in), optional :: inv_rho(ixo^s)
 
    if (present(inv_rho)) then
       ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) * inv_rho
    else
       ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) / w(ixo^s, rho_)
    end if
  end function hd_kin_en
 
  function hd_inv_rho(w, ixI^L, ixO^L) result(inv_rho)
    use mod_global_parameters, only: nw, ndim
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s, nw)
    double precision              :: inv_rho(ixo^s)
 
    ! Can make this more robust
    inv_rho = 1.0d0 / w(ixo^s, rho_)
  end function hd_inv_rho
 
  subroutine hd_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
    ! handles hydro (density,pressure,velocity) bootstrapping
    ! any negative dust density is flagged as well (and throws an error)
    ! small_values_method=replace also for dust
    use mod_global_parameters
    use mod_small_values
    use mod_dust, only: dust_n_species, dust_mom, dust_rho
    logical, intent(in)             :: primitive
    integer, intent(in)             :: ixi^l,ixo^l
    double precision, intent(inout) :: w(ixi^s,1:nw)
    double precision, intent(in)    :: x(ixi^s,1:ndim)
    character(len=*), intent(in)    :: subname
 
    integer :: n,idir
    logical :: flag(ixi^s,1:nw)
 
    call hd_check_w(primitive, ixi^l, ixo^l, w, flag)
 
    if (any(flag)) then
      select case (small_values_method)
      case ("replace")
        where(flag(ixo^s,rho_)) w(ixo^s,rho_) = small_density
        do idir = 1, ndir
          if(small_values_fix_iw(mom(idir))) then
            where(flag(ixo^s,rho_)) w(ixo^s, mom(idir)) = 0.0d0
          end if
        end do
        if(hd_radiation_fld)then
          if (small_values_fix_iw(r_e)) then
            where(flag(ixo^s,r_e)) w(ixo^s,r_e) = small_r_e
          end if
        end if
        if(hd_energy)then
          if(small_values_fix_iw(e_)) then
            if(primitive) then
              where(flag(ixo^s,rho_)) w(ixo^s, p_) = small_pressure
            else
              where(flag(ixo^s,rho_)) w(ixo^s, e_) = small_e + hd_kin_en(w,ixi^l,ixo^l)
            endif
          end if
        endif
 
        if(hd_energy) then
          if(primitive) then
            where(flag(ixo^s,e_)) w(ixo^s,p_) = small_pressure
          else
            where(flag(ixo^s,e_))
              ! Add kinetic energy
              w(ixo^s,e_) = small_e + hd_kin_en(w,ixi^l,ixo^l)
            end where
          end if
        end if
 
        if(hd_dust)then
           do n=1,dust_n_species
              where(flag(ixo^s,dust_rho(n))) w(ixo^s,dust_rho(n)) = 0.0d0
              do idir = 1, ndir
                  where(flag(ixo^s,dust_rho(n))) w(ixo^s,dust_mom(idir,n)) = 0.0d0
              enddo
           enddo
        endif
      case ("average")
        if(primitive)then
           ! averaging for all primitive fields, including dust
           call small_values_average(ixi^l, ixo^l, w, x, flag)
        else
           ! do averaging of density
           call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
           if(hd_energy) then
             ! do averaging of pressure
             w(ixi^s,p_)=(hd_gamma-1.d0)*(w(ixi^s,e_) &
              -0.5d0*sum(w(ixi^s, mom(:))**2, dim=ndim+1)/w(ixi^s,rho_))
             call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
             w(ixi^s,e_)=w(ixi^s,p_)/(hd_gamma-1.d0) &
               +0.5d0*sum(w(ixi^s, mom(:))**2, dim=ndim+1)/w(ixi^s,rho_)
           end if
           if(hd_radiation_fld) then
              ! do averaging of radiative energy density
              call small_values_average(ixi^l, ixo^l, w, x, flag, r_e)
           endif
           if(hd_dust)then
              do n=1,dust_n_species
                 where(flag(ixo^s,dust_rho(n))) w(ixo^s,dust_rho(n)) = 0.0d0
                 do idir = 1, ndir
                    where(flag(ixo^s,dust_rho(n))) w(ixo^s,dust_mom(idir,n)) = 0.0d0
                 enddo
              enddo
          endif
        endif
      case default
        if(.not.primitive) then
          !convert w to primitive
          ! Calculate pressure = (gamma-1) * (e-ek)
          if(hd_energy) then
            w(ixo^s,p_)=(hd_gamma-1.d0)*(w(ixo^s,e_)-hd_kin_en(w,ixi^l,ixo^l))
          end if
          ! Convert gas momentum to velocity
          do idir = 1, ndir
            w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))/w(ixo^s,rho_)
          end do
        end if
        ! NOTE: dust entries may still have conserved values here
        call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
      end select
    end if
  end subroutine hd_handle_small_values
 
  subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
    use mod_global_parameters
    use mod_ionization_degree
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s,1:nw)
    double precision, intent(in) :: x(ixi^s,1:ndim)
    double precision, intent(out):: rfactor(ixi^s)
 
    double precision :: iz_h(ixo^s),iz_he(ixo^s)
 
    call ionization_degree_from_temperature(ixi^l,ixo^l,w(ixi^s,te_),iz_h,iz_he)
    ! assume the first and second ionization of Helium have the same degree
    rfactor(ixo^s)=(1.d0+iz_h(ixo^s)+0.1d0*(1.d0+iz_he(ixo^s)*(1.d0+iz_he(ixo^s))))/2.3d0
 
  end subroutine rfactor_from_temperature_ionization
 
  subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
    use mod_global_parameters
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s,1:nw)
    double precision, intent(in) :: x(ixi^s,1:ndim)
    double precision, intent(out):: rfactor(ixi^s)
 
    rfactor(ixo^s)=rr
 
  end subroutine rfactor_from_constant_ionization
 
  subroutine hd_update_temperature(ixI^L,ixO^L,wCT,w,x)
    use mod_global_parameters
    use mod_ionization_degree
 
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: wct(ixi^s,1:nw), x(ixi^s,1:ndim)
    double precision, intent(inout) :: w(ixi^s,1:nw)
 
    double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
 
    call ionization_degree_from_temperature(ixi^l,ixo^l,wct(ixi^s,te_),iz_h,iz_he)
 
    call hd_get_pthermal(w,x,ixi^l,ixo^l,pth)
 
    w(ixo^s,te_)=(2.d0+3.d0*he_abundance)*pth(ixo^s)/(w(ixo^s,rho_)*(1.d0+iz_h(ixo^s)+&
     he_abundance*(iz_he(ixo^s)*(iz_he(ixo^s)+1.d0)+1.d0)))
 
  end subroutine hd_update_temperature
 
end module mod_hd_phys