mod_ffhd_phys.t

Go to the documentation of this file.

!> Frozen-field hydrodynamics module
module mod_ffhd_phys
 
#include "amrvac.h"
 
  use mod_global_parameters, only: std_len, const_c
  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              :: ffhd_energy = .true.
 
  !> Whether thermal conduction is used
  logical, public, protected              :: ffhd_thermal_conduction = .false.
  !> Whether hyperbolic type thermal conduction is used
  logical, public, protected              :: ffhd_hyperbolic_thermal_conduction = .false.
  !> type of fluid for thermal conduction
  type(tc_fluid), public, allocatable     :: tc_fl
  !> type of fluid for thermal emission synthesis
  type(te_fluid), public, allocatable     :: te_fl_ffhd
 
  !> Whether radiative cooling is added
  logical, public, protected              :: ffhd_radiative_cooling = .false.
  !> type of fluid for radiative cooling
  type(rc_fluid), public, allocatable     :: rc_fl
 
  !> Whether viscosity is added
  logical, public, protected              :: ffhd_viscosity = .false.
 
  !> Whether gravity is added
  logical, public, protected              :: ffhd_gravity = .false.
 
  !> Whether TRAC method is used
  logical, public, protected              :: ffhd_trac = .false.
 
  !> Which TRAC method is used
  integer, public, protected              :: ffhd_trac_type=1
 
  !> Height of the mask used in the TRAC method
  double precision, public, protected     :: ffhd_trac_mask = 0.d0
 
  !> Distance between two adjacent traced magnetic field lines (in finest cell size)
  integer, public, protected              :: ffhd_trac_finegrid=4
 
  !> Whether plasma is partially ionized
  logical, public, protected              :: ffhd_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(:)
 
  !> 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_
 
  !> Indices of temperature
  integer, public, protected :: te_
 
  !> Index of the cutoff temperature for the TRAC method
  integer, public, protected              :: tcoff_
  integer, public, protected              :: tweight_
  integer, public, protected              :: q_
 
  !> The adiabatic index
  double precision, public                :: ffhd_gamma = 5.d0/3.0d0
 
  !> The adiabatic constant
  double precision, public                :: ffhd_adiab = 1.0d0
 
  !> The small_est allowed energy
  double precision, protected             :: small_e
 
  !> The thermal conductivity kappa in hyperbolic thermal conduction
  double precision, public                :: hypertc_kappa
 
  !> 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.
 
  !> gamma minus one and its inverse
  double precision :: gamma_1, inv_gamma_1
 
  !define the subroutine interface for the ambipolar mask
  abstract interface
 
    function fun_kin_en(w, ixI^L, ixO^L, inv_rho) result(ke)
      use mod_global_parameters, only: nw, ndim,block
      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)
    end function fun_kin_en
 
  end interface
 
  procedure(sub_convert), pointer      :: ffhd_to_primitive        => null()
  procedure(sub_convert), pointer      :: ffhd_to_conserved        => null()
  procedure(sub_small_values), pointer :: ffhd_handle_small_values => null()
  procedure(sub_get_pthermal), pointer :: ffhd_get_pthermal        => null()
  procedure(sub_get_pthermal), pointer :: ffhd_get_rfactor         => null()
  procedure(sub_get_pthermal), pointer :: ffhd_get_temperature     => null()
  procedure(sub_get_v), pointer        :: ffhd_get_v               => null()
  procedure(fun_kin_en), pointer       :: ffhd_kin_en              => null()
  ! Public methods
  public :: ffhd_phys_init
  public :: ffhd_kin_en
  public :: ffhd_get_pthermal
  public :: ffhd_get_temperature
  public :: ffhd_get_v
  public :: ffhd_get_rho
  public :: ffhd_get_v_idim
  public :: ffhd_to_conserved
  public :: ffhd_to_primitive
  public :: ffhd_get_csound2
  public :: ffhd_e_to_ei
  public :: ffhd_ei_to_e
 
contains
 
  subroutine ffhd_read_params(files)
    use mod_global_parameters
    use mod_particles, only: particles_eta, particles_etah
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /ffhd_list/ ffhd_energy, ffhd_gamma, ffhd_adiab,&
      ffhd_thermal_conduction, ffhd_hyperbolic_thermal_conduction, ffhd_radiative_cooling, ffhd_gravity,&
      ffhd_viscosity, he_abundance, h_ion_fr, he_ion_fr, he_ion_fr2, eq_state_units, si_unit,&
      b0field, busr, ffhd_partial_ionization, ffhd_trac, ffhd_trac_type, ffhd_trac_mask, ffhd_trac_finegrid
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status="old")
       read(unitpar, ffhd_list, end=111)
111    close(unitpar)
    end do
  end subroutine ffhd_read_params
 
  !> Write this module's parameters to a snapsoht
  subroutine ffhd_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) = ffhd_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 ffhd_write_info
 
  subroutine ffhd_phys_init()
    use mod_global_parameters
    use mod_thermal_conduction
    use mod_radiative_cooling
    use mod_viscosity, only: viscosity_init
    use mod_gravity, only: gravity_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
    integer :: itr, idir
 
    call ffhd_read_params(par_files)
 
    if(.not. ffhd_energy) then
      if(ffhd_thermal_conduction) then
        ffhd_thermal_conduction=.false.
        if(mype==0) write(*,*) 'WARNING: set ffhd_thermal_conduction=F when ffhd_energy=F'
      end if
      if(ffhd_thermal_conduction) then
        ffhd_hyperbolic_thermal_conduction=.false.
        if(mype==0) write(*,*) 'WARNING: set ffhd_hyperbolic_thermal_conduction=F when ffhd_energy=F'
      end if
      if(ffhd_radiative_cooling) then
        ffhd_radiative_cooling=.false.
        if(mype==0) write(*,*) 'WARNING: set ffhd_radiative_cooling=F when ffhd_energy=F'
      end if
      if(ffhd_trac) then
        ffhd_trac=.false.
        if(mype==0) write(*,*) 'WARNING: set ffhd_trac=F when ffhd_energy=F'
      end if
      if(ffhd_partial_ionization) then
        ffhd_partial_ionization=.false.
        if(mype==0) write(*,*) 'WARNING: set ffhd_partial_ionization=F when ffhd_energy=F'
      end if
    end if
    if(.not.eq_state_units) then
      if(ffhd_partial_ionization) then
        ffhd_partial_ionization=.false.
        if(mype==0) write(*,*) 'WARNING: set ffhd_partial_ionization=F when eq_state_units=F'
      end if
    end if
 
    if(ffhd_hyperbolic_thermal_conduction) then
      ffhd_thermal_conduction=.false.
      if(mype==0) write(*,*) 'WARNING: turn off parabolic TC when using hyperbolic TC'
    end if
 
    physics_type = "ffhd"
    phys_energy=ffhd_energy
    phys_internal_e=.false.
    phys_trac=ffhd_trac
    phys_trac_type=ffhd_trac_type
    phys_partial_ionization=ffhd_partial_ionization
 
    phys_gamma = ffhd_gamma
    phys_total_energy=ffhd_energy
    phys_trac_finegrid=ffhd_trac_finegrid
 
    {^ifoned
    if(ffhd_trac .and. ffhd_trac_type .gt. 2) then
      ffhd_trac_type=1
      if(mype==0) write(*,*) 'WARNING: reset ffhd_trac_type=1 for 1D simulation'
    end if
    }
    if(ffhd_trac .and. ffhd_trac_type .le. 4) then
      ffhd_trac_mask=bigdouble
      if(mype==0) write(*,*) 'WARNING: set ffhd_trac_mask==bigdouble for global TRAC method'
    end if
    phys_trac_mask=ffhd_trac_mask
 
    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(1))
    mom(:) = var_set_momentum(1)
 
    ! Set index of energy variable
    if(ffhd_energy) then
      e_     = var_set_energy() ! energy density
      p_     = e_               ! gas pressure
    else
      e_     = -1
      p_     = -1
    end if
 
    if(ffhd_hyperbolic_thermal_conduction) then
      q_ = var_set_q()
      need_global_cs2max=.true.
    else
      q_=-1
    end if
 
    ! set number of variables which need update ghostcells
    nwgc=nwflux
 
    ! set the index of the last flux variable for species 1
    stop_indices(1)=nwflux
 
    !  set temperature as an auxiliary variable to get ionization degree
    if(ffhd_partial_ionization) then
      te_ = var_set_auxvar('Te','Te')
    else
      te_ = -1
    end if
 
    ! set cutoff temperature when using the TRAC method, as well as an auxiliary weight
    tweight_ = -1
    if(ffhd_trac) then
      tcoff_ = var_set_wextra()
      iw_tcoff=tcoff_
      if(ffhd_trac_type .ge. 3) then
        tweight_ = var_set_wextra()
        iw_tweight=tweight_
      end if
    else
      tcoff_ = -1
    end if
 
    nvector      = 0 ! No. vector vars
 
    ! Check whether custom flux types have been defined
    if(.not. allocated(flux_type)) then
       allocate(flux_type(ndir, nwflux))
       flux_type = flux_default
    else if(any(shape(flux_type) /= [ndir, nwflux])) then
       call mpistop("phys_check error: flux_type has wrong shape")
    end if
 
    phys_get_dt              => ffhd_get_dt
    phys_get_cmax            => ffhd_get_cmax_origin
    phys_get_a2max           => ffhd_get_a2max
    phys_get_cs2max          => ffhd_get_cs2max
    phys_get_tcutoff         => ffhd_get_tcutoff
    phys_get_cbounds         => ffhd_get_cbounds
    phys_to_primitive        => ffhd_to_primitive_origin
    ffhd_to_primitive        => ffhd_to_primitive_origin
    phys_to_conserved        => ffhd_to_conserved_origin
    ffhd_to_conserved        => ffhd_to_conserved_origin
    phys_get_flux            => ffhd_get_flux
    phys_get_v               => ffhd_get_v_origin
    ffhd_get_v               => ffhd_get_v_origin
    phys_get_rho             => ffhd_get_rho
    ffhd_kin_en              => ffhd_kin_en_origin
    !> need to check source geom here
    !phys_add_source_geom     => ffhd_add_source_geom
    phys_add_source          => ffhd_add_source
    phys_check_params        => ffhd_check_params
    phys_write_info          => ffhd_write_info
    phys_handle_small_values => ffhd_handle_small_values_origin
    ffhd_handle_small_values => ffhd_handle_small_values_origin
    phys_check_w             => ffhd_check_w_origin
 
    if(.not.ffhd_energy) then
      phys_get_pthermal      => ffhd_get_pthermal_iso
      ffhd_get_pthermal      => ffhd_get_pthermal_iso
    else
      phys_get_pthermal      => ffhd_get_pthermal_origin
      ffhd_get_pthermal      => ffhd_get_pthermal_origin
    end if
 
    ! choose Rfactor in ideal gas law
    if(ffhd_partial_ionization) then
      ffhd_get_rfactor=>rfactor_from_temperature_ionization
      phys_update_temperature => ffhd_update_temperature
    else if(associated(usr_rfactor)) then
      ffhd_get_rfactor=>usr_rfactor
    else
      ffhd_get_rfactor=>rfactor_from_constant_ionization
    end if
 
    if(ffhd_partial_ionization) then
      ffhd_get_temperature => ffhd_get_temperature_from_te
    else
      ffhd_get_temperature => ffhd_get_temperature_from_etot
    end if
 
    ! derive units from basic units
    call ffhd_physical_units()
 
    if(ffhd_hyperbolic_thermal_conduction) then
      hypertc_kappa=8.d-7*unit_temperature**3.5d0/unit_length/unit_density/unit_velocity**3
    end if
    if(.not. ffhd_energy .and. ffhd_thermal_conduction) then
      call mpistop("thermal conduction needs ffhd_energy=T")
    end if
    if(.not. ffhd_energy .and. ffhd_hyperbolic_thermal_conduction) then
      call mpistop("hyperbolic thermal conduction needs ffhd_energy=T")
    end if
    if(.not. ffhd_energy .and. ffhd_radiative_cooling) then
      call mpistop("radiative cooling needs ffhd_energy=T")
    end if
 
    ! initialize thermal conduction module
    if(ffhd_thermal_conduction) then
      phys_req_diagonal = .true.
 
      call sts_init()
      call tc_init_params(ffhd_gamma)
 
      allocate(tc_fl)
      call tc_get_ffhd_params(tc_fl,tc_params_read_ffhd)
      call add_sts_method(ffhd_get_tc_dt_ffhd,ffhd_sts_set_source_tc_ffhd,e_,1,e_,1,.false.)
      tc_fl%get_temperature_from_conserved => ffhd_get_temperature_from_etot
      tc_fl%get_temperature_from_eint => ffhd_get_temperature_from_eint
      call set_conversion_methods_to_head(ffhd_e_to_ei, ffhd_ei_to_e)
      call set_error_handling_to_head(ffhd_tc_handle_small_e)
      tc_fl%get_rho => ffhd_get_rho
      tc_fl%e_ = e_
      tc_fl%Tcoff_ = tcoff_
    end if
 
    ! Initialize radiative cooling module
    if(ffhd_radiative_cooling) then
      call radiative_cooling_init_params(ffhd_gamma,he_abundance)
      allocate(rc_fl)
      call radiative_cooling_init(rc_fl,rc_params_read)
      rc_fl%get_rho => ffhd_get_rho
      rc_fl%get_pthermal => ffhd_get_pthermal
      rc_fl%get_var_Rfactor => ffhd_get_rfactor
      rc_fl%e_ = e_
      rc_fl%Tcoff_ = tcoff_
      rc_fl%has_equi = .false.
    end if
    allocate(te_fl_ffhd)
    te_fl_ffhd%get_rho=> ffhd_get_rho
    te_fl_ffhd%get_pthermal=> ffhd_get_pthermal
    te_fl_ffhd%get_var_Rfactor => ffhd_get_rfactor
{^ifthreed
    phys_te_images => ffhd_te_images
}
    ! Initialize viscosity module
    if(ffhd_viscosity) call viscosity_init(phys_wider_stencil,phys_req_diagonal)
 
    ! Initialize gravity module
    if(ffhd_gravity) then
      call gravity_init()
    end if
 
    ! initialize ionization degree table
    if(ffhd_partial_ionization) call ionization_degree_init()
  end subroutine ffhd_phys_init
 
{^ifthreed
  subroutine ffhd_te_images
    use mod_global_parameters
    use mod_thermal_emission
 
    select case(convert_type)
      case('EIvtiCCmpi','EIvtuCCmpi')
        call get_euv_image(unitconvert,te_fl_ffhd)
      case('ESvtiCCmpi','ESvtuCCmpi')
        call get_euv_spectrum(unitconvert,te_fl_ffhd)
      case('SIvtiCCmpi','SIvtuCCmpi')
        call get_sxr_image(unitconvert,te_fl_ffhd)
      case('WIvtiCCmpi','WIvtuCCmpi')
        call get_whitelight_image(unitconvert,te_fl_ffhd)
      case default
        call mpistop("Error in synthesize emission: Unknown convert_type")
      end select
  end subroutine ffhd_te_images
}
 
  subroutine  ffhd_sts_set_source_tc_ffhd(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_ffhd
    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_ffhd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl)
  end subroutine ffhd_sts_set_source_tc_ffhd
 
  function ffhd_get_tc_dt_ffhd(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
    !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
    !where                      tc_k_para_i=tc_k_para*B_i**2/B**2
    !and                        T=p/rho
    use mod_global_parameters
    use mod_thermal_conduction, only: get_tc_dt_ffhd
 
    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_ffhd(w,ixi^l,ixo^l,dx^d,x,tc_fl)
  end function ffhd_get_tc_dt_ffhd
 
  subroutine ffhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
    use mod_global_parameters
 
    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
    character(len=140) :: error_msg
 
    write(error_msg,"(a,i3)") "Thermal conduction step ", step
    call ffhd_handle_small_ei(w,x,ixi^l,ixo^l,e_,error_msg)
  end subroutine ffhd_tc_handle_small_e
 
  subroutine tc_params_read_ffhd(fl)
    use mod_global_parameters, only: unitpar,par_files
    type(tc_fluid), intent(inout) :: fl
    integer                      :: n
    ! list parameters
    logical :: tc_saturate=.false.
    double precision :: tc_k_para=0d0
    character(len=std_len)  :: tc_slope_limiter="MC"
 
    namelist /tc_list/ tc_saturate, tc_slope_limiter, 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
    select case(tc_slope_limiter)
     case ('no','none')
       fl%tc_slope_limiter = 0
     case ('MC')
       ! montonized central limiter Woodward and Collela limiter (eq.3.51h), a factor of 2 is pulled out
       fl%tc_slope_limiter = 1
     case('minmod')
       ! minmod limiter
       fl%tc_slope_limiter = 2
     case ('superbee')
       ! Roes superbee limiter (eq.3.51i)
       fl%tc_slope_limiter = 3
     case ('koren')
       ! Barry Koren Right variant
       fl%tc_slope_limiter = 4
     case default
       call mpistop("Unknown tc_slope_limiter, choose MC, minmod")
    end select
  end subroutine tc_params_read_ffhd
 
  subroutine rc_params_read(fl)
    use mod_global_parameters, only: unitpar,par_files
    use mod_constants, only: bigdouble
    type(rc_fluid), intent(inout) :: fl
    integer                      :: n
    integer :: ncool = 4000
    double precision :: cfrac=0.1d0
  
    !> Name of cooling curve
    character(len=std_len)  :: coolcurve='JCcorona'
  
    !> Name of cooling method
    character(len=std_len)  :: coolmethod='exact'
  
    !> 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.
    logical    :: rad_cut=.false.
    double precision :: rad_cut_hgt=0.5d0
    double precision :: rad_cut_dey=0.15d0
 
    namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split, rad_cut, rad_cut_hgt, rad_cut_dey
 
    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%coolmethod=coolmethod
    fl%tlow=tlow
    fl%Tfix=tfix
    fl%rc_split=rc_split
    fl%cfrac=cfrac
    fl%rad_cut=rad_cut
    fl%rad_cut_hgt=rad_cut_hgt
    fl%rad_cut_dey=rad_cut_dey
  end subroutine rc_params_read
 
  subroutine ffhd_check_params
    use mod_global_parameters
    use mod_usr_methods
    use mod_convert, only: add_convert_method
 
    gamma_1=ffhd_gamma-1.d0
    if (.not. ffhd_energy) then
       if (ffhd_gamma <= 0.0d0) call mpistop ("Error: ffhd_gamma <= 0")
       if (ffhd_adiab < 0.0d0) call mpistop ("Error: ffhd_adiab < 0")
       small_pressure = ffhd_adiab*small_density**ffhd_gamma
    else
       if (ffhd_gamma <= 0.0d0 .or. ffhd_gamma == 1.0d0) &
            call mpistop ("Error: ffhd_gamma <= 0 or ffhd_gamma == 1")
       inv_gamma_1=1.d0/gamma_1
       small_e = small_pressure * inv_gamma_1
    end if
 
    if (number_equi_vars > 0 .and. .not. associated(usr_set_equi_vars)) then
      call mpistop("usr_set_equi_vars has to be implemented in the user file")
    end if
  end subroutine ffhd_check_params
 
  subroutine ffhd_physical_units()
    use mod_global_parameters
    double precision :: mp,kB
    double precision :: a,b
 
    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(ffhd_partial_ionization) then
        b = 2+.3d0
      else
        b = 1d0 + h_ion_fr + he_abundance*(he_ion_fr*(he_ion_fr2 + 1d0)+1d0)
      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) then
      unit_numberdensity=unit_density/(a*mp)
    else
      unit_density=a*mp*unit_numberdensity
    end if
    if(unit_velocity/=1.d0) then
      unit_pressure=unit_density*unit_velocity**2
      unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
    else if(unit_pressure/=1.d0) then
      unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
      unit_velocity=sqrt(unit_pressure/unit_density)
    else
      unit_pressure=b*unit_numberdensity*kb*unit_temperature
      unit_velocity=sqrt(unit_pressure/unit_density)
    end if
    if(unit_time/=1.d0) then
      unit_length=unit_time*unit_velocity
    else
      unit_time=unit_length/unit_velocity
    end if
    unit_mass=unit_density*unit_length**3
  end subroutine ffhd_physical_units
 
  subroutine ffhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
    use mod_global_parameters
    logical, intent(in) :: primitive
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S,nw)
    double precision :: tmp(ixI^S)
    logical, intent(inout) :: flag(ixI^S,1:nw)
 
    flag=.false.
    where(w(ixo^s,rho_) < small_density) flag(ixo^s,rho_) = .true.
 
    if(ffhd_energy) then
      if(primitive) then
        where(w(ixo^s,e_) < small_pressure) flag(ixo^s,e_) = .true.
      else
        tmp(ixo^s)=w(ixo^s,e_)-ffhd_kin_en(w,ixi^l,ixo^l)
        where(tmp(ixo^s) < small_e) flag(ixo^s,e_) = .true.
      end if
    end if
  end subroutine ffhd_check_w_origin
 
  subroutine ffhd_to_conserved_origin(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)
    double precision :: inv_gamma2(ixO^S)
    integer                         :: idir
 
    if(ffhd_energy) then
      w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1+half*w(ixo^s,mom(1))**2*w(ixo^s,rho_)
    end if
    w(ixo^s,mom(1))=w(ixo^s,rho_)*w(ixo^s,mom(1))
  end subroutine ffhd_to_conserved_origin
 
  subroutine ffhd_to_primitive_origin(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)
    double precision                :: inv_rho(ixO^S), gamma2(ixO^S)
 
    if(fix_small_values) then
      !> fix small values preventing NaN numbers in the following converting
      call ffhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'ffhd_to_primitive_origin')
    end if
 
    w(ixo^s,mom(1)) = w(ixo^s,mom(1))/w(ixo^s,rho_)
    if(ffhd_energy) then
      w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)-half*w(ixo^s,rho_)*w(ixo^s,mom(1))**2)
    end if
  end subroutine ffhd_to_primitive_origin
 
  subroutine ffhd_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)
 
    w(ixi^s,e_)=w(ixi^s,e_)+ffhd_kin_en(w,ixi^l,ixi^l)
  end subroutine ffhd_ei_to_e
 
  subroutine ffhd_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)
 
    w(ixi^s,e_)=w(ixi^s,e_)-ffhd_kin_en(w,ixi^l,ixi^l)
    if(fix_small_values) then
      call ffhd_handle_small_ei(w,x,ixi^l,ixi^l,e_,'ffhd_e_to_ei')
    end if
  end subroutine ffhd_e_to_ei
 
  subroutine ffhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
    use mod_global_parameters
    use mod_small_values
    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
 
    logical :: flag(ixI^S,1:nw)
    double precision :: tmp2(ixI^S)
 
    call phys_check_w(primitive, ixi^l, ixi^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
        if(small_values_fix_iw(mom(1))) then
          where(flag(ixo^s,rho_)) w(ixo^s, mom(1)) = 0.0d0
        end if
        if(ffhd_energy) then
          if(primitive) then
            where(flag(ixo^s,e_)) w(ixo^s,p_) = small_pressure
          else
            where(flag(ixo^s,e_))
              w(ixo^s,e_) = small_e+ffhd_kin_en(w,ixi^l,ixo^l)
            end where
          end if
        end if
      case ("average")
        call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
        if(ffhd_energy) then
          if(primitive) then
            call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
          else
            w(ixi^s,e_)=w(ixi^s,e_)-ffhd_kin_en(w,ixi^l,ixi^l)
            call small_values_average(ixi^l, ixo^l, w, x, flag, e_)
            w(ixi^s,e_)=w(ixi^s,e_)+ffhd_kin_en(w,ixi^l,ixi^l)
          end if
        end if
      case default
        if(.not.primitive) then
          if(ffhd_energy) then
            w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)-ffhd_kin_en(w,ixi^l,ixo^l))
          end if
          w(ixo^s,mom(1))=w(ixo^s,mom(1))/w(ixo^s,rho_)
        end if
        call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
      end select
    end if
  end subroutine ffhd_handle_small_values_origin
 
  subroutine ffhd_get_v_origin(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:ndim)
    double precision, intent(out) :: v(ixI^S,ndir)
    double precision :: rho(ixI^S)
    integer :: idir
 
    call ffhd_get_rho(w,x,ixi^l,ixo^l,rho)
    rho(ixo^s)=1.d0/rho(ixo^s)
    do idir=1,ndir
      v(ixo^s,ndir) = w(ixo^s,mom(1))*block%B0(ixo^s,idir,0)*rho(ixo^s)
    end do
  end subroutine ffhd_get_v_origin
 
  subroutine ffhd_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)
    double precision              :: rho(ixi^s)
 
    call ffhd_get_rho(w,x,ixi^l,ixo^l,rho)
    v(ixo^s) = (w(ixo^s, mom(1))*block%B0(ixo^s,idim,0)) / rho(ixo^s)
  end subroutine ffhd_get_v_idim
 
  subroutine ffhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
    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(inout) :: cmax(ixI^S)
    double precision :: vel(ixI^S)
 
    call ffhd_get_csound(w,x,ixi^l,ixo^l,idim,cmax)
    call ffhd_get_v_idim(w,x,ixi^l,ixo^l,idim,vel)
 
    cmax(ixo^s)=abs(vel(ixo^s))+cmax(ixo^s)
  end subroutine ffhd_get_cmax_origin
 
  subroutine ffhd_get_cs2max(w,x,ixI^L,ixO^L,cs2max)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: cs2max
    double precision :: cs2(ixI^S)
 
    call ffhd_get_csound2(w,x,ixi^l,ixo^l,cs2)
    cs2max=maxval(cs2(ixo^s))
  end subroutine ffhd_get_cs2max
 
  subroutine ffhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: a2max(ndim)
    double precision :: a2(ixI^S,ndim,nw)
    integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
 
    a2=zero
    do i = 1,ndim
      !> 4th order
      hxo^l=ixo^l-kr(i,^d);
      gxo^l=hxo^l-kr(i,^d);
      jxo^l=ixo^l+kr(i,^d);
      kxo^l=jxo^l+kr(i,^d);
      a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
         -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
      a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/dxlevel(i)**2
    end do
  end subroutine ffhd_get_a2max
 
  subroutine ffhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in) :: ixI^L,ixO^L
    double precision, intent(in) :: x(ixI^S,1:ndim)
    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, dimension(ixI^S,1:ndir) :: bunitvec
    double precision, dimension(ixI^S,1:ndim) :: gradT
    double precision :: Bdir(ndim)
    double precision :: ltrc,ltrp,altr(ixI^S)
    integer :: idims,jxO^L,hxO^L,ixA^D,ixB^D
    integer :: jxP^L,hxP^L,ixP^L,ixQ^L
    logical :: lrlt(ixI^S)
 
    call ffhd_get_temperature(w,x,ixi^l,ixi^l,te)
    tco_local=zero
    tmax_local=maxval(te(ixo^s))
 
    {^ifoned
    select case(ffhd_trac_type)
    case(0)
      !> test case, fixed cutoff temperature
      block%wextra(ixi^s,tcoff_)=2.5d5/unit_temperature
    case(1)
      hxo^l=ixo^l-1;
      jxo^l=ixo^l+1;
      lts(ixo^s)=0.5d0*abs(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=4.d0
      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)
      lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
      lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
      block%wextra(ixo^s,tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
    case default
      call mpistop("ffhd_trac_type not allowed for 1D simulation")
    end select
    }
    {^nooned
    select case(ffhd_trac_type)
    case(0)
      !> test case, fixed cutoff temperature
      if(slab_uniform) then
        !> assume cgs units
        block%wextra(ixi^s,tcoff_)=max(min(3.d5/unit_temperature,6.d5/unit_temperature-3.d-4/unit_temperature*unit_length*x(ixi^s,ndim)),zero)
      else
        block%wextra(ixi^s,tcoff_)=2.5d5/unit_temperature
      end if
    case(1,4,6)
      do idims=1,ndim
        call gradient(te,ixi^l,ixo^l,idims,tmp1)
        gradt(ixo^s,idims)=tmp1(ixo^s)
      end do
      bunitvec(ixo^s,:)=block%B0(ixo^s,:,0)
      if(ffhd_trac_type .gt. 1) then
        ! B direction at cell center
        bdir=zero
        {do ixa^d=0,1\}
          ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
          bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
        {end do\}
        if(sum(bdir(:)**2) .gt. zero) then
          bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
        end if
        block%special_values(3:ndim+2)=bdir(1:ndim)
      end if
      tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
      where(tmp1(ixo^s)/=0.d0)
        tmp1(ixo^s)=1.d0/tmp1(ixo^s)
      else where
        tmp1(ixo^s)=bigdouble
      end where
      ! b unit vector: magnetic field direction vector
      do idims=1,ndim
        bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
      end do
      ! temperature length scale inversed
      lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
      ! fraction of cells size to temperature length scale
      if(slab_uniform) then
        lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
      else
        lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
      end if
      lrlt=.false.
      where(lts(ixo^s) > trac_delta)
        lrlt(ixo^s)=.true.
      end where
      if(any(lrlt(ixo^s))) then
        block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
      else
        block%special_values(1)=zero
      end if
      block%special_values(2)=tmax_local
    case(2)
      !> iijima et al. 2021, LTRAC method
      ltrc=1.5d0
      ltrp=4.d0
      ixp^l=ixo^l^ladd2;
      do idims=1,ndim
        ixq^l=ixp^l;
        hxp^l=ixp^l;
        jxp^l=ixp^l;
        select case(idims)
          {case(^d)
           ixqmin^d=ixqmin^d+1
           ixqmax^d=ixqmax^d-1
           hxpmax^d=ixpmin^d
           jxpmin^d=ixpmax^d
          \}
        end select
        call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
        call gradientx(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),.false.)
        call gradientq(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims))
      end do
      bunitvec(ixp^s,:)=block%B0(ixp^s,:,0)
      lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
      if(slab_uniform) then
        lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
      else
        lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
      end if
      lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
  
      altr=zero
      ixp^l=ixo^l^ladd1;
      do idims=1,ndim
        hxo^l=ixp^l-kr(idims,^d);
        jxo^l=ixp^l+kr(idims,^d);
        altr(ixp^s)=altr(ixp^s)+0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
      end do
      block%wextra(ixp^s,tcoff_)=te(ixp^s)*altr(ixp^s)**0.4d0
    case(3,5)
      !> do nothing here
    case default
      call mpistop("unknown ffhd_trac_type")
    end select
    }
  end subroutine ffhd_get_tcutoff
 
  subroutine ffhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L, idim
    double precision, intent(in)    :: wLC(ixI^S, nw), wRC(ixI^S, nw)
    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)=sqrt(wlp(ixo^s,rho_))
      tmp2(ixo^s)=sqrt(wrp(ixo^s,rho_))
      tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
      umean(ixo^s)=(wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim)*tmp1(ixo^s)&
                   +wrp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim)*tmp2(ixo^s))*tmp3(ixo^s)
      call ffhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
      call ffhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
      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(1))*block%B0(ixo^s,idim,idim)-wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,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)
      else
        cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
      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(1))*block%B0(ixo^s,idim,idim)/wmean(ixo^s,rho_)
      call ffhd_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
      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)
      else
        cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
      end if
    case (3)
      ! Miyoshi 2005 JCP 208, 315 equation (67)
      call ffhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
      call ffhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
      csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
      if(present(cmin)) then
        cmin(ixo^s,1)=min(wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim),&
                          wrp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim))-csoundl(ixo^s)
        cmax(ixo^s,1)=max(wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim),&
                          wrp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim))+csoundl(ixo^s)
      else
        cmax(ixo^s,1)=max(wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim),&
                          wrp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim))+csoundl(ixo^s)
      end if
    end select
  end subroutine ffhd_get_cbounds
 
  subroutine ffhd_get_csound(w,x,ixI^L,ixO^L,idim,csound)
    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):: csound(ixI^S)
 
    call ffhd_get_csound2(w,x,ixi^l,ixo^l,csound)
    csound(ixo^s) = dsqrt(csound(ixo^s))*abs(block%B0(ixo^s,idim,idim))
  end subroutine ffhd_get_csound
 
  !> Calculate fast magnetosonic wave speed
  subroutine ffhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
    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):: csound(ixI^S)
 
    if(ffhd_energy) then
      csound(ixo^s)=ffhd_gamma*w(ixo^s,e_)/w(ixo^s,rho_)
    else
      csound(ixo^s)=ffhd_gamma*ffhd_adiab*w(ixo^s,rho_)**gamma_1
    end if
    csound(ixo^s) = dsqrt(csound(ixo^s))
  end subroutine ffhd_get_csound_prim
 
  subroutine ffhd_get_pthermal_iso(w,x,ixI^L,ixO^L,pth)
    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):: pth(ixI^S)
 
    call ffhd_get_rho(w,x,ixi^l,ixo^l,pth)
    pth(ixo^s)=ffhd_adiab*pth(ixo^s)**ffhd_gamma
  end subroutine ffhd_get_pthermal_iso
 
  subroutine ffhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
    use mod_global_parameters
    use mod_small_values, only: trace_small_values
    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):: pth(ixI^S)
    integer                      :: iw, ix^D
 
    pth(ixo^s)=gamma_1*(w(ixo^s,e_)-ffhd_kin_en(w,ixi^l,ixo^l))
    if (fix_small_values) then
      {do ix^db= ixo^lim^db\}
         if(pth(ix^d)<small_pressure) then
            pth(ix^d)=small_pressure
         end if
      {end do^D&\}
    elseif(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 ffhd_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(*,*) sqrt(pth(ix^d)-bigdouble)
           write(*,*) "Saving status at the previous time step"
           crash=.true.
         end if
      {end do^D&\}
    end if
  end subroutine ffhd_get_pthermal_origin
 
  subroutine ffhd_get_temperature_from_te(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)
 
    res(ixo^s) = w(ixo^s, te_)
  end subroutine ffhd_get_temperature_from_te
 
  subroutine ffhd_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 ffhd_get_rfactor(w,x,ixi^l,ixo^l,r)
    res(ixo^s) = gamma_1 * w(ixo^s, e_)/(w(ixo^s,rho_)*r(ixo^s))
  end subroutine ffhd_get_temperature_from_eint
 
  subroutine ffhd_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 ffhd_get_rfactor(w,x,ixi^l,ixo^l,r)
    call ffhd_get_pthermal(w,x,ixi^l,ixo^l,res)
    res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,rho_))
  end subroutine ffhd_get_temperature_from_etot
 
  subroutine ffhd_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)
    double precision    :: rho(ixi^s)
    
    call ffhd_get_rho(w,x,ixi^l,ixo^l,rho)
    if(ffhd_energy) then
      call ffhd_get_pthermal(w,x,ixi^l,ixo^l,csound2)
      csound2(ixo^s)=ffhd_gamma*csound2(ixo^s)/rho(ixo^s)
    else
      csound2(ixo^s)=ffhd_gamma*ffhd_adiab*rho(ixo^s)**gamma_1
    end if
  end subroutine ffhd_get_csound2
 
  subroutine ffhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in)          :: ixI^L, ixO^L, idim
    ! conservative w
    double precision, intent(in) :: wC(ixI^S,nw)
    ! primitive w
    double precision, intent(in) :: w(ixI^S,nw)
    double precision, intent(in) :: x(ixI^S,1:ndim)
    double precision,intent(out) :: f(ixI^S,nwflux)
    double precision             :: ptotal(ixO^S)
    double precision             :: tmp(ixI^S)
    integer                      :: idirmin, iw, idir, jdir, kdir
    double precision, dimension(ixI^S) :: Te,tau,sigT
 
    f(ixo^s,rho_)=w(ixo^s,mom(1))*w(ixo^s,rho_)*block%B0(ixo^s,idim,idim)
 
    if(ffhd_energy) then
      ptotal(ixo^s)=w(ixo^s,p_)
    else
      ptotal(ixo^s)=ffhd_adiab*w(ixo^s,rho_)**ffhd_gamma
    end if
 
    ! Get flux of momentum
    f(ixo^s,mom(1))=(wc(ixo^s,mom(1))*w(ixo^s,mom(1))+ptotal(ixo^s))*block%B0(ixo^s,idim,idim)
 
    ! Get flux of energy
    if(ffhd_energy) then
      f(ixo^s,e_)=w(ixo^s,mom(1))*(wc(ixo^s,e_)+ptotal(ixo^s))*block%B0(ixo^s,idim,idim)
      if(ffhd_hyperbolic_thermal_conduction) then
        f(ixo^s,e_)=f(ixo^s,e_)+w(ixo^s,q_)*block%B0(ixo^s,idim,idim)
        f(ixo^s,q_)=zero
      end if
    end if
  end subroutine ffhd_get_flux
 
  subroutine ffhd_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_viscosity, only: viscosity_add_source
    use mod_gravity, only: gravity_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
 
    if (.not. qsourcesplit) then
      active = .true.
      call add_punitb(qdt,ixi^l,ixo^l,wct,w,x,wctprim)
      if(ffhd_hyperbolic_thermal_conduction) then
        call add_hypertc_source(qdt,ixi^l,ixo^l,wct,w,x,wctprim)
      end if
    end if
 
    if(ffhd_radiative_cooling) then
      call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,wctprim,&
           w,x,qsourcesplit,active, rc_fl)
    end if
 
    if(ffhd_viscosity) then
      call viscosity_add_source(qdt,ixi^l,ixo^l,wct,&
           w,x,ffhd_energy,qsourcesplit,active)
    end if
 
    if(ffhd_gravity) then
      call gravity_add_source(qdt,ixi^l,ixo^l,wct,wctprim,&
           w,x,ffhd_energy,.false.,qsourcesplit,active)
    end if
 
    ! update temperature from new pressure, density, and old ionization degree
    if(ffhd_partial_ionization) then
      if(.not.qsourcesplit) then
        active = .true.
        call ffhd_update_temperature(ixi^l,ixo^l,wct,w,x)
      end if
    end if
  end subroutine ffhd_add_source
 
  subroutine add_punitb(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in)             :: ixI^L,ixO^L
    double precision, intent(in)    :: qdt
    double precision, intent(in)    :: wCT(ixI^S,1:nw),x(ixI^S,1:ndim)
    double precision, intent(in)    :: wCTprim(ixI^S,1:nw)
    double precision, intent(inout) :: w(ixI^S,1:nw)
 
    integer                         :: idims,hxO^L
    double precision                :: divb(ixI^S)
 
    divb=zero
    if(slab_uniform) then
      do idims=1,ndim
        hxo^l=ixo^l-kr(idims,^d);
        divb(ixo^s)=divb(ixo^s)+(block%B0(ixo^s,idims,idims)-block%B0(hxo^s,idims,idims))/dxlevel(idims)
      end do
    else
      call divvector(block%B0(ixi^s,1:ndir,0),ixi^l,ixo^l,divb)
    end if
    w(ixo^s,mom(1))=w(ixo^s,mom(1))+qdt*wctprim(ixo^s,p_)*divb(ixo^s)
  end subroutine add_punitb
 
  subroutine ffhd_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 ffhd_get_rho
 
  subroutine ffhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
    use mod_global_parameters
    use mod_small_values
    integer, intent(in)             :: ixI^L,ixO^L, ie
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    character(len=*), intent(in)    :: subname
    integer :: idir
    logical :: flag(ixI^S,1:nw)
    double precision              :: rho(ixI^S)
 
    flag=.false.
    where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
    if(any(flag(ixo^s,ie))) then
      select case (small_values_method)
      case ("replace")
        where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
      case ("average")
        call small_values_average(ixi^l, ixo^l, w, x, flag, ie)
      case default
        w(ixo^s,e_)=w(ixo^s,e_)*gamma_1
        call ffhd_get_rho(w,x,ixi^l,ixo^l,rho)
        w(ixo^s,mom(1)) = w(ixo^s,mom(1))/rho(ixo^s)
        call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
      end select
    end if
  end subroutine ffhd_handle_small_ei
 
  subroutine ffhd_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 ffhd_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 ffhd_update_temperature
 
  subroutine ffhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
    use mod_global_parameters
    use mod_usr_methods
    use mod_radiative_cooling, only: cooling_get_dt
    use mod_viscosity, only: viscosity_get_dt
    use mod_gravity, only: gravity_get_dt
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(inout) :: dtnew
    double precision, intent(in)    :: dx^D
    double precision, intent(in)    :: w(ixI^S,1:nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    integer                       :: idirmin,idim
    double precision              :: dxarr(ndim)
    double precision              :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
 
    dtnew = bigdouble
 
    if(ffhd_radiative_cooling) then
      call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x,rc_fl)
    end if
 
    if(ffhd_viscosity) then
      call viscosity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
    end if
 
    if(ffhd_gravity) then
      call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
    end if
  end subroutine ffhd_get_dt
 
!  subroutine ffhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
!    use mod_global_parameters
!    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), w(ixI^S,1:nw)
!
!    integer          :: iw,idir, h1x^L{^NOONED, h2x^L}
!    double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),invrho(ixO^S),invr(ixO^S)
!
!    integer :: mr_,mphi_ ! Polar var. names
!    integer :: br_,bphi_
!
!    mr_=mom(1); mphi_=mom(1)-1+phi_  ! Polar var. names
!    br_=mag(1); bphi_=mag(1)-1+phi_
!
!    ! 1/rho
!    invrho(ixO^S)=1.d0/wCT(ixO^S,rho_)
!    ! include dt in invr, invr is always used with qdt
!    if(local_timestep) then
!      invr(ixO^S) = block%dt(ixO^S) * dtfactor/x(ixO^S,1)
!    else
!      invr(ixO^S) = qdt/x(ixO^S,1)
!    end if  
!
!
!    select case (coordinate)
!    case (cylindrical)
!      call ffhd_get_p_total(wCT,x,ixI^L,ixO^L,tmp)
!      if(phi_>0) then
!        w(ixO^S,mr_)=w(ixO^S,mr_)+invr(ixO^S)*(tmp(ixO^S)-&
!                  wCT(ixO^S,bphi_)**2+wCT(ixO^S,mphi_)**2*invrho(ixO^S))
!        w(ixO^S,mphi_)=w(ixO^S,mphi_)+invr(ixO^S)*(&
!                 -wCT(ixO^S,mphi_)*wCT(ixO^S,mr_)*invrho(ixO^S) &
!                 +wCT(ixO^S,bphi_)*wCT(ixO^S,br_))
!        if(.not.stagger_grid) then
!          w(ixO^S,bphi_)=w(ixO^S,bphi_)+invr(ixO^S)*&
!                   (wCT(ixO^S,bphi_)*wCT(ixO^S,mr_) &
!                   -wCT(ixO^S,br_)*wCT(ixO^S,mphi_)) &
!                   *invrho(ixO^S)
!        end if
!      else
!        w(ixO^S,mr_)=w(ixO^S,mr_)+invr(ixO^S)*tmp(ixO^S)
!      end if
!      if(ffhd_glm) w(ixO^S,br_)=w(ixO^S,br_)+wCT(ixO^S,psi_)*invr(ixO^S)
!    case (spherical)
!       h1x^L=ixO^L-kr(1,^D); {^NOONED h2x^L=ixO^L-kr(2,^D);}
!       call ffhd_get_p_total(wCT,x,ixI^L,ixO^L,tmp1)
!       ! m1
!       tmp(ixO^S)=tmp1(ixO^S)*x(ixO^S,1) &
!                  *(block%surfaceC(ixO^S,1)-block%surfaceC(h1x^S,1))/block%dvolume(ixO^S)
!       do idir=2,ndir
!         tmp(ixO^S)=tmp(ixO^S)+wCT(ixO^S,mom(idir))**2*invrho(ixO^S)-wCT(ixO^S,mag(idir))**2
!       end do
!       w(ixO^S,mom(1))=w(ixO^S,mom(1))+tmp(ixO^S)*invr(ixO^S)
!       ! b1
!       if(ffhd_glm) then
!         w(ixO^S,mag(1))=w(ixO^S,mag(1))+invr(ixO^S)*2.0d0*wCT(ixO^S,psi_)
!       end if
!
!       {^NOONED
!       ! m2
!       ! This will make hydrostatic p=const an exact solution
!       if(local_timestep) then
!          tmp(ixO^S) = block%dt(ixO^S) * tmp1(ixO^S)
!       else
!          tmp(ixO^S) = qdt * tmp1(ixO^S)
!       end if  
!       w(ixO^S,mom(2))=w(ixO^S,mom(2))+tmp(ixO^S) &
!            *(block%surfaceC(ixO^S,2)-block%surfaceC(h2x^S,2)) &
!            /block%dvolume(ixO^S)
!       tmp(ixO^S)=-(wCT(ixO^S,mom(1))*wCT(ixO^S,mom(2))*invrho(ixO^S) &
!            -wCT(ixO^S,mag(1))*wCT(ixO^S,mag(2)))
!       if(ndir==3) then
!         tmp(ixO^S)=tmp(ixO^S)+(wCT(ixO^S,mom(3))**2*invrho(ixO^S) &
!              -wCT(ixO^S,mag(3))**2)*dcos(x(ixO^S,2))/dsin(x(ixO^S,2))
!       end if
!       w(ixO^S,mom(2))=w(ixO^S,mom(2))+tmp(ixO^S)*invr(ixO^S)
!       ! b2
!       if(.not.stagger_grid) then
!         tmp(ixO^S)=(wCT(ixO^S,mom(1))*wCT(ixO^S,mag(2)) &
!              -wCT(ixO^S,mom(2))*wCT(ixO^S,mag(1)))*invrho(ixO^S)
!         if(ffhd_glm) then
!           tmp(ixO^S)=tmp(ixO^S) &
!                + dcos(x(ixO^S,2))/dsin(x(ixO^S,2))*wCT(ixO^S,psi_)
!         end if
!         w(ixO^S,mag(2))=w(ixO^S,mag(2))+tmp(ixO^S)*invr(ixO^S)
!       end if
!       }
!
!       if(ndir==3) then
!         ! m3
!         tmp(ixO^S)=-(wCT(ixO^S,mom(3))*wCT(ixO^S,mom(1))*invrho(ixO^S) &
!              -wCT(ixO^S,mag(3))*wCT(ixO^S,mag(1))) {^NOONED &
!              -(wCT(ixO^S,mom(2))*wCT(ixO^S,mom(3))*invrho(ixO^S) &
!              -wCT(ixO^S,mag(2))*wCT(ixO^S,mag(3))) &
!              *dcos(x(ixO^S,2))/dsin(x(ixO^S,2)) }
!         w(ixO^S,mom(3))=w(ixO^S,mom(3))+tmp(ixO^S)*invr(ixO^S)
!         ! b3
!         if(.not.stagger_grid) then
!           tmp(ixO^S)=(wCT(ixO^S,mom(1))*wCT(ixO^S,mag(3)) &
!                -wCT(ixO^S,mom(3))*wCT(ixO^S,mag(1)))*invrho(ixO^S) {^NOONED &
!                -(wCT(ixO^S,mom(3))*wCT(ixO^S,mag(2)) &
!                -wCT(ixO^S,mom(2))*wCT(ixO^S,mag(3)))*dcos(x(ixO^S,2)) &
!                /(wCT(ixO^S,rho_)*dsin(x(ixO^S,2))) }
!           w(ixO^S,mag(3))=w(ixO^S,mag(3))+tmp(ixO^S)*invr(ixO^S)
!         end if
!       end if
!    end select
!  end subroutine ffhd_add_source_geom
 
  function ffhd_kin_en_origin(w, ixI^L, ixO^L, inv_rho) result(ke)
    use mod_global_parameters, only: nw, ndim,block
    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(ixo^s)=0.5d0*w(ixo^s,mom(1))**2*inv_rho(ixo^s)
    else
      ke(ixo^s)=0.5d0*w(ixo^s,mom(1))**2/w(ixo^s,rho_)
    end if
  end function ffhd_kin_en_origin
 
  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)
    rfactor(ixo^s)=(1.d0+iz_h(ixo^s)+0.1d0*(1.d0+iz_he(ixo^s)*(1.d0+iz_he(ixo^s))))/(2.d0+3.d0*he_abundance)
  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 get_tau(ixI^L,ixO^L,w,Te,tau,sigT5)
    use mod_global_parameters
    integer, intent(in) :: ixI^L, ixO^L
    double precision, dimension(ixI^S,1:nw), intent(in) :: w
    double precision, dimension(ixI^S), intent(in) :: Te
    double precision, dimension(ixI^S), intent(out) :: tau,sigT5
    integer :: ix^D
    double precision :: dxmin,taumin
    double precision, dimension(ixI^S) :: sigT7,eint
 
    taumin=4.d0
    !> w supposed to be wCTprim here
    if(ffhd_trac) then
      where(te(ixo^s) .lt. block%wextra(ixo^s,tcoff_))
        sigt5(ixo^s)=hypertc_kappa*sqrt(block%wextra(ixo^s,tcoff_)**5)
        sigt7(ixo^s)=sigt5(ixo^s)*block%wextra(ixo^s,tcoff_)
      else where
        sigt5(ixo^s)=hypertc_kappa*sqrt(te(ixo^s)**5)
        sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
      end where
    else
      sigt5(ixo^s)=hypertc_kappa*sqrt(te(ixo^s)**5)
      sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
    end if
    eint(ixo^s)=w(ixo^s,p_)/(ffhd_gamma-one)
    tau(ixo^s)=max(taumin*dt,sigt7(ixo^s)/eint(ixo^s)/cs2max_global)
  end subroutine get_tau
 
  subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
    use mod_global_parameters
    integer, intent(in) :: ixI^L,ixO^L
    double precision, intent(in) :: qdt
    double precision, dimension(ixI^S,1:ndim), intent(in) :: x
    double precision, dimension(ixI^S,1:nw), intent(in) :: wCT,wCTprim
    double precision, dimension(ixI^S,1:nw), intent(inout) :: w
    integer :: idims
    integer :: hxC^L,hxO^L,ixC^L,jxC^L,jxO^L,kxC^L
    double precision :: invdx
    double precision, dimension(ixI^S) :: Te,tau,sigT,htc_qsrc,Tface
    double precision, dimension(ixI^S) :: htc_esrc
 
    te(ixi^s)=wctprim(ixi^s,p_)/wct(ixi^s,rho_)
    call get_tau(ixi^l,ixo^l,wctprim,te,tau,sigt)
    htc_qsrc=zero
    do idims=1,ndim
      invdx=1.d0/dxlevel(idims)
      ixc^l=ixo^l;
      ixcmin^d=ixomin^d-kr(idims,^d);ixcmax^d=ixomax^d;
      jxc^l=ixc^l+kr(idims,^d);
      kxc^l=jxc^l+kr(idims,^d);
      hxc^l=ixc^l-kr(idims,^d);
      hxo^l=ixo^l-kr(idims,^d);
      tface(ixc^s)=(7.d0*(te(ixc^s)+te(jxc^s))-(te(hxc^s)+te(kxc^s)))/12.d0
      htc_qsrc(ixo^s)=htc_qsrc(ixo^s)+sigt(ixo^s)*block%B0(ixo^s,idims,0)*(tface(ixo^s)-tface(hxo^s))*invdx
    end do
    htc_qsrc(ixo^s)=(htc_qsrc(ixo^s)+wct(ixo^s,q_))/tau(ixo^s)
    w(ixo^s,q_)=w(ixo^s,q_)-qdt*htc_qsrc(ixo^s)
  end subroutine add_hypertc_source
end module mod_ffhd_phys