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
  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 has_equi = .true. get_temperature_equi and get_rho_equi have to be set
    logical :: has_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.
    ! 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()
    procedure(get_var_subr), pointer,nopass :: get_temperature_from_conserved => null()
    procedure(get_var_subr), pointer,nopass :: get_temperature_equi => 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
 
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)
    !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:ndim),te(ixi^s),rho(ixi^s),gradt(ixi^s)
    double precision :: tmp(ixo^s),hfs(ixo^s),blocal(1:ndim)
    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\}
          ^d&blocal(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
         {^iftwod
          if(blocal(1)/=0.d0) then
            mf(ix^d,1)=sign(1.d0,blocal(1))/dsqrt(1.d0+(blocal(2)/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)
          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=ixomin^db,ixomax^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+(w(ix^d,iw_mag(2))/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)
          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
      ! 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)
 
    !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
        ! 3^1.5*kB^2/(4*sqrt(pi)*e^4)
        ! l_mfpe=3.d0**1.5d0*kB_cgs**2/(4.d0*sqrt(dpi)*e_cgs**4*37.d0)=7093.9239487765044d0
        if(si_unit) call mpistop("needs adjusting in get_tc_dt_mhd")
        tmp(ixo^s)=te(ixo^s)**2/rho(ixo^s)*7093.9239487765044d0*unit_temperature**2/(unit_numberdensity*unit_length)
        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*dsqrt(te(ixo^s)**5)/(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*dsqrt(te(ixo^s)**5)
      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)
    double precision, intent(in) ::   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)
    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 fl%get_temperature_from_eint(w, x, ixi^l, ixi^l, te)  !calculate Te in whole domain (+ghosts)
    call fl%get_rho(w, x, ixi^l, ixi^l, rho)  !calculate rho in whole domain (+ghosts)
    if(slab_uniform) then
      call set_source_tc_mhd(ixi^l,ixo^l,w,x,fl,qvec,rho,te,alpha)
    else
      call set_source_tc_mhd_geo(ixi^l,ixo^l,w,x,fl,qvec,rho,te,alpha)
    end if
    if(fl%has_equi) then
      allocate(qvec_equi(ixi^s,1:ndim))
      call fl%get_temperature_equi(w, x, ixi^l, ixi^l, te)  !calculate Te in whole domain (+ghosts)
      call fl%get_rho_equi(w, x, ixi^l, ixi^l, rho)  !calculate rho in whole domain (+ghosts)
      if(slab_uniform) then
        call set_source_tc_mhd(ixi^l,ixo^l,w,x,fl,qvec_equi,rho,te,alpha)
      else
        call set_source_tc_mhd_geo(ixi^l,ixo^l,w,x,fl,qvec_equi,rho,te,alpha)
      end if
      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,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)
    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(ndim)
    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\}
          ^d&blocal(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
         {^iftwod
          if(blocal(1)/=0.d0) then
            mf(ix^d,1)=sign(1.d0,blocal(1))/dsqrt(1.d0+(blocal(2)/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)
          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+(w(ix^d,iw_mag(2))/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)
          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/((^d&(w(ix^s,iw_mag(^d))+block%B0(ix^s,^d,0))**2+)*dsqrt(te(ix^s))+smalldouble)
        else
          qdd(ix^s)=fl%tc_k_perp*rho(ix^s)**2/((^d&w(ix^s,iw_mag(^d))**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
    if(fl%tc_slope_limiter==0) then
      ! calculate thermal conduction flux with symmetric scheme
      do idims=1,ndim
        !qdd corner values
        qdd=0.d0
        {do ix^db=0,1 \}
           if({ ix^d==0 .and. ^d==idims | .or.}) then
             ixbmin^d=ixcmin^d+ix^d;
             ixbmax^d=ixcmax^d+ix^d;
             qdd(ixc^s)=qdd(ixc^s)+gradt(ixb^s,idims)
           end if
        {end do\}
        ! temperature gradient at cell corner
        qvec(ixc^s,idims)=qdd(ixc^s)*0.5d0**(ndim-1)
      end do
      ! b grad T at cell corner
      qdd(ixc^s)=sum(qvec(ixc^s,1:ndim)*bc(ixc^s,1:ndim),dim=ndim+1)
      do idims=1,ndim
        ! TC flux at cell corner
        gradt(ixc^s,idims)=ka(ixc^s)*bc(ixc^s,idims)*qdd(ixc^s)
        if(fl%tc_perpendicular) gradt(ixc^s,idims)=gradt(ixc^s,idims)+ke(ixc^s)*qvec(ixc^s,idims)
      end do
      ! TC flux at cell face
      qvec=0.d0
      do idims=1,ndim
        ixb^l=ixo^l-kr(idims,^d);
        ixamax^d=ixomax^d; ixamin^d=ixbmin^d;
        {do ix^db=0,1 \}
           if({ ix^d==0 .and. ^d==idims | .or.}) then
             ixbmin^d=ixamin^d-ix^d;
             ixbmax^d=ixamax^d-ix^d;
             qvec(ixa^s,idims)=qvec(ixa^s,idims)+gradt(ixb^s,idims)
           end if
        {end do\}
        qvec(ixa^s,idims)=qvec(ixa^s,idims)*0.5d0**(ndim-1)
        if(fl%tc_saturate) then
          ! consider saturation (Cowie and Mckee 1977 ApJ, 211, 135)
          ! unsigned saturated TC flux = 5 phi rho c**3, c is isothermal sound speed
          bcf=0.d0
          {do ix^db=0,1 \}
             if({ ix^d==0 .and. ^d==idims | .or.}) then
               ixbmin^d=ixamin^d-ix^d;
               ixbmax^d=ixamax^d-ix^d;
               bcf(ixa^s,idims)=bcf(ixa^s,idims)+bc(ixb^s,idims)
             end if
          {end do\}
          ! averaged b at face centers
          bcf(ixa^s,idims)=bcf(ixa^s,idims)*0.5d0**(ndim-1)
          ixb^l=ixa^l+kr(idims,^d);
          qdd(ixa^s)=2.75d0*(rho(ixa^s)+rho(ixb^s))*dsqrt(0.5d0*(te(ixa^s)+te(ixb^s)))**3*dabs(bcf(ixa^s,idims))
         {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
    else
      ! 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)
          ! unsigned saturated TC flux = 5 phi rho c**3, c=sqrt(p/rho) is isothermal sound speed, phi=1.1
          ixb^l=ixa^l+kr(idims,^d);
          qdd(ixa^s)=2.75d0*(rho(ixa^s)+rho(ixb^s))*dsqrt(0.5d0*(te(ixa^s)+te(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 if
  end subroutine set_source_tc_mhd
 
  subroutine set_source_tc_mhd_geo(ixI^L,ixO^L,w,x,fl,qvec,rho,Te,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)
    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(ndim)
    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\}
          ^d&blocal(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
         {^iftwod
          if(blocal(1)/=0.d0) then
            mf(ix^d,1)=sign(1.d0,blocal(1))/dsqrt(1.d0+(blocal(2)/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)
          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+(w(ix^d,iw_mag(2))/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)
          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/((^d&(w(ix^s,iw_mag(^d))+block%B0(ix^s,^d,0))**2+)*dsqrt(te(ix^s))+smalldouble)
        else
          qdd(ix^s)=fl%tc_k_perp*rho(ix^s)**2/((^d&w(ix^s,iw_mag(^d))**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
    if(fl%tc_slope_limiter==0) then
      ! calculate thermal conduction flux with symmetric scheme
      do idims=1,ndim
        !qdd corner values
        qdd=0.d0
        {do ix^db=0,1 \}
           if({ ix^d==0 .and. ^d==idims | .or.}) then
             ixbmin^d=ixcmin^d+ix^d;
             ixbmax^d=ixcmax^d+ix^d;
             qdd(ixc^s)=qdd(ixc^s)+gradt(ixb^s,idims)
           end if
        {end do\}
        ! temperature gradient at cell corner
        qvec(ixc^s,idims)=qdd(ixc^s)*0.5d0**(ndim-1)
      end do
      ! b grad T at cell corner
      qdd(ixc^s)=sum(qvec(ixc^s,1:ndim)*bc(ixc^s,1:ndim),dim=ndim+1)
      do idims=1,ndim
        ! TC flux at cell corner
        gradt(ixc^s,idims)=ka(ixc^s)*bc(ixc^s,idims)*qdd(ixc^s)
        if(fl%tc_perpendicular) gradt(ixc^s,idims)=gradt(ixc^s,idims)+ke(ixc^s)*qvec(ixc^s,idims)
      end do
      ! TC flux at cell face
      qvec=0.d0
      do idims=1,ndim
        ixb^l=ixo^l-kr(idims,^d);
        ixamax^d=ixomax^d; ixamin^d=ixbmin^d;
        {do ix^db=0,1 \}
           if({ ix^d==0 .and. ^d==idims | .or.}) then
             ixbmin^d=ixamin^d-ix^d;
             ixbmax^d=ixamax^d-ix^d;
             qvec(ixa^s,idims)=qvec(ixa^s,idims)+gradt(ixb^s,idims)
           end if
        {end do\}
        qvec(ixa^s,idims)=qvec(ixa^s,idims)*0.5d0**(ndim-1)
        if(fl%tc_saturate) then
          ! consider saturation (Cowie and Mckee 1977 ApJ, 211, 135)
          ! unsigned saturated TC flux = 5 phi rho c**3, c is isothermal sound speed
          bcf=0.d0
          {do ix^db=0,1 \}
             if({ ix^d==0 .and. ^d==idims | .or.}) then
               ixbmin^d=ixamin^d-ix^d;
               ixbmax^d=ixamax^d-ix^d;
               bcf(ixa^s,idims)=bcf(ixa^s,idims)+bc(ixb^s,idims)
             end if
          {end do\}
          ! averaged b at face centers
          bcf(ixa^s,idims)=bcf(ixa^s,idims)*0.5d0**(ndim-1)
          ixb^l=ixa^l+kr(idims,^d);
          qdd(ixa^s)=2.75d0*(rho(ixa^s)+rho(ixb^s))*dsqrt(0.5d0*(te(ixa^s)+te(ixb^s)))**3*dabs(bcf(ixa^s,idims))
         {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
    else
      ! 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))
         {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)+1.d-300)
         {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)+1.d-300)
         {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)
          ! unsigned saturated TC flux = 5 phi rho c**3, c=sqrt(p/rho) is isothermal sound speed, phi=1.1
          ixb^l=ixa^l+kr(idims,^d);
          qdd(ixa^s)=2.75d0*(rho(ixa^s)+rho(ixb^s))*dsqrt(0.5d0*(te(ixa^s)+te(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 if
  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)
    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 :: 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)
 
    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
        ! 3^1.5*kB^2/(4*sqrt(pi)*e^4)
        ! l_mfpe=3.d0**1.5d0*kB_cgs**2/(4.d0*sqrt(dpi)*e_cgs**4*37.d0)=7093.9239487765044d0
        if(si_unit) call mpistop("needs adjusting in get_tc_dt_hd")
        tmp2(ixo^s)=te(ixo^s)**2/rho(ixo^s)*7093.9239487765044d0*unit_temperature**2/(unit_numberdensity*unit_length)
        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*dsqrt((te(ixo^s))**5)/(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*dsqrt((te(ixo^s))**5)/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
 
  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)
    double precision, intent(in) ::  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)
    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,ixa^l
 
    dxinv=1.d0/dxlevel
 
    !calculate Te in whole domain (+ghosts)
    call fl%get_temperature_from_eint(w, x, ixi^l, ixi^l, te)
    call fl%get_rho(w, x, ixi^l, ixi^l, rho)
    if(slab_uniform) then
      call set_source_tc_hd(ixi^l,ixo^l,w,x,fl,qvec,rho,te)
    else
      call set_source_tc_hd_geo(ixi^l,ixo^l,w,x,fl,qvec,rho,te)
    end if
    if(fl%has_equi) then
      allocate(qvec_equi(ixi^s,1:ndim))
      call fl%get_temperature_equi(w, x, ixi^l, ixi^l, te)  !calculate Te in whole domain (+ghosts)
      call fl%get_rho_equi(w, x, ixi^l, ixi^l, rho)  !calculate rho in whole domain (+ghosts)
      if(slab_uniform) then
        call set_source_tc_hd(ixi^l,ixo^l,w,x,fl,qvec_equi,rho,te)
      else
        call set_source_tc_hd_geo(ixi^l,ixo^l,w,x,fl,qvec_equi,rho,te)
      end if
      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)
    endif
 
    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_hd
 
  !> get HD thermal conduction source and flux in uniform Cartesian grid
  subroutine set_source_tc_hd(ixI^L,ixO^L,w,x,fl,qvec,rho,Te)
    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)
    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)
          ! unsigned saturated TC flux = 5 phi rho c**3, c=sqrt(p/rho) is isothermal sound speed, phi=1.1
          ixb^l=ixa^l+kr(idims,^d);
          qd(ixa^s)=2.75d0*(rho(ixa^s)+rho(ixb^s))*dsqrt(0.5d0*(te(ixa^s)+te(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
 
  !> get HD thermal conduction source and flux in non-uniform geometry grid
  subroutine set_source_tc_hd_geo(ixI^L,ixO^L,w,x,fl,qvec,rho,Te)
    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)
    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))
       {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)+1.d-300)
          ! 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)
          ! unsigned saturated TC flux = 5 phi rho c**3, c=sqrt(p/rho) is isothermal sound speed, phi=1.1
          ixb^l=ixa^l+kr(idims,^d);
          qd(ixa^s)=2.75d0*(rho(ixa^s)+rho(ixb^s))*dsqrt(0.5d0*(te(ixa^s)+te(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
 
end module mod_thermal_conduction