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_eos
  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_tc = .false.
  !> Whether saturation is considered for hyperbolic TC
  logical, public, protected              :: ffhd_hyperbolic_tc_sat = .false.
  !> Whether the perpendicular hyperbolic-TC channel is enabled
  logical, public, protected              :: ffhd_hyperbolic_tc_use_perp = .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 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
 
  !> 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 and electron number density (LTE stored aux state)
  integer, public, protected :: te_
  integer, public, protected :: ne_
 
  !> Index of the cutoff temperature for the TRAC method
  integer, public, protected              :: tcoff_
  integer, public, protected              :: tweight_
  integer, public, protected              :: q_
 
  !> The adiabatic index (now owned by eos%; use eos%gamma)
 
  !> The adiabatic constant
  double precision, public                :: ffhd_adiab = 1.0d0
 
  !> The thermal conductivity kappa in hyperbolic thermal conduction
  double precision, public                :: hyperbolic_tc_kappa
 
  !> Helium abundance over Hydrogen (now owned by eos%; use eos%He_abundance)
  !> 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.
 
 
  !define the function interface for the kinetic energy
  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, public :: 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_ei
  public :: ffhd_get_pthermal
  public :: ffhd_get_rfactor
  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
  ! concrete routines referenced by the mod_ffhd_eos seam
  public :: ffhd_to_conserved_origin
  public :: ffhd_to_primitive_origin
  public :: ffhd_get_pthermal_origin
  public :: ffhd_get_temperature_from_etot
  public :: ffhd_get_temperature_from_eint
  public :: ffhd_get_temperature_from_te
  public :: rfactor_from_constant_ionization
 
contains
 
  subroutine ffhd_read_params(files)
    use mod_global_parameters
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    ! gamma and eos%He_abundance migrated to &eos_list (eos%gamma / eos%He_abundance)
    namelist /ffhd_list/ ffhd_energy, ffhd_adiab, &
      ffhd_thermal_conduction, ffhd_hyperbolic_tc, ffhd_hyperbolic_tc_sat, ffhd_radiative_cooling, ffhd_gravity,&
      h_ion_fr, he_ion_fr, he_ion_fr2, si_unit,&
      b0field, busr, 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) = eos%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_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_eos_pi_tables
    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_hyperbolic_tc) then
        ffhd_hyperbolic_tc=.false.
        if(mype==0) write(*,*) 'WARNING: set ffhd_hyperbolic_tc=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
      ! PI/LTE carry a thermodynamic state and so need the energy equation; only
      ! FI is meaningful without it (also guarded in ffhd_link_eos).
      if (eos%eos_type /= 'FI') &
        call mpistop("eos_type "//trim(eos%eos_type)//" requires ffhd_energy=T")
    end if
 
    if(ffhd_hyperbolic_tc) 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
    ! eos_type='PI' is the single partial-ionisation selector (no legacy flag).
 
    phys_gamma = eos%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))
    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_tc) then
      q_ = var_set_fluxvar('q', 'q', need_bc=.false.)
      need_global_cmax=.true.
    else
      q_=-1
    end if
 
    !> LTE stores Ne_/Te_ as advection-free extra vars (repopulated each step by
    !> eos%update_eos); PI keeps Te_ as an auxiliary var for the ionisation
    !> degree. Mirrors mhd: register here (before stop_indices) so FI leaves
    !> nwaux=0 and the var layout byte-identical.
    if (eos%eos_type == 'LTE') then
      ne_ = var_set_ne()
      te_ = var_set_te()
    else if (eos%eos_type == 'PI') then !  PI stores Te via var_set_te (sets iw_te) so the generic mod_eos_pi getters address it like LTE
      ne_ = -1
      te_ = var_set_te()
    else
      ne_ = -1
      te_ = -1
    end if
 
    ! set number of variables which need update ghostcells
    ! set number of variables which need update ghostcells.
    ! The EoS-derived slots Ne/Te are derived state that must stay consistent with (rho,e)
    ! wherever the conserved state is valid, so they have to travel with it. var_set_ne /
    ! var_set_te bump nw but neither nwflux nor nwaux, so the historic nwflux+nwaux silently
    ! excluded them: they were never communicated, only derived, and every ghost cell held
    ! zero until something derived it. Extend the window to whichever of them exist (they are
    ! registered contiguously just above); FI leaves both at -1 and the window is unchanged.
    ! bc_phys is unaffected -- it iterates nwflux+nwaux independently.
    nwgc=nwflux+nwaux
    if (iw_ne > 0) nwgc = max(nwgc, iw_ne)
    if (iw_te > 0) nwgc = max(nwgc, iw_te)
 
    ! set the index of the last flux variable for species 1
    stop_indices(1)=nwflux
 
    ! 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_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
    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
 
    ! Rfactor / temperature getter dispatch now lives in ffhd_link_eos (per eos_type)
 
    ! derive units from basic units
    call ffhd_physical_units()
 
    if(ffhd_hyperbolic_tc) then
      if(si_unit)then
        ! parallel conduction Spitzer
        hyperbolic_tc_kappa=8.d-12*unit_temperature**3.5d0/unit_length/unit_density/unit_velocity**3
      else
        ! in cgs
        hyperbolic_tc_kappa=8.d-7*unit_temperature**3.5d0/unit_length/unit_density/unit_velocity**3
      endif
    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_tc) 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
      call sts_init()
      call tc_init_params(eos%gamma)
 
      allocate(tc_fl)
      call tc_get_hd_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.)
      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%e_ = e_
      tc_fl%Tcoff_ = tcoff_
      ! get_temperature_*/get_rho/get_ne_nH/get_var_Rfactor/scalars wired in ffhd_bind_eos_to_source
    end if
 
    ! Initialize radiative cooling module
    if(ffhd_radiative_cooling) then
      call radiative_cooling_init_params(eos%gamma,eos%He_abundance)
      allocate(rc_fl)
      call radiative_cooling_init(rc_fl,rc_params_read)
      rc_fl%e_ = e_
      rc_fl%Tcoff_ = tcoff_
      rc_fl%subtract_equi = .false.
      ! get_rho/get_pthermal/get_var_Rfactor/get_Te/get_ne_nH + EoS scalars/kernels
      ! wired in ffhd_bind_eos_to_source
    end if
 
{^ifthreed
    allocate(te_fl_ffhd)
    phys_te_images => ffhd_te_images
    ! te_fl_ffhd%get_rho/get_pthermal/get_var_Rfactor/get_ne_nH wired in ffhd_bind_eos_to_source
}
 
    ! Initialize gravity module
    if(ffhd_gravity) then
      call gravity_init()
    end if
 
    ! The PI ionisation backend (both ionE modes) is initialised centrally in
    ! eos_finalise_PI (mod_eos_pi) -- mirrors hd/mhd, which no longer call
    ! ionization_degree_init from the phys module (one source of truth; calling
    ! it here too triggers "ionization_degree_init called more than once").
  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_mhd
    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_mhd(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_mhd
 
    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_mhd(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., tc_frozen_ion=.true.
    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, tc_frozen_ion
 
    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%frozen_ion = tc_frozen_ion
    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
  
    !> 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.
    logical    :: rad_damp=.false.
    double precision :: rad_damp_height=0.5d0
    double precision :: rad_damp_scale=0.15d0
 
    namelist /rc_list/ coolcurve, ncool, tlow, tfix, rc_split, rad_damp, rad_damp_height, rad_damp_scale
 
    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
    fl%rad_damp=rad_damp
    fl%rad_damp_height=rad_damp_height
    fl%rad_damp_scale=rad_damp_scale
  end subroutine rc_params_read
 
  subroutine ffhd_check_params
    use mod_global_parameters
    use mod_usr_methods
    use mod_convert, only: add_convert_method
    use mod_geometry, only: coordinate
 
 
    if (.not. ffhd_energy) then
       if (eos%gamma <= 0.0d0) call mpistop ("Error: eos%gamma <= 0")
       if (ffhd_adiab < 0.0d0) call mpistop ("Error: ffhd_adiab < 0")
       small_pressure = ffhd_adiab*small_density**eos%gamma
    else
       if (eos%gamma <= 0.0d0 .or. eos%gamma == 1.0d0) &
            call mpistop ("Error: eos%gamma <= 0 or eos%gamma == 1")
       small_e = small_pressure * eos%inv_gamma_minus_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
 
 
    if(mype==0)then
           write(*,*)'====FFHD run with settings===================='
           write(*,*)'Using mod_ffhd_phys with settings:'
           write(*,*)'SI_unit=',si_unit
           write(*,*)'Dimensionality   :',ndim
           write(*,*)'vector components:',ndir
           write(*,*)'coordinate set to type,slab:',coordinate,slab
           write(*,*)'number of variables          nw=',nw
           write(*,*)'    start index         iwstart=',iwstart
           write(*,*)'number of      vector variables=',nvector
           write(*,*)'number of stagger variables nws=',nws
           write(*,*)'number of    variables with BCs=',nwgc
           write(*,*)'number of      vars with fluxes=',nwflux
           write(*,*)'number of   vars with flux + BC=',nwfluxbc
           write(*,*)'number of   auxiliary variables=',nwaux
           write(*,*)'number of extra vars without flux=',nwextra
           write(*,*)'number of extra vars   for wextra=',nw_extra
           write(*,*)'number of auxiliary I/O variables=',nwauxio
           write(*,*)'    ffhd_energy=',ffhd_energy
           write(*,*)'    ffhd_gravity=',ffhd_gravity
           write(*,*)'    ffhd_radiative_cooling=',ffhd_radiative_cooling
           write(*,*)'    ffhd_hyperbolic_tc=',ffhd_hyperbolic_tc
           write(*,*)'    ffhd_trac=',ffhd_trac
           write(*,*)'number of             ghostcells=',nghostcells
           write(*,*)'number due to phys_wider_stencil=',phys_wider_stencil
           write(*,*)'==========================================='
    endif
  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
    
    ! Normalisation dispatch keyed solely on eos%eos_type (FI is the default, so
    ! legacy parfiles that set neither eos_type nor any flag land in the FI/PI
    ! absorbed-(a,b), RR=1 branch -- the historical eq_state_units=.true. result).
    if (eos%eos_type == 'LTE') then
      !> Remove the assumed FI normalisation from the units; handle in EoS.
      a=1d0
      b=1d0
      eos%nH2rhoFactor = 1d0+4d0*eos%He_abundance
      rr=(2d0+3d0*eos%He_abundance)/(1d0+4d0*eos%He_abundance)
    else
      !> FI / PI: absorbed-(a,b), RR=1. The (1+4He) factor is absorbed into
      !> unit_density, so nH2rhoFactor stays at the eos_init default of 1 (else
      !> get_ne_nH would double-divide rho and under-cool by (1+4He)^2).
      a=1d0+4d0*eos%He_abundance
      if(eos%eos_type=='PI') then
        b=1d0+h_ion_fr+eos%He_abundance*(he_ion_fr*(he_ion_fr2+1d0)+1d0)
      else
        b=2d0+3d0*eos%He_abundance
      end if
      rr=1d0
    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=dsqrt(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=dsqrt(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=dsqrt(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
  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)
 
    if(ffhd_energy) then
      w(ixo^s,e_)=w(ixo^s,p_)*eos%inv_gamma_minus_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)
 
    if(fix_small_values) then
      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_)=eos%gamma_minus_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(ixo^s,e_)=w(ixo^s,e_)+ffhd_kin_en(w,ixi^l,ixo^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)
 
    ! Restrict to ixO^L (mirrors hd/mhd): update_eos_LTE calls phys_e_to_ei
    ! over the mesh interior while ixI^S still spans unfilled corner ghosts
    ! (rho=0 -> SIGFPE in ffhd_kin_en). FI is unaffected (update_eos_FI no-op;
    ! the TC STS heads pass the range they need).
    w(ixo^s,e_)=w(ixo^s,e_)-ffhd_kin_en(w,ixi^l,ixo^l)
    if(fix_small_values) then
      call ffhd_handle_small_ei(w,x,ixi^l,ixo^l,e_,'ffhd_e_to_ei')
    end if
  end subroutine ffhd_e_to_ei
 
  !> Internal energy eint = E_total - E_kinetic (single field-aligned momentum).
  !> Wired to phys_get_ei; the LTE+ionE radiative cooling uses it to recover the
  !> gas internal energy from the conserved state.
  function ffhd_get_ei(w, ixI^L, ixO^L) result(ei)
    use mod_global_parameters
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, nw)
    double precision             :: ei(ixo^s)
 
    ei(ixo^s) = w(ixo^s,e_) - ffhd_kin_en(w,ixi^l,ixo^l)
  end function ffhd_get_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_)=eos%gamma_minus_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(wprim,x,ixI^L,ixO^L,idim,cmax)
    use mod_global_parameters
    integer, intent(in)          :: ixi^l, ixo^l, idim
    ! w in primitive form
    double precision, intent(in) :: wprim(ixi^s, nw), x(ixi^s,1:ndim)
    double precision, intent(inout) :: cmax(ixi^s)
 
    if(ffhd_energy) then
      if (eos%ionE) then
        ! LTE/PI energy: Gamma_1(p,rho)*p/rho from the EoS (wprim is primitive)
        call eos%get_csound2(wprim,x,ixi^l,ixo^l,cmax)
        cmax(ixo^s)=dsqrt(cmax(ixo^s))
      else
        cmax(ixo^s)=dsqrt(eos%gamma*wprim(ixo^s,p_)/wprim(ixo^s,rho_))
      end if
    else
      cmax(ixo^s)=dsqrt(eos%gamma*ffhd_adiab*wprim(ixo^s,rho_)**eos%gamma_minus_1)
    end if
    cmax(ixo^s)=dabs(wprim(ixo^s,mom(1))*block%B0(ixo^s,idim,0))+cmax(ixo^s)
 
  end subroutine ffhd_get_cmax_origin
 
  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 :: te(ixi^s),lts(ixi^s)
    double precision, dimension(ixI^S,1:ndim) :: gradt
    double precision :: bdir(ndim)
    double precision :: ltrc,ltrp,altr
    integer :: idims,ix^d,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*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=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*dabs(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,gradt(ixi^s,idims))
      end do
      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)+block%B0(ixb^d,1:ndim,0)
        {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
     {do ix^db=ixomin^db,ixomax^db\}
        {^iftwod
        ! temperature length scale inversed
        lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2))*&
          abs(^d&gradt({ix^d},^d)*block%B0({ix^d},^d,0)+)/te(ix^d)
        }
        {^ifthreed
        ! temperature length scale inversed
        lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2),block%ds(ix^d,3))*&
          abs(^d&gradt({ix^d},^d)*block%B0({ix^d},^d,0)+)/te(ix^d)
        }
        if(lts(ix^d)>trac_delta) then
          block%special_values(1)=max(block%special_values(1),te(ix^d))
        end if
     {end do\}
      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 gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
        call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
      end do
     {do ix^db=ixpmin^db,ixpmax^db\}
        ! temperature length scale inversed
        lts(ix^d)=abs(^d&gradt({ix^d},^d)*block%B0({ix^d},^d,0)+)/te(ix^d)
        ! fraction of cells size to temperature length scale
        lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
        lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
     {end do\}
  
      ! need one ghost layer for thermal conductivity
      ixp^l=ixo^l^ladd1;
     {do ix^db=ixpmin^db,ixpmax^db\}
        {^iftwod
        altr=0.25d0*((lts(ix1-1,ix2)+two*lts(ix^d)+lts(ix1+1,ix2))*block%B0(ix^d,1,0)**2+&
                     (lts(ix1,ix2-1)+two*lts(ix^d)+lts(ix1,ix2+1))*block%B0(ix^d,2,0)**2)
        block%wextra(ix^d,tcoff_)=te(ix^d)*altr**0.4d0
        }
        {^ifthreed
        altr=0.25d0*((lts(ix1-1,ix2,ix3)+two*lts(ix^d)+lts(ix1+1,ix2,ix3))*block%B0(ix^d,1,0)**2+&
                     (lts(ix1,ix2-1,ix3)+two*lts(ix^d)+lts(ix1,ix2+1,ix3))*block%B0(ix^d,2,0)**2+&
                     (lts(ix1,ix2,ix3-1)+two*lts(ix^d)+lts(ix1,ix2,ix3+1))*block%B0(ix^d,3,0)**2)
        block%wextra(ix^d,tcoff_)=te(ix^d)*altr**0.4d0
        }
     {end do\}
    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), wmeanp(ixi^s,nw)
    double precision, dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
 
    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(1))*tmp1(ixo^s)&
                   +wrp(ixo^s,mom(1))*tmp2(ixo^s))*tmp3(ixo^s)
      umean(ixo^s)=umean(ixo^s)*block%B0(ixo^s,idim,idim)
      call ffhd_csound2_cbounds(wlc,wlp,x,ixi^l,ixo^l,csoundl)
      call ffhd_csound2_cbounds(wrc,wrp,x,ixi^l,ixo^l,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))-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)=dabs(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_)
      if (eos%ionE) then
        wmeanp(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
        call eos%get_csound2(wmeanp,x,ixi^l,ixo^l,csoundr)
      else
        call ffhd_get_csound2(wmean,x,ixi^l,ixo^l,csoundr)
      end if
      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)
      else
        cmax(ixo^s,1)=dabs(tmp1(ixo^s))+csoundr(ixo^s)
      end if
    case (3)
      ! Miyoshi 2005 JCP 208, 315 equation (67)
      call ffhd_csound2_cbounds(wlc,wlp,x,ixi^l,ixo^l,csoundl)
      call ffhd_csound2_cbounds(wrc,wrp,x,ixi^l,ixo^l,csoundr)
      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(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
 
  !> Sound speed^2 at a Riemann interface. FI uses the conserved-state pthermal
  !> path (bit-identical to the legacy ffhd_get_csound2(wC) call); LTE uses the
  !> primitive Gamma_1*p/rho via eos% because the stored Te_/Ne_ are NOT
  !> reconstructed onto interface states (they are advection-free extra vars).
  subroutine ffhd_csound2_cbounds(wC,wp,x,ixI^L,ixO^L,cs2)
    use mod_global_parameters
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: wc(ixi^s,nw), wp(ixi^s,nw), x(ixi^s,1:ndim)
    double precision, intent(out):: cs2(ixi^s)
 
    if (eos%ionE) then
      call eos%get_csound2(wp,x,ixi^l,ixo^l,cs2)
    else
      call ffhd_get_csound2(wc,x,ixi^l,ixo^l,cs2)
    end if
  end subroutine ffhd_csound2_cbounds
 
  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)**eos%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)=eos%gamma_minus_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(*,*) dsqrt(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) = eos%gamma_minus_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)=eos%gamma*csound2(ixo^s)/rho(ixo^s)
    else
      csound2(ixo^s)=eos%gamma*ffhd_adiab*rho(ixo^s)**eos%gamma_minus_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)
 
    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_)**eos%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_tc) 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_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_tc) then
        call add_hyperbolic_tc_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_gravity) then
      call gravity_add_source(qdt,ixi^l,ixo^l,wct,wctprim,&
           w,x,ffhd_energy,qsourcesplit,active)
    end if
 
    ! update temperature from new pressure, density, and old ionization degree
    if(eos%eos_type == 'PI') then
      if(.not.qsourcesplit) then
        active = .true.
        call eos%update_eos(ixi^l,ixo^l,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)
    double precision                :: rhovpar(ixi^s),gradrhov(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_)*eos%gamma_minus_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_get_dt(wprim,ixI^L,ixO^L,dtnew,dx^D,x)
    use mod_global_parameters
    use mod_usr_methods
    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)    :: wprim(ixi^s,1:nw)
    double precision, intent(in)    :: x(ixi^s,1:ndim)
 
    dtnew = bigdouble
 
    if(ffhd_gravity) then
      call gravity_get_dt(wprim,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,wprim,w,x)
    use mod_global_parameters
    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)
 
    ! no geometric source terms needed for ffhd, no divergences of tensors
  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_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 add_hyperbolic_tc_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
 
    double precision, dimension(ixI^S) :: te, r
    double precision :: sigma_t5,sigma_t7,sigmat5_bgradt,f_sat,tau
    integer :: ix^d
 
    ! EoS-aware temperature: T = p/(R*rho). For FI R=RR (=1 under eq_state_units)
    ! recovers p/rho; for LTE R=Rfactor_from_LTE varies, so LTE+HTC is physically
    ! consistent. Computed over ixI^L (the +-2 stencil reaches ghost layers;
    ! their R uses the update_eos_4_bc-refreshed aux state, as in mhd).
    call ffhd_get_rfactor(wct,x,ixi^l,ixi^l,r)
    te(ixi^s)=wctprim(ixi^s,p_)/(r(ixi^s)*wct(ixi^s,rho_))
    ! temperature on face T_(i+1/2)=(7(T_i+T_(i+1))-(T_(i-1)+T_(i+2)))/12
    ! T_(i+1/2)-T_(i-1/2)=(8(T_(i+1)-T_(i-1))-T_(i+2)+T_(i-2))/12
   {^iftwod
    do ix2=ixomin2,ixomax2
      do ix1=ixomin1,ixomax1
        if(ffhd_trac) then
          if(te(ix^d)<block%wextra(ix^d,tcoff_)) then
            sigma_t5=hyperbolic_tc_kappa*sqrt(block%wextra(ix^d,tcoff_)**5)
            sigma_t7=sigma_t5*block%wextra(ix^d,tcoff_)
          else
            sigma_t5=hyperbolic_tc_kappa*sqrt(te(ix^d)**5)
            sigma_t7=sigma_t5*te(ix^d)
          end if
        else
          sigma_t5=hyperbolic_tc_kappa*sqrt(te(ix^d)**5)
          sigma_t7=sigma_t5*te(ix^d)
        end if
        sigmat5_bgradt=sigma_t5*(&
           block%B0(ix^d,1,0)*((8.d0*(te(ix1+1,ix2)-te(ix1-1,ix2))-te(ix1+2,ix2)+te(ix1-2,ix2))/12.d0)/block%ds(ix^d,1)&
          +block%B0(ix^d,2,0)*((8.d0*(te(ix1,ix2+1)-te(ix1,ix2-1))-te(ix1,ix2+2)+te(ix1,ix2-2))/12.d0)/block%ds(ix^d,2))
        if(ffhd_hyperbolic_tc_sat) then
          ! 5 phi rho c^3, phi=0.3, c=sqrt(p/rho) isothermal sound speed
          f_sat=one/(one+dabs(sigmat5_bgradt)/(1.5d0*wct(ix^d,rho_)*(wctprim(ix^d,p_)/wct(ix^d,rho_))**1.5d0))
          tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(wctprim(ix^d,p_)*eos%inv_gamma_minus_1*cmax_global**2))
          w(ix^d,q_)=w(ix^d,q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,q_))/tau
        else
          w(ix^d,q_)=w(ix^d,q_)-qdt*(sigmat5_bgradt+wct(ix^d,q_))/&
           max(4.d0*dt, sigma_t7*courantpar**2/(wctprim(ix^d,p_)*eos%inv_gamma_minus_1*cmax_global**2))
        end if
      end do
    end do
   }
   {^ifthreed
    do ix3=ixomin3,ixomax3
      do ix2=ixomin2,ixomax2
        do ix1=ixomin1,ixomax1
          if(ffhd_trac) then
            if(te(ix^d)<block%wextra(ix^d,tcoff_)) then
              sigma_t5=hyperbolic_tc_kappa*sqrt(block%wextra(ix^d,tcoff_)**5)
              sigma_t7=sigma_t5*block%wextra(ix^d,tcoff_)
            else
              sigma_t5=hyperbolic_tc_kappa*sqrt(te(ix^d)**5)
              sigma_t7=sigma_t5*te(ix^d)
            end if
          else
            sigma_t5=hyperbolic_tc_kappa*sqrt(te(ix^d)**5)
            sigma_t7=sigma_t5*te(ix^d)
          end if
          sigmat5_bgradt=sigma_t5*(&
             block%B0(ix^d,1,0)*((8.d0*(te(ix1+1,ix2,ix3)-te(ix1-1,ix2,ix3))-te(ix1+2,ix2,ix3)+te(ix1-2,ix2,ix3))/12.d0)/block%ds(ix^d,1)&
            +block%B0(ix^d,2,0)*((8.d0*(te(ix1,ix2+1,ix3)-te(ix1,ix2-1,ix3))-te(ix1,ix2+2,ix3)+te(ix1,ix2-2,ix3))/12.d0)/block%ds(ix^d,2)&
            +block%B0(ix^d,3,0)*((8.d0*(te(ix1,ix2,ix3+1)-te(ix1,ix2,ix3-1))-te(ix1,ix2,ix3+2)+te(ix1,ix2,ix3-2))/12.d0)/block%ds(ix^d,3))
          if(ffhd_hyperbolic_tc_sat) then
            ! 5 phi rho c^3, phi=0.3, c=sqrt(p/rho) isothermal sound speed
            f_sat=one/(one+dabs(sigmat5_bgradt)/(1.5d0*wct(ix^d,rho_)*(wctprim(ix^d,p_)/wct(ix^d,rho_))**1.5d0))
            tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(wctprim(ix^d,p_)*eos%inv_gamma_minus_1*cmax_global**2))
            w(ix^d,q_)=w(ix^d,q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,q_))/tau
          else
            w(ix^d,q_)=w(ix^d,q_)-qdt*(sigmat5_bgradt+wct(ix^d,q_))/&
             max(4.d0*dt, sigma_t7*courantpar**2/(wctprim(ix^d,p_)*eos%inv_gamma_minus_1*cmax_global**2))
          end if
        end do
      end do
    end do
    }
  end subroutine add_hyperbolic_tc_source
 
end module mod_ffhd_phys