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)
!>
!> 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 :: tc_get_ffhd_params
  public :: get_tc_dt_mhd
  public :: get_tc_dt_hd
  public :: get_tc_dt_ffhd
  public :: sts_set_source_tc_mhd
  public :: sts_set_source_tc_hd
  public :: sts_set_source_tc_ffhd
 
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.
    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
    end if
 
  end subroutine tc_get_hd_params
 
  subroutine tc_get_ffhd_params(fl,read_ffhd_params)
    use mod_global_parameters
 
    interface
      subroutine read_ffhd_params(fl)
        use mod_global_parameters, only: unitpar,par_files
        import tc_fluid
        type(tc_fluid), intent(inout) :: fl
      end subroutine read_ffhd_params
    end interface
 
    type(tc_fluid), intent(inout) :: fl
 
    fl%tc_k_para=0.d0
    !> Read tc parameters from par file: ffHD case
    call read_ffhd_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
      endif
      if(mype .eq. 0) print*, "Spitzer ffHD par: ",fl%tc_k_para
    else
      fl%tc_constant=.true.
    endif
  end subroutine tc_get_ffhd_params
 
  !> Get the explicut timestep for the TC (mhd implementation)
  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)
    double precision :: dtdiff_tcond,maxtmp2
    integer          :: idims,ix^d
 
    !temperature
    call fl%get_temperature_from_conserved(w,x,ixi^l,ixi^l,te)
 
    ! B
   {do ix^db=ixomin^db,ixomax^db\}
      if(b0field) then
        ^d&mf({ix^d},^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
      else
        ^d&mf({ix^d},^d)=w({ix^d},iw_mag(^d))\
      end if
      ! Bsquared
      tmp(ix^d)=(^d&mf({ix^d},^d)**2+)
      ! B_i**2/B**2
      if(tmp(ix^d)/=0.d0) then
        ^d&mf({ix^d},^d)=mf({ix^d},^d)**2/tmp({ix^d})\
      else
        ^d&mf({ix^d},^d)=1.d0\
      end if
   {end do\}
 
    dtnew=bigdouble
    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
        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)*sqrt(mf(ixo^s,idims))
          else
            hfs(ixo^s)=hfs(ixo^s)+gradt(ixo^s)*sqrt(mf(ixo^s,idims))
          end if
        end do
        ! kappa=kappa_Spizer/(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_Spizer
        tmp(ixo^s)=fl%tc_k_para*dsqrt(te(ixo^s)**5)
      end if
    end if
 
    if(slab_uniform) then
      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)/(rho(ixo^s)*dxlevel(idims)))
        ! dt< dx_idim**2/((gamma-1)*tc_k_para_i/rho*B_i**2/B**2)
        dtdiff_tcond=dxlevel(idims)/(tc_gamma_1*maxtmp2+smalldouble)
        ! limit the time step
        dtnew=min(dtnew,dtdiff_tcond)
      end do
    else
      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)/(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
    end if
    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)
    call set_source_tc_mhd(ixi^l,ixo^l,w,x,fl,qvec,rho,te,alpha)
    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)
      call set_source_tc_mhd(ixi^l,ixo^l,w,x,fl,qvec_equi,rho,te,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,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))
    integer :: idims,idir,ix^D,ix^L,ixC^L,ixA^L,ixB^L
 
    ix^l=ixo^l^ladd1;
 
    ! T gradient at cell faces
    ! B vector
    if(b0field) then
     {do ix^db=ixmin^db,ixmax^db\}
        ^d&mf({ix^d},^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
        ! |B|
        bnorm(ix^d)=dsqrt(^d&mf({ix^d},^d)**2+)
        if(bnorm(ix^d)/=0.d0) then
          bnorm(ix^d)=1.d0/bnorm(ix^d)
        else
          bnorm(ix^d)=bigdouble
        end if
        ! b unit vector: magnetic field direction vector
        ^d&mf({ix^d},^d)=mf({ix^d},^d)*bnorm({ix^d})\
     {end do\}
    else
     {do ix^db=ixmin^db,ixmax^db\}
        ! |B|
        bnorm(ix^d)=dsqrt(^d&w({ix^d},iw_mag(^d))**2+)
        if(bnorm(ix^d)/=0.d0) then
          bnorm(ix^d)=1.d0/bnorm(ix^d)
        else
          bnorm(ix^d)=bigdouble
        end if
        ! b unit vector: magnetic field direction vector
        ^d&mf({ix^d},^d)=w({ix^d},iw_mag(^d))*bnorm({ix^d})\
     {end do\}
    end if
    ! 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 gradientc(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*sqrt(block%wextra(ix^d,fl%Tcoff_)**5)
          else
            qdd(ix^d)=fl%tc_k_para*sqrt(te(ix^d)**5)
          end if
       {end do\}
      else
        qdd(ix^s)=fl%tc_k_para*sqrt(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
        qdd(ix^s)=fl%tc_k_perp*rho(ix^s)**2*bnorm(ix^s)**2/dsqrt(te(ix^s))
       {^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
 
  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 :: 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' montonized 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' Roes 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 default
       call mpistop("Unknown slope limiter for thermal conduction")
    end select
  end function slope_limiter
 
  !> Calculate gradient of a scalar q at cell interfaces in direction idir
  subroutine gradientc(q,x,ixI^L,ixO^L,idir,gradq)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L, idir
    double precision, intent(in)    :: q(ixI^S),x(ixI^S,1:ndim)
    double precision, intent(inout) :: gradq(ixI^S)
    integer                         :: jxO^L
 
    jxo^l=ixo^l+kr(idir,^d);
    select case(coordinate)
    case(cartesian)
      gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/dxlevel(idir)
    case(cartesian_stretched)
      gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/(x(jxo^s,idir)-x(ixo^s,idir))
    case(spherical)
      select case(idir)
      case(1)
        gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/(x(jxo^s,1)-x(ixo^s,1))
        {^nooned
      case(2)
        gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/( (x(jxo^s,2)-x(ixo^s,2))*x(ixo^s,1) )
        }
        {^ifthreed
      case(3)
        gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/( (x(jxo^s,3)-x(ixo^s,3))*x(ixo^s,1)*dsin(x(ixo^s,2)) )
        }
      end select
    case(cylindrical)
      if(idir==phi_) then
        gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/((x(jxo^s,phi_)-x(ixo^s,phi_))*x(ixo^s,r_))
      else
        gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/(x(jxo^s,idir)-x(ixo^s,idir))
      end if
    case default
      call mpistop('Unknown geometry')
    end select
 
  end subroutine gradientc
 
  !> Get the explicit timestep for the TC (hd implementation)
  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_i/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_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
      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_Spizer/(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)*sqrt(hfs(ixo^s))/te(ixo^s)))
    else
      tmp(ixo^s)=fl%tc_k_para*dsqrt((te(ixo^s))**5)/rho(ixo^s)
    end if
    dtnew = bigdouble
 
    if(slab_uniform) then
      do idim=1,ndim
        ! approximate thermal conduction flux: tc_k_para_i/rho/dx
        tmp2(ixo^s)=tmp(ixo^s)/dxlevel(idim)
        maxtmp2=maxval(tmp2(ixo^s))
        ! dt< dx_idim**2/((gamma-1)*tc_k_para_i/rho)
        dtdiff_tcond=dxlevel(idim)/(tc_gamma_1*maxtmp2+smalldouble)
        ! limit the time step
        dtnew=min(dtnew,dtdiff_tcond)
      end do
    else
      do idim=1,ndim
        ! approximate thermal conduction flux: tc_k_para_i/rho/dx
        tmp2(ixo^s)=tmp(ixo^s)/block%ds(ixo^s,idim)
        maxtmp2=maxval(tmp2(ixo^s)/block%ds(ixo^s,idim))
        ! 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
    end if
    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, allocatable, dimension(:^D&,:,:) :: fluxall
 
    double precision :: dxinv(ndim)
    integer :: idims,ix^l,ixb^l
 
    ix^l=ixo^l^ladd1;
 
    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)
    call set_source_tc_hd(ixi^l,ixo^l,w,x,fl,qvec,rho,te)
    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)
      call set_source_tc_hd(ixi^l,ixo^l,w,x,fl,qvec_equi,rho,te)
      do idims=1,ndim
        qvec(ix^s,idims)=qvec(ix^s,idims) - qvec_equi(ix^s,idims)
      end do
      deallocate(qvec_equi)
    endif
 
    qd=0.d0
    if(slab_uniform) then
      do idims=1,ndim
        qvec(ix^s,idims)=dxinv(idims)*qvec(ix^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
        qvec(ix^s,idims)=qvec(ix^s,idims)*block%surfaceC(ix^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
      allocate(fluxall(ixi^s,1,1:ndim))
      fluxall(ixi^s,1,1:ndim)=my_dt*qvec(ixi^s,1:ndim)
      call store_flux(igrid,fluxall,1,ndim,nflux)
      deallocate(fluxall)
    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)
    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,ixD^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 gradientc(te,x,ixi^l,ixb^l,idims,ke)
      qd=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;
          qd(ixc^s)=qd(ixc^s)+ke(ixb^s)
        end if
     {end do\}
      ! temperature gradient at cell corner
      qvec(ixc^s,idims)=qd(ixc^s)*0.5d0**(ndim-1)
    end do
    ! conductivity at cell center
    if(phys_trac) then
      ! transition region adaptive conduction
      where(te(ix^s) < block%wextra(ix^s,fl%Tcoff_))
        qd(ix^s)=fl%tc_k_para*dsqrt(block%wextra(ix^s,fl%Tcoff_))**5
      else where
        qd(ix^s)=fl%tc_k_para*dsqrt(te(ix^s))**5
      end where
    else
      qd(ix^s)=fl%tc_k_para*dsqrt(te(ix^s))**5
    end if
    ke=0.d0
    {do ix^db=0,1\}
      ixbmin^d=ixcmin^d+ix^d;
      ixbmax^d=ixcmax^d+ix^d;
      ke(ixc^s)=ke(ixc^s)+qd(ixb^s)
    {end do\}
    ! cell corner conductivity
    ke(ixc^s)=0.5d0**ndim*ke(ixc^s)
    ! cell corner conduction flux
    do idims=1,ndim
      gradt(ixc^s,idims)=ke(ixc^s)*qvec(ixc^s,idims)
    end do
 
    if(fl%tc_saturate) then
      ! consider saturation with unsigned saturated TC flux = 5 phi rho c**3
      ! saturation flux at cell center
      qd(ix^s)=5.5d0*rho(ix^s)*dsqrt(te(ix^s)**3)
      !cell corner values of qd in ke
      ke=0.d0
      {do ix^db=0,1\}
        ixbmin^d=ixcmin^d+ix^d;
        ixbmax^d=ixcmax^d+ix^d;
        ke(ixc^s)=ke(ixc^s)+qd(ixb^s)
      {end do\}
      ! cell corner saturation flux
      ke(ixc^s)=0.5d0**ndim*ke(ixc^s)
      ! magnitude of cell corner conduction flux
      qd(ixc^s)=norm2(gradt(ixc^s,:),dim=ndim+1)
      {do ix^db=ixcmin^db,ixcmax^db\}
        if(qd(ix^d)>ke(ix^d)) then
          ke(ix^d)=ke(ix^d)/(qd(ix^d)+smalldouble)
          ^d&gradt({ix^d},^d)=ke({ix^d})*gradt({ix^d},^d)\
        end if
      {end do\}
    end if
 
    ! conduction 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)
    end do
  end subroutine set_source_tc_hd
 
  !> Get the explicut timestep for the TC (ffhd implementation)
  function get_tc_dt_ffhd(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 :: tmp2(ixo^s),tmp(ixo^s),hfs(ixo^s)
    double precision :: dtdiff_tcond,maxtmp2
    integer          :: idims
 
    !temperature
    call fl%get_temperature_from_conserved(w,x,ixi^l,ixi^l,te)
    mf(ixo^s,1:ndim)=(block%B0(ixo^s,1:ndim,0))**2
 
    dtnew=bigdouble
    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
        tmp2(ixo^s)=te(ixo^s)**2/rho(ixo^s)*7093.9239487765044d0*unit_temperature**2/(unit_numberdensity*unit_length)
        hfs=0.d0
        do idims=1,ndim
          call gradient(te,ixi^l,ixo^l,idims,gradt)
          hfs(ixo^s)=hfs(ixo^s)+gradt(ixo^s)*sqrt(mf(ixo^s,idims))
        end do
        ! kappa=kappa_Spizer/(1+4.2*l_mfpe/(T/|gradT.b|))
        tmp(ixo^s)=fl%tc_k_para*dsqrt(te(ixo^s)**5)/(1.d0+4.2d0*tmp2(ixo^s)*dabs(hfs(ixo^s))/te(ixo^s))
      else
        ! kappa=kappa_Spizer
        tmp(ixo^s)=fl%tc_k_para*dsqrt(te(ixo^s)**5)
      end if
    end if
 
    if(slab_uniform) then
      do idims=1,ndim
        ! approximate thermal conduction flux: tc_k_para_i/rho/dx*B_i**2/B**2
        tmp2(ixo^s)=tmp(ixo^s)*mf(ixo^s,idims)/(rho(ixo^s)*dxlevel(idims))
        maxtmp2=maxval(tmp2(ixo^s))
        ! dt< dx_idim**2/((gamma-1)*tc_k_para_i/rho*B_i**2/B**2)
        dtdiff_tcond=dxlevel(idims)/(tc_gamma_1*maxtmp2+smalldouble)
        ! limit the time step
        dtnew=min(dtnew,dtdiff_tcond)
      end do
    else
      do idims=1,ndim
        ! approximate thermal conduction flux: tc_k_para_i/rho/dx*B_i**2/B**2
        tmp2(ixo^s)=tmp(ixo^s)*mf(ixo^s,idims)/(rho(ixo^s)*block%ds(ixo^s,idims))
        maxtmp2=maxval(tmp2(ixo^s)/block%ds(ixo^s,idims))
        ! 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
    end if
    dtnew=dtnew/dble(ndim)
  end function get_tc_dt_ffhd
 
  !> anisotropic thermal conduction with slope limited symmetric scheme
  !> Sharma 2007 Journal of Computational Physics 227, 123
  subroutine sts_set_source_tc_ffhd(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(ixi^s)
    double precision :: rho(ixi^s),te(ixi^s)
    double precision :: qvec(ixi^s,1:ndim)
    double precision, allocatable, dimension(:^D&,:) :: qvec_equi
    double precision, allocatable, dimension(:^D&,:,:) :: fluxall
    double precision :: alpha,dxinv(ndim)
    integer :: idims,idir,ix^d,ix^l,ixc^l,ixa^l,ixb^l
 
    ! coefficient of limiting on normal component
    if(ndim<3) then
      alpha=0.75d0
    else
      alpha=0.85d0
    end if
    ix^l=ixo^l^ladd1;
    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)
    call set_source_tc_ffhd(ixi^l,ixo^l,w,x,fl,qvec,rho,te,alpha)
 
    qd=0.d0
    if(slab_uniform) then
      do idims=1,ndim
        qvec(ix^s,idims)=dxinv(idims)*qvec(ix^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
        qvec(ix^s,idims)=qvec(ix^s,idims)*block%surfaceC(ix^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
      allocate(fluxall(ixi^s,1,1:ndim))
      fluxall(ixi^s,1,1:ndim)=my_dt*qvec(ixi^s,1:ndim)
      call store_flux(igrid,fluxall,1,ndim,nflux)
      deallocate(fluxall)
    end if
    wres(ixo^s,fl%e_)=qd(ixo^s)
  end subroutine sts_set_source_tc_ffhd
 
  subroutine set_source_tc_ffhd(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)
    !! qd store the heat conduction energy changing rate
    double precision :: qd(ixI^S)
    double precision, dimension(ixI^S,1:ndim) :: mf,Bc,Bcf
    double precision, dimension(ixI^S,1:ndim) :: gradT
    double precision, dimension(ixI^S) :: ka,kaf,ke,kef,qdd,qe,Binv,minq,maxq,Bnorm
    double precision, allocatable, dimension(:^D&,:,:) :: fluxall
    integer, dimension(ndim) :: lowindex
    integer :: idims,idir,ix^D,ix^L,ixC^L,ixA^L,ixB^L,ixA^D,ixB^D
 
    ix^l=ixo^l^ladd1;
 
    ! T gradient at cell faces
    ! B vector
    mf(ixi^s,1:ndim)=block%B0(ixi^s,1:ndim,0)
    ! ixC is cell-corner index
    ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
    ! b unit vector at cell corner
    bc=0.d0
    {do ix^db=0,1\}
      ixamin^d=ixcmin^d+ix^d;
      ixamax^d=ixcmax^d+ix^d;
      bc(ixc^s,1:ndim)=bc(ixc^s,1:ndim)+mf(ixa^s,1:ndim)
    {end do\}
    bc(ixc^s,1:ndim)=bc(ixc^s,1:ndim)*0.5d0**ndim
    ! T gradient at cell faces
    do idims=1,ndim
      ixbmin^d=ixmin^d;
      ixbmax^d=ixmax^d-kr(idims,^d);
      call gradientc(te,x,ixi^l,ixb^l,idims,minq)
      gradt(ixb^s,idims)=minq(ixb^s)
    end do
    if(fl%tc_constant) then
      ka(ixc^s)=fl%tc_k_para
    else
      ! conductivity at cell center
      if(phys_trac) then
        minq(ix^s)=te(ix^s)
       {do ix^db=ixmin^db,ixmax^db\}
          if(minq(ix^d) < block%wextra(ix^d,fl%Tcoff_)) then
            minq(ix^d)=block%wextra(ix^d,fl%Tcoff_)
          end if
       {end do\}
        minq(ix^s)=fl%tc_k_para*sqrt(minq(ix^s)**5)
      else
        minq(ix^s)=fl%tc_k_para*sqrt(te(ix^s)**5)
      end if
      ka=0.d0
      {do ix^db=0,1\}
        ixbmin^d=ixcmin^d+ix^d;
        ixbmax^d=ixcmax^d+ix^d;
        ka(ixc^s)=ka(ixc^s)+minq(ixb^s)
      {end do\}
      ! cell corner conductivity
      ka(ixc^s)=0.5d0**ndim*ka(ixc^s)
    end if
    if(fl%tc_slope_limiter==0) then
      ! calculate thermal conduction flux with symmetric scheme
      do idims=1,ndim
        !qd corner values
        qd=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;
             qd(ixc^s)=qd(ixc^s)+gradt(ixb^s,idims)
           end if
        {end do\}
        ! temperature gradient at cell corner
        qvec(ixc^s,idims)=qd(ixc^s)*0.5d0**(ndim-1)
      end do
      ! b grad T at cell corner
      qd(ixc^s)=0.d0
      do idims=1,ndim
       qd(ixc^s)=qvec(ixc^s,idims)*bc(ixc^s,idims)+qd(ixc^s)
      end do
      do idims=1,ndim
        ! TC flux at cell corner
        gradt(ixc^s,idims)=ka(ixc^s)*bc(ixc^s,idims)*qd(ixc^s)
      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);
          qd(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))>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
    else
      ! calculate thermal conduction flux with slope-limited symmetric scheme
      qvec=0.d0
      do idims=1,ndim
        ixb^l=ixo^l-kr(idims,^d);
        ixamax^d=ixomax^d; ixamin^d=ixbmin^d;
        ! 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))
        bcf=0.d0
        kaf=0.d0
        kef=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,1:ndim)=bcf(ixa^s,1:ndim)+bc(ixb^s,1:ndim)
             kaf(ixa^s)=kaf(ixa^s)+ka(ixb^s)
           end if
        {end do\}
        ! averaged b at face centers
        bcf(ixa^s,1:ndim)=bcf(ixa^s,1:ndim)*0.5d0**(ndim-1)
        ! averaged thermal conductivity at face centers
        kaf(ixa^s)=kaf(ixa^s)*0.5d0**(ndim-1)
        ! limited normal component
        minq(ixa^s)=min(alpha*gradt(ixa^s,idims),gradt(ixa^s,idims)/alpha)
        maxq(ixa^s)=max(alpha*gradt(ixa^s,idims),gradt(ixa^s,idims)/alpha)
        ! eq (19)
        !corner values of gradT
        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
        qdd(ixc^s)=qdd(ixc^s)*0.5d0**(ndim-1)
        ! eq (21)
        qe=0.d0
        {do ix^db=0,1 \}
           qd(ixc^s)=qdd(ixc^s)
           if({ ix^d==0 .and. ^d==idims | .or.}) then
             ixbmin^d=ixamin^d-ix^d;
             ixbmax^d=ixamax^d-ix^d;
            {do ixa^db=ixamin^db,ixamax^db
               ixb^db=ixa^db-ix^db\}
               if(qd(ixb^d)<=minq(ixa^d)) then
                 qd(ixb^d)=minq(ixa^d)
               else if(qd(ixb^d)>=maxq(ixa^d)) then
                 qd(ixb^d)=maxq(ixa^d)
               end if 
            {end do\}
             qvec(ixa^s,idims)=qvec(ixa^s,idims)+bc(ixb^s,idims)**2*qd(ixb^s)
           end if
        {end do\}
        qvec(ixa^s,idims)=kaf(ixa^s)*qvec(ixa^s,idims)*0.5d0**(ndim-1)
        ! limited transverse component, eq (17)
        ixbmin^d=ixamin^d;
        ixbmax^d=ixamax^d+kr(idims,^d);
        do idir=1,ndim
          if(idir==idims) cycle
          qd(ixi^s)=slope_limiter(gradt(ixi^s,idir),ixi^l,ixb^l,idir,-1,fl%tc_slope_limiter)
          qd(ixi^s)=slope_limiter(qd,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)*qd(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);
          qd(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))>qd(ix^d)) then
              !write(*,*) 'it',it,qvec(ix^D,idims),qd(ix^D),' TC saturated at ',&
              !x(ix^D,:),' rho',rho(ix^D),' Te',Te(ix^D)
              qvec(ix^d,idims)=sign(1.d0,qvec(ix^d,idims))*qd(ix^d)
            end if
         {end do\}
        end if
      end do
    end if
  end subroutine set_source_tc_ffhd
end module mod_thermal_conduction