mod_thermal_conduction.t

Go to the documentation of this file.

!> Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module
!> Adaptation of mod_thermal_conduction for the mod_supertimestepping
!>
!> The TC is set by calling 
!> tc_init_params()
!>
!> Organized such that it can call either isotropic (HD) or anisotropic (MHD) variants
!> it adds a heat conduction source to each energy equation
!> and can be recycled within a multi-fluid context (such as plasma-neutral twofl module)
!>
!>
!> 10.07.2011 developed by Chun Xia and Rony Keppens
!> 01.09.2012 moved to modules folder by Oliver Porth
!> 13.10.2013 optimized further by Chun Xia
!> 12.03.2014 implemented RKL2 super timestepping scheme to reduce iterations
!> and improve stability and accuracy up to second order in time by Chun Xia.
!> 23.08.2014 implemented saturation and perpendicular TC by Chun Xia
!> 12.01.2017 modulized by Chun Xia
!>            adapted by Beatrice Popescu to twofluid settings
!> 06.09.2024 cleaned up for use in rhd and rmhd modules (Nishant Narechania and Rony Keppens)
!> 30.11.2025 Minor cleanup (for consistency between hd and mhd)
!>
!> PURPOSE:
!> IN MHD ADD THE HEAT CONDUCTION SOURCE TO THE ENERGY EQUATION
!> S=DIV(KAPPA_i,j . GRAD_j T)
!> where KAPPA_i,j = tc_k_para b_i b_j + tc_k_perp (I - b_i b_j)
!> b_i b_j = B_i B_j / B**2, I is the unit matrix, and i, j= 1, 2, 3 for 3D
!> IN HD ADD THE HEAT CONDUCTION SOURCE TO THE ENERGY EQUATION
!> S=DIV(tc_k_para . GRAD T)
!> USAGE:
!> 1. in mod_usr.t -> subroutine usr_init(), add
!>        unit_length=your length unit
!>        unit_numberdensity=your number density unit
!>        unit_velocity=your velocity unit
!>        unit_temperature=your temperature unit
!>    before call (m)hd_activate()
!> 2. to switch on thermal conduction in the (r)(m)hd_list of amrvac.par add:
!>    (r)(m)hd_thermal_conduction=.true.
!> 3. in the tc_list of amrvac.par :
!>    tc_perpendicular=.true.  ! (default .false.) turn on thermal conduction perpendicular to magnetic field
!>    tc_saturate=.true.  ! (default .false. ) turn on thermal conduction saturate effect
!>    tc_slope_limiter='MC' ! choose limiter for slope-limited anisotropic thermal conduction in MHD
!> note: twofl_list incorporates instances for charges and neutrals
 
module mod_thermal_conduction
  use mod_global_parameters, only: std_len
  use mod_geometry
  use mod_comm_lib, only: mpistop
  implicit none
 
    !> The adiabatic index
    double precision :: tc_gamma_1
 
  abstract interface
    subroutine get_var_subr(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,nw)
      double precision, intent(in) :: x(ixI^S,1:ndim)
      double precision, intent(out):: res(ixI^S)
    end subroutine get_var_subr
 
    !> temperature from the energy variable and, if present, ciso2 = p/rho
    subroutine get_tciso2_subr(w,x,ixI^L,ixO^L,res,ciso2)
      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):: res(ixI^S)
      double precision, intent(out), optional :: ciso2(ixI^S)
    end subroutine get_tciso2_subr
 
    subroutine get_2var_subr(ixI^L, ixO^L, w, x, ne, nH)
      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):: ne(ixI^S), nH(ixI^S)
    end subroutine get_2var_subr
 
    !> Scalar EoS inverse, e.g. eint_nH_from_T(log_nH, log_T)
    double precision function eos_scalar2_func(a, b)
      double precision, intent(in) :: a, b
    end function eos_scalar2_func
  end interface
 
  type tc_fluid
 
    ! BEGIN the following are read from param file or set in tc_read_hd_params or tc_read_mhd_params
    !> Coefficient of thermal conductivity (parallel to magnetic field)
    double precision :: tc_k_para
 
    !> Coefficient of thermal conductivity perpendicular to magnetic field
    double precision :: tc_k_perp
 
     !> Indices of the variables
    integer :: e_=-1
    !> Index of cut off temperature for TRAC
    integer :: tcoff_
    !> Name of slope limiter for transverse component of thermal flux
    integer :: tc_slope_limiter
 
    ! if subtract_equi = .true. get_temperature_equi and get_rho_equi have to be set
    logical :: subtract_equi=.false.
 
    !> Logical switch for test constant conductivity
    logical :: tc_constant=.false.
 
    !> Calculate thermal conduction perpendicular to magnetic field (.true.) or not (.false.)
    logical :: tc_perpendicular=.false.
 
    !> Consider thermal conduction saturation effect (.true.) or not (.false.)
    logical :: tc_saturate=.false.
 
    !> Patch cells where e_int < 0 during STS Chebyshev substeps by
    !> neighbor-averaging the temperature. Prevents NaN from sqrt(T<0)
    !> in the Spitzer conductivity when the RKL2 polynomial overshoots.
    logical :: tc_patch_eint=.false.
 
    !> LTE with ionisation energy: hold the ionisation state fixed through each conduction
    !> step. T and R = p/(rho T) come from the full EoS at the first STS stage; the later
    !> stages map the internal energy to T with the frozen-composition heat capacity
    logical :: frozen_ion=.true.
 
    !> true when R = p/(rho T) is not 1, so saturation takes ciso2 = R T from get_var_Rfactor
    logical :: sat_r_variable=.false.
 
    !> Minimum temperature (code units) below which TRAC does not modify conductivity.
    !> Below this T, the energy balance is dominated by optically thick radiation
    !> and recombination, not optically thin cooling + Spitzer conduction, so TRAC's
    !> broadening assumption does not apply. Read in Kelvin via tc_list (trac_T_floor),
    !> converted to code units at init. Default 0 (no floor). Only the HD TRAC cutoff (types 0
    !> and 7) reads it, with Tcoff = Te below it, so conduction, losses and heating are left
    !> unmodified there
    double precision :: trac_t_floor=0.d0
    ! END the following are read from param file or set in tc_read_hd_params or tc_read_mhd_params
    procedure(get_var_subr), pointer, nopass :: get_rho => null()
    procedure(get_var_subr), pointer, nopass :: get_rho_equi => null()
    procedure(get_var_subr), pointer,nopass :: get_temperature_from_eint => null()
    !> temperature getter that also returns ciso2 = p/rho, bound instead of
    !> get_temperature_from_eint for LTE with ionisation energy
    procedure(get_tciso2_subr), pointer,nopass :: get_temperature_ciso2_from_eint => null()
    procedure(get_var_subr), pointer,nopass :: get_temperature_from_conserved => null()
    procedure(get_var_subr), pointer,nopass :: get_temperature_equi => null()
    procedure(get_2var_subr), pointer,nopass :: get_ne_nh => null()
    procedure(get_var_subr), pointer,nopass :: get_var_rfactor => null()
    !> EoS snapshots + inverse accessor (set in bind_eos_to_source); let TC reach
    !> thermodynamics only through this object, never mod_eos directly.
    double precision :: inv_gamma_minus_1
    double precision :: nh2rhofactor
    double precision :: log_t_floor
    procedure(eos_scalar2_func), pointer, nopass :: eint_from_t => null()
  end type tc_fluid
 
  public :: tc_get_mhd_params
  public :: tc_get_hd_params
  public :: get_tc_dt_mhd
  public :: get_tc_dt_hd
  public :: sts_set_source_tc_mhd
  public :: sts_set_source_tc_hd
  public :: tc_patch_negative_eint
 
  !> Per-block anchor of the frozen-ionisation conduction step: T, R and the internal
  !> energy at the step start, and the inverse frozen-composition heat capacity per volume
  type tc_anchor_t
    integer :: sts_id = -1               !< STS application the anchor belongs to
    double precision, allocatable :: t0(:^d&), e0(:^D&), R0(:^D&), cvinv(:^D&)
  end type tc_anchor_t
  type(tc_anchor_t), allocatable :: tc_anchor(:)
 
contains
 
  subroutine tc_init_params(phys_gamma)
    use mod_global_parameters
    double precision, intent(in) :: phys_gamma
 
    tc_gamma_1=phys_gamma-1d0
  end subroutine tc_init_params
 
  !> Init TC coefficients: MHD case
  subroutine tc_get_mhd_params(fl,read_mhd_params)
    use mod_global_parameters
 
    interface
      subroutine read_mhd_params(fl)
        use mod_global_parameters, only: unitpar,par_files
        import tc_fluid
        type(tc_fluid), intent(inout) :: fl
 
      end subroutine read_mhd_params
    end interface
    type(tc_fluid), intent(inout) :: fl
 
    fl%tc_slope_limiter=1
    fl%tc_k_para=0.d0
    fl%tc_k_perp=0.d0
 
    !> Read tc module parameters from par file: MHD case
    call read_mhd_params(fl)
 
    if(fl%tc_k_para==0.d0 .and. fl%tc_k_perp==0.d0) then
      if(si_unit) then
        ! Spitzer thermal conductivity with SI units
        fl%tc_k_para=8.d-12*unit_temperature**3.5d0/unit_length/unit_density/unit_velocity**3
        ! thermal conductivity perpendicular to magnetic field
        fl%tc_k_perp=4.d-30*unit_numberdensity**2/unit_magneticfield**2/unit_temperature**3*fl%tc_k_para
      else
        ! Spitzer thermal conductivity with cgs units
        fl%tc_k_para=8.d-7*unit_temperature**3.5d0/unit_length/unit_density/unit_velocity**3
        ! thermal conductivity perpendicular to magnetic field
        fl%tc_k_perp=4.d-10*unit_numberdensity**2/unit_magneticfield**2/unit_temperature**3*fl%tc_k_para
      end if
      if(mype .eq. 0) print*, "Spitzer MHD: par: ",fl%tc_k_para, &
          " ,perp: ",fl%tc_k_perp
    else
      fl%tc_constant=.true.
      if(mype .eq. 0) print*, "Constant thermal conduction coefficients with values: ",fl%tc_k_para,fl%tc_k_perp
    end if
  end subroutine tc_get_mhd_params
 
  !> Init  TC coefficients: HD case
  subroutine tc_get_hd_params(fl,read_hd_params)
    use mod_global_parameters
 
    interface
      subroutine read_hd_params(fl)
        use mod_global_parameters, only: unitpar,par_files
        import tc_fluid
        type(tc_fluid), intent(inout) :: fl
 
      end subroutine read_hd_params
    end interface
    type(tc_fluid), intent(inout) :: fl
 
    fl%tc_k_para=0.d0
 
    !> Read tc parameters from par file: HD case
    call read_hd_params(fl)
 
    if(fl%tc_k_para==0.d0) then
      if(si_unit) then
        ! Spitzer thermal conductivity with SI units
        fl%tc_k_para=8.d-12*unit_temperature**3.5d0/unit_length/unit_density/unit_velocity**3
      else
        ! Spitzer thermal conductivity with cgs units
        fl%tc_k_para=8.d-7*unit_temperature**3.5d0/unit_length/unit_density/unit_velocity**3
      end if
      if(mype .eq. 0) print*, "Spitzer HD par: ",fl%tc_k_para
    else
      fl%tc_constant=.true.
      if(mype .eq. 0) print*, "Constant thermal conduction coefficient at value: ",fl%tc_k_para
    end if
 
  end subroutine tc_get_hd_params
 
  !> Get the explicit timestep for the TC (mhd implementation)
  !> Note: for multi-D MHD (1D MHD will use HD fall-back)
  function get_tc_dt_mhd(w,ixI^L,ixO^L,dx^D,x,fl) result(dtnew)
    !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/(rho R))
    !where                      tc_k_para_i=tc_k_para*B_i**2/B**2
    !and                        T=p/rho
    use mod_global_parameters
 
    type(tc_fluid), intent(in)  ::  fl
    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
 
    double precision :: mf(ixo^s,1:ndir),te(ixi^s),rho(ixi^s),gradt(ixi^s)
    double precision :: ne(ixi^s), nh_arr(ixi^s), rf(ixi^s)
    double precision :: tmp(ixo^s),hfs(ixo^s),blocal(1:ndir),bmag
    double precision :: dtdiff_tcond,maxtmp2
    integer          :: idims,ix^d
 
 
    ! B
    if(allocated(iw_mag)) then
      if(b0field) then
       {do ix^db=ixomin^db,ixomax^db\}
          ^c&blocal(^c)=w({ix^d},iw_mag(^c))+block%B0({ix^d},^c,0)\
          bmag=dsqrt(^c&blocal(^c)**2+)+smalldouble
          ^c&mf(ix^d,^c)=blocal(^c)/bmag\
       {end do\}
      else
       {do ix^db=ixomin^db,ixomax^db\}
          bmag=dsqrt(^c&w(ix^d,iw_mag(^c))**2+)+smalldouble
          ^c&mf(ix^d,^c)=w(ix^d,iw_mag(^c))/bmag\
       {end do\}
      end if
    else
      ! this if for ffhd or other physics modules without B components 
      mf(ixo^s,1:ndim)=block%B0(ixo^s,1:ndim,0)
    end if
 
    !temperature
    call fl%get_temperature_from_conserved(w,x,ixi^l,ixi^l,te)
    call fl%get_rho(w,x,ixi^l,ixo^l,rho)
    call fl%get_ne_nH(ixi^l, ixo^l, w, x, ne, nh_arr)
    ! the diffusivity is kappa/(rho c_v) with c_v = R/(gamma-1): weight rho by R
    if(associated(fl%get_var_Rfactor)) then
      call fl%get_var_Rfactor(w, x, ixi^l, ixo^l, rf)
      rho(ixo^s)=rho(ixo^s)*rf(ixo^s)
    end if
 
    !tc_k_para_i
    if(fl%tc_constant) then
      tmp(ixo^s)=fl%tc_k_para
    else
      if(fl%tc_saturate) then
        ! Kannan 2016 MN 458, 410
        ! l_mfpe = 3^1.5*kB^2/(4*sqrt(pi)*e^4*lnLambda) * T^2/n_e
        if(si_unit) then
          ! 5.730205638843984e27 = 3^1.5*kB_SI^2/(4*sqrt(pi)*e_SI^4*37)
          tmp(ixo^s)=te(ixo^s)**2/ne(ixo^s)*5.730205638843984d27*unit_temperature**2/(unit_numberdensity*unit_length)
        else
          ! 7093.9239487765044 = 3^1.5*kB_cgs^2/(4*sqrt(pi)*e_cgs^4*37)
          tmp(ixo^s)=te(ixo^s)**2/ne(ixo^s)*7093.9239487765044d0*unit_temperature**2/(unit_numberdensity*unit_length)
        end if
        do idims=1,ndim
          call gradient(te,ixi^l,ixo^l,idims,gradt)
          if(idims==1) then
            hfs(ixo^s)=gradt(ixo^s)*mf(ixo^s,idims)
          else
            hfs(ixo^s)=hfs(ixo^s)+gradt(ixo^s)*mf(ixo^s,idims)
          end if
        end do
        ! kappa=kappa_Spitzer/(1+4.2*l_mfpe/(T/|gradT.b|))
        tmp(ixo^s)=fl%tc_k_para*te(ixo^s)*te(ixo^s)*dsqrt(te(ixo^s))/(1.d0+4.2d0*tmp(ixo^s)*dabs(hfs(ixo^s))/te(ixo^s))
      else
        ! kappa=kappa_Spitzer
        tmp(ixo^s)=fl%tc_k_para*te(ixo^s)*te(ixo^s)*dsqrt(te(ixo^s))
      end if
    end if
 
    dtnew=bigdouble
    do idims=1,ndim
      ! approximate thermal conduction flux: tc_k_para_i/rho/dx*B_i**2/B**2
      maxtmp2=maxval(tmp(ixo^s)*mf(ixo^s,idims)**2/(rho(ixo^s)*block%ds(ixo^s,idims)**2))
      ! dt< dx_idim**2/((gamma-1)*tc_k_para_i/rho*B_i**2/B**2)
      dtdiff_tcond=1.d0/(tc_gamma_1*maxtmp2+smalldouble)
      ! limit the time step
      dtnew=min(dtnew,dtdiff_tcond)
    end do
    dtnew=dtnew/dble(ndim)
  end function get_tc_dt_mhd
 
  !> anisotropic thermal conduction with slope limited symmetric scheme
  !> Sharma 2007 Journal of Computational Physics 227, 123
  subroutine sts_set_source_tc_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,fl)
    use mod_global_parameters
    use mod_fix_conserve
    integer, intent(in) :: ixi^l, ixo^l, igrid, nflux
    double precision, intent(in) ::  x(ixi^s,1:ndim)
    ! intent(inout) so tc_patch_negative_eint can repair w(:, ie) in place
    ! when tc_patch_eint = .true.; legacy path (patch off) leaves w untouched.
    double precision, intent(inout) ::   w(ixi^s,1:nw)
    double precision, intent(inout) ::  wres(ixi^s,1:nw)
    double precision, intent(in) :: my_dt
    logical, intent(in) :: fix_conserve_at_step
    type(tc_fluid), intent(in) :: fl
 
    !! qd store the heat conduction energy changing rate
    double precision :: qd(ixo^s)
    double precision :: rho(ixi^s),te(ixi^s),ciso2(ixi^s)
    double precision :: qvec(ixi^s,1:ndim)
    double precision :: fluxall(ixi^s,1,1:ndim)
    double precision :: alpha,dxinv(ndim)
    double precision, allocatable, dimension(:^D&,:) :: qvec_equi
    integer :: idims,ixa^l
 
    ! coefficient of limiting on normal component
    if(ndim<3) then
      alpha=0.75d0
    else
      alpha=0.85d0
    end if
 
    dxinv=1.d0/dxlevel
 
    call tc_get_te_ciso2(ixi^l,w,x,fl,igrid,te,ciso2)
    call fl%get_rho(w, x, ixi^l, ixi^l, rho)
    call tc_flux_mhd(ixi^l,ixo^l,w,x,fl,qvec,rho,te,ciso2,alpha)
    if(fl%subtract_equi) then
      allocate(qvec_equi(ixi^s,1:ndim))
      call tc_equi_state(ixi^l, w, x, fl, te, rho, ciso2)
      call tc_flux_mhd(ixi^l,ixo^l,w,x,fl,qvec_equi,rho,te,ciso2,alpha)
      do idims=1,ndim
        ixamax^d=ixomax^d; ixamin^d=ixomin^d-kr(idims,^d);
        qvec(ixa^s,idims)=qvec(ixa^s,idims)-qvec_equi(ixa^s,idims)
      end do
      deallocate(qvec_equi)
    end if
 
    if(slab_uniform) then
      do idims=1,ndim
        ixamax^d=ixomax^d; ixamin^d=ixomin^d-kr(idims,^d);
        qvec(ixa^s,idims)=dxinv(idims)*qvec(ixa^s,idims)
        ixa^l=ixo^l-kr(idims,^d);
        if(idims==1) then
          qd(ixo^s)=qvec(ixo^s,idims)-qvec(ixa^s,idims)
        else
          qd(ixo^s)=qd(ixo^s)+qvec(ixo^s,idims)-qvec(ixa^s,idims)
        end if
      end do
    else
      do idims=1,ndim
        ixamax^d=ixomax^d; ixamin^d=ixomin^d-kr(idims,^d);
        qvec(ixa^s,idims)=qvec(ixa^s,idims)*block%surfaceC(ixa^s,idims)
        ixa^l=ixo^l-kr(idims,^d);
        if(idims==1) then
          qd(ixo^s)=qvec(ixo^s,idims)-qvec(ixa^s,idims)
        else
          qd(ixo^s)=qd(ixo^s)+qvec(ixo^s,idims)-qvec(ixa^s,idims)
        end if
      end do
      qd(ixo^s)=qd(ixo^s)/block%dvolume(ixo^s)
    end if
 
    if(fix_conserve_at_step) then
      do idims=1,ndim
        ixamax^d=ixomax^d; ixamin^d=ixomin^d-kr(idims,^d);
        fluxall(ixa^s,1,idims)=my_dt*qvec(ixa^s,idims)
      end do
      call store_flux(igrid,fluxall,1,ndim,nflux)
    end if
 
    wres(ixo^s,fl%e_)=qd(ixo^s)
  end subroutine sts_set_source_tc_mhd
 
  subroutine set_source_tc_mhd(ixI^L,ixO^L,w,x,fl,qvec,rho,Te,ciso2,alpha)
    use mod_global_parameters
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) ::  x(ixI^S,1:ndim)
    double precision, intent(in) ::  w(ixI^S,1:nw)
    type(tc_fluid), intent(in) :: fl
    double precision, intent(in) :: rho(ixI^S),Te(ixI^S),ciso2(ixI^S)
    double precision, intent(in) :: alpha
    double precision, intent(out) :: qvec(ixI^S,1:ndim)
 
    !! qdd store the heat conduction energy changing rate
    double precision, dimension(ixI^S,1:ndim) :: mf,Bc,Bcf,gradT
    double precision, dimension(ixI^S) :: ka,kaf,ke,kef,qdd,Bnorm
    double precision :: minq,maxq,qd(ixI^S,2**(ndim-1)), blocal(ndir)
    integer :: idims,idir,ix^D,ix^L,ixC^L,ixA^L,ixB^L
 
    ix^l=ixo^l^ladd1;
 
    ! T gradient at cell faces
    ! b unit vector mf: magnetic field direction vector
    if(allocated(iw_mag)) then
      if(b0field) then
       {do ix^db=ixmin^db,ixmax^db\}
          ^c&blocal(^c)=w({ix^d},iw_mag(^c))+block%B0({ix^d},^c,0)\
         {^iftwod
          if(blocal(1)/=0.d0) then
            mf(ix^d,1)=sign(1.d0,blocal(1))/dsqrt(1.d0+(^ce&(blocal(^ce)/blocal(1))**2+))
          else
            mf(ix^d,1)=0.d0
          end if
          if(blocal(2)/=0.d0) then
            mf(ix^d,2)=sign(1.d0,blocal(2))/dsqrt(1.d0+(^cf&(blocal(^cf)/blocal(2))**2+))
          else
            mf(ix^d,2)=0.d0
          end if
         }
         {^ifthreed
          if(blocal(1)/=0.d0) then
            mf(ix^d,1)=sign(1.d0,blocal(1))/dsqrt(1.d0+(blocal(2)/blocal(1))**2+(blocal(3)/blocal(1))**2)
          else
            mf(ix^d,1)=0.d0
          end if
          if(blocal(2)/=0.d0) then
            mf(ix^d,2)=sign(1.d0,blocal(2))/dsqrt(1.d0+(blocal(1)/blocal(2))**2+(blocal(3)/blocal(2))**2)
          else
            mf(ix^d,2)=0.d0
          end if
          if(blocal(3)/=0.d0) then
            mf(ix^d,3)=sign(1.d0,blocal(3))/dsqrt(1.d0+(blocal(1)/blocal(3))**2+(blocal(2)/blocal(3))**2)
          else
            mf(ix^d,3)=0.d0
          end if
         }
       {end do\}
      else
       {do ix^db=ixmin^db,ixmax^db\}
         {^iftwod
          if(w(ix^d,iw_mag(1))/=0.d0) then
            mf(ix^d,1)=sign(1.d0,w(ix^d,iw_mag(1)))/dsqrt(1.d0+(^ce&(w(ix^d,iw_mag(^ce))/w(ix^d,iw_mag(1)))**2+))
          else
            mf(ix^d,1)=0.d0
          end if
          if(w(ix^d,iw_mag(2))/=0.d0) then
            mf(ix^d,2)=sign(1.d0,w(ix^d,iw_mag(2)))/dsqrt(1.d0+(^cf&(w(ix^d,iw_mag(^cf))/w(ix^d,iw_mag(2)))**2+))
          else
            mf(ix^d,2)=0.d0
          end if
         }
         {^ifthreed
          if(w(ix^d,iw_mag(1))/=0.d0) then
            mf(ix^d,1)=sign(1.d0,w(ix^d,iw_mag(1)))/dsqrt(1.d0+(w(ix^d,iw_mag(2))/w(ix^d,iw_mag(1)))**2+&
              (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(1)))**2)
          else
            mf(ix^d,1)=0.d0
          end if
          if(w(ix^d,iw_mag(2))/=0.d0) then
            mf(ix^d,2)=sign(1.d0,w(ix^d,iw_mag(2)))/dsqrt(1.d0+(w(ix^d,iw_mag(1))/w(ix^d,iw_mag(2)))**2+&
              (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(2)))**2)
          else
            mf(ix^d,2)=0.d0
          end if
          if(w(ix^d,iw_mag(3))/=0.d0) then
            mf(ix^d,3)=sign(1.d0,w(ix^d,iw_mag(3)))/dsqrt(1.d0+(w(ix^d,iw_mag(1))/w(ix^d,iw_mag(3)))**2+&
              (w(ix^d,iw_mag(2))/w(ix^d,iw_mag(3)))**2)
          else
            mf(ix^d,3)=0.d0
          end if
         }
       {end do\}
      end if
    else
      mf(ix^s,1:ndim)=block%B0(ix^s,1:ndim,0)
    endif
    ! ixC is cell-corner index
    ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
    ! b unit vector at cell corner
   {^ifthreed
    do idims=1,3
   {do ix^db=ixcmin^db,ixcmax^db\}
      bc(ix^d,idims)=0.125d0*(mf(ix1,ix2,ix3,idims)+mf(ix1+1,ix2,ix3,idims)&
                     +mf(ix1,ix2+1,ix3,idims)+mf(ix1+1,ix2+1,ix3,idims)&
                     +mf(ix1,ix2,ix3+1,idims)+mf(ix1+1,ix2,ix3+1,idims)&
                     +mf(ix1,ix2+1,ix3+1,idims)+mf(ix1+1,ix2+1,ix3+1,idims))
   {end do\}
    end do
   }
   {^iftwod
    do idims=1,2
   {do ix^db=ixcmin^db,ixcmax^db\}
      bc(ix^d,idims)=0.25d0*(mf(ix1,ix2,idims)+mf(ix1+1,ix2,idims)&
                     +mf(ix1,ix2+1,idims)+mf(ix1+1,ix2+1,idims))
   {end do\}
    end do
   }
    ! T gradient at cell faces
    do idims=1,ndim
      ixbmin^d=ixmin^d;
      ixbmax^d=ixmax^d-kr(idims,^d);
      call gradientf(te,x,ixi^l,ixb^l,idims,gradt(ixi^s,idims))
    end do
    if(fl%tc_constant) then
      if(fl%tc_perpendicular) then
        ka(ixc^s)=fl%tc_k_para-fl%tc_k_perp
        ke(ixc^s)=fl%tc_k_perp
      else
        ka(ixc^s)=fl%tc_k_para
      end if
    else
      ! conductivity at cell center
      if(phys_trac) then
       {do ix^db=ixmin^db,ixmax^db\}
          if(te(ix^d) < block%wextra(ix^d,fl%Tcoff_)) then
            qdd(ix^d)=fl%tc_k_para*dsqrt(block%wextra(ix^d,fl%Tcoff_)**5)
          else
            qdd(ix^d)=fl%tc_k_para*dsqrt(te(ix^d)**5)
          end if
       {end do\}
      else
        qdd(ix^s)=fl%tc_k_para*dsqrt(te(ix^s)**5)
      end if
     ! cell corner parallel conductivity in ka
     {^ifthreed
     {do ix^db=ixcmin^db,ixcmax^db\}
        ka(ix^d)=0.125d0*(qdd(ix1,ix2,ix3)+qdd(ix1+1,ix2,ix3)&
                       +qdd(ix1,ix2+1,ix3)+qdd(ix1+1,ix2+1,ix3)&
                       +qdd(ix1,ix2,ix3+1)+qdd(ix1+1,ix2,ix3+1)&
                       +qdd(ix1,ix2+1,ix3+1)+qdd(ix1+1,ix2+1,ix3+1))
     {end do\}
     }
     {^iftwod
     {do ix^db=ixcmin^db,ixcmax^db\}
        ka(ix^d)=0.25d0*(qdd(ix1,ix2)+qdd(ix1+1,ix2)&
                       +qdd(ix1,ix2+1)+qdd(ix1+1,ix2+1))
     {end do\}
     }
      ! compensate with perpendicular conductivity
      if(fl%tc_perpendicular) then
        if(b0field) then
          qdd(ix^s)=fl%tc_k_perp*rho(ix^s)**2/((^c&(w(ix^s,iw_mag(^c))+block%B0(ix^s,^c,0))**2+)*dsqrt(te(ix^s))+smalldouble)
        else
          qdd(ix^s)=fl%tc_k_perp*rho(ix^s)**2/((^c&w(ix^s,iw_mag(^c))**2+)*dsqrt(te(ix^s))+smalldouble)
        end if
       {^ifthreed
       {do ix^db=ixcmin^db,ixcmax^db\}
          ke(ix^d)=0.125d0*(qdd(ix1,ix2,ix3)+qdd(ix1+1,ix2,ix3)&
                         +qdd(ix1,ix2+1,ix3)+qdd(ix1+1,ix2+1,ix3)&
                         +qdd(ix1,ix2,ix3+1)+qdd(ix1+1,ix2,ix3+1)&
                         +qdd(ix1,ix2+1,ix3+1)+qdd(ix1+1,ix2+1,ix3+1))
          if(ke(ix^d)<ka(ix^d)) then
            ka(ix^d)=ka(ix^d)-ke(ix^d)
          else
            ke(ix^d)=ka(ix^d)
            ka(ix^d)=0.d0
          end if
       {end do\}
       }
       {^iftwod
       {do ix^db=ixcmin^db,ixcmax^db\}
          ke(ix^d)=0.25d0*(qdd(ix1,ix2)+qdd(ix1+1,ix2)&
                         +qdd(ix1,ix2+1)+qdd(ix1+1,ix2+1))
          if(ke(ix^d)<ka(ix^d)) then
            ka(ix^d)=ka(ix^d)-ke(ix^d)
          else
            ke(ix^d)=ka(ix^d)
            ka(ix^d)=0.d0
          end if
       {end do\}
       }
      end if
    end if
    ! calculate thermal conduction flux with slope-limited symmetric scheme
    do idims=1,ndim
      ixamax^d=ixomax^d; ixamin^d=ixomin^d-kr(idims,^d);
     {^ifthreed
      if(idims==1) then
       {do ix^db=ixamin^db,ixamax^db\}
          ! averaged b at face centers
          ^d&bcf({ix^d},^d)=0.25d0*(bc({ix^d},^d)+bc(ix1,ix2-1,ix3,^d)&
                         +bc(ix1,ix2,ix3-1,^d)+bc(ix1,ix2-1,ix3-1,^d))\
          kaf(ix^d)=0.25d0*(ka(ix1,ix2,ix3)+ka(ix1,ix2-1,ix3)&
                         +ka(ix1,ix2,ix3-1)+ka(ix1,ix2-1,ix3-1))
          ! averaged thermal conductivity at face centers
          if(fl%tc_perpendicular) &
          kef(ix^d)=0.25d0*(ke(ix1,ix2,ix3)+ke(ix1,ix2-1,ix3)&
                         +ke(ix1,ix2,ix3-1)+ke(ix1,ix2-1,ix3-1))
       {end do\}
      else if(idims==2) then
       {do ix^db=ixamin^db,ixamax^db\}
          ^d&bcf({ix^d},^d)=0.25d0*(bc({ix^d},^d)+bc(ix1-1,ix2,ix3,^d)&
                         +bc(ix1,ix2,ix3-1,^d)+bc(ix1-1,ix2,ix3-1,^d))\
          kaf(ix^d)=0.25d0*(ka(ix1,ix2,ix3)+ka(ix1-1,ix2,ix3)&
                         +ka(ix1,ix2,ix3-1)+ka(ix1-1,ix2,ix3-1))
          if(fl%tc_perpendicular) &
          kef(ix^d)=0.25d0*(ke(ix1,ix2,ix3)+ke(ix1-1,ix2,ix3)&
                         +ke(ix1,ix2,ix3-1)+ke(ix1-1,ix2,ix3-1))
       {end do\}
      else
       {do ix^db=ixamin^db,ixamax^db\}
          ^d&bcf({ix^d},^d)=0.25d0*(bc({ix^d},^d)+bc(ix1,ix2-1,ix3,^d)&
                         +bc(ix1-1,ix2,ix3,^d)+bc(ix1-1,ix2-1,ix3,^d))\
          kaf(ix^d)=0.25d0*(ka(ix1,ix2,ix3)+ka(ix1,ix2-1,ix3)&
                         +ka(ix1-1,ix2,ix3)+ka(ix1-1,ix2-1,ix3))
          if(fl%tc_perpendicular) &
          kef(ix^d)=0.25d0*(ke(ix1,ix2,ix3)+ke(ix1,ix2-1,ix3)&
                         +ke(ix1-1,ix2,ix3)+ke(ix1-1,ix2-1,ix3))
       {end do\}
      end if
     }
     {^iftwod
      if(idims==1) then
       {do ix^db=ixamin^db,ixamax^db\}
          ^d&bcf({ix^d},^d)=0.5d0*(bc(ix1,ix2,^d)+bc(ix1,ix2-1,^d))\
          kaf(ix^d)=0.5d0*(ka(ix1,ix2)+ka(ix1,ix2-1))
          if(fl%tc_perpendicular) &
          kef(ix^d)=0.5d0*(ke(ix1,ix2)+ke(ix1,ix2-1))
       {end do\}
      else
       {do ix^db=ixamin^db,ixamax^db\}
          ^d&bcf({ix^d},^d)=0.5d0*(bc(ix1,ix2,^d)+bc(ix1-1,ix2,^d))\
          kaf(ix^d)=0.5d0*(ka(ix1,ix2)+ka(ix1-1,ix2))
          if(fl%tc_perpendicular) &
          kef(ix^d)=0.5d0*(ke(ix1,ix2)+ke(ix1-1,ix2))
       {end do\}
      end if
     }
      ! eq (19)
      ! temperature gradient at cell corner
     {^ifthreed
      if(idims==1) then
       {do ix^db=ixcmin^db,ixcmax^db\}
          qdd(ix^d)=0.25d0*(gradt(ix1,ix2,ix3,idims)+gradt(ix1,ix2+1,ix3,idims)&
                         +gradt(ix1,ix2,ix3+1,idims)+gradt(ix1,ix2+1,ix3+1,idims))
       {end do\}
      else if(idims==2) then
       {do ix^db=ixcmin^db,ixcmax^db\}
          qdd(ix^d)=0.25d0*(gradt(ix1,ix2,ix3,idims)+gradt(ix1+1,ix2,ix3,idims)&
                         +gradt(ix1,ix2,ix3+1,idims)+gradt(ix1+1,ix2,ix3+1,idims))
       {end do\}
      else
       {do ix^db=ixcmin^db,ixcmax^db\}
          qdd(ix^d)=0.25d0*(gradt(ix1,ix2,ix3,idims)+gradt(ix1+1,ix2,ix3,idims)&
                         +gradt(ix1,ix2+1,ix3,idims)+gradt(ix1+1,ix2+1,ix3,idims))
       {end do\}
      end if
     }
     {^iftwod
      if(idims==1) then
       {do ix^db=ixcmin^db,ixcmax^db\}
          qdd(ix^d)=0.5d0*(gradt(ix1,ix2,idims)+gradt(ix1,ix2+1,idims))
       {end do\}
      else
       {do ix^db=ixcmin^db,ixcmax^db\}
          qdd(ix^d)=0.5d0*(gradt(ix1,ix2,idims)+gradt(ix1+1,ix2,idims))
       {end do\}
      end if
     }
      ! eq (21)
     {^ifthreed
      if(idims==1) then
       {do ix^db=ixamin^db,ixamax^db\}
          minq=min(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          maxq=max(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          if(qdd(ix^d)<minq) then
            qd(ix^d,1)=minq
          else if(qdd(ix^d)>maxq) then
            qd(ix^d,1)=maxq
          else
            qd(ix^d,1)=qdd(ix^d)
          end if
          if(qdd(ix1,ix2-1,ix3)<minq) then
            qd(ix^d,2)=minq
          else if(qdd(ix1,ix2-1,ix3)>maxq) then
            qd(ix^d,2)=maxq
          else
            qd(ix^d,2)=qdd(ix1,ix2-1,ix3)
          end if
          if(qdd(ix1,ix2,ix3-1)<minq) then
            qd(ix^d,3)=minq
          else if(qdd(ix1,ix2,ix3-1)>maxq) then
            qd(ix^d,3)=maxq
          else
            qd(ix^d,3)=qdd(ix1,ix2,ix3-1)
          end if
          if(qdd(ix1,ix2-1,ix3-1)<minq) then
            qd(ix^d,4)=minq
          else if(qdd(ix1,ix2-1,ix3-1)>maxq) then
            qd(ix^d,4)=maxq
          else
            qd(ix^d,4)=qdd(ix1,ix2-1,ix3-1)
          end if
          qvec(ix^d,idims)=kaf(ix^d)*0.25d0*(bc(ix^d,idims)**2*qd(ix^d,1)+bc(ix1,ix2-1,ix3,idims)**2*qd(ix^d,2)&
                         +bc(ix1,ix2,ix3-1,idims)**2*qd(ix^d,3)+bc(ix1,ix2-1,ix3-1,idims)**2*qd(ix^d,4))
          if(fl%tc_perpendicular) &
          qvec(ix^d,idims)=qvec(ix^d,idims)+kef(ix^d)*0.25d0*(qd(ix^d,1)+qd(ix^d,2)+qd(ix^d,3)+qd(ix^d,4))
       {end do\}
      else if(idims==2) then
       {do ix^db=ixamin^db,ixamax^db\}
          minq=min(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          maxq=max(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          if(qdd(ix^d)<minq) then
            qd(ix^d,1)=minq
          else if(qdd(ix^d)>maxq) then
            qd(ix^d,1)=maxq
          else
            qd(ix^d,1)=qdd(ix^d)
          end if
          if(qdd(ix1-1,ix2,ix3)<minq) then
            qd(ix^d,2)=minq
          else if(qdd(ix1-1,ix2,ix3)>maxq) then
            qd(ix^d,2)=maxq
          else
            qd(ix^d,2)=qdd(ix1-1,ix2,ix3)
          end if
          if(qdd(ix1,ix2,ix3-1)<minq) then
            qd(ix^d,3)=minq
          else if(qdd(ix1,ix2,ix3-1)>maxq) then
            qd(ix^d,3)=maxq
          else
            qd(ix^d,3)=qdd(ix1,ix2,ix3-1)
          end if
          if(qdd(ix1-1,ix2,ix3-1)<minq) then
            qd(ix^d,4)=minq
          else if(qdd(ix1-1,ix2,ix3-1)>maxq) then
            qd(ix^d,4)=maxq
          else
            qd(ix^d,4)=qdd(ix1-1,ix2,ix3-1)
          end if
          qvec(ix^d,idims)=kaf(ix^d)*0.25d0*(bc(ix^d,idims)**2*qd(ix^d,1)+bc(ix1-1,ix2,ix3,idims)**2*qd(ix^d,2)&
                         +bc(ix1,ix2,ix3-1,idims)**2*qd(ix^d,3)+bc(ix1-1,ix2,ix3-1,idims)**2*qd(ix^d,4))
          if(fl%tc_perpendicular) &
          qvec(ix^d,idims)=qvec(ix^d,idims)+kef(ix^d)*0.25d0*(qd(ix^d,1)+qd(ix^d,2)+qd(ix^d,3)+qd(ix^d,4))
       {end do\}
      else
       {do ix^db=ixamin^db,ixamax^db\}
          minq=min(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          maxq=max(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          if(qdd(ix^d)<minq) then
            qd(ix^d,1)=minq
          else if(qdd(ix^d)>maxq) then
            qd(ix^d,1)=maxq
          else
            qd(ix^d,1)=qdd(ix^d)
          end if
          if(qdd(ix1-1,ix2,ix3)<minq) then
            qd(ix^d,2)=minq
          else if(qdd(ix1-1,ix2,ix3)>maxq) then
            qd(ix^d,2)=maxq
          else
            qd(ix^d,2)=qdd(ix1-1,ix2,ix3)
          end if
          if(qdd(ix1,ix2-1,ix3)<minq) then
            qd(ix^d,3)=minq
          else if(qdd(ix1,ix2-1,ix3)>maxq) then
            qd(ix^d,3)=maxq
          else
            qd(ix^d,3)=qdd(ix1,ix2-1,ix3)
          end if
          if(qdd(ix1-1,ix2-1,ix3)<minq) then
            qd(ix^d,4)=minq
          else if(qdd(ix1-1,ix2-1,ix3)>maxq) then
            qd(ix^d,4)=maxq
          else
            qd(ix^d,4)=qdd(ix1-1,ix2-1,ix3)
          end if
          qvec(ix^d,idims)=kaf(ix^d)*0.25d0*(bc(ix^d,idims)**2*qd(ix^d,1)+bc(ix1-1,ix2,ix3,idims)**2*qd(ix^d,2)&
                         +bc(ix1,ix2-1,ix3,idims)**2*qd(ix^d,3)+bc(ix1-1,ix2-1,ix3,idims)**2*qd(ix^d,4))
          if(fl%tc_perpendicular) &
          qvec(ix^d,idims)=qvec(ix^d,idims)+kef(ix^d)*0.25d0*(qd(ix^d,1)+qd(ix^d,2)+qd(ix^d,3)+qd(ix^d,4))
       {end do\}
      end if
     }
     {^iftwod
      if(idims==1) then
       {do ix^db=ixamin^db,ixamax^db\}
          minq=min(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          maxq=max(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          if(qdd(ix^d)<minq) then
            qd(ix^d,1)=minq
          else if(qdd(ix^d)>maxq) then
            qd(ix^d,1)=maxq
          else
            qd(ix^d,1)=qdd(ix^d)
          end if
          if(qdd(ix1,ix2-1)<minq) then
            qd(ix^d,2)=minq
          else if(qdd(ix1,ix2-1)>maxq) then
            qd(ix^d,2)=maxq
          else
            qd(ix^d,2)=qdd(ix1,ix2-1)
          end if
          qvec(ix^d,idims)=kaf(ix^d)*0.5d0*(bc(ix1,ix2,idims)**2*qd(ix^d,1)+bc(ix1,ix2-1,idims)**2*qd(ix^d,2))
          if(fl%tc_perpendicular) &
          qvec(ix^d,idims)=qvec(ix^d,idims)+kef(ix^d)*0.5d0*(qd(ix^d,1)+qd(ix^d,2))
       {end do\}
      else
       {do ix^db=ixamin^db,ixamax^db\}
          minq=min(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          maxq=max(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          if(qdd(ix^d)<minq) then
            qd(ix^d,1)=minq
          else if(qdd(ix^d)>maxq) then
            qd(ix^d,1)=maxq
          else
            qd(ix^d,1)=qdd(ix^d)
          end if
          if(qdd(ix1-1,ix2)<minq) then
            qd(ix^d,2)=minq
          else if(qdd(ix1-1,ix2)>maxq) then
            qd(ix^d,2)=maxq
          else
            qd(ix^d,2)=qdd(ix1-1,ix2)
          end if
          qvec(ix^d,idims)=kaf(ix^d)*0.5d0*(bc(ix1,ix2,idims)**2*qd(ix^d,1)+bc(ix1-1,ix2,idims)**2*qd(ix^d,2))
          if(fl%tc_perpendicular) &
          qvec(ix^d,idims)=qvec(ix^d,idims)+kef(ix^d)*0.5d0*(qd(ix^d,1)+qd(ix^d,2))
       {end do\}
      end if
     }
      ! calculate normal of magnetic field
      ixb^l=ixa^l+kr(idims,^d);
      bnorm(ixa^s)=0.5d0*(mf(ixa^s,idims)+mf(ixb^s,idims))
      ! limited transverse component, eq (17)
      ixbmin^d=ixamin^d;
      ixbmax^d=ixamax^d+kr(idims,^d);
      do idir=1,ndim
        if(idir==idims) cycle
        qdd(ixi^s)=slope_limiter(gradt(ixi^s,idir),ixi^l,ixb^l,idir,-1,fl%tc_slope_limiter)
        qdd(ixi^s)=slope_limiter(qdd,ixi^l,ixa^l,idims,1,fl%tc_slope_limiter)
        qvec(ixa^s,idims)=qvec(ixa^s,idims)+kaf(ixa^s)*bnorm(ixa^s)*bcf(ixa^s,idir)*qdd(ixa^s)
      end do
      if(fl%tc_saturate) then
        ! consider saturation (Cowie and Mckee 1977 ApJ, 211, 135: phi=1.1, Balbus and Mckee 1982 ApJ, 252, 529: phi=0.3)
        ! unsigned saturated TC flux = 5 phi rho c**3, c=sqrt(p/rho) is isothermal sound speed, phi=0.3
        ixb^l=ixa^l+kr(idims,^d);
        qdd(ixa^s)=0.75d0*(rho(ixa^s)+rho(ixb^s))*dsqrt(0.5d0*(ciso2(ixa^s)+ciso2(ixb^s)))**3*dabs(bnorm(ixa^s))
       {do ix^db=ixamin^db,ixamax^db\}
          if(dabs(qvec(ix^d,idims))>qdd(ix^d)) then
            qvec(ix^d,idims)=sign(1.d0,qvec(ix^d,idims))*qdd(ix^d)
          end if
       {end do\}
      end if
    end do
  end subroutine set_source_tc_mhd
  subroutine set_source_tc_mhd_geo(ixI^L,ixO^L,w,x,fl,qvec,rho,Te,ciso2,alpha)
    use mod_global_parameters
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) ::  x(ixI^S,1:ndim)
    double precision, intent(in) ::  w(ixI^S,1:nw)
    type(tc_fluid), intent(in) :: fl
    double precision, intent(in) :: rho(ixI^S),Te(ixI^S),ciso2(ixI^S)
    double precision, intent(in) :: alpha
    double precision, intent(out) :: qvec(ixI^S,1:ndim)
 
    !! qdd store the heat conduction energy changing rate
    double precision, dimension(ixI^S,1:ndim) :: mf,Bc,Bcf,gradT
    double precision, dimension(ixI^S) :: ka,kaf,ke,kef,qdd,Bnorm
    double precision :: minq,maxq,qd(ixI^S,2**(ndim-1)), blocal(ndir)
    integer :: idims,idir,ix^D,ix^L,ixC^L,ixA^L,ixB^L
 
    ix^l=ixo^l^ladd1;
 
    ! T gradient at cell faces
    ! b unit vector mf: magnetic field direction vector
    if(allocated(iw_mag)) then
      if(b0field) then
       {do ix^db=ixmin^db,ixmax^db\}
          ^c&blocal(^c)=w({ix^d},iw_mag(^c))+block%B0({ix^d},^c,0)\
         {^iftwod
          if(blocal(1)/=0.d0) then
            mf(ix^d,1)=sign(1.d0,blocal(1))/dsqrt(1.d0+(^ce&(blocal(^ce)/blocal(1))**2+))
          else
            mf(ix^d,1)=0.d0
          end if
          if(blocal(2)/=0.d0) then
            mf(ix^d,2)=sign(1.d0,blocal(2))/dsqrt(1.d0+(^cf&(blocal(^cf)/blocal(2))**2+))
          else
            mf(ix^d,2)=0.d0
          end if
         }
         {^ifthreed
          if(blocal(1)/=0.d0) then
            mf(ix^d,1)=sign(1.d0,blocal(1))/dsqrt(1.d0+(blocal(2)/blocal(1))**2+(blocal(3)/blocal(1))**2)
          else
            mf(ix^d,1)=0.d0
          end if
          if(blocal(2)/=0.d0) then
            mf(ix^d,2)=sign(1.d0,blocal(2))/dsqrt(1.d0+(blocal(1)/blocal(2))**2+(blocal(3)/blocal(2))**2)
          else
            mf(ix^d,2)=0.d0
          end if
          if(blocal(3)/=0.d0) then
            mf(ix^d,3)=sign(1.d0,blocal(3))/dsqrt(1.d0+(blocal(1)/blocal(3))**2+(blocal(2)/blocal(3))**2)
          else
            mf(ix^d,3)=0.d0
          end if
         }
       {end do\}
      else
       {do ix^db=ixmin^db,ixmax^db\}
         {^iftwod
          if(w(ix^d,iw_mag(1))/=0.d0) then
            mf(ix^d,1)=sign(1.d0,w(ix^d,iw_mag(1)))/dsqrt(1.d0+(^ce&(w(ix^d,iw_mag(^ce))/w(ix^d,iw_mag(1)))**2+))
          else
            mf(ix^d,1)=0.d0
          end if
          if(w(ix^d,iw_mag(2))/=0.d0) then
            mf(ix^d,2)=sign(1.d0,w(ix^d,iw_mag(2)))/dsqrt(1.d0+(^cf&(w(ix^d,iw_mag(^cf))/w(ix^d,iw_mag(2)))**2+))
          else
            mf(ix^d,2)=0.d0
          end if
         }
         {^ifthreed
          if(w(ix^d,iw_mag(1))/=0.d0) then
            mf(ix^d,1)=sign(1.d0,w(ix^d,iw_mag(1)))/dsqrt(1.d0+(w(ix^d,iw_mag(2))/w(ix^d,iw_mag(1)))**2+&
              (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(1)))**2)
          else
            mf(ix^d,1)=0.d0
          end if
          if(w(ix^d,iw_mag(2))/=0.d0) then
            mf(ix^d,2)=sign(1.d0,w(ix^d,iw_mag(2)))/dsqrt(1.d0+(w(ix^d,iw_mag(1))/w(ix^d,iw_mag(2)))**2+&
              (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(2)))**2)
          else
            mf(ix^d,2)=0.d0
          end if
          if(w(ix^d,iw_mag(3))/=0.d0) then
            mf(ix^d,3)=sign(1.d0,w(ix^d,iw_mag(3)))/dsqrt(1.d0+(w(ix^d,iw_mag(1))/w(ix^d,iw_mag(3)))**2+&
              (w(ix^d,iw_mag(2))/w(ix^d,iw_mag(3)))**2)
          else
            mf(ix^d,3)=0.d0
          end if
         }
       {end do\}
      end if
    else
      mf(ix^s,1:ndim)=block%B0(ix^s,1:ndim,0)
    endif
    ! ixC is cell-corner index
    ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
    ! b unit vector at cell corner
   {^ifthreed
    do idims=1,3
   {do ix^db=ixcmin^db,ixcmax^db\}
      bc(ix^d,idims)=(mf(ix1,ix2,ix3,idims)*block%dvolume(ix1,ix2,ix3)+mf(ix1+1,ix2,ix3,idims)*block%dvolume(ix1+1,ix2,ix3)&
                     +mf(ix1,ix2+1,ix3,idims)*block%dvolume(ix1,ix2+1,ix3)+mf(ix1+1,ix2+1,ix3,idims)*block%dvolume(ix1+1,ix2+1,ix3)&
                     +mf(ix1,ix2,ix3+1,idims)*block%dvolume(ix1,ix2,ix3+1)+mf(ix1+1,ix2,ix3+1,idims)*block%dvolume(ix1+1,ix2,ix3+1)&
                     +mf(ix1,ix2+1,ix3+1,idims)*block%dvolume(ix1,ix2+1,ix3+1)+mf(ix1+1,ix2+1,ix3+1,idims)*block%dvolume(ix1+1,ix2+1,ix3+1))&
            /(block%dvolume(ix1,ix2,ix3)+block%dvolume(ix1+1,ix2,ix3)+block%dvolume(ix1,ix2+1,ix3)+block%dvolume(ix1+1,ix2+1,ix3)&
             +block%dvolume(ix1,ix2,ix3+1)+block%dvolume(ix1+1,ix2,ix3+1)+block%dvolume(ix1,ix2+1,ix3+1)+block%dvolume(ix1+1,ix2+1,ix3+1))
   {end do\}
    end do
   }
   {^iftwod
    do idims=1,2
   {do ix^db=ixcmin^db,ixcmax^db\}
      bc(ix^d,idims)=(mf(ix1,ix2,idims)*block%dvolume(ix1,ix2)+mf(ix1+1,ix2,idims)*block%dvolume(ix1+1,ix2)&
                     +mf(ix1,ix2+1,idims)*block%dvolume(ix1,ix2+1)+mf(ix1+1,ix2+1,idims)*block%dvolume(ix1+1,ix2+1))&
                 /(block%dvolume(ix1,ix2)+block%dvolume(ix1+1,ix2)+block%dvolume(ix1,ix2+1)+block%dvolume(ix1+1,ix2+1))
   {end do\}
    end do
   }
    ! T gradient at cell faces
    do idims=1,ndim
      ixbmin^d=ixmin^d;
      ixbmax^d=ixmax^d-kr(idims,^d);
      call gradientf(te,x,ixi^l,ixb^l,idims,gradt(ixi^s,idims))
    end do
    if(fl%tc_constant) then
      if(fl%tc_perpendicular) then
        ka(ixc^s)=fl%tc_k_para-fl%tc_k_perp
        ke(ixc^s)=fl%tc_k_perp
      else
        ka(ixc^s)=fl%tc_k_para
      end if
    else
      ! conductivity at cell center
      if(phys_trac) then
       {do ix^db=ixmin^db,ixmax^db\}
          if(te(ix^d) < block%wextra(ix^d,fl%Tcoff_)) then
            qdd(ix^d)=fl%tc_k_para*dsqrt(block%wextra(ix^d,fl%Tcoff_)**5)
          else
            qdd(ix^d)=fl%tc_k_para*dsqrt(te(ix^d)**5)
          end if
       {end do\}
      else
        qdd(ix^s)=fl%tc_k_para*dsqrt(te(ix^s)**5)
      end if
     ! cell corner parallel conductivity in ka
     {^ifthreed
     {do ix^db=ixcmin^db,ixcmax^db\}
        ka(ix^d)=(qdd(ix1,ix2,ix3)*block%dvolume(ix1,ix2,ix3)+qdd(ix1+1,ix2,ix3)*block%dvolume(ix1+1,ix2,ix3)&
           +qdd(ix1,ix2+1,ix3)*block%dvolume(ix1,ix2+1,ix3)+qdd(ix1+1,ix2+1,ix3)*block%dvolume(ix1+1,ix2+1,ix3)&
           +qdd(ix1,ix2,ix3+1)*block%dvolume(ix1,ix2,ix3+1)+qdd(ix1+1,ix2,ix3+1)*block%dvolume(ix1+1,ix2,ix3+1)&
           +qdd(ix1,ix2+1,ix3+1)*block%dvolume(ix1,ix2+1,ix3+1)+qdd(ix1+1,ix2+1,ix3+1)*block%dvolume(ix1+1,ix2+1,ix3+1))&
          /(block%dvolume(ix1,ix2,ix3)+block%dvolume(ix1+1,ix2,ix3)+block%dvolume(ix1,ix2+1,ix3)+block%dvolume(ix1+1,ix2+1,ix3)&
           +block%dvolume(ix1,ix2,ix3+1)+block%dvolume(ix1+1,ix2,ix3+1)+block%dvolume(ix1,ix2+1,ix3+1)+block%dvolume(ix1+1,ix2+1,ix3+1))
     {end do\}
     }
     {^iftwod
     {do ix^db=ixcmin^db,ixcmax^db\}
        ka(ix^d)=(qdd(ix1,ix2)*block%dvolume(ix1,ix2)+qdd(ix1+1,ix2)*block%dvolume(ix1+1,ix2)&
           +qdd(ix1,ix2+1)*block%dvolume(ix1,ix2+1)+qdd(ix1+1,ix2+1)*block%dvolume(ix1+1,ix2+1))&
         /(block%dvolume(ix1,ix2)+block%dvolume(ix1+1,ix2)+block%dvolume(ix1,ix2+1)+block%dvolume(ix1+1,ix2+1))
     {end do\}
     }
      ! compensate with perpendicular conductivity
      if(fl%tc_perpendicular) then
        if(b0field) then
          qdd(ix^s)=fl%tc_k_perp*rho(ix^s)**2/((^c&(w(ix^s,iw_mag(^c))+block%B0(ix^s,^c,0))**2+)*dsqrt(te(ix^s))+smalldouble)
        else
          qdd(ix^s)=fl%tc_k_perp*rho(ix^s)**2/((^c&w(ix^s,iw_mag(^c))**2+)*dsqrt(te(ix^s))+smalldouble)
        end if
       {^ifthreed
       {do ix^db=ixcmin^db,ixcmax^db\}
          ke(ix^d)=(qdd(ix1,ix2,ix3)*block%dvolume(ix1,ix2,ix3)+qdd(ix1+1,ix2,ix3)*block%dvolume(ix1+1,ix2,ix3)&
             +qdd(ix1,ix2+1,ix3)*block%dvolume(ix1,ix2+1,ix3)+qdd(ix1+1,ix2+1,ix3)*block%dvolume(ix1+1,ix2+1,ix3)&
             +qdd(ix1,ix2,ix3+1)*block%dvolume(ix1,ix2,ix3+1)+qdd(ix1+1,ix2,ix3+1)*block%dvolume(ix1+1,ix2,ix3+1)&
             +qdd(ix1,ix2+1,ix3+1)*block%dvolume(ix1,ix2+1,ix3+1)+qdd(ix1+1,ix2+1,ix3+1)*block%dvolume(ix1+1,ix2+1,ix3+1))&
            /(block%dvolume(ix1,ix2,ix3)+block%dvolume(ix1+1,ix2,ix3)+block%dvolume(ix1,ix2+1,ix3)+block%dvolume(ix1+1,ix2+1,ix3)&
             +block%dvolume(ix1,ix2,ix3+1)+block%dvolume(ix1+1,ix2,ix3+1)+block%dvolume(ix1,ix2+1,ix3+1)+block%dvolume(ix1+1,ix2+1,ix3+1))
          if(ke(ix^d)<ka(ix^d)) then
            ka(ix^d)=ka(ix^d)-ke(ix^d)
          else
            ke(ix^d)=ka(ix^d)
            ka(ix^d)=0.d0
          end if
       {end do\}
       }
       {^iftwod
       {do ix^db=ixcmin^db,ixcmax^db\}
          ke(ix^d)=(qdd(ix1,ix2)*block%dvolume(ix1,ix2)+qdd(ix1+1,ix2)*block%dvolume(ix1+1,ix2)&
             +qdd(ix1,ix2+1)*block%dvolume(ix1,ix2+1)+qdd(ix1+1,ix2+1)*block%dvolume(ix1+1,ix2+1))&
           /(block%dvolume(ix1,ix2)+block%dvolume(ix1+1,ix2)+block%dvolume(ix1,ix2+1)+block%dvolume(ix1+1,ix2+1))
          if(ke(ix^d)<ka(ix^d)) then
            ka(ix^d)=ka(ix^d)-ke(ix^d)
          else
            ke(ix^d)=ka(ix^d)
            ka(ix^d)=0.d0
          end if
       {end do\}
       }
      end if
    end if
    ! calculate thermal conduction flux with slope-limited symmetric scheme
    do idims=1,ndim
      ixamax^d=ixomax^d; ixamin^d=ixomin^d-kr(idims,^d);
     {^ifthreed
      if(idims==1) then
       {do ix^db=ixamin^db,ixamax^db\}
          ! averaged b at face centers
          ^d&bcf({ix^d},^d)=0.25d0*(bc({ix^d},^d)+bc(ix1,ix2-1,ix3,^d)&
                         +bc(ix1,ix2,ix3-1,^d)+bc(ix1,ix2-1,ix3-1,^d))\
          kaf(ix^d)=0.25d0*(ka(ix1,ix2,ix3)+ka(ix1,ix2-1,ix3)&
                         +ka(ix1,ix2,ix3-1)+ka(ix1,ix2-1,ix3-1))
          ! averaged thermal conductivity at face centers
          if(fl%tc_perpendicular) &
          kef(ix^d)=0.25d0*(ke(ix1,ix2,ix3)+ke(ix1,ix2-1,ix3)&
                         +ke(ix1,ix2,ix3-1)+ke(ix1,ix2-1,ix3-1))
       {end do\}
      else if(idims==2) then
       {do ix^db=ixamin^db,ixamax^db\}
          ^d&bcf({ix^d},^d)=0.25d0*(bc({ix^d},^d)+bc(ix1-1,ix2,ix3,^d)&
                         +bc(ix1,ix2,ix3-1,^d)+bc(ix1-1,ix2,ix3-1,^d))\
          kaf(ix^d)=0.25d0*(ka(ix1,ix2,ix3)+ka(ix1-1,ix2,ix3)&
                         +ka(ix1,ix2,ix3-1)+ka(ix1-1,ix2,ix3-1))
          if(fl%tc_perpendicular) &
          kef(ix^d)=0.25d0*(ke(ix1,ix2,ix3)+ke(ix1-1,ix2,ix3)&
                         +ke(ix1,ix2,ix3-1)+ke(ix1-1,ix2,ix3-1))
       {end do\}
      else
       {do ix^db=ixamin^db,ixamax^db\}
          ^d&bcf({ix^d},^d)=0.25d0*(bc({ix^d},^d)+bc(ix1,ix2-1,ix3,^d)&
                         +bc(ix1-1,ix2,ix3,^d)+bc(ix1-1,ix2-1,ix3,^d))\
          kaf(ix^d)=0.25d0*(ka(ix1,ix2,ix3)+ka(ix1,ix2-1,ix3)&
                         +ka(ix1-1,ix2,ix3)+ka(ix1-1,ix2-1,ix3))
          if(fl%tc_perpendicular) &
          kef(ix^d)=0.25d0*(ke(ix1,ix2,ix3)+ke(ix1,ix2-1,ix3)&
                         +ke(ix1-1,ix2,ix3)+ke(ix1-1,ix2-1,ix3))
       {end do\}
      end if
     }
     {^iftwod
      if(idims==1) then
       {do ix^db=ixamin^db,ixamax^db\}
          ^d&bcf({ix^d},^d)=0.5d0*(bc(ix1,ix2,^d)+bc(ix1,ix2-1,^d))\
          kaf(ix^d)=0.5d0*(ka(ix1,ix2)+ka(ix1,ix2-1))
          if(fl%tc_perpendicular) &
          kef(ix^d)=0.5d0*(ke(ix1,ix2)+ke(ix1,ix2-1))
       {end do\}
      else
       {do ix^db=ixamin^db,ixamax^db\}
          ^d&bcf({ix^d},^d)=0.5d0*(bc(ix1,ix2,^d)+bc(ix1-1,ix2,^d))\
          kaf(ix^d)=0.5d0*(ka(ix1,ix2)+ka(ix1-1,ix2))
          if(fl%tc_perpendicular) &
          kef(ix^d)=0.5d0*(ke(ix1,ix2)+ke(ix1-1,ix2))
       {end do\}
      end if
     }
      ! eq (19)
      ! temperature gradient at cell corner
     {^ifthreed
      if(idims==1) then
       {do ix^db=ixcmin^db,ixcmax^db\}
          qdd(ix^d)=(gradt(ix1,ix2,ix3,idims)*block%surfaceC(ix1,ix2,ix3,1)&
                  +gradt(ix1,ix2+1,ix3,idims)*block%surfaceC(ix1,ix2+1,ix3,1)&
                  +gradt(ix1,ix2,ix3+1,idims)*block%surfaceC(ix1,ix2,ix3+1,1)&
                +gradt(ix1,ix2+1,ix3+1,idims)*block%surfaceC(ix1,ix2+1,ix3+1,1))/&
              (block%surfaceC(ix1,ix2,ix3,1)+block%surfaceC(ix1,ix2+1,ix3,1)&
              +block%surfaceC(ix1,ix2,ix3+1,1)+block%surfaceC(ix1,ix2+1,ix3+1,1)+smalldouble**2)
       {end do\}
      else if(idims==2) then
       {do ix^db=ixcmin^db,ixcmax^db\}
          qdd(ix^d)=(gradt(ix1,ix2,ix3,idims)*block%surfaceC(ix1,ix2,ix3,2)&
                  +gradt(ix1+1,ix2,ix3,idims)*block%surfaceC(ix1+1,ix2,ix3,2)&
                  +gradt(ix1,ix2,ix3+1,idims)*block%surfaceC(ix1,ix2,ix3+1,2)&
                +gradt(ix1+1,ix2,ix3+1,idims)*block%surfaceC(ix1+1,ix2,ix3+1,2))/&
            (block%surfaceC(ix1,ix2,ix3,2)+block%surfaceC(ix1+1,ix2,ix3,2)&
          +block%surfaceC(ix1,ix2,ix3+1,2)+block%surfaceC(ix1+1,ix2,ix3+1,2)+smalldouble**2)
       {end do\}
      else
       {do ix^db=ixcmin^db,ixcmax^db\}
          qdd(ix^d)=(gradt(ix1,ix2,ix3,idims)*block%surfaceC(ix1,ix2,ix3,3)&
                  +gradt(ix1+1,ix2,ix3,idims)*block%surfaceC(ix1+1,ix2,ix3,3)&
                  +gradt(ix1,ix2+1,ix3,idims)*block%surfaceC(ix1,ix2+1,ix3,3)&
                +gradt(ix1+1,ix2+1,ix3,idims)*block%surfaceC(ix1+1,ix2+1,ix3,3))/&
               (block%surfaceC(ix1,ix2,ix3,3)+block%surfaceC(ix1+1,ix2,ix3,3)&
             +block%surfaceC(ix1,ix2+1,ix3,3)+block%surfaceC(ix1+1,ix2+1,ix3,3))
       {end do\}
      end if
     }
     {^iftwod
      if(idims==1) then
       {do ix^db=ixcmin^db,ixcmax^db\}
          qdd(ix^d)=(gradt(ix1,ix2,idims)*block%surfaceC(ix1,ix2,1)&
                  +gradt(ix1,ix2+1,idims)*block%surfaceC(ix1,ix2+1,1))/&
               (block%surfaceC(ix1,ix2,1)+block%surfaceC(ix1,ix2+1,1))
       {end do\}
      else
       {do ix^db=ixcmin^db,ixcmax^db\}
          qdd(ix^d)=(gradt(ix1,ix2,idims)*block%surfaceC(ix1,ix2,2)&
                  +gradt(ix1+1,ix2,idims)*block%surfaceC(ix1+1,ix2,2))/&
               (block%surfaceC(ix1,ix2,2)+block%surfaceC(ix1+1,ix2,2)+smalldouble)
       {end do\}
      end if
     }
      ! eq (21)
     {^ifthreed
      if(idims==1) then
       {do ix^db=ixamin^db,ixamax^db\}
          minq=min(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          maxq=max(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          if(qdd(ix^d)<minq) then
            qd(ix^d,1)=minq
          else if(qdd(ix^d)>maxq) then
            qd(ix^d,1)=maxq
          else
            qd(ix^d,1)=qdd(ix^d)
          end if
          if(qdd(ix1,ix2-1,ix3)<minq) then
            qd(ix^d,2)=minq
          else if(qdd(ix1,ix2-1,ix3)>maxq) then
            qd(ix^d,2)=maxq
          else
            qd(ix^d,2)=qdd(ix1,ix2-1,ix3)
          end if
          if(qdd(ix1,ix2,ix3-1)<minq) then
            qd(ix^d,3)=minq
          else if(qdd(ix1,ix2,ix3-1)>maxq) then
            qd(ix^d,3)=maxq
          else
            qd(ix^d,3)=qdd(ix1,ix2,ix3-1)
          end if
          if(qdd(ix1,ix2-1,ix3-1)<minq) then
            qd(ix^d,4)=minq
          else if(qdd(ix1,ix2-1,ix3-1)>maxq) then
            qd(ix^d,4)=maxq
          else
            qd(ix^d,4)=qdd(ix1,ix2-1,ix3-1)
          end if
          qvec(ix^d,idims)=kaf(ix^d)*0.25d0*(bc(ix^d,idims)**2*qd(ix^d,1)+bc(ix1,ix2-1,ix3,idims)**2*qd(ix^d,2)&
                         +bc(ix1,ix2,ix3-1,idims)**2*qd(ix^d,3)+bc(ix1,ix2-1,ix3-1,idims)**2*qd(ix^d,4))
          if(fl%tc_perpendicular) &
          qvec(ix^d,idims)=qvec(ix^d,idims)+kef(ix^d)*0.25d0*(qd(ix^d,1)+qd(ix^d,2)+qd(ix^d,3)+qd(ix^d,4))
       {end do\}
      else if(idims==2) then
       {do ix^db=ixamin^db,ixamax^db\}
          minq=min(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          maxq=max(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          if(qdd(ix^d)<minq) then
            qd(ix^d,1)=minq
          else if(qdd(ix^d)>maxq) then
            qd(ix^d,1)=maxq
          else
            qd(ix^d,1)=qdd(ix^d)
          end if
          if(qdd(ix1-1,ix2,ix3)<minq) then
            qd(ix^d,2)=minq
          else if(qdd(ix1-1,ix2,ix3)>maxq) then
            qd(ix^d,2)=maxq
          else
            qd(ix^d,2)=qdd(ix1-1,ix2,ix3)
          end if
          if(qdd(ix1,ix2,ix3-1)<minq) then
            qd(ix^d,3)=minq
          else if(qdd(ix1,ix2,ix3-1)>maxq) then
            qd(ix^d,3)=maxq
          else
            qd(ix^d,3)=qdd(ix1,ix2,ix3-1)
          end if
          if(qdd(ix1-1,ix2,ix3-1)<minq) then
            qd(ix^d,4)=minq
          else if(qdd(ix1-1,ix2,ix3-1)>maxq) then
            qd(ix^d,4)=maxq
          else
            qd(ix^d,4)=qdd(ix1-1,ix2,ix3-1)
          end if
          qvec(ix^d,idims)=kaf(ix^d)*0.25d0*(bc(ix^d,idims)**2*qd(ix^d,1)+bc(ix1-1,ix2,ix3,idims)**2*qd(ix^d,2)&
                         +bc(ix1,ix2,ix3-1,idims)**2*qd(ix^d,3)+bc(ix1-1,ix2,ix3-1,idims)**2*qd(ix^d,4))
          if(fl%tc_perpendicular) &
          qvec(ix^d,idims)=qvec(ix^d,idims)+kef(ix^d)*0.25d0*(qd(ix^d,1)+qd(ix^d,2)+qd(ix^d,3)+qd(ix^d,4))
       {end do\}
      else
       {do ix^db=ixamin^db,ixamax^db\}
          minq=min(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          maxq=max(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          if(qdd(ix^d)<minq) then
            qd(ix^d,1)=minq
          else if(qdd(ix^d)>maxq) then
            qd(ix^d,1)=maxq
          else
            qd(ix^d,1)=qdd(ix^d)
          end if
          if(qdd(ix1-1,ix2,ix3)<minq) then
            qd(ix^d,2)=minq
          else if(qdd(ix1-1,ix2,ix3)>maxq) then
            qd(ix^d,2)=maxq
          else
            qd(ix^d,2)=qdd(ix1-1,ix2,ix3)
          end if
          if(qdd(ix1,ix2-1,ix3)<minq) then
            qd(ix^d,3)=minq
          else if(qdd(ix1,ix2-1,ix3)>maxq) then
            qd(ix^d,3)=maxq
          else
            qd(ix^d,3)=qdd(ix1,ix2-1,ix3)
          end if
          if(qdd(ix1-1,ix2-1,ix3)<minq) then
            qd(ix^d,4)=minq
          else if(qdd(ix1-1,ix2-1,ix3)>maxq) then
            qd(ix^d,4)=maxq
          else
            qd(ix^d,4)=qdd(ix1-1,ix2-1,ix3)
          end if
          qvec(ix^d,idims)=kaf(ix^d)*0.25d0*(bc(ix^d,idims)**2*qd(ix^d,1)+bc(ix1-1,ix2,ix3,idims)**2*qd(ix^d,2)&
                         +bc(ix1,ix2-1,ix3,idims)**2*qd(ix^d,3)+bc(ix1-1,ix2-1,ix3,idims)**2*qd(ix^d,4))
          if(fl%tc_perpendicular) &
          qvec(ix^d,idims)=qvec(ix^d,idims)+kef(ix^d)*0.25d0*(qd(ix^d,1)+qd(ix^d,2)+qd(ix^d,3)+qd(ix^d,4))
       {end do\}
      end if
     }
     {^iftwod
      if(idims==1) then
       {do ix^db=ixamin^db,ixamax^db\}
          minq=min(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          maxq=max(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          if(qdd(ix^d)<minq) then
            qd(ix^d,1)=minq
          else if(qdd(ix^d)>maxq) then
            qd(ix^d,1)=maxq
          else
            qd(ix^d,1)=qdd(ix^d)
          end if
          if(qdd(ix1,ix2-1)<minq) then
            qd(ix^d,2)=minq
          else if(qdd(ix1,ix2-1)>maxq) then
            qd(ix^d,2)=maxq
          else
            qd(ix^d,2)=qdd(ix1,ix2-1)
          end if
          qvec(ix^d,idims)=kaf(ix^d)*0.5d0*(bc(ix1,ix2,idims)**2*qd(ix^d,1)+bc(ix1,ix2-1,idims)**2*qd(ix^d,2))
          if(fl%tc_perpendicular) &
          qvec(ix^d,idims)=qvec(ix^d,idims)+kef(ix^d)*0.5d0*(qd(ix^d,1)+qd(ix^d,2))
       {end do\}
      else
       {do ix^db=ixamin^db,ixamax^db\}
          minq=min(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          maxq=max(alpha*gradt(ix^d,idims),gradt(ix^d,idims)/alpha)
          if(qdd(ix^d)<minq) then
            qd(ix^d,1)=minq
          else if(qdd(ix^d)>maxq) then
            qd(ix^d,1)=maxq
          else
            qd(ix^d,1)=qdd(ix^d)
          end if
          if(qdd(ix1-1,ix2)<minq) then
            qd(ix^d,2)=minq
          else if(qdd(ix1-1,ix2)>maxq) then
            qd(ix^d,2)=maxq
          else
            qd(ix^d,2)=qdd(ix1-1,ix2)
          end if
          qvec(ix^d,idims)=kaf(ix^d)*0.5d0*(bc(ix1,ix2,idims)**2*qd(ix^d,1)+bc(ix1-1,ix2,idims)**2*qd(ix^d,2))
          if(fl%tc_perpendicular) &
          qvec(ix^d,idims)=qvec(ix^d,idims)+kef(ix^d)*0.5d0*(qd(ix^d,1)+qd(ix^d,2))
       {end do\}
      end if
     }
      ! calculate normal of magnetic field
      ixb^l=ixa^l+kr(idims,^d);
      bnorm(ixa^s)=0.5d0*(mf(ixa^s,idims)+mf(ixb^s,idims))
      ! limited transverse component, eq (17)
      ixbmin^d=ixamin^d;
      ixbmax^d=ixamax^d+kr(idims,^d);
      do idir=1,ndim
        if(idir==idims) cycle
        qdd(ixi^s)=slope_limiter(gradt(ixi^s,idir),ixi^l,ixb^l,idir,-1,fl%tc_slope_limiter)
        qdd(ixi^s)=slope_limiter(qdd,ixi^l,ixa^l,idims,1,fl%tc_slope_limiter)
        qvec(ixa^s,idims)=qvec(ixa^s,idims)+kaf(ixa^s)*bnorm(ixa^s)*bcf(ixa^s,idir)*qdd(ixa^s)
      end do
      if(fl%tc_saturate) then
        ! consider saturation (Cowie and Mckee 1977 ApJ, 211, 135: phi=1.1, Balbus and Mckee 1982 ApJ, 252, 529: phi=0.3)
        ! unsigned saturated TC flux = 5 phi rho c**3, c=sqrt(p/rho) is isothermal sound speed, phi=0.3
        ixb^l=ixa^l+kr(idims,^d);
        qdd(ixa^s)=0.75d0*(rho(ixa^s)+rho(ixb^s))*dsqrt(0.5d0*(ciso2(ixa^s)+ciso2(ixb^s)))**3*dabs(bnorm(ixa^s))
       {do ix^db=ixamin^db,ixamax^db\}
          if(dabs(qvec(ix^d,idims))>qdd(ix^d)) then
            qvec(ix^d,idims)=sign(1.d0,qvec(ix^d,idims))*qdd(ix^d)
          end if
       {end do\}
      end if
    end do
  end subroutine set_source_tc_mhd_geo
 
  function slope_limiter(f,ixI^L,ixO^L,idims,pm,tc_slope_limiter) result(lf)
    use mod_global_parameters
    integer, intent(in) :: ixi^l, ixo^l, idims, pm
    double precision, intent(in) :: f(ixi^s)
    double precision :: lf(ixi^s)
    integer, intent(in)  :: tc_slope_limiter
 
    double precision, parameter :: qsmall=1.d-12
    double precision :: signf(ixi^s)
    integer :: ixb^l
 
    ixb^l=ixo^l+pm*kr(idims,^d);
    signf(ixo^s)=sign(1.d0,f(ixo^s))
    select case(tc_slope_limiter)
     case(1)
       ! 'MC' monotonized central limiter Woodward and Collela limiter (eq.3.51h), a factor of 2 is pulled out
       lf(ixo^s)=two*signf(ixo^s)* &
            max(zero,min(dabs(f(ixo^s)),signf(ixo^s)*f(ixb^s),&
            signf(ixo^s)*quarter*(f(ixb^s)+f(ixo^s))))
     case(2)
       ! 'minmod' limiter
       lf(ixo^s)=signf(ixo^s)*max(0.d0,min(abs(f(ixo^s)),signf(ixo^s)*f(ixb^s)))
     case(3)
       ! 'superbee' Roe superbee limiter (eq.3.51i)
       lf(ixo^s)=signf(ixo^s)* &
            max(zero,min(two*dabs(f(ixo^s)),signf(ixo^s)*f(ixb^s)),&
            min(dabs(f(ixo^s)),two*signf(ixo^s)*f(ixb^s)))
     case(4)
       ! 'koren' Barry Koren Right variant
       lf(ixo^s)=signf(ixo^s)* &
            max(zero,min(two*dabs(f(ixo^s)),two*signf(ixo^s)*f(ixb^s),&
            (two*f(ixb^s)*signf(ixo^s)+dabs(f(ixo^s)))*third))
     case(5)
       ! van Leer limiter
       lf(ixo^s)=two*max(f(ixb^s)*f(ixo^s),zero)/(f(ixo^s)+f(ixb^s)+qsmall)
     case default
       call mpistop("Unknown slope limiter for thermal conduction")
    end select
  end function slope_limiter
 
  !> Get the explicit timestep for the TC (hd implementation)
  !> Note: also used in 1D MHD (or for neutrals in twofl)
  function get_tc_dt_hd(w,ixI^L,ixO^L,dx^D,x,fl)  result(dtnew)
    ! Check diffusion time limit dt < dx_i**2 / ((gamma-1)*tc_k_para/(rho R))
    use mod_global_parameters
 
    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)
    type(tc_fluid), intent(in) :: fl
    double precision :: dtnew
 
    double precision :: tmp(ixo^s),tmp2(ixo^s),te(ixi^s),rho(ixi^s),hfs(ixo^s),gradt(ixi^s)
    double precision :: ne(ixi^s), nh_arr(ixi^s), rf(ixi^s)
    double precision :: dtdiff_tcond,maxtmp2
    integer          :: idim
 
    call fl%get_temperature_from_conserved(w,x,ixi^l,ixi^l,te)
    call fl%get_rho(w,x,ixi^l,ixo^l,rho)
    call fl%get_ne_nH(ixi^l, ixo^l, w, x, ne, nh_arr)
    ! the diffusivity is kappa/(rho c_v) with c_v = R/(gamma-1): weight rho by R
    if(associated(fl%get_var_Rfactor)) then
      call fl%get_var_Rfactor(w, x, ixi^l, ixo^l, rf)
      rho(ixo^s)=rho(ixo^s)*rf(ixo^s)
    end if
 
    if(fl%tc_constant) then
      tmp(ixo^s)=fl%tc_k_para/rho(ixo^s)
    else
      if(fl%tc_saturate) then
        ! Kannan 2016 MN 458, 410
        ! l_mfpe = 3^1.5*kB^2/(4*sqrt(pi)*e^4*lnLambda) * T^2/n_e
        if(si_unit) then
          ! 5.730205638843984e27 = 3^1.5*kB_SI^2/(4*sqrt(pi)*e_SI^4*37)
          tmp2(ixo^s)=te(ixo^s)**2/ne(ixo^s)*5.730205638843984d27*unit_temperature**2/(unit_numberdensity*unit_length)
        else
          ! 7093.9239487765044 = 3^1.5*kB_cgs^2/(4*sqrt(pi)*e_cgs^4*37)
          tmp2(ixo^s)=te(ixo^s)**2/ne(ixo^s)*7093.9239487765044d0*unit_temperature**2/(unit_numberdensity*unit_length)
        end if
        hfs=0.d0
        do idim=1,ndim
          call gradient(te,ixi^l,ixo^l,idim,gradt)
          hfs(ixo^s)=hfs(ixo^s)+gradt(ixo^s)**2
        end do
        ! kappa=kappa_Spitzer/(1+4.2*l_mfpe/(T/|gradT|))
        tmp(ixo^s)=fl%tc_k_para*te(ixo^s)*te(ixo^s)*dsqrt(te(ixo^s))/(rho(ixo^s)*(1.d0+4.2d0*tmp2(ixo^s)*dsqrt(hfs(ixo^s))/te(ixo^s)))
      else
        tmp(ixo^s)=fl%tc_k_para*te(ixo^s)*te(ixo^s)*dsqrt(te(ixo^s))/rho(ixo^s)
      end if
    end if
 
    dtnew = bigdouble
    do idim=1,ndim
      ! approximate thermal conduction flux: tc_k_para/rho/dx**2
      maxtmp2=maxval(tmp(ixo^s)/(block%ds(ixo^s,idim)**2))
      ! dt< dx_idim**2/((gamma-1)*tc_k_para/rho)
      dtdiff_tcond=1.d0/(tc_gamma_1*maxtmp2+smalldouble)
      ! limit the time step
      dtnew=min(dtnew,dtdiff_tcond)
    end do
    dtnew=dtnew/dble(ndim)
  end function get_tc_dt_hd
 
  !> Te over the block (+ghosts) and ciso2 = p/rho for the saturation limit: the frozen-
  !> ionisation map within a conduction step with frozen_ion, else the full EoS with the
  !> optional negative-eint repair (ciso2 = Te where R = 1 or without saturation)
  subroutine tc_get_te_ciso2(ixI^L,w,x,fl,igrid,Te,ciso2)
    use mod_global_parameters
    use mod_supertimestepping, only: sts_app_id
    integer, intent(in) :: ixI^L, igrid
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision, intent(in) :: x(ixI^S,1:ndim)
    type(tc_fluid), intent(in) :: fl
    double precision, intent(out) :: Te(ixI^S), ciso2(ixI^S)
 
    double precision :: Tdummy(ixI^S)
    logical :: patched, want_ciso2
 
    if(fl%frozen_ion .and. associated(fl%get_temperature_ciso2_from_eint)) then
      ! The anchors are allocated once, by the first thread to arrive
      !$OMP CRITICAL (tc_anchor_alloc)
      if(.not. allocated(tc_anchor)) allocate(tc_anchor(max_blocks))
      !$OMP END CRITICAL (tc_anchor_alloc)
      ! the first stage of each conduction step anchors on the full EoS
      if(tc_anchor(igrid)%sts_id /= sts_app_id) then
        tc_anchor(igrid)%sts_id = sts_app_id
        call tc_frozen_anchor(ixi^l, w, x, fl, igrid)
      end if
      associate(a => tc_anchor(igrid))
      te(ixi^s) = a%T0(ixi^s) + (w(ixi^s,fl%e_) - a%e0(ixi^s))*a%cvinv(ixi^s)
      if(all(te(ixi^s) > 0.d0)) then
        ciso2(ixi^s) = a%R0(ixi^s)*te(ixi^s)
        return
      end if
      ! a non-positive frozen T: this block takes the full EoS for this stage
      end associate
    end if
 
    want_ciso2 = fl%tc_saturate .and. associated(fl%get_temperature_ciso2_from_eint)
    if(want_ciso2) then
      call fl%get_temperature_ciso2_from_eint(w, x, ixi^l, ixi^l, te, ciso2)
    else if(associated(fl%get_temperature_ciso2_from_eint)) then
      call fl%get_temperature_ciso2_from_eint(w, x, ixi^l, ixi^l, te)
    else
      call fl%get_temperature_from_eint(w, x, ixi^l, ixi^l, te)
    end if
    patched = .false.
    if(fl%tc_patch_eint) call tc_patch_negative_eint(w, x, ixi^l, ixi^l, te, fl%e_, fl, patched)
 
    if(want_ciso2) then
      ! ciso2 from the repaired energy, keeping the repaired Te
      if(patched) call fl%get_temperature_ciso2_from_eint(w, x, ixi^l, ixi^l, tdummy, ciso2)
    else if(fl%tc_saturate .and. fl%sat_R_variable .and. associated(fl%get_var_Rfactor)) then
      call fl%get_var_Rfactor(w, x, ixi^l, ixi^l, ciso2)
      ciso2(ixi^s) = ciso2(ixi^s)*te(ixi^s)
    else
      ciso2(ixi^s) = te(ixi^s)
    end if
  end subroutine tc_get_te_ciso2
 
  !> Anchor of the frozen-ionisation conduction step for block igrid, on the step-start
  !> state (internal energy in w(:,e_)): T and R = ciso2/T from the full EoS, the energy,
  !> and the inverse frozen-composition heat capacity. Within the step
  !> T = T0 + (e - e0)/(R0 rho/(gamma-1)), exact where the gas is fully ionised or neutral
  subroutine tc_frozen_anchor(ixI^L, w, x, fl, igrid)
    use mod_global_parameters
    integer, intent(in) :: ixI^L, igrid
    double precision, intent(in) :: w(ixI^S,1:nw), x(ixI^S,1:ndim)
    type(tc_fluid), intent(in) :: fl
    double precision :: rho(ixI^S), ciso2(ixI^S)
 
    associate(a => tc_anchor(igrid))
    if(.not. allocated(a%T0)) allocate(a%T0(ixi^s), a%e0(ixi^s), a%R0(ixi^s), a%cvinv(ixi^s))
    call fl%get_temperature_ciso2_from_eint(w, x, ixi^l, ixi^l, a%T0, ciso2)
    call fl%get_rho(w, x, ixi^l, ixi^l, rho)
    a%R0(ixi^s) = ciso2(ixi^s)/a%T0(ixi^s)
    a%e0(ixi^s) = w(ixi^s,fl%e_)
    a%cvinv(ixi^s) = 1.d0/(fl%inv_gamma_minus_1*a%R0(ixi^s)*rho(ixi^s))
    end associate
  end subroutine tc_frozen_anchor
 
  !> Background Te and rho for the equilibrium flux. ciso2 keeps R = ciso2/Te of the full state
  subroutine tc_equi_state(ixI^L, w, x, fl, Te, rho, ciso2)
    use mod_global_parameters
    integer, intent(in) :: ixI^L
    double precision, intent(in) :: w(ixI^S,1:nw), x(ixI^S,1:ndim)
    type(tc_fluid), intent(in) :: fl
    double precision, intent(inout) :: Te(ixI^S), rho(ixI^S), ciso2(ixI^S)
    ciso2(ixi^s)=ciso2(ixi^s)/te(ixi^s)
    call fl%get_temperature_equi(w, x, ixi^l, ixi^l, te)
    call fl%get_rho_equi(w, x, ixi^l, ixi^l, rho)
    ciso2(ixi^s)=ciso2(ixi^s)*te(ixi^s)
  end subroutine tc_equi_state
 
  !> Face fluxes for the HD STS source, uniform or general geometry
  subroutine tc_flux_hd(ixI^L,ixO^L,w,x,fl,qvec,rho,Te,ciso2)
    use mod_global_parameters
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: x(ixI^S,1:ndim), w(ixI^S,1:nw)
    type(tc_fluid), intent(in) :: fl
    double precision, intent(in) :: rho(ixI^S), Te(ixI^S), ciso2(ixI^S)
    double precision, intent(out) :: qvec(ixI^S,1:ndim)
 
    if(slab_uniform) then
      call set_source_tc_hd(ixi^l,ixo^l,w,x,fl,qvec,rho,te,ciso2)
    else
      call set_source_tc_hd_geo(ixi^l,ixo^l,w,x,fl,qvec,rho,te,ciso2)
    end if
  end subroutine tc_flux_hd
 
  !> Face fluxes for the MHD STS source, uniform or general geometry
  subroutine tc_flux_mhd(ixI^L,ixO^L,w,x,fl,qvec,rho,Te,ciso2,alpha)
    use mod_global_parameters
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: x(ixI^S,1:ndim), w(ixI^S,1:nw)
    type(tc_fluid), intent(in) :: fl
    double precision, intent(in) :: rho(ixI^S), Te(ixI^S), ciso2(ixI^S), alpha
    double precision, intent(out) :: qvec(ixI^S,1:ndim)
 
    if(slab_uniform) then
      call set_source_tc_mhd(ixi^l,ixo^l,w,x,fl,qvec,rho,te,ciso2,alpha)
    else
      call set_source_tc_mhd_geo(ixi^l,ixo^l,w,x,fl,qvec,rho,te,ciso2,alpha)
    end if
  end subroutine tc_flux_mhd
 
  subroutine sts_set_source_tc_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,fl)
    use mod_global_parameters
    use mod_fix_conserve
 
    integer, intent(in) :: ixi^l, ixo^l, igrid, nflux
    double precision, intent(in) ::  x(ixi^s,1:ndim)
    ! intent(inout) so tc_patch_negative_eint can repair w(:, ie) in place
    ! when tc_patch_eint = .true.; legacy path (patch off) leaves w untouched.
    double precision, intent(inout) ::  w(ixi^s,1:nw)
    double precision, intent(inout) ::  wres(ixi^s,1:nw)
    double precision, intent(in) :: my_dt
    logical, intent(in) :: fix_conserve_at_step
    type(tc_fluid), intent(in)    :: fl
 
    double precision :: te(ixi^s),rho(ixi^s),ciso2(ixi^s)
    double precision :: qvec(ixi^s,1:ndim),qd(ixi^s)
    double precision, allocatable, dimension(:^D&,:) :: qvec_equi
    double precision :: fluxall(ixi^s,1,1:ndim)
 
    double precision :: dxinv(ndim)
    integer :: idims,ix^l,ixb^l,ixa^l
 
    ix^l=ixo^l^ladd1;
 
    dxinv=1.d0/dxlevel
 
    call tc_get_te_ciso2(ixi^l,w,x,fl,igrid,te,ciso2)
    call fl%get_rho(w, x, ixi^l, ixi^l, rho)
    call tc_flux_hd(ixi^l,ixo^l,w,x,fl,qvec,rho,te,ciso2)
    if(fl%subtract_equi) then
      allocate(qvec_equi(ixi^s,1:ndim))
      call tc_equi_state(ixi^l, w, x, fl, te, rho, ciso2)
      call tc_flux_hd(ixi^l,ixo^l,w,x,fl,qvec_equi,rho,te,ciso2)
      do idims=1,ndim
        ! operate only on the per-idims face range that set_source_tc_hd fills
        ! (ixOmin-kr .. ixOmax); ix^S (=ixO+1) over-reaches into the
        ! uninitialised ixOmax+1 layer -> snan/Inf -> SIGFPE. Matches tc_mhd.
        ixamax^d=ixomax^d; ixamin^d=ixomin^d-kr(idims,^d);
        qvec(ixa^s,idims)=qvec(ixa^s,idims) - qvec_equi(ixa^s,idims)
      end do
      deallocate(qvec_equi)
    endif
 
    qd=0.d0
    if(slab_uniform) then
      do idims=1,ndim
        ixamax^d=ixomax^d; ixamin^d=ixomin^d-kr(idims,^d);
        qvec(ixa^s,idims)=dxinv(idims)*qvec(ixa^s,idims)
        ixb^l=ixo^l-kr(idims,^d);
        qd(ixo^s)=qd(ixo^s)+qvec(ixo^s,idims)-qvec(ixb^s,idims)
      end do
    else
      do idims=1,ndim
        ixamax^d=ixomax^d; ixamin^d=ixomin^d-kr(idims,^d);
        qvec(ixa^s,idims)=qvec(ixa^s,idims)*block%surfaceC(ixa^s,idims)
        ixb^l=ixo^l-kr(idims,^d);
        qd(ixo^s)=qd(ixo^s)+qvec(ixo^s,idims)-qvec(ixb^s,idims)
      end do
      qd(ixo^s)=qd(ixo^s)/block%dvolume(ixo^s)
    end if
 
    if(fix_conserve_at_step) then
      do idims=1,ndim
        ixamax^d=ixomax^d; ixamin^d=ixomin^d-kr(idims,^d);
        fluxall(ixa^s,1,idims)=my_dt*qvec(ixa^s,idims)
      end do
      call store_flux(igrid,fluxall,1,ndim,nflux)
    end if
 
    wres(ixo^s,fl%e_)=qd(ixo^s)
  end subroutine sts_set_source_tc_hd
 
  subroutine set_source_tc_hd(ixI^L,ixO^L,w,x,fl,qvec,rho,Te,ciso2)
    use mod_global_parameters
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) ::  x(ixI^S,1:ndim)
    double precision, intent(in) ::  w(ixI^S,1:nw)
    type(tc_fluid), intent(in)    :: fl
    double precision, intent(in) :: Te(ixI^S),rho(ixI^S),ciso2(ixI^S)
    double precision, intent(out) :: qvec(ixI^S,1:ndim)
    double precision :: gradT(ixI^S,1:ndim),ke(ixI^S),qd(ixI^S)
    integer :: idims,ix^D,ix^L,ixC^L,ixA^L,ixB^L
 
    ix^l=ixo^l^ladd1;
    ! ixC is cell-corner index
    ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
 
    ! calculate thermal conduction flux with symmetric scheme
    ! T gradient (central difference) at cell corners
    do idims=1,ndim
      ixbmin^d=ixmin^d;
      ixbmax^d=ixmax^d-kr(idims,^d);
      call gradientf(te,x,ixi^l,ixb^l,idims,ke)
     {^ifthreed
      if(idims==1) then
       {do ix^db=ixcmin^db,ixcmax^db\}
          qvec(ix^d,idims)=0.25d0*(ke(ix1,ix2,ix3)+ke(ix1,ix2+1,ix3)&
                         +ke(ix1,ix2,ix3+1)+ke(ix1,ix2+1,ix3+1))
       {end do\}
      else if(idims==2) then
       {do ix^db=ixcmin^db,ixcmax^db\}
          qvec(ix^d,idims)=0.25d0*(ke(ix1,ix2,ix3)+ke(ix1+1,ix2,ix3)&
                         +ke(ix1,ix2,ix3+1)+ke(ix1+1,ix2,ix3+1))
       {end do\}
      else
       {do ix^db=ixcmin^db,ixcmax^db\}
          qvec(ix^d,idims)=0.25d0*(ke(ix1,ix2,ix3)+ke(ix1+1,ix2,ix3)&
                         +ke(ix1,ix2+1,ix3)+ke(ix1+1,ix2+1,ix3))
       {end do\}
      end if
     }
     {^iftwod
      if(idims==1) then
       {do ix^db=ixcmin^db,ixcmax^db\}
          qvec(ix^d,idims)=0.5d0*(ke(ix1,ix2)+ke(ix1,ix2+1))
       {end do\}
      else
       {do ix^db=ixcmin^db,ixcmax^db\}
          qvec(ix^d,idims)=0.5d0*(ke(ix1,ix2)+ke(ix1+1,ix2))
       {end do\}
      end if
     }
     {^ifoned
     do ix1=ixcmin1,ixcmax1
       qvec(ix1,idims)=ke(ix1)
     end do
     }
    end do
    ! conductivity at cell center
    if(fl%tc_constant) then
      qd(ix^s)=fl%tc_k_para
    else
      if(phys_trac) then
       {do ix^db=ixmin^db,ixmax^db\}
          if(te(ix^d) < block%wextra(ix^d,fl%Tcoff_)) then
            qd(ix^d)=fl%tc_k_para*dsqrt(block%wextra(ix^d,fl%Tcoff_)**5)
          else
            qd(ix^d)=fl%tc_k_para*dsqrt(te(ix^d)**5)
          end if
       {end do\}
      else
        qd(ix^s)=fl%tc_k_para*dsqrt(te(ix^s)**5)
      end if
    end if
    ! conductivity Ke at cell corner
    ! cell corner conduction flux gradT
    {^ifthreed
    {do ix^db=ixcmin^db,ixcmax^db\}
       ke(ix^d)=0.125d0*(qd(ix1,ix2,ix3)+qd(ix1+1,ix2,ix3)&
                      +qd(ix1,ix2+1,ix3)+qd(ix1+1,ix2+1,ix3)&
                      +qd(ix1,ix2,ix3+1)+qd(ix1+1,ix2,ix3+1)&
                      +qd(ix1,ix2+1,ix3+1)+qd(ix1+1,ix2+1,ix3+1))
       gradt(ix^d,1)=ke(ix^d)*qvec(ix^d,1)
       gradt(ix^d,2)=ke(ix^d)*qvec(ix^d,2)
       gradt(ix^d,3)=ke(ix^d)*qvec(ix^d,3)
    {end do\}
    }
    {^iftwod
    {do ix^db=ixcmin^db,ixcmax^db\}
       ke(ix^d)=0.25d0*(qd(ix1,ix2)+qd(ix1+1,ix2)+qd(ix1,ix2+1)+qd(ix1+1,ix2+1))
       gradt(ix^d,1)=ke(ix^d)*qvec(ix^d,1)
       gradt(ix^d,2)=ke(ix^d)*qvec(ix^d,2)
    {end do\}
    }
    {^ifoned
     do ix1=ixcmin1,ixcmax1
       gradt(ix^d,1)=0.5d0*(qd(ix1)+qd(ix1+1))*qvec(ix^d,1)
     end do
    }
 
    ! conduction flux qvec at cell face
    do idims=1,ndim
      ixamax^d=ixomax^d; ixamin^d=ixomin^d-kr(idims,^d);
     {^ifthreed
      if(idims==1) then
       {do ix^db=ixamin^db,ixamax^db\}
          qvec(ix^d,idims)=0.25d0*(gradt(ix1,ix2,ix3,idims)+gradt(ix1,ix2-1,ix3,idims)&
                                +gradt(ix1,ix2,ix3-1,idims)+gradt(ix1,ix2-1,ix3-1,idims))
       {end do\}
      else if(idims==2) then
       {do ix^db=ixamin^db,ixamax^db\}
          qvec(ix^d,idims)=0.25d0*(gradt(ix1,ix2,ix3,idims)+gradt(ix1-1,ix2,ix3,idims)&
                                +gradt(ix1,ix2,ix3-1,idims)+gradt(ix1-1,ix2,ix3-1,idims))
       {end do\}
      else
       {do ix^db=ixamin^db,ixamax^db\}
          qvec(ix^d,idims)=0.25d0*(gradt(ix1,ix2,ix3,idims)+gradt(ix1,ix2-1,ix3,idims)&
                                +gradt(ix1-1,ix2,ix3,idims)+gradt(ix1-1,ix2-1,ix3,idims))
       {end do\}
      end if
     }
     {^iftwod
      if(idims==1) then
       {do ix^db=ixamin^db,ixamax^db\}
          qvec(ix^d,idims)=0.5d0*(gradt(ix1,ix2,idims)+gradt(ix1,ix2-1,idims))
       {end do\}
      else
       {do ix^db=ixamin^db,ixamax^db\}
          qvec(ix^d,idims)=0.5d0*(gradt(ix1,ix2,idims)+gradt(ix1-1,ix2,idims))
       {end do\}
      end if
     }
     {^ifoned
      do ix1=ixamin1,ixamax1
        qvec(ix1,idims)=gradt(ix1,idims)
      end do
     }
      if(fl%tc_saturate) then
        ! consider saturation (Cowie and Mckee 1977 ApJ, 211, 135: phi=1.1, Balbus and Mckee 1982 ApJ, 252, 529: phi=0.3)
        ! unsigned saturated TC flux = 5 phi rho c**3, c=sqrt(p/rho) is isothermal sound speed, phi=0.3
        ixb^l=ixa^l+kr(idims,^d);
        qd(ixa^s)=0.75d0*(rho(ixa^s)+rho(ixb^s))*dsqrt(0.5d0*(ciso2(ixa^s)+ciso2(ixb^s)))**3
       {do ix^db=ixamin^db,ixamax^db\}
          if(dabs(qvec(ix^d,idims))>qd(ix^d)) then
            qvec(ix^d,idims)=sign(1.d0,qvec(ix^d,idims))*qd(ix^d)
          end if
       {end do\}
      end if
    end do
  end subroutine set_source_tc_hd
  subroutine set_source_tc_hd_geo(ixI^L,ixO^L,w,x,fl,qvec,rho,Te,ciso2)
    use mod_global_parameters
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) ::  x(ixI^S,1:ndim)
    double precision, intent(in) ::  w(ixI^S,1:nw)
    type(tc_fluid), intent(in)    :: fl
    double precision, intent(in) :: Te(ixI^S),rho(ixI^S),ciso2(ixI^S)
    double precision, intent(out) :: qvec(ixI^S,1:ndim)
    double precision :: gradT(ixI^S,1:ndim),ke(ixI^S),qd(ixI^S)
    integer :: idims,ix^D,ix^L,ixC^L,ixA^L,ixB^L
 
    ix^l=ixo^l^ladd1;
    ! ixC is cell-corner index
    ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
 
    ! calculate thermal conduction flux with symmetric scheme
    ! T gradient from face centers to cell corners: surface weighted average
    do idims=1,ndim
      ixbmin^d=ixmin^d;
      ixbmax^d=ixmax^d-kr(idims,^d);
      call gradientf(te,x,ixi^l,ixb^l,idims,ke)
     {^ifthreed
      if(idims==1) then
       {do ix^db=ixcmin^db,ixcmax^db\}
          qvec(ix^d,1)=(ke(ix1,ix2,ix3)*block%surfaceC(ix1,ix2,ix3,1)&
                     +ke(ix1,ix2+1,ix3)*block%surfaceC(ix1,ix2+1,ix3,1)&
                     +ke(ix1,ix2,ix3+1)*block%surfaceC(ix1,ix2,ix3+1,1)&
                   +ke(ix1,ix2+1,ix3+1)*block%surfaceC(ix1,ix2+1,ix3+1,1))/&
            (block%surfaceC(ix1,ix2,ix3,1)+block%surfaceC(ix1,ix2+1,ix3,1)&
          +block%surfaceC(ix1,ix2,ix3+1,1)+block%surfaceC(ix1,ix2+1,ix3+1,1)+smalldouble**2)
       {end do\}
      else if(idims==2) then
       {do ix^db=ixcmin^db,ixcmax^db\}
          qvec(ix^d,2)=(ke(ix1,ix2,ix3)*block%surfaceC(ix1,ix2,ix3,2)&
                     +ke(ix1+1,ix2,ix3)*block%surfaceC(ix1+1,ix2,ix3,2)&
                     +ke(ix1,ix2,ix3+1)*block%surfaceC(ix1,ix2,ix3+1,2)&
                   +ke(ix1+1,ix2,ix3+1)*block%surfaceC(ix1+1,ix2,ix3+1,2))/&
            (block%surfaceC(ix1,ix2,ix3,2)+block%surfaceC(ix1+1,ix2,ix3,2)&
          +block%surfaceC(ix1,ix2,ix3+1,2)+block%surfaceC(ix1+1,ix2,ix3+1,2)+smalldouble**2)
          ! zero theta-normal surface area at pole axis
       {end do\}
      else
       {do ix^db=ixcmin^db,ixcmax^db\}
          qvec(ix^d,3)=(ke(ix1,ix2,ix3)*block%surfaceC(ix1,ix2,ix3,3)&
                     +ke(ix1+1,ix2,ix3)*block%surfaceC(ix1+1,ix2,ix3,3)&
                     +ke(ix1,ix2+1,ix3)*block%surfaceC(ix1,ix2+1,ix3,3)&
                   +ke(ix1+1,ix2+1,ix3)*block%surfaceC(ix1+1,ix2+1,ix3,3))/&
            (block%surfaceC(ix1,ix2,ix3,3)+block%surfaceC(ix1+1,ix2,ix3,3)&
          +block%surfaceC(ix1,ix2+1,ix3,3)+block%surfaceC(ix1+1,ix2+1,ix3,3))
       {end do\}
      end if
     }
     {^iftwod
      if(idims==1) then
       {do ix^db=ixcmin^db,ixcmax^db\}
          qvec(ix^d,1)=(ke(ix1,ix2)*block%surfaceC(ix1,ix2,1)+ke(ix1,ix2+1)*block%surfaceC(ix1,ix2+1,1))&
                      /(block%surfaceC(ix1,ix2,1)+block%surfaceC(ix1,ix2+1,1))
       {end do\}
      else
       {do ix^db=ixcmin^db,ixcmax^db\}
          qvec(ix^d,2)=(ke(ix1,ix2)*block%surfaceC(ix1,ix2,2)+ke(ix1+1,ix2)*block%surfaceC(ix1+1,ix2,2))&
                      /(block%surfaceC(ix1,ix2,2)+block%surfaceC(ix1+1,ix2,2))
       {end do\}
      end if
     }
     {^ifoned
     do ix1=ixcmin1,ixcmax1
       qvec(ix1,idims)=ke(ix1)
     end do
     }
    end do
    ! conductivity at cell center
    if(fl%tc_constant) then
      qd(ix^s)=fl%tc_k_para
    else
      if(phys_trac) then
       {do ix^db=ixmin^db,ixmax^db\}
          if(te(ix^d) < block%wextra(ix^d,fl%Tcoff_)) then
            qd(ix^d)=fl%tc_k_para*dsqrt(block%wextra(ix^d,fl%Tcoff_)**5)
          else
            qd(ix^d)=fl%tc_k_para*dsqrt(te(ix^d)**5)
          end if
       {end do\}
      else
        qd(ix^s)=fl%tc_k_para*dsqrt(te(ix^s)**5)
      end if
    end if
    ! conductivity Ke at cell corner
    ! cell corner conduction flux gradT
    {^ifthreed
    {do ix^db=ixcmin^db,ixcmax^db\}
       ke(ix^d)=(qd(ix1,ix2,ix3)*block%dvolume(ix1,ix2,ix3)+qd(ix1+1,ix2,ix3)*block%dvolume(ix1+1,ix2,ix3)&
              +qd(ix1,ix2+1,ix3)*block%dvolume(ix1,ix2+1,ix3)+qd(ix1+1,ix2+1,ix3)*block%dvolume(ix1+1,ix2+1,ix3)&
              +qd(ix1,ix2,ix3+1)*block%dvolume(ix1,ix2,ix3+1)+qd(ix1+1,ix2,ix3+1)*block%dvolume(ix1+1,ix2,ix3+1)&
            +qd(ix1,ix2+1,ix3+1)*block%dvolume(ix1,ix2+1,ix3+1)+qd(ix1+1,ix2+1,ix3+1)*block%dvolume(ix1+1,ix2+1,ix3+1))&
           /(block%dvolume(ix1,ix2,ix3)+block%dvolume(ix1+1,ix2,ix3)+block%dvolume(ix1,ix2+1,ix3)+block%dvolume(ix1+1,ix2+1,ix3)&
        +block%dvolume(ix1,ix2,ix3+1)+block%dvolume(ix1+1,ix2,ix3+1)+block%dvolume(ix1,ix2+1,ix3+1)+block%dvolume(ix1+1,ix2+1,ix3+1))
       gradt(ix^d,1)=ke(ix^d)*qvec(ix^d,1)
       gradt(ix^d,2)=ke(ix^d)*qvec(ix^d,2)
       gradt(ix^d,3)=ke(ix^d)*qvec(ix^d,3)
    {end do\}
    }
    {^iftwod
    {do ix^db=ixcmin^db,ixcmax^db\}
       ke(ix^d)=(qd(ix1,ix2)*block%dvolume(ix1,ix2)+qd(ix1+1,ix2)*block%dvolume(ix1+1,ix2)&
              +qd(ix1,ix2+1)*block%dvolume(ix1,ix2+1)+qd(ix1+1,ix2+1)*block%dvolume(ix1+1,ix2+1))&
             /(block%dvolume(ix1,ix2)+block%dvolume(ix1+1,ix2)+block%dvolume(ix1,ix2+1)+block%dvolume(ix1+1,ix2+1))
       gradt(ix^d,1)=ke(ix^d)*qvec(ix^d,1)
       gradt(ix^d,2)=ke(ix^d)*qvec(ix^d,2)
    {end do\}
    }
    {^ifoned
     do ix1=ixcmin1,ixcmax1
       gradt(ix^d,1)=(qd(ix1)*block%dvolume(ix1)+qd(ix1+1)*block%dvolume(ix1+1))/(block%dvolume(ix1)+block%dvolume(ix1+1))*qvec(ix^d,1)
     end do
    }
 
    ! conduction flux qvec at cell face
    do idims=1,ndim
      ixamax^d=ixomax^d; ixamin^d=ixomin^d-kr(idims,^d);
     {^ifthreed
      if(idims==1) then
       {do ix^db=ixamin^db,ixamax^db\}
          qvec(ix^d,idims)=0.25d0*(gradt(ix1,ix2,ix3,idims)+gradt(ix1,ix2-1,ix3,idims)&
                                +gradt(ix1,ix2,ix3-1,idims)+gradt(ix1,ix2-1,ix3-1,idims))
       {end do\}
      else if(idims==2) then
       {do ix^db=ixamin^db,ixamax^db\}
          qvec(ix^d,idims)=0.25d0*(gradt(ix1,ix2,ix3,idims)+gradt(ix1-1,ix2,ix3,idims)&
                                +gradt(ix1,ix2,ix3-1,idims)+gradt(ix1-1,ix2,ix3-1,idims))
       {end do\}
      else
       {do ix^db=ixamin^db,ixamax^db\}
          qvec(ix^d,idims)=0.25d0*(gradt(ix1,ix2,ix3,idims)+gradt(ix1,ix2-1,ix3,idims)&
                                +gradt(ix1-1,ix2,ix3,idims)+gradt(ix1-1,ix2-1,ix3,idims))
       {end do\}
      end if
     }
     {^iftwod
      if(idims==1) then
       {do ix^db=ixamin^db,ixamax^db\}
          qvec(ix^d,idims)=0.5d0*(gradt(ix1,ix2,idims)+gradt(ix1,ix2-1,idims))
       {end do\}
      else
       {do ix^db=ixamin^db,ixamax^db\}
          qvec(ix^d,idims)=0.5d0*(gradt(ix1,ix2,idims)+gradt(ix1-1,ix2,idims))
       {end do\}
      end if
     }
     {^ifoned
      do ix1=ixamin1,ixamax1
        qvec(ix1,idims)=gradt(ix1,idims)
      end do
     }
      if(fl%tc_saturate) then
        ! consider saturation (Cowie and Mckee 1977 ApJ, 211, 135: phi=1.1, Balbus and Mckee 1982 ApJ, 252, 529: phi=0.3)
        ! unsigned saturated TC flux = 5 phi rho c**3, c=sqrt(p/rho) is isothermal sound speed, phi=0.3
        ixb^l=ixa^l+kr(idims,^d);
        qd(ixa^s)=0.75d0*(rho(ixa^s)+rho(ixb^s))*dsqrt(0.5d0*(ciso2(ixa^s)+ciso2(ixb^s)))**3
       {do ix^db=ixamin^db,ixamax^db\}
          if(dabs(qvec(ix^d,idims))>qd(ix^d)) then
            qvec(ix^d,idims)=sign(1.d0,qvec(ix^d,idims))*qd(ix^d)
          end if
       {end do\}
      end if
    end do
  end subroutine set_source_tc_hd_geo
 
  !> Patch cells where e_int <= 0 by neighbor-averaging the temperature
  !> AND repairing the conserved internal energy w(:, ie) in place.
  !> Called after get_temperature_from_eint in the TC source routines when
  !> tc_patch_eint is .true.  During STS RKL2 Chebyshev substeps the
  !> polynomial can overshoot e_int to negative values at low-density
  !> coronal cells.  A negative T would produce NaN via sqrt(T) in the
  !> Spitzer conductivity.  This routine replaces those cells' T with
  !> the average of valid (e_int > 0) neighbors and writes the matching
  !> e_int back to w(ix, ie) via the EoS inverse helper eint_nH_from_T,
  !> so the next STS substep does not re-read negative w and amplify it.
  subroutine tc_patch_negative_eint(w, x, ixI^L, ixO^L, Te, ie, fl, patched)
    use mod_global_parameters
    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)
    double precision, intent(inout) :: te(ixi^s)
    type(tc_fluid), intent(in)      :: fl
    !> true when at least one cell was repaired
    logical, intent(out)            :: patched
 
    integer :: ix^d, count
    integer :: ipatch(ixi^s)
    double precision :: t_avg, t_use, rho_cell, nh_cell, log_nh, log_t
    double precision :: eint_new, small_e_local, log_t_min, t_floor
    double precision :: rho(ixi^s)
 
    ipatch(ixi^s) = 0
 
    ! Mark cells with negative e_int
    {do ix^db=ixomin^db,ixomax^db\}
        if (w(ix^d, ie) <= 0.0d0) ipatch(ix^d) = 1
    {end do\}
 
    ! Quick exit if no cells need patching
    patched = any(ipatch(ixo^s) /= 0)
    if (.not. patched) return
 
    ! Need rho per-cell to back out e_int from a prescribed T via the EoS.
    call fl%get_rho(w, x, ixi^l, ixi^l, rho)
 
    ! Floor on log_T for the (rho, T) inverse table, from the EoS via the port
    ! (eos_get_log_T_floor encapsulates the method->container choice; picking the
    ! wrong container would leave var2_min=0 and clobber cold cells).
    log_t_min = fl%log_T_floor
 
    ! Local minimum e_int (code units): equivalent to small_pressure / (gamma-1).
    small_e_local = small_pressure * fl%inv_gamma_minus_1
 
    ! Fallback T (code units) for cells that have no valid neighbour.
    t_floor = 1.0d0 / unit_temperature
 
    ! Replace marked cells with neighbor-averaged T, then back out e_int
    ! and floor w(:, ie) so subsequent STS substeps see a positive state.
    {do ix^db=ixomin^db+1,ixomax^db-1\}
        if (ipatch(ix^d) == 1) then
            t_avg = 0.0d0
            count = 0
            {^ifoned
            if (ipatch(ix1-1)==0) then; t_avg=t_avg+te(ix1-1); count=count+1; end if
            if (ipatch(ix1+1)==0) then; t_avg=t_avg+te(ix1+1); count=count+1; end if
            }
            {^iftwod
            if (ipatch(ix1-1,ix2)==0) then; t_avg=t_avg+te(ix1-1,ix2); count=count+1; end if
            if (ipatch(ix1+1,ix2)==0) then; t_avg=t_avg+te(ix1+1,ix2); count=count+1; end if
            if (ipatch(ix1,ix2-1)==0) then; t_avg=t_avg+te(ix1,ix2-1); count=count+1; end if
            if (ipatch(ix1,ix2+1)==0) then; t_avg=t_avg+te(ix1,ix2+1); count=count+1; end if
            }
            {^ifthreed
            if (ipatch(ix1-1,ix2,ix3)==0) then; t_avg=t_avg+te(ix1-1,ix2,ix3); count=count+1; end if
            if (ipatch(ix1+1,ix2,ix3)==0) then; t_avg=t_avg+te(ix1+1,ix2,ix3); count=count+1; end if
            if (ipatch(ix1,ix2-1,ix3)==0) then; t_avg=t_avg+te(ix1,ix2-1,ix3); count=count+1; end if
            if (ipatch(ix1,ix2+1,ix3)==0) then; t_avg=t_avg+te(ix1,ix2+1,ix3); count=count+1; end if
            if (ipatch(ix1,ix2,ix3-1)==0) then; t_avg=t_avg+te(ix1,ix2,ix3-1); count=count+1; end if
            if (ipatch(ix1,ix2,ix3+1)==0) then; t_avg=t_avg+te(ix1,ix2,ix3+1); count=count+1; end if
            }
            if (count > 0) then
                t_use = t_avg / dble(count)
            else
                t_use = t_floor
            end if
 
            ! Repair the conserved internal energy in place. Use the same
            ! EoS-aware inverse used by hd/mhd_from_prolong_LTE: e_int =
            ! nH * eint_nH_from_T(log_nH, log_T).  Floor with small_e.
            rho_cell = max(rho(ix^d), small_density)
            nh_cell  = rho_cell / fl%nH2rhoFactor
            log_nh   = dlog10(max(nh_cell, smalldouble))
            log_t    = dlog10(max(t_use, 10.0d0**log_t_min))
            eint_new = nh_cell * fl%eint_from_T(log_nh, log_t)
            w(ix^d, ie) = max(small_e_local, eint_new)
 
            ! Keep Te consistent with the floored w(:, ie): if small_e
            ! kicked in, the implied T is slightly different from T_use.
            ! Use the prescribed T_use here (Te is only consumed within
            ! this TC substep; the next call recomputes it from w).
            te(ix^d) = t_use
 
            write(*,*) ' WARNING: tc_patch_eint it=', it, &
                ' e_int=', w(ix^d, ie), &
                ' pe=', mype, ' x=', x(ix^d, 1:ndim)
            ipatch(ix^d) = 2  ! mark as handled
        end if
    {end do\}
 
    ! Cells at the block boundary that the interior loop above skipped
    ! (ixOmin/ixOmax edges).  Fall back to T = 1/unit_temperature and
    ! repair w(:, ie) consistently so STS does not re-read negative e.
    {do ix^db=ixomin^db,ixomax^db\}
        if (ipatch(ix^d) == 1) then
            t_use    = t_floor
            rho_cell = max(rho(ix^d), small_density)
            nh_cell  = rho_cell / fl%nH2rhoFactor
            log_nh   = dlog10(max(nh_cell, smalldouble))
            log_t    = dlog10(max(t_use, 10.0d0**log_t_min))
            eint_new = nh_cell * fl%eint_from_T(log_nh, log_t)
            w(ix^d, ie) = max(small_e_local, eint_new)
            te(ix^d) = t_use
            write(*,*) ' WARNING: tc_patch_eint it=', it, &
                ' e_int=', w(ix^d, ie), &
                ' pe=', mype, ' x=', x(ix^d, 1:ndim), ' (edge)'
        end if
    {end do\}
 
  end subroutine tc_patch_negative_eint
 
end module mod_thermal_conduction