mod_twofl_phys.t

Go to the documentation of this file.

!> Magneto-hydrodynamics module
module mod_twofl_phys
 
#include "amrvac.h"
 
  use mod_physics
  use mod_global_parameters, only: std_len
  use mod_thermal_conduction, only: tc_fluid
  use mod_radiative_cooling, only: rc_fluid
  use mod_thermal_emission, only: te_fluid
  use mod_functions_bfield, only: get_divb,mag
  use mod_comm_lib, only: mpistop
 
  implicit none
  private
  !! E_c = E_kin + E_mag + E_int
  !! E_n = E_kin + E_int
  integer, public, parameter              :: eq_energy_tot=2
  !! E_c = E_int
  !! E_n = E_int
  integer, public, parameter              :: eq_energy_int=1
  !! E_n, E_c are calculated from density as c_adiab rho^gamma
  !! No energy equation => no variable assigned for it
  integer, public, parameter              :: eq_energy_none=0
  !! E_c = E_kin + E_int
  !! E_n = E_kin + E_int
  integer, public, parameter              :: eq_energy_ki=3
 
  integer, public, protected              :: twofl_eq_energy = eq_energy_tot
 
  !> Whether hyperdiffusivity is used
  logical, public, protected              :: twofl_hyperdiffusivity = .false.
  logical, public, protected              :: twofl_dump_hyperdiffusivity_coef = .false.
  double precision, public, protected, allocatable :: c_shk(:)
  double precision, public, protected, allocatable :: c_hyp(:)
 
  !> Whether thermal conduction is used
  logical, public, protected              :: twofl_thermal_conduction_c = .false.
  !> type of TC used: 1: adapted module (mhd implementation), 2: adapted module (hd implementation)
  integer, parameter, private             :: mhd_tc =1
  integer, parameter, private             :: hd_tc =2
  integer, protected                      :: use_twofl_tc_c = mhd_tc
 
  !> Whether radiative cooling is added
  logical, public, protected              :: twofl_radiative_cooling_c = .false.
  type(rc_fluid), public, allocatable :: rc_fl_c
 
  !> Whether viscosity is added
  logical, public, protected              :: twofl_viscosity = .false.
 
  !> Whether gravity is added: common flag for charges and neutrals
  logical, public, protected              :: twofl_gravity = .false.
 
  !> whether dump full variables (when splitting is used) in a separate dat file
  logical, public, protected              :: twofl_dump_full_vars = .false.
 
  !> Whether Hall-MHD is used
  logical, public, protected              :: twofl_hall = .false.
 
  type(tc_fluid), public, allocatable :: tc_fl_c
  type(te_fluid), public, allocatable :: te_fl_c
 
  type(tc_fluid), allocatable :: tc_fl_n
  logical, public, protected              :: twofl_thermal_conduction_n = .false.
  logical, public, protected              :: twofl_radiative_cooling_n = .false.
  type(rc_fluid), allocatable :: rc_fl_n
 
  !> Whether TRAC method is used
  logical, public, protected              :: twofl_trac = .false.
 
  !> Whether GLM-MHD is used
  logical, public, protected              :: twofl_glm = .false.
 
  !> Which TRAC method is used            
  integer, public, protected              :: twofl_trac_type=1
 
  !> Height of the mask used in the TRAC method
  double precision, public, protected              :: twofl_trac_mask = 0.d0
 
  !> Whether divB cleaning sources are added splitting from fluid solver
  logical, public, protected              :: source_split_divb = .false.
 
  !> GLM-MHD parameter: ratio of the diffusive and advective time scales for div b
  !> taking values within [0, 1]
  double precision, public                :: twofl_glm_alpha = 0.5d0
 
  !> MHD fourth order
  logical, public, protected              :: twofl_4th_order = .false.
 
  !> Index of the density (in the w array)
  integer, public             :: rho_c_
 
  !> Indices of the momentum density
  integer, allocatable, public :: mom_c(:)
 
  !> Index of the energy density (-1 if not present)
  integer, public             :: e_c_=-1
 
  !> Index of the cutoff temperature for the TRAC method
  integer, public              :: tcoff_c_
  integer, public              :: tweight_c_
 
  !> Indices of the GLM psi
  integer, public, protected :: psi_
 
  !> equi vars flags
  logical, public :: has_equi_rho_c0 = .false.  
  logical, public :: has_equi_pe_c0 = .false.  
 
  !> equi vars indices in the state%equi_vars array
  integer, public :: equi_rho_c0_ = -1
  integer, public :: equi_pe_c0_ = -1
  logical, public                         :: twofl_equi_thermal_c = .false.
 
  logical, public                         :: twofl_equi_thermal = .false.
  !neutrals:
 
  integer, public              :: rho_n_
  integer, allocatable, public :: mom_n(:)
  integer, public              :: e_n_
  integer, public              :: tcoff_n_
  integer, public              :: tweight_n_
  logical, public :: has_equi_rho_n0 = .false. 
  logical, public :: has_equi_pe_n0 = .false.  
  integer, public :: equi_rho_n0_ = -1
  integer, public :: equi_pe_n0_ = -1
 
  ! related to collisions:
  !> collisional alpha
  double precision, public                :: twofl_alpha_coll = 0d0
  logical, public                         :: twofl_alpha_coll_constant = .true.
  !> whether include thermal exchange collisional terms
  logical, public                         :: twofl_coll_inc_te = .true.
  !> whether include ionization/recombination inelastic collisional terms
  logical, public                         :: twofl_coll_inc_ionrec = .false.
  logical, public                         :: twofl_equi_thermal_n = .false.
  double precision, public                :: dtcollpar = -1d0 !negative value does not impose restriction on the timestep
  !> whether dump collisional terms in a separte dat file
  logical, public, protected              :: twofl_dump_coll_terms = .false.
 
  ! TODO Helium abundance not used, radiative cooling init uses it
  ! not in parameters list anymore 
  double precision, public, protected  :: he_abundance = 0d0
  ! two fluid is only H plasma
  double precision, public, protected  :: rc = 2d0
  double precision, public, protected  :: rn = 1d0
 
  !> The adiabatic index
  double precision, public                :: twofl_gamma = 5.d0/3.0d0
 
  !> The adiabatic constant
  double precision, public                :: twofl_adiab = 1.0d0
 
  !> The MHD resistivity
  double precision, public                :: twofl_eta = 0.0d0
 
  !> The MHD hyper-resistivity
  double precision, public                :: twofl_eta_hyper = 0.0d0
 
  !> The MHD Hall coefficient
  double precision, public                :: twofl_etah = 0.0d0
 
  !> The small_est allowed energy
  double precision, protected             :: small_e
 
  !> Method type to clean divergence of B
  character(len=std_len), public, protected :: typedivbfix  = 'linde'
 
  !> Method type of constrained transport
  character(len=std_len), public, protected :: type_ct  = 'uct_contact'
 
  !> Whether divB is computed with a fourth order approximation
  logical, public, protected :: twofl_divb_4thorder = .false.
 
  !> Method type in a integer for good performance
  integer :: type_divb
 
  !> Coefficient of diffusive divB cleaning
  double precision :: divbdiff     = 0.8d0
 
  !> Update all equations due to divB cleaning
  character(len=std_len) ::    typedivbdiff = 'all'
 
  !> clean initial divB
  logical, public :: clean_initial_divb     = .false.
 
  !> Add divB wave in Roe solver
  logical, public :: divbwave     = .true.
 
  !> To control divB=0 fix for boundary
  logical, public, protected :: boundary_divbfix(2*^nd)=.true.
 
  !> To skip * layer of ghost cells during divB=0 fix for boundary
  integer, public, protected :: boundary_divbfix_skip(2*^nd)=0
 
  !> B0 field is force-free
  logical, public, protected :: b0field_forcefree=.true.
 
  logical :: twofl_cbounds_species = .true.
 
  !> added from modules: gravity
  !> source split or not
  logical :: grav_split= .false.
 
  !> gamma minus one and its inverse
  double precision :: gamma_1, inv_gamma_1
 
  ! DivB cleaning methods
  integer, parameter :: divb_none          = 0
  integer, parameter :: divb_multigrid     = -1
  integer, parameter :: divb_glm           = 1
  integer, parameter :: divb_powel         = 2
  integer, parameter :: divb_janhunen      = 3
  integer, parameter :: divb_linde         = 4
  integer, parameter :: divb_lindejanhunen = 5
  integer, parameter :: divb_lindepowel    = 6
  integer, parameter :: divb_lindeglm      = 7
  integer, parameter :: divb_ct            = 8
 
  ! Public methods
  public :: twofl_phys_init
  public :: twofl_to_conserved
  public :: twofl_to_primitive
  public :: get_divb
  public :: get_rhoc_tot
  public :: twofl_get_v_c_idim
  ! TODO needed for the roe, see if can be used for n
  public :: twofl_get_csound2_c_from_conserved
  public :: get_rhon_tot
  public :: get_alpha_coll
  public :: get_gamma_ion_rec
  public :: twofl_get_v_n_idim
  public :: get_current
  public :: twofl_get_pthermal_c
  public :: twofl_get_pthermal_n
  public :: twofl_face_to_center
  public :: get_normalized_divb
  public :: b_from_vector_potential
  public :: usr_mask_gamma_ion_rec
  public :: usr_mask_alpha
 
  {^nooned
  public :: twofl_clean_divb_multigrid
  }
 
  abstract interface
 
    subroutine implicit_mult_factor_subroutine(ixI^L, ixO^L, step_dt, JJ, res)
      integer, intent(in)                :: ixi^l, ixo^l
      double precision, intent(in) :: step_dt
      double precision, intent(in) :: jj(ixi^s)
      double precision, intent(out) :: res(ixi^s)
 
    end subroutine implicit_mult_factor_subroutine
 
    subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
      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)
      double precision, intent(inout) :: res(ixi^s)
    end subroutine mask_subroutine
 
    subroutine mask_subroutine2(ixI^L,ixO^L,w,x,res1, res2)
      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)
      double precision, intent(inout) :: res1(ixI^S),res2(ixI^S)
    end subroutine mask_subroutine2
 
  end interface
 
   procedure(implicit_mult_factor_subroutine), pointer :: calc_mult_factor => null()
   integer, protected ::  twofl_implicit_calc_mult_method = 1
  procedure(mask_subroutine), pointer  :: usr_mask_alpha => null()
  procedure(mask_subroutine2), pointer  :: usr_mask_gamma_ion_rec => null()
 
contains
 
  !> Read this module"s parameters from a file
  subroutine twofl_read_params(files)
    use mod_global_parameters
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /twofl_list/ twofl_eq_energy, twofl_gamma, twofl_adiab,&
      twofl_eta, twofl_eta_hyper, twofl_etah, twofl_glm_alpha,& 
      twofl_thermal_conduction_c, use_twofl_tc_c, twofl_radiative_cooling_c, twofl_hall, twofl_gravity,&
      twofl_viscosity, twofl_4th_order, typedivbfix, source_split_divb, divbdiff,&
      typedivbdiff, type_ct, divbwave, si_unit, b0field,&
      b0field_forcefree, bdip, bquad, boct, busr,twofl_equi_thermal_c,twofl_equi_thermal,&
      twofl_dump_full_vars, has_equi_rho_c0, has_equi_pe_c0, twofl_hyperdiffusivity,twofl_dump_hyperdiffusivity_coef,&
      has_equi_pe_n0, has_equi_rho_n0, twofl_thermal_conduction_n, twofl_radiative_cooling_n,  &
      twofl_alpha_coll,twofl_alpha_coll_constant,&
      twofl_coll_inc_te, twofl_coll_inc_ionrec,&
      twofl_equi_thermal_n,dtcollpar,&
      twofl_dump_coll_terms,twofl_implicit_calc_mult_method,&
      boundary_divbfix, boundary_divbfix_skip, twofl_divb_4thorder, &
      clean_initial_divb,  &
      twofl_trac, twofl_trac_type, twofl_trac_mask,twofl_cbounds_species 
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status="old")
       read(unitpar, twofl_list, end=111)
111    close(unitpar)
    end do
 
  end subroutine twofl_read_params
 
  subroutine twofl_init_hyper(files)
    use mod_global_parameters
    use mod_hyperdiffusivity, only: hyperdiffusivity_init
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /hyperdiffusivity_list/ c_shk, c_hyp
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status="old")
       read(unitpar, hyperdiffusivity_list, end=113)
113    close(unitpar)
    end do
    
    call hyperdiffusivity_init()
 
    !!DEBUG
    if(mype .eq. 0) then
      print*, "Using Hyperdiffusivity"
      print*, "C_SHK ", c_shk(:)
      print*, "C_HYP ", c_hyp(:)
    endif
 
  end subroutine twofl_init_hyper
 
  !> Write this module's parameters to a snapsoht
  subroutine twofl_write_info(fh)
    use mod_global_parameters
    integer, intent(in)                 :: fh
    integer, parameter                  :: n_par = 1
    double precision                    :: values(n_par)
    character(len=name_len)             :: names(n_par)
    integer, dimension(MPI_STATUS_SIZE) :: st
    integer                             :: er
 
    call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
 
    names(1) = "gamma"
    values(1) = twofl_gamma
    call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
    call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
  end subroutine twofl_write_info
 
  subroutine twofl_phys_init()
    use mod_global_parameters
    use mod_thermal_conduction
    use mod_radiative_cooling
    use mod_viscosity, only: viscosity_init
    !use mod_gravity, only: gravity_init
    use mod_supertimestepping, only: sts_init, add_sts_method,&
            set_conversion_methods_to_head, set_error_handling_to_head
    {^nooned
    use mod_multigrid_coupling
    }
    integer :: itr, idir
 
    call twofl_read_params(par_files)
    physics_type = "twofl"
    if (twofl_cbounds_species) then
      number_species = 2
    endif
    phys_energy=.true.
  !> Solve total energy equation or not
  ! for the two fluid the true value means 
  ! E_charges = E_mag + E_kin_charges + E_int_charges
  ! E_neutrals =  E_kin_neutrals + E_int_neutrals
    phys_total_energy=.false.
 
  !> Solve internal energy instead of total energy
  ! for the two fluid the true value means 
  ! E_charges = E_int_charges
  ! E_neutrals = E_int_neutrals
    phys_internal_e=.false.
 
  ! For the two fluid phys_energy=.true. and phys_internal_e=.false. and phys_total_energy = .false. means
  ! E_charges = E_kin_charges + E_int_charges
  ! E_neutrals =  E_kin_neutrals + E_int_neutrals
    phys_gamma = twofl_gamma
 
    if(twofl_eq_energy == eq_energy_int) then
      phys_internal_e = .true.
    elseif(twofl_eq_energy == eq_energy_tot) then
      phys_total_energy = .true.
    elseif(twofl_eq_energy == eq_energy_none) then
      phys_energy = .false.
    endif
 
    phys_trac=twofl_trac
    phys_trac_type=twofl_trac_type
 
    if(.not. phys_energy) then
      if(twofl_thermal_conduction_n) then
        twofl_thermal_conduction_n=.false.
        if(mype==0) write(*,*) 'WARNING: set twofl_thermal_conduction_n=F when twofl_energy=F'
      end if
      if(twofl_radiative_cooling_n) then
        twofl_radiative_cooling_n=.false.
        if(mype==0) write(*,*) 'WARNING: set twofl_radiative_cooling_n=F when twofl_energy=F'
      end if
      if(twofl_thermal_conduction_c) then
        twofl_thermal_conduction_c=.false.
        if(mype==0) write(*,*) 'WARNING: set twofl_thermal_conduction_c=F when twofl_energy=F'
      end if
      if(twofl_radiative_cooling_c) then
        twofl_radiative_cooling_c=.false.
        if(mype==0) write(*,*) 'WARNING: set twofl_radiative_cooling_c=F when twofl_energy=F'
      end if
      if(twofl_trac) then
        twofl_trac=.false.
        if(mype==0) write(*,*) 'WARNING: set twofl_trac=F when twofl_energy=F'
      end if
    end if
    {^ifoned
      if(twofl_trac .and. twofl_trac_type .gt. 1) then
        twofl_trac_type=1
        if(mype==0) write(*,*) 'WARNING: set twofl_trac_type=1 for 1D simulation'
      end if
    }
    if(twofl_trac .and. twofl_trac_type .le. 3) then
      twofl_trac_mask=bigdouble
      if(mype==0) write(*,*) 'WARNING: set twofl_trac_mask==bigdouble for global TRAC method'
    end if
    phys_trac_mask=twofl_trac_mask
 
    ! set default gamma for polytropic/isothermal process
    if(ndim==1) typedivbfix='none'
    select case (typedivbfix)
    case ('none')
       type_divb = divb_none
    {^nooned
    case ('multigrid')
       type_divb = divb_multigrid
       use_multigrid = .true.
       mg%operator_type = mg_laplacian
       phys_global_source_after => twofl_clean_divb_multigrid
    }
    case ('glm')
      twofl_glm          = .true.
      need_global_cmax = .true.
      type_divb        = divb_glm
    case ('powel', 'powell')
      type_divb = divb_powel
    case ('janhunen')
      type_divb = divb_janhunen
    case ('linde')
      type_divb = divb_linde
    case ('lindejanhunen')
      type_divb = divb_lindejanhunen
    case ('lindepowel')
      type_divb = divb_lindepowel
    case ('lindeglm')
      twofl_glm          = .true.
      need_global_cmax = .true.
      type_divb        = divb_lindeglm
    case ('ct')
      type_divb = divb_ct
      stagger_grid = .true.
    case default
      call mpistop('Unknown divB fix')
    end select
 
    allocate(start_indices(number_species))
    allocate(stop_indices(number_species))
    start_indices(1)=1
    !allocate charges first and the same order as in mhd module
    rho_c_ = var_set_fluxvar("rho_c", "rho_c")
    !set variables from mod_variables to point to charges vars
    iw_rho = rho_c_
 
    allocate(mom_c(ndir))
    do idir=1,ndir
      mom_c(idir) = var_set_fluxvar("m_c","v_c",idir)
    enddo
 
    allocate(iw_mom(ndir))
    iw_mom(1:ndir) = mom_c(1:ndir)
 
    ! Set index of energy variable
    if (phys_energy) then
      e_c_ = var_set_fluxvar("e_c", "p_c")
      iw_e = e_c_
    else
      e_c_ = -1
    end if
 
  ! ambipolar sts assumes mag and energy charges are continuous
    allocate(mag(ndir))
    mag(:) = var_set_bfield(ndir)
 
    if (twofl_glm) then
      psi_ = var_set_fluxvar('psi', 'psi', need_bc=.false.)
    else
      psi_ = -1
    end if
 
    ! set cutoff temperature when using the TRAC method, as well as an auxiliary weight
    tweight_c_ = -1
    if(twofl_trac) then
      tcoff_c_ = var_set_wextra()
      iw_tcoff = tcoff_c_
      if(twofl_trac_type > 2) then
        tweight_c_ = var_set_wextra()
      endif
    else
      tcoff_c_ = -1
    end if
 
  !now allocate neutrals
 
    ! TODO so far number_species is only used to treat them differently
    ! in the solvers (different cbounds)
    if (twofl_cbounds_species) then
      stop_indices(1)=nwflux
      start_indices(2)=nwflux+1
    endif
 
    ! Determine flux variables
    rho_n_ = var_set_fluxvar("rho_n", "rho_n")
    allocate(mom_n(ndir))
    do idir=1,ndir
      mom_n(idir) = var_set_fluxvar("m_n","v_n",idir)
    enddo
    if (phys_energy) then
      e_n_ = var_set_fluxvar("e_n", "p_n")
    else
      e_n_     = -1
    end if
 
    tweight_n_ = -1
    if(twofl_trac) then
      tcoff_n_ = var_set_wextra()
      if(twofl_trac_type > 2) then
        tweight_n_ = var_set_wextra()
      endif
    else
      tcoff_n_ = -1
    end if
 
    stop_indices(number_species)=nwflux
 
    ! set indices of equi vars and update number_equi_vars
    number_equi_vars = 0
    if(has_equi_rho_n0) then
      number_equi_vars = number_equi_vars + 1
      equi_rho_n0_ = number_equi_vars
    endif  
    if(has_equi_pe_n0) then
      number_equi_vars = number_equi_vars + 1
      equi_pe_n0_ = number_equi_vars
      phys_equi_pe=.true.
    endif  
    if(has_equi_rho_c0) then
      number_equi_vars = number_equi_vars + 1
      equi_rho_c0_ = number_equi_vars
      iw_equi_rho = equi_rho_c0_
    endif  
    if(has_equi_pe_c0) then
      number_equi_vars = number_equi_vars + 1
      equi_pe_c0_ = number_equi_vars
      iw_equi_p = equi_pe_c0_
      phys_equi_pe=.true.
    endif  
 
    ! set number of variables which need update ghostcells
    nwgc=nwflux+nwaux
 
    ! determine number of stagger variables
    nws=ndim
 
    ! Check whether custom flux types have been defined
    if (.not. allocated(flux_type)) then
       allocate(flux_type(ndir, nw))
       flux_type = flux_default
    else if (any(shape(flux_type) /= [ndir, nw])) then
       call mpistop("phys_check error: flux_type has wrong shape")
    end if
 
    if(ndim>1) then
      if(twofl_glm) then
        flux_type(:,psi_)=flux_special
        do idir=1,ndir
          flux_type(idir,mag(idir))=flux_special
        end do
      else
        do idir=1,ndir
          flux_type(idir,mag(idir))=flux_tvdlf
        end do
      end if
    end if
 
    phys_get_dt              => twofl_get_dt
    phys_get_cmax            => twofl_get_cmax
    phys_get_a2max           => twofl_get_a2max
    !phys_get_tcutoff         => twofl_get_tcutoff_c
    if(twofl_cbounds_species) then
      if (mype .eq. 0) print*, "Using different cbounds for each species nspecies = ", number_species
      phys_get_cbounds         => twofl_get_cbounds_species
      phys_get_h_speed         => twofl_get_h_speed_species
    else
      if (mype .eq. 0) print*, "Using same cbounds for all species"
      phys_get_cbounds         => twofl_get_cbounds_one
      phys_get_h_speed         => twofl_get_h_speed_one
    endif
    phys_get_flux            => twofl_get_flux
    phys_add_source_geom     => twofl_add_source_geom
    phys_add_source          => twofl_add_source
    phys_to_conserved        => twofl_to_conserved
    phys_to_primitive        => twofl_to_primitive
    phys_check_params        => twofl_check_params
    phys_check_w             => twofl_check_w
    phys_write_info          => twofl_write_info
    phys_handle_small_values => twofl_handle_small_values
    !set equilibrium variables for the new grid
    if(number_equi_vars>0) then
      phys_set_equi_vars => set_equi_vars_grid
    endif
    ! convert_type is not known here, so associate the corresp. subroutine in check_params
    if(type_divb==divb_glm) then
      phys_modify_wlr => twofl_modify_wlr
    end if
 
    ! if using ct stagger grid, boundary divb=0 is not done here
    if(stagger_grid) then
      phys_get_ct_velocity => twofl_get_ct_velocity
      phys_update_faces => twofl_update_faces
      phys_face_to_center => twofl_face_to_center
      phys_modify_wlr => twofl_modify_wlr
    else if(ndim>1) then
      phys_boundary_adjust => twofl_boundary_adjust
    end if
 
    {^nooned
    ! clean initial divb
    if(clean_initial_divb) phys_clean_divb => twofl_clean_divb_multigrid
    }
 
    ! Whether diagonal ghost cells are required for the physics
    if(type_divb < divb_linde) phys_req_diagonal = .false.
 
    ! derive units from basic units
    call twofl_physical_units()
 
    if(.not. phys_energy .and. (twofl_thermal_conduction_c& 
        .or. twofl_thermal_conduction_n)) then
      call mpistop("thermal conduction needs twofl_energy=T")
    end if
 
    ! initialize thermal conduction module
    if (twofl_thermal_conduction_c &
        .or. twofl_thermal_conduction_n) then
      phys_req_diagonal = .true.
      call sts_init()
      call tc_init_params(twofl_gamma)
    endif
    if (twofl_thermal_conduction_c) then
      allocate(tc_fl_c)
      if(has_equi_pe_c0 .and. has_equi_rho_c0) then
        tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c_with_equi
        if(phys_internal_e) then
          tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c_with_equi
        else
            if(twofl_eq_energy == eq_energy_ki) then
              tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c_with_equi
            else
              tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c_with_equi
            endif
        endif
        if(twofl_equi_thermal_c) then
          tc_fl_c%has_equi = .true.
          tc_fl_c%get_temperature_equi => twofl_get_temperature_c_equi
          tc_fl_c%get_rho_equi => twofl_get_rho_c_equi
        else  
          tc_fl_c%has_equi = .false.
        endif
      else
        if(phys_internal_e) then
          tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c
         else
            if(twofl_eq_energy == eq_energy_ki) then
              tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c
            else
              tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c
          endif  
         endif  
        tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c
      endif
      if(use_twofl_tc_c .eq. mhd_tc) then
        call tc_get_mhd_params(tc_fl_c,tc_c_params_read_mhd)
        call add_sts_method(twofl_get_tc_dt_mhd_c,twofl_sts_set_source_tc_c_mhd,e_c_,1,e_c_,1,.false.)
      else if(use_twofl_tc_c .eq. hd_tc) then
        call tc_get_hd_params(tc_fl_c,tc_c_params_read_hd)
        call add_sts_method(twofl_get_tc_dt_hd_c,twofl_sts_set_source_tc_c_hd,e_c_,1,e_c_,1,.false.)
      endif
      if(.not. phys_internal_e) then
        call set_conversion_methods_to_head(twofl_e_to_ei_c, twofl_ei_to_e_c)
      endif
      call set_error_handling_to_head(twofl_tc_handle_small_e_c)
      tc_fl_c%get_rho => get_rhoc_tot
      tc_fl_c%e_ = e_c_
      tc_fl_c%Tcoff_ = tcoff_c_
    end if
    if (twofl_thermal_conduction_n) then
      allocate(tc_fl_n)
      call tc_get_hd_params(tc_fl_n,tc_n_params_read_hd)
      if(has_equi_pe_n0 .and. has_equi_rho_n0) then
        tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n_with_equi
        if(twofl_equi_thermal_n) then
          tc_fl_n%has_equi = .true.
          tc_fl_n%get_temperature_equi => twofl_get_temperature_n_equi
          tc_fl_n%get_rho_equi => twofl_get_rho_n_equi
        else  
          tc_fl_n%has_equi = .false.
        endif
      else
        tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n
      endif
      if(phys_internal_e) then
        if(has_equi_pe_n0 .and. has_equi_rho_n0) then
          tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n_with_equi
        else
          tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n
        endif 
        call add_sts_method(twofl_get_tc_dt_hd_n,twofl_sts_set_source_tc_n_hd,e_n_,1,e_n_,1,.false.)
      else
        if(has_equi_pe_n0 .and. has_equi_rho_n0) then
          tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n_with_equi
        else
          tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n
        endif 
        call add_sts_method(twofl_get_tc_dt_hd_n,twofl_sts_set_source_tc_n_hd,e_n_,1,e_n_,1,.false.)
        call set_conversion_methods_to_head(twofl_e_to_ei_n, twofl_ei_to_e_n)
      endif
      call set_error_handling_to_head(twofl_tc_handle_small_e_n)
      tc_fl_n%get_rho => get_rhon_tot
      tc_fl_n%e_ = e_n_
      tc_fl_n%Tcoff_ = tcoff_n_
    end if
 
 
    if(.not. phys_energy .and. (twofl_radiative_cooling_c& 
        .or. twofl_radiative_cooling_n)) then
      call mpistop("radiative cooling needs twofl_energy=T")
    end if
 
    if(twofl_equi_thermal .and. (.not. has_equi_pe_c0 .or. .not.  has_equi_pe_n0)) then
      call mpistop("twofl_equi_thermal=T has_equi_pe_n0 and has _equi_pe_c0=T")
    endif
 
    ! initialize thermal conduction module
    if (twofl_radiative_cooling_c &
        .or. twofl_radiative_cooling_n) then
    ! Initialize radiative cooling module
      call radiative_cooling_init_params(twofl_gamma,he_abundance)
      if(twofl_radiative_cooling_c) then
        allocate(rc_fl_c)
        call radiative_cooling_init(rc_fl_c,rc_params_read_c)
        rc_fl_c%get_rho => get_rhoc_tot
        rc_fl_c%get_pthermal => twofl_get_pthermal_c
        rc_fl_c%get_var_Rfactor => rfactor_c
        rc_fl_c%e_ = e_c_
        rc_fl_c%Tcoff_ = tcoff_c_
        if(has_equi_pe_c0 .and. has_equi_rho_c0 .and. twofl_equi_thermal_c) then
          rc_fl_c%has_equi = .true.
          rc_fl_c%get_rho_equi => twofl_get_rho_c_equi
          rc_fl_c%get_pthermal_equi => twofl_get_pe_c_equi
        else
          rc_fl_c%has_equi = .false.
        end if
      end if
    end if
    allocate(te_fl_c)
    te_fl_c%get_rho=> get_rhoc_tot
    te_fl_c%get_pthermal=> twofl_get_pthermal_c
    te_fl_c%get_var_Rfactor => rfactor_c
{^ifthreed
    phys_te_images => twofl_te_images
}
 
    ! Initialize viscosity module
    !!TODO
    !if (twofl_viscosity) call viscosity_init(phys_wider_stencil,phys_req_diagonal)
 
    ! Initialize gravity module
    if(twofl_gravity) then
    !  call gravity_init()
       call grav_params_read(par_files)
    end if
 
    ! Initialize particles module
    ! For Hall, we need one more reconstructed layer since currents are computed
    ! in getflux: assuming one additional ghost layer (two for FOURTHORDER) was
    ! added in nghostcells.
    if (twofl_hall) then
       phys_req_diagonal = .true.
       if (twofl_4th_order) then
          phys_wider_stencil = 2
       else
          phys_wider_stencil = 1
       end if
    end if
 
    if(twofl_hyperdiffusivity) then
      allocate(c_shk(1:nwflux))
      allocate(c_hyp(1:nwflux))
      call twofl_init_hyper(par_files)
    end if
 
  end subroutine twofl_phys_init
 
{^ifthreed
  subroutine twofl_te_images
    use mod_global_parameters
    use mod_thermal_emission
 
    select case(convert_type)
      case('EIvtiCCmpi','EIvtuCCmpi')
        call get_euv_image(unitconvert,te_fl_c)
      case('ESvtiCCmpi','ESvtuCCmpi')
        call get_euv_spectrum(unitconvert,te_fl_c)
      case('SIvtiCCmpi','SIvtuCCmpi')
        call get_sxr_image(unitconvert,te_fl_c)
      case('WIvtiCCmpi','WIvtuCCmpi')
        call get_whitelight_image(unitconvert,te_fl_c)
      case default
        call mpistop("Error in synthesize emission: Unknown convert_type")
      end select
  end subroutine twofl_te_images
}  
 
  ! wrappers for STS functions in thermal_conductivity module
  ! which take as argument the tc_fluid (defined in the physics module)
  subroutine  twofl_sts_set_source_tc_c_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
    use mod_global_parameters
    use mod_fix_conserve
    use mod_thermal_conduction, only: sts_set_source_tc_mhd
    integer, intent(in) :: ixI^L, ixO^L, igrid, nflux
    double precision, intent(in) ::  x(ixI^S,1:ndim)
    double precision, intent(inout) ::  wres(ixI^S,1:nw), w(ixI^S,1:nw)
    double precision, intent(in) :: my_dt
    logical, intent(in) :: fix_conserve_at_step
    call sts_set_source_tc_mhd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl_c)
  end subroutine twofl_sts_set_source_tc_c_mhd
 
  subroutine  twofl_sts_set_source_tc_c_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
    use mod_global_parameters
    use mod_fix_conserve
    use mod_thermal_conduction, only: sts_set_source_tc_hd
    integer, intent(in) :: ixI^L, ixO^L, igrid, nflux
    double precision, intent(in) ::  x(ixI^S,1:ndim)
    double precision, intent(inout) ::  wres(ixI^S,1:nw), w(ixI^S,1:nw)
    double precision, intent(in) :: my_dt
    logical, intent(in) :: fix_conserve_at_step
    call sts_set_source_tc_hd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl_c)
  end subroutine twofl_sts_set_source_tc_c_hd
 
  function twofl_get_tc_dt_mhd_c(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
    !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
    !where                      tc_k_para_i=tc_k_para*B_i**2/B**2
    !and                        T=p/rho
    use mod_global_parameters
    use mod_thermal_conduction, only: get_tc_dt_mhd
 
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
    double precision, intent(in) :: w(ixi^s,1:nw)
    double precision :: dtnew
 
    dtnew=get_tc_dt_mhd(w,ixi^l,ixo^l,dx^d,x,tc_fl_c) 
  end function twofl_get_tc_dt_mhd_c
 
  function twofl_get_tc_dt_hd_c(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
    !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
    !where                      tc_k_para_i=tc_k_para*B_i**2/B**2
    !and                        T=p/rho
    use mod_global_parameters
    use mod_thermal_conduction, only: get_tc_dt_hd
 
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
    double precision, intent(in) :: w(ixi^s,1:nw)
    double precision :: dtnew
 
    dtnew=get_tc_dt_hd(w,ixi^l,ixo^l,dx^d,x,tc_fl_c) 
  end function twofl_get_tc_dt_hd_c
 
  subroutine twofl_tc_handle_small_e_c(w, x, ixI^L, ixO^L, step)
    use mod_global_parameters
    use mod_small_values
 
    integer, intent(in)             :: ixI^L,ixO^L
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    integer, intent(in)    :: step
 
    character(len=140) :: error_msg
 
    write(error_msg,"(a,i3)") "Charges thermal conduction step ", step
    call twofl_handle_small_ei_c(w,x,ixi^l,ixo^l,e_c_,error_msg)
  end subroutine twofl_tc_handle_small_e_c
 
  subroutine  twofl_sts_set_source_tc_n_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
    use mod_global_parameters
    use mod_fix_conserve
    use mod_thermal_conduction, only: sts_set_source_tc_hd
    integer, intent(in) :: ixI^L, ixO^L, igrid, nflux
    double precision, intent(in) ::  x(ixI^S,1:ndim)
    double precision, intent(inout) ::  wres(ixI^S,1:nw), w(ixI^S,1:nw)
    double precision, intent(in) :: my_dt
    logical, intent(in) :: fix_conserve_at_step
    call sts_set_source_tc_hd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl_n)
  end subroutine twofl_sts_set_source_tc_n_hd
 
  subroutine twofl_tc_handle_small_e_n(w, x, ixI^L, ixO^L, step)
    use mod_global_parameters
 
    integer, intent(in)             :: ixI^L,ixO^L
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    integer, intent(in)    :: step
 
    character(len=140) :: error_msg
 
    write(error_msg,"(a,i3)") "Neutral thermal conduction step ", step
    call twofl_handle_small_ei_n(w,x,ixi^l,ixo^l,e_n_,error_msg)
  end subroutine twofl_tc_handle_small_e_n
 
  function twofl_get_tc_dt_hd_n(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
    !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
    !where                      tc_k_para_i=tc_k_para*B_i**2/B**2
    !and                        T=p/rho
    use mod_global_parameters
    use mod_thermal_conduction, only: get_tc_dt_hd
 
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
    double precision, intent(in) :: w(ixi^s,1:nw)
    double precision :: dtnew
 
    dtnew=get_tc_dt_hd(w,ixi^l,ixo^l,dx^d,x,tc_fl_n) 
  end function twofl_get_tc_dt_hd_n
 
  subroutine tc_n_params_read_hd(fl)
    use mod_global_parameters, only: unitpar,par_files
    use mod_global_parameters, only: unitpar
    type(tc_fluid), intent(inout) :: fl
    integer                      :: n
    logical :: tc_saturate=.false.
    double precision :: tc_k_para=0d0
 
    namelist /tc_n_list/ tc_saturate, tc_k_para
 
    do n = 1, size(par_files)
       open(unitpar, file=trim(par_files(n)), status="old")
       read(unitpar, tc_n_list, end=111)
111      close(unitpar)
    end do
    fl%tc_saturate = tc_saturate
    fl%tc_k_para = tc_k_para
 
  end subroutine tc_n_params_read_hd
 
  subroutine rc_params_read_n(fl)
    use mod_global_parameters, only: unitpar,par_files
    use mod_constants, only: bigdouble
    type(rc_fluid), intent(inout) :: fl
    integer                      :: n
    ! list parameters
    integer :: ncool = 4000
    double precision :: cfrac=0.1d0
  
    !> Name of cooling curve
    character(len=std_len)  :: coolcurve='JCorona'
  
    !> Name of cooling method
    character(len=std_len)  :: coolmethod='exact'
  
    !> Fixed temperature not lower than tlow
    logical    :: Tfix=.false.
  
    !> Lower limit of temperature
    double precision   :: tlow=bigdouble
  
    !> Add cooling source in a split way (.true.) or un-split way (.false.)
    logical    :: rc_split=.false.
 
    namelist /rc_list_n/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
 
    do n = 1, size(par_files)
      open(unitpar, file=trim(par_files(n)), status="old")
      read(unitpar, rc_list_n, end=111)
111     close(unitpar)
    end do
 
    fl%ncool=ncool
    fl%coolcurve=coolcurve
    fl%coolmethod=coolmethod
    fl%tlow=tlow
    fl%Tfix=tfix
    fl%rc_split=rc_split
    fl%cfrac=cfrac
  end subroutine rc_params_read_n
 
  !end wrappers
 
  ! fill in tc_fluid fields from namelist
  subroutine tc_c_params_read_mhd(fl)
    use mod_global_parameters, only: unitpar,par_files
    type(tc_fluid), intent(inout) :: fl
 
    integer                      :: n
 
    ! list parameters
    logical :: tc_perpendicular=.false.
    logical :: tc_saturate=.false.
    double precision :: tc_k_para=0d0
    double precision :: tc_k_perp=0d0
    character(len=std_len)  :: tc_slope_limiter="MC"
 
    namelist /tc_c_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
    do n = 1, size(par_files)
      open(unitpar, file=trim(par_files(n)), status="old")
      read(unitpar, tc_c_list, end=111)
111     close(unitpar)
    end do
 
    fl%tc_perpendicular = tc_perpendicular
    fl%tc_saturate = tc_saturate
    fl%tc_k_para = tc_k_para
    fl%tc_k_perp = tc_k_perp
    select case(tc_slope_limiter)
     case ('no','none')
       fl%tc_slope_limiter = 0
     case ('MC')
       ! montonized central limiter Woodward and Collela limiter (eq.3.51h), a factor of 2 is pulled out
       fl%tc_slope_limiter = 1
     case('minmod')
       ! minmod limiter
       fl%tc_slope_limiter = 2
     case ('superbee')
       ! Roes superbee limiter (eq.3.51i)
       fl%tc_slope_limiter = 3
     case ('koren')
       ! Barry Koren Right variant
       fl%tc_slope_limiter = 4
     case default
       call mpistop("Unknown tc_slope_limiter, choose MC, minmod")
    end select
  end subroutine tc_c_params_read_mhd
 
  subroutine tc_c_params_read_hd(fl)
    use mod_global_parameters, only: unitpar,par_files
    use mod_global_parameters, only: unitpar
    type(tc_fluid), intent(inout) :: fl
    integer                      :: n
    logical :: tc_saturate=.false.
    double precision :: tc_k_para=0d0
 
    namelist /tc_c_list/ tc_saturate, tc_k_para
 
    do n = 1, size(par_files)
       open(unitpar, file=trim(par_files(n)), status="old")
       read(unitpar, tc_c_list, end=111)
111      close(unitpar)
    end do
    fl%tc_saturate = tc_saturate
    fl%tc_k_para = tc_k_para
 
  end subroutine tc_c_params_read_hd
 
!! end th cond
 
!!rad cool
  subroutine rc_params_read_c(fl)
    use mod_global_parameters, only: unitpar,par_files
    use mod_constants, only: bigdouble
    type(rc_fluid), intent(inout) :: fl
    integer                      :: n
    ! list parameters
    integer :: ncool = 4000
    double precision :: cfrac=0.1d0
  
    !> Name of cooling curve
    character(len=std_len)  :: coolcurve='JCcorona'
  
    !> Name of cooling method
    character(len=std_len)  :: coolmethod='exact'
  
    !> Fixed temperature not lower than tlow
    logical    :: Tfix=.false.
  
    !> Lower limit of temperature
    double precision   :: tlow=bigdouble
  
    !> Add cooling source in a split way (.true.) or un-split way (.false.)
    logical    :: rc_split=.false.
 
 
    namelist /rc_list_c/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
 
    do n = 1, size(par_files)
      open(unitpar, file=trim(par_files(n)), status="old")
      read(unitpar, rc_list_c, end=111)
111     close(unitpar)
    end do
 
    fl%ncool=ncool
    fl%coolcurve=coolcurve
    fl%coolmethod=coolmethod
    fl%tlow=tlow
    fl%Tfix=tfix
    fl%rc_split=rc_split
    fl%cfrac=cfrac
  end subroutine rc_params_read_c
 
!! end rad cool
 
  !> sets the equilibrium variables
  subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
    use mod_global_parameters
    use mod_global_parameters
    use mod_usr_methods
    integer, intent(in) :: igrid, ixI^L, ixO^L
    double precision, intent(in) :: x(ixI^S,1:ndim)
 
    double precision :: delx(ixI^S,1:ndim)
    double precision :: xC(ixI^S,1:ndim),xshift^D
    integer :: idims, ixC^L, hxO^L, ix, idims2
 
    if(slab_uniform)then
      ^d&delx(ixi^s,^d)=rnode(rpdx^d_,igrid)\
    else
      ! for all non-cartesian and stretched cartesian coordinates
      delx(ixi^s,1:ndim)=ps(igrid)%dx(ixi^s,1:ndim)
    endif
  
  
    do idims=1,ndim
      hxo^l=ixo^l-kr(idims,^d);
      if(stagger_grid) then
        ! ct needs all transverse cells
        ixcmax^d=ixomax^d+nghostcells-nghostcells*kr(idims,^d); ixcmin^d=hxomin^d-nghostcells+nghostcells*kr(idims,^d);
      else
        ! ixC is centered index in the idims direction from ixOmin-1/2 to ixOmax+1/2
        ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
      end if
      ! always xshift=0 or 1/2
      xshift^d=half*(one-kr(^d,idims));
      do idims2=1,ndim
        select case(idims2)
        {case(^d)
          do ix = ixc^lim^d
            ! xshift=half: this is the cell center coordinate
            ! xshift=0: this is the cell edge i+1/2 coordinate
            xc(ix^d%ixC^s,^d)=x(ix^d%ixC^s,^d)+(half-xshift^d)*delx(ix^d%ixC^s,^d)
          end do\}
        end select
      end do
      call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
    end do
 
  end subroutine set_equi_vars_grid_faces
 
  !> sets the equilibrium variables
  subroutine set_equi_vars_grid(igrid)
    use mod_global_parameters
    use mod_usr_methods
 
    integer, intent(in) :: igrid
 
    !values at the center
    call usr_set_equi_vars(ixg^ll,ixg^ll,ps(igrid)%x,ps(igrid)%equi_vars(ixg^t,1:number_equi_vars,0))
 
    !values at the interfaces
    call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^ll,ixm^ll)
 
  end subroutine set_equi_vars_grid
 
  ! w, wnew conserved
  function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc) result(wnew)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l,ixo^l, nwc
    double precision, intent(in)    :: w(ixi^s, 1:nw)
    double precision, intent(in)    :: x(ixi^s,1:ndim) 
    double precision   :: wnew(ixo^s, 1:nwc)
    double precision   :: rho(ixi^s)
 
    call  get_rhon_tot(w,x,ixi^l,ixo^l,rho(ixi^s))
    wnew(ixo^s,rho_n_) = rho(ixo^s)
    wnew(ixo^s,mom_n(:)) =  w(ixo^s,mom_n(:))
    call  get_rhoc_tot(w,x,ixi^l,ixo^l,rho(ixi^s))
    wnew(ixo^s,rho_c_) = rho(ixo^s)
    wnew(ixo^s,mom_c(:)) =  w(ixo^s,mom_c(:))
 
    if (b0field) then
      ! add background magnetic field B0 to B
      wnew(ixo^s,mag(:))=w(ixo^s,mag(:))+block%B0(ixo^s,:,0)
    else 
      wnew(ixo^s,mag(:))=w(ixo^s,mag(:))
    end if
 
    if(phys_energy) then
      wnew(ixo^s,e_n_) = w(ixo^s,e_n_)
      if(has_equi_pe_n0) then
        wnew(ixo^s,e_n_) = wnew(ixo^s,e_n_) + block%equi_vars(ixo^s,equi_pe_n0_,0)* inv_gamma_1  
      endif
      wnew(ixo^s,e_c_) = w(ixo^s,e_c_)
      if(has_equi_pe_c0) then
        wnew(ixo^s,e_c_) = wnew(ixo^s,e_c_) + block%equi_vars(ixo^s,equi_pe_c0_,0)* inv_gamma_1  
      endif
      if(b0field .and. phys_total_energy) then
          wnew(ixo^s,e_c_)=wnew(ixo^s,e_c_)+0.5d0*sum(block%B0(ixo^s,:,0)**2,dim=ndim+1) &
              + sum(w(ixo^s,mag(:))*block%B0(ixo^s,:,0),dim=ndim+1)
      endif
    endif
 
  end function convert_vars_splitting
 
  !> copied from mod_gravity
  subroutine grav_params_read(files)
    use mod_global_parameters, only: unitpar
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /grav_list/ grav_split
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status="old")
       read(unitpar, grav_list, end=111)
111    close(unitpar)
    end do
 
  end subroutine grav_params_read
 
  subroutine associate_dump_hyper()
    use mod_global_parameters
    use mod_convert, only: add_convert_method
    integer :: ii
    do ii = 1,ndim 
      if(ii==1) then
        call add_convert_method(dump_hyperdiffusivity_coef_x, nw, cons_wnames(1:nw), "hyper_x")
      elseif(ii==2) then
        call add_convert_method(dump_hyperdiffusivity_coef_y, nw, cons_wnames(1:nw), "hyper_y")
      else
        call add_convert_method(dump_hyperdiffusivity_coef_z, nw, cons_wnames(1:nw), "hyper_z")
      endif
    enddo
  end subroutine associate_dump_hyper
 
  subroutine twofl_check_params
    use mod_global_parameters
    use mod_usr_methods
    use mod_convert, only: add_convert_method
 
    ! after user parameter setting
    gamma_1=twofl_gamma-1.d0
    if (.not. phys_energy) then
       if (twofl_gamma <= 0.0d0) call mpistop ("Error: twofl_gamma <= 0")
       if (twofl_adiab < 0.0d0) call mpistop ("Error: twofl_adiab < 0")
       small_pressure = twofl_adiab*small_density**twofl_gamma
    else
       if (twofl_gamma <= 0.0d0 .or. twofl_gamma == 1.0d0) &
            call mpistop ("Error: twofl_gamma <= 0 or twofl_gamma == 1")
       inv_gamma_1=1.d0/gamma_1
       small_e = small_pressure * inv_gamma_1
    end if
 
    ! this has to be done here as use_imex_scheme is not set in init subroutine, 
    ! but here it is
    if(use_imex_scheme) then
      if(has_collisions()) then
        ! implicit collisional terms update
        phys_implicit_update => twofl_implicit_coll_terms_update
        phys_evaluate_implicit => twofl_evaluate_implicit
        if(mype .eq. 1) then
            print*, "IMPLICIT UPDATE with calc_mult_factor", twofl_implicit_calc_mult_method
        endif
        if(twofl_implicit_calc_mult_method == 1) then
          calc_mult_factor => calc_mult_factor1
        else
          calc_mult_factor => calc_mult_factor2
        endif
      endif
    else
      ! check dtcoll par for explicit implementation of the coll. terms
      if(dtcollpar .le. 0d0 .or. dtcollpar .ge. 1d0) then
        if (mype .eq. 0) print*, "Explicit update of coll terms requires 0<dtcollpar<1, dtcollpar set to 0.8."
        dtcollpar = 0.8
      endif 
        
    endif
!    if(H_ion_fr == 0d0 .and. He_ion_fr == 0d0) then
!      call mpistop("H_ion_fr or He_ion_fr must be > 0 or use hd module")
!    endif
!    if(H_ion_fr == 1d0 .and. He_ion_fr == 1d0) then
!      call mpistop("H_ion_fr or He_ion_fr must be < 1 or use mhd module")
!    endif
    if (number_equi_vars > 0 .and. .not. associated(usr_set_equi_vars)) then
      call mpistop("usr_set_equi_vars has to be implemented in the user file")
    endif
    if(convert .or. autoconvert) then
      if(convert_type .eq. 'dat_generic_mpi') then
        if(twofl_dump_full_vars) then
          if(mype .eq. 0) print*, " add conversion method: split -> full "
          call add_convert_method(convert_vars_splitting, nw, cons_wnames, "new")
        endif
        if(twofl_dump_coll_terms) then
          if(mype .eq. 0) print*, " add conversion method: dump coll terms "
          call add_convert_method(dump_coll_terms, 3, (/"alpha    ", "gamma_rec", "gamma_ion"/), "_coll")
        endif
        if(twofl_hyperdiffusivity .and. twofl_dump_hyperdiffusivity_coef) then
          if(mype .eq. 0) print*, " add conversion method: dump hyperdiffusivity coeff. "
          call associate_dump_hyper()
        endif
      endif
    endif
  end subroutine twofl_check_params
 
  subroutine twofl_physical_units()
    use mod_global_parameters
    double precision :: mp,kB,miu0,c_lightspeed
    !double precision :: a,b,c,d
    double precision :: a,b
    ! Derive scaling units
    if(si_unit) then
      mp=mp_si
      kb=kb_si
      miu0=miu0_si
      c_lightspeed=c_si
    else
      mp=mp_cgs
      kb=kb_cgs
      miu0=4.d0*dpi
      c_lightspeed=const_c
    end if
 
    a=1d0  
    b=1d0
    rc=2d0
    rn=1d0  
 
    !now the unit choice:
    !unit 1 from number density or density -> mH
    !unit 2 from 
 
    if(unit_density/=1.d0) then
      unit_numberdensity=unit_density/(a*mp)
    else
      ! unit of numberdensity is independent by default
      unit_density=a*mp*unit_numberdensity
    end if
    if(unit_velocity/=1.d0) then
      unit_pressure=unit_density*unit_velocity**2
      unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
      unit_magneticfield=sqrt(miu0*unit_pressure)
    else if(unit_pressure/=1.d0) then
      unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
      unit_velocity=sqrt(unit_pressure/unit_density)
      unit_magneticfield=sqrt(miu0*unit_pressure)
    else if(unit_magneticfield/=1.d0) then
      unit_pressure=unit_magneticfield**2/miu0
      unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
      unit_velocity=sqrt(unit_pressure/unit_density)
    else if(unit_temperature/=1.d0) then
      unit_pressure=b*unit_numberdensity*kb*unit_temperature
      unit_velocity=sqrt(unit_pressure/unit_density)
      unit_magneticfield=sqrt(miu0*unit_pressure)
    end if
    if(unit_time/=1.d0) then
      unit_length=unit_time*unit_velocity
    else
      ! unit of length is independent by default
      unit_time=unit_length/unit_velocity
    end if
    ! Additional units needed for the particles
    c_norm=c_lightspeed/unit_velocity
    unit_charge=unit_magneticfield*unit_length**2/unit_velocity/miu0
    if (.not. si_unit) unit_charge = unit_charge*const_c
    unit_mass=unit_density*unit_length**3
  end subroutine twofl_physical_units
 
  subroutine twofl_check_w(primitive,ixI^L,ixO^L,w,flag)
    use mod_global_parameters
 
    logical, intent(in) :: primitive
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S,nw)
    double precision :: tmp(ixI^S)
    logical, intent(inout) :: flag(ixI^S,1:nw)
 
    flag=.false.
        
    if(has_equi_rho_n0) then
      tmp(ixo^s) = w(ixo^s,rho_n_) + block%equi_vars(ixo^s,equi_rho_n0_,0)
    else  
      tmp(ixo^s) = w(ixo^s,rho_n_) 
    endif
    where(tmp(ixo^s) < small_density) flag(ixo^s,rho_n_) = .true.
    if(has_equi_rho_c0) then
      tmp(ixo^s) = w(ixo^s,rho_c_) + block%equi_vars(ixo^s,equi_rho_c0_,0)
    else  
      tmp(ixo^s) = w(ixo^s,rho_c_) 
    endif
    where(tmp(ixo^s) < small_density) flag(ixo^s,rho_c_) = .true.
    if(phys_energy) then
      if(primitive) then
        tmp(ixo^s) = w(ixo^s,e_n_)
        if(has_equi_pe_n0) then
          tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe_n0_,0)
        endif
        where(tmp(ixo^s) < small_pressure) flag(ixo^s,e_n_) = .true.
        tmp(ixo^s) = w(ixo^s,e_c_)
        if(has_equi_pe_c0) then
          tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe_c0_,0)
        endif
        where(tmp(ixo^s) < small_pressure) flag(ixo^s,e_c_) = .true.
      else
        if(phys_internal_e) then
          tmp(ixo^s)=w(ixo^s,e_n_)
          if(has_equi_pe_n0) then
            tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe_n0_,0)*inv_gamma_1
          endif
          where(tmp(ixo^s) < small_e) flag(ixo^s,e_n_) = .true.
          tmp(ixo^s)=w(ixo^s,e_c_)
          if(has_equi_pe_c0) then
            tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe_c0_,0)*inv_gamma_1
          endif
          where(tmp(ixo^s) < small_e) flag(ixo^s,e_c_) = .true.
        else
          !neutrals
          tmp(ixo^s)=w(ixo^s,e_n_)-&
                twofl_kin_en_n(w,ixi^l,ixo^l)
          if(has_equi_pe_n0) then
            tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe_n0_,0)*inv_gamma_1
          endif
          where(tmp(ixo^s) < small_e) flag(ixo^s,e_n_) = .true.
          if(phys_total_energy) then
            tmp(ixo^s)=w(ixo^s,e_c_)-&
                twofl_kin_en_c(w,ixi^l,ixo^l)-twofl_mag_en(w,ixi^l,ixo^l)
          else
            tmp(ixo^s)=w(ixo^s,e_c_)-&
                twofl_kin_en_c(w,ixi^l,ixo^l)
          end if
          if(has_equi_pe_c0) then
            tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe_c0_,0)*inv_gamma_1
          endif
          where(tmp(ixo^s) < small_e) flag(ixo^s,e_c_) = .true.
        end if
      endif
    end if
 
  end subroutine twofl_check_w
 
  !> Transform primitive variables into conservative ones
  subroutine twofl_to_conserved(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
    integer                         :: idir
    double precision                :: rhoc(ixi^s)
    double precision                :: rhon(ixi^s)
 
    !if (fix_small_values) then
    !  call twofl_handle_small_values(.true., w, x, ixI^L, ixO^L, 'twofl_to_conserved')
    !end if
 
    call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
    call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
 
    ! Calculate total energy from pressure, kinetic and magnetic energy
    if(phys_energy) then
      if(phys_internal_e) then
        w(ixo^s,e_n_)=w(ixo^s,e_n_)*inv_gamma_1
        w(ixo^s,e_c_)=w(ixo^s,e_c_)*inv_gamma_1
      else
        w(ixo^s,e_n_)=w(ixo^s,e_n_)*inv_gamma_1&
                   +half*sum(w(ixo^s,mom_n(:))**2,dim=ndim+1)*rhon(ixo^s)
        if(phys_total_energy) then
          w(ixo^s,e_c_)=w(ixo^s,e_c_)*inv_gamma_1&
                   +half*sum(w(ixo^s,mom_c(:))**2,dim=ndim+1)*rhoc(ixo^s)&
                   +twofl_mag_en(w, ixi^l, ixo^l)
        else
          ! kinetic energy + internal energy is evolved   
          w(ixo^s,e_c_)=w(ixo^s,e_c_)*inv_gamma_1&
                   +half*sum(w(ixo^s,mom_c(:))**2,dim=ndim+1)*rhoc(ixo^s)
        end if
      end if
    end if
 
    ! Convert velocity to momentum
    do idir = 1, ndir
       w(ixo^s, mom_n(idir)) = rhon(ixo^s) * w(ixo^s, mom_n(idir))
       w(ixo^s, mom_c(idir)) = rhoc(ixo^s) * w(ixo^s, mom_c(idir))
    end do
  end subroutine twofl_to_conserved
 
  !> Transform conservative variables into primitive ones
  subroutine twofl_to_primitive(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
    integer                         :: idir
    double precision                :: rhoc(ixi^s)
    double precision                :: rhon(ixi^s)
 
    if (fix_small_values) then
      call twofl_handle_small_values(.false., w, x, ixi^l, ixo^l, 'twofl_to_primitive')
    end if
 
    call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
    call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
 
    if(phys_energy) then
      if(phys_internal_e) then
        w(ixo^s,e_n_)=w(ixo^s,e_n_)*gamma_1
        w(ixo^s,e_c_)=w(ixo^s,e_c_)*gamma_1
      else
        ! neutrals evolved energy = ke + e_int 
        w(ixo^s,e_n_)=gamma_1*(w(ixo^s,e_n_)&
                    -twofl_kin_en_n(w,ixi^l,ixo^l))
        ! charges
        if(phys_total_energy) then
         ! evolved energy = ke + e_int + e_mag 
          w(ixo^s,e_c_)=gamma_1*(w(ixo^s,e_c_)&
                    -twofl_kin_en_c(w,ixi^l,ixo^l)&
                    -twofl_mag_en(w,ixi^l,ixo^l))
        else
         ! evolved energy = ke + e_int 
          w(ixo^s,e_c_)=gamma_1*(w(ixo^s,e_c_)&
                    -twofl_kin_en_c(w,ixi^l,ixo^l))
        end if
      end if
    end if
 
    ! Convert momentum to velocity
    do idir = 1, ndir
       w(ixo^s, mom_c(idir)) = w(ixo^s, mom_c(idir))/rhoc(ixo^s)
       w(ixo^s, mom_n(idir)) = w(ixo^s, mom_n(idir))/rhon(ixo^s)
    end do
 
  end subroutine twofl_to_primitive
 
!!USED IN TC
  !> Transform internal energy to total energy
  subroutine twofl_ei_to_e_c(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(inout) :: w(ixI^S, nw)
    double precision, intent(in)    :: x(ixI^S, 1:ndim)
 
    ! Calculate total energy from internal, kinetic and magnetic energy
    if(twofl_eq_energy == eq_energy_ki) then
      w(ixo^s,e_c_)=w(ixo^s,e_c_)&
                +twofl_kin_en_c(w,ixi^l,ixo^l)
    else
      w(ixo^s,e_c_)=w(ixo^s,e_c_)&
               +twofl_kin_en_c(w,ixi^l,ixo^l)&
               +twofl_mag_en(w,ixi^l,ixo^l)
    endif
  end subroutine twofl_ei_to_e_c
 
  !> Transform total energy to internal energy
  subroutine twofl_e_to_ei_c(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(inout) :: w(ixI^S, nw)
    double precision, intent(in)    :: x(ixI^S, 1:ndim)
 
    if(twofl_eq_energy == eq_energy_ki) then
      w(ixo^s,e_c_)=w(ixo^s,e_c_)&
                -twofl_kin_en_c(w,ixi^l,ixo^l)
    else
    ! Calculate ei = e - ek - eb
      w(ixo^s,e_c_)=w(ixo^s,e_c_)&
                -twofl_kin_en_c(w,ixi^l,ixo^l)&
                -twofl_mag_en(w,ixi^l,ixo^l)
    endif
  end subroutine twofl_e_to_ei_c
 
   !Neutrals 
  subroutine twofl_ei_to_e_n(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(inout) :: w(ixI^S, nw)
    double precision, intent(in)    :: x(ixI^S, 1:ndim)
 
    ! Calculate total energy from internal and kinetic  energy
 
    w(ixo^s,e_n_)=w(ixo^s,e_n_)+twofl_kin_en_n(w,ixi^l,ixo^l)
 
  end subroutine twofl_ei_to_e_n
 
  !> Transform total energy to internal energy
  subroutine twofl_e_to_ei_n(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(inout) :: w(ixI^S, nw)
    double precision, intent(in)    :: x(ixI^S, 1:ndim)
 
    ! Calculate ei = e - ek 
    w(ixo^s,e_n_)=w(ixo^s,e_n_)-twofl_kin_en_n(w,ixi^l,ixo^l)
 
    call twofl_handle_small_ei_n(w,x,ixi^l,ixo^l,e_n_,"e_to_ei_n")
  end subroutine twofl_e_to_ei_n
 
  subroutine twofl_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
    use mod_global_parameters
    use mod_small_values
    logical, intent(in)             :: primitive
    integer, intent(in)             :: ixI^L,ixO^L
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    character(len=*), intent(in)    :: subname
 
    integer :: idir
    logical :: flag(ixI^S,1:nw)
    double precision :: tmp2(ixI^S)
    double precision :: tmp1(ixI^S)
 
    call twofl_check_w(primitive, ixi^l, ixo^l, w, flag)
 
    if(any(flag)) then
      select case (small_values_method)
      case ("replace")
        if(has_equi_rho_c0) then
          where(flag(ixo^s,rho_c_)) w(ixo^s,rho_c_) = &
                    small_density-block%equi_vars(ixo^s,equi_rho_c0_,0)
        else
          where(flag(ixo^s,rho_c_)) w(ixo^s,rho_c_) = small_density
        end if
        if(has_equi_rho_n0) then
          where(flag(ixo^s,rho_n_)) w(ixo^s,rho_n_) = &
                  small_density-block%equi_vars(ixo^s,equi_rho_n0_,0)
        else
          where(flag(ixo^s,rho_n_)) w(ixo^s,rho_n_) = small_density
        end if
        do idir = 1, ndir
          if(small_values_fix_iw(mom_n(idir))) then
            where(flag(ixo^s,rho_n_)) w(ixo^s, mom_n(idir)) = 0.0d0
          end if
          if(small_values_fix_iw(mom_c(idir))) then
            where(flag(ixo^s,rho_c_)) w(ixo^s, mom_c(idir)) = 0.0d0
          end if
        end do
 
        if(phys_energy) then
          if(primitive) then
            if(has_equi_pe_n0) then 
             tmp1(ixo^s) = small_pressure - &
               block%equi_vars(ixo^s,equi_pe_n0_,0)
            else
             tmp1(ixo^s) = small_pressure
            end if  
            if(has_equi_pe_c0) then 
             tmp2(ixo^s) = small_e - &
               block%equi_vars(ixo^s,equi_pe_c0_,0)
            else
             tmp2(ixo^s) = small_pressure
            end if
          else
            ! conserved  
            if(has_equi_pe_n0) then 
              tmp1(ixo^s) = small_e - &
              block%equi_vars(ixo^s,equi_pe_n0_,0)*inv_gamma_1 
            else
              tmp1(ixo^s) = small_e
            end if  
            if(has_equi_pe_c0) then 
              tmp2(ixo^s) = small_e - &
                block%equi_vars(ixo^s,equi_pe_c0_,0)*inv_gamma_1 
            else
              tmp2(ixo^s) = small_e
            end if  
            if(phys_internal_e) then
              where(flag(ixo^s,e_n_))
                w(ixo^s,e_n_)=tmp1(ixo^s)
              end where
              where(flag(ixo^s,e_c_))
                w(ixo^s,e_c_)=tmp2(ixo^s)
              end where
            else
              where(flag(ixo^s,e_n_))
                w(ixo^s,e_n_) = tmp1(ixo^s)+&
                 twofl_kin_en_n(w,ixi^l,ixo^l)
              end where
              if(phys_total_energy) then
                where(flag(ixo^s,e_c_))
                  w(ixo^s,e_c_) = tmp2(ixo^s)+&
                   twofl_kin_en_c(w,ixi^l,ixo^l)+&
                   twofl_mag_en(w,ixi^l,ixo^l)
                end where
              else
                where(flag(ixo^s,e_c_))
                  w(ixo^s,e_c_) = tmp2(ixo^s)+&
                   twofl_kin_en_c(w,ixi^l,ixo^l)
                end where
              end if
            end if
          end if
        end if
      case ("average")
        call small_values_average(ixi^l, ixo^l, w, x, flag)
      case default
        if(.not.primitive) then
          !convert w to primitive
          ! Calculate pressure = (gamma-1) * (e-ek-eb)
          if(phys_energy) then
            if(phys_internal_e) then
              w(ixo^s,e_c_)=w(ixo^s,e_c_)*gamma_1
              w(ixo^s,e_n_)=w(ixo^s,e_n_)*gamma_1
            else
              w(ixo^s,e_n_)=gamma_1*(w(ixo^s,e_n_)&
                         -twofl_kin_en_n(w,ixi^l,ixo^l))
              if(phys_total_energy) then
                w(ixo^s,e_c_)=gamma_1*(w(ixo^s,e_c_)&
                          -twofl_kin_en_c(w,ixi^l,ixo^l)&
                          -twofl_mag_en(w,ixi^l,ixo^l))
               else
                 w(ixo^s,e_c_)=gamma_1*(w(ixo^s,e_c_)&
                          -twofl_kin_en_c(w,ixi^l,ixo^l))
 
               end if  
            end if
          end if
          ! Convert momentum to velocity
          if(has_equi_rho_n0) then
            tmp1(ixo^s) = w(ixo^s,rho_n_) + block%equi_vars(ixo^s,equi_rho_n0_,0)
          else  
            tmp1(ixo^s) = w(ixo^s,rho_n_) 
          end if
    
          if(has_equi_rho_c0) then
            tmp2(ixo^s) = w(ixo^s,rho_c_) + block%equi_vars(ixo^s,equi_rho_c0_,0)
          else  
            tmp2(ixo^s) = w(ixo^s,rho_c_) 
          end if
          do idir = 1, ndir
             w(ixo^s, mom_n(idir)) = w(ixo^s, mom_n(idir))/tmp1(ixo^s)
             w(ixo^s, mom_c(idir)) = w(ixo^s, mom_c(idir))/tmp2(ixo^s)
          end do
        end if
        call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
      end select
    end if
  end subroutine twofl_handle_small_values
 
  !> Calculate cmax_idim=csound+abs(v_idim) within ixO^L
  subroutine twofl_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
    use mod_global_parameters
 
    integer, intent(in)          :: ixI^L, ixO^L, idim
    double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: cmax(ixI^S)
    double precision                :: vc(ixI^S)
    double precision                :: cmax2(ixI^S)
    double precision                :: vn(ixI^S)
 
    call twofl_get_csound_c_idim(w,x,ixi^l,ixo^l,idim,cmax)
    call twofl_get_v_c_idim(w,x,ixi^l,ixo^l,idim,vc)
    call twofl_get_v_n_idim(w,x,ixi^l,ixo^l,idim,vn)
    call twofl_get_csound_n(w,x,ixi^l,ixo^l,cmax2)
    cmax(ixo^s)=max(abs(vn(ixo^s))+cmax2(ixo^s),& 
            abs(vc(ixo^s))+cmax(ixo^s))
        
  end subroutine twofl_get_cmax
 
  subroutine twofl_get_a2max(w,x,ixI^L,ixO^L,a2max)
    use mod_global_parameters
 
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: a2max(ndim)
    double precision :: a2(ixI^S,ndim,nw)
    integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
 
    a2=zero
    do i = 1,ndim
      !> 4th order
      hxo^l=ixo^l-kr(i,^d);
      gxo^l=hxo^l-kr(i,^d);
      jxo^l=ixo^l+kr(i,^d);
      kxo^l=jxo^l+kr(i,^d);
      a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
         -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
      a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/dxlevel(i)**2
    end do
  end subroutine twofl_get_a2max
 
  ! COPIED from hd/moh_hd_phys
  !> get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
  subroutine twofl_get_tcutoff_n(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
    use mod_global_parameters
    integer, intent(in) :: ixI^L,ixO^L
    double precision, intent(in) :: x(ixI^S,1:ndim),w(ixI^S,1:nw)
    double precision, intent(out) :: tco_local, Tmax_local
 
    double precision, parameter :: delta=0.25d0
    double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
    integer :: jxO^L,hxO^L
    logical :: lrlt(ixI^S)
 
    {^ifoned
    ! reuse lts as rhon
    call get_rhon_tot(w,x,ixi^l,ixi^l,lts)
    tmp1(ixi^s)=w(ixi^s,e_n_)-0.5d0*sum(w(ixi^s,mom_n(:))**2,dim=ndim+1)/lts(ixi^s)
    te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(twofl_gamma-1.d0)
 
    tmax_local=maxval(te(ixo^s))
 
    hxo^l=ixo^l-1;
    jxo^l=ixo^l+1;
    lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
    lrlt=.false.
    where(lts(ixo^s) > delta)
      lrlt(ixo^s)=.true.
    end where
    tco_local=zero
    if(any(lrlt(ixo^s))) then
      tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
    end if
    }
  end subroutine twofl_get_tcutoff_n
 
  !> get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
  subroutine twofl_get_tcutoff_c(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in) :: ixI^L,ixO^L
    double precision, intent(in) :: x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision, intent(out) :: Tco_local,Tmax_local
 
    double precision, parameter :: trac_delta=0.25d0
    double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
    double precision, dimension(ixI^S,1:ndir) :: bunitvec
    double precision, dimension(ixI^S,1:ndim) :: gradT
    double precision :: Bdir(ndim)
    double precision :: ltr(ixI^S),ltrc,ltrp,altr(ixI^S)
    integer :: idims,jxO^L,hxO^L,ixA^D,ixB^D
    integer :: jxP^L,hxP^L,ixP^L
    logical :: lrlt(ixI^S)
 
    ! reuse lts as rhoc
    call get_rhoc_tot(w,x,ixi^l,ixi^l,lts)
    if(phys_internal_e) then
      tmp1(ixi^s)=w(ixi^s,e_c_)
    else
      tmp1(ixi^s)=w(ixi^s,e_c_)-0.5d0*(sum(w(ixi^s,mom_c(:))**2,dim=ndim+1)/&
                       lts(ixi^s)+sum(w(ixi^s,mag(:))**2,dim=ndim+1))
    end if
    te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(twofl_gamma-1.d0)
    tmax_local=maxval(te(ixo^s))
 
    {^ifoned
    select case(twofl_trac_type)
    case(0)
      !> test case, fixed cutoff temperature
      w(ixi^s,tcoff_c_)=2.5d5/unit_temperature
    case(1)
      hxo^l=ixo^l-1;
      jxo^l=ixo^l+1;
      lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
      lrlt=.false.
      where(lts(ixo^s) > trac_delta)
        lrlt(ixo^s)=.true.
      end where
      if(any(lrlt(ixo^s))) then
        tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
      end if
    case(2)
      !> iijima et al. 2021, LTRAC method
      ltrc=1.5d0
      ltrp=2.5d0
      ixp^l=ixo^l^ladd1;
      hxo^l=ixo^l-1;
      jxo^l=ixo^l+1;
      hxp^l=ixp^l-1;
      jxp^l=ixp^l+1;
      lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
      ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
      w(ixo^s,tcoff_c_)=te(ixo^s)*&
        (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
    case default
      call mpistop("twofl_trac_type not allowed for 1D simulation")
    end select
    }
    {^nooned
    select case(twofl_trac_type)
    case(0)
      !> test case, fixed cutoff temperature
      w(ixi^s,tcoff_c_)=2.5d5/unit_temperature
    case(1,4,6)
      ! temperature gradient at cell centers
      do idims=1,ndim
        call gradient(te,ixi^l,ixo^l,idims,tmp1)
        gradt(ixo^s,idims)=tmp1(ixo^s)
      end do
      ! B vector
      if(b0field) then
        bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+block%B0(ixo^s,:,0)
      else
        bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
      end if
      if(twofl_trac_type .gt. 1) then
        ! B direction at cell center
        bdir=zero
        {do ixa^d=0,1\}
          ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
          bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
        {end do\}
        if(sum(bdir(:)**2) .gt. zero) then
          bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
        end if
        block%special_values(3:ndim+2)=bdir(1:ndim)
      end if
      tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
      where(tmp1(ixo^s)/=0.d0)
        tmp1(ixo^s)=1.d0/tmp1(ixo^s)
      elsewhere
        tmp1(ixo^s)=bigdouble
      end where
      ! b unit vector: magnetic field direction vector
      do idims=1,ndim
        bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
      end do
      ! temperature length scale inversed
      lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
      ! fraction of cells size to temperature length scale
      if(slab_uniform) then
        lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
      else
        lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
      end if
      lrlt=.false.
      where(lts(ixo^s) > trac_delta)
        lrlt(ixo^s)=.true.
      end where
      if(any(lrlt(ixo^s))) then
        block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
      else
        block%special_values(1)=zero
      end if
      block%special_values(2)=tmax_local
    case(2)
      !> iijima et al. 2021, LTRAC method
      ltrc=1.5d0
      ltrp=4.d0
      ixp^l=ixo^l^ladd1;
      ! temperature gradient at cell centers
      do idims=1,ndim
        call gradient(te,ixi^l,ixp^l,idims,tmp1)
        gradt(ixp^s,idims)=tmp1(ixp^s)
      end do
      ! B vector
      if(b0field) then
        bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
      else
        bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
      end if
      tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
      where(tmp1(ixp^s)/=0.d0)
        tmp1(ixp^s)=1.d0/tmp1(ixp^s)
      elsewhere
        tmp1(ixp^s)=bigdouble
      end where
      ! b unit vector: magnetic field direction vector
      do idims=1,ndim
        bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
      end do
      ! temperature length scale inversed
      lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
      ! fraction of cells size to temperature length scale
      if(slab_uniform) then
        lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
      else
        lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
      end if
      ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
  
      altr(ixi^s)=zero  
      do idims=1,ndim 
        hxo^l=ixo^l-kr(idims,^d);
        jxo^l=ixo^l+kr(idims,^d);
        altr(ixo^s)=altr(ixo^s) &
          +0.25*(ltr(hxo^s)+two*ltr(ixo^s)+ltr(jxo^s))*bunitvec(ixo^s,idims)**2
        w(ixo^s,tcoff_c_)=te(ixo^s)*altr(ixo^s)**(0.4*ltrp)
      end do
    case(3,5)
      !> do nothing here
    case default
      call mpistop("unknown twofl_trac_type")
    end select 
    }
  end subroutine twofl_get_tcutoff_c
 
  !> get H speed for H-correction to fix the carbuncle problem at grid-aligned shock front
  subroutine twofl_get_h_speed_one(wprim,x,ixI^L,ixO^L,idim,Hspeed) 
    use mod_global_parameters
 
    integer, intent(in)             :: ixI^L, ixO^L, idim
    double precision, intent(in)    :: wprim(ixI^S, nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    double precision, intent(out)   :: Hspeed(ixI^S,1:number_species)
 
    double precision :: csound(ixI^S,ndim),tmp(ixI^S)
    integer :: jxC^L, ixC^L, ixA^L, id, ix^D
 
    hspeed=0.d0
    ixa^l=ixo^l^ladd1;
    do id=1,ndim
      call twofl_get_csound_prim(wprim,x,ixi^l,ixa^l,id,tmp)
      csound(ixa^s,id)=tmp(ixa^s)
    end do
    ixcmax^d=ixomax^d;
    ixcmin^d=ixomin^d+kr(idim,^d)-1;
    jxcmax^d=ixcmax^d+kr(idim,^d);
    jxcmin^d=ixcmin^d+kr(idim,^d);
    hspeed(ixc^s,1)=0.5d0*abs(&
      0.5d0 * (wprim(jxc^s,mom_c(idim))+ wprim(jxc^s,mom_n(idim))) &
      +csound(jxc^s,idim)- &
      0.5d0 * (wprim(ixc^s,mom_c(idim)) + wprim(ixc^s,mom_n(idim)))&
      +csound(ixc^s,idim))
 
    do id=1,ndim
      if(id==idim) cycle
      ixamax^d=ixcmax^d+kr(id,^d);
      ixamin^d=ixcmin^d+kr(id,^d);
      hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
        0.5d0 * (wprim(ixa^s,mom_c(id)) + wprim(ixa^s,mom_n(id)))&
        +csound(ixa^s,id)-&
        0.5d0 * (wprim(ixc^s,mom_c(id)) + wprim(ixc^s,mom_n(id)))&
        +csound(ixc^s,id)))
 
 
      ixamax^d=ixcmax^d-kr(id,^d);
      ixamin^d=ixcmin^d-kr(id,^d);
      hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
        0.5d0 * (wprim(ixc^s,mom_c(id)) + wprim(ixc^s,mom_n(id)))&
        +csound(ixc^s,id)-&
        0.5d0 * (wprim(ixa^s,mom_c(id)) + wprim(ixa^s,mom_n(id)))&
        +csound(ixa^s,id)))
 
    end do
 
    do id=1,ndim
      if(id==idim) cycle
      ixamax^d=jxcmax^d+kr(id,^d);
      ixamin^d=jxcmin^d+kr(id,^d);
      hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
        0.5d0 * (wprim(ixa^s,mom_c(id)) + wprim(ixa^s,mom_n(id)))&
        +csound(ixa^s,id)-&
        0.5d0 * (wprim(jxc^s,mom_c(id)) + wprim(jxc^s,mom_n(id)))&
        +csound(jxc^s,id)))
      ixamax^d=jxcmax^d-kr(id,^d);
      ixamin^d=jxcmin^d-kr(id,^d);
      hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
        0.5d0 * (wprim(jxc^s,mom_c(id)) + wprim(jxc^s,mom_n(id)))&
        +csound(jxc^s,id)-&
        0.5d0 * (wprim(ixa^s,mom_c(id)) + wprim(ixa^s,mom_n(id)))&
        +csound(ixa^s,id)))
    end do
 
  end subroutine twofl_get_h_speed_one
 
  !> get H speed for H-correction to fix the carbuncle problem at grid-aligned shock front
  subroutine twofl_get_h_speed_species(wprim,x,ixI^L,ixO^L,idim,Hspeed) 
    use mod_global_parameters
 
    integer, intent(in)             :: ixI^L, ixO^L, idim
    double precision, intent(in)    :: wprim(ixI^S, nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    double precision, intent(out)   :: Hspeed(ixI^S,1:number_species)
 
    double precision :: csound(ixI^S,ndim),tmp(ixI^S)
    integer :: jxC^L, ixC^L, ixA^L, id, ix^D
 
    hspeed=0.d0
    ! charges
    ixa^l=ixo^l^ladd1;
    do id=1,ndim
      call twofl_get_csound_prim_c(wprim,x,ixi^l,ixa^l,id,tmp)
      csound(ixa^s,id)=tmp(ixa^s)
    end do
    ixcmax^d=ixomax^d;
    ixcmin^d=ixomin^d+kr(idim,^d)-1;
    jxcmax^d=ixcmax^d+kr(idim,^d);
    jxcmin^d=ixcmin^d+kr(idim,^d);
    hspeed(ixc^s,1)=0.5d0*abs(wprim(jxc^s,mom_c(idim))+csound(jxc^s,idim)-wprim(ixc^s,mom_c(idim))+csound(ixc^s,idim))
 
    do id=1,ndim
      if(id==idim) cycle
      ixamax^d=ixcmax^d+kr(id,^d);
      ixamin^d=ixcmin^d+kr(id,^d);
      hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,mom_c(id))+csound(ixa^s,id)-wprim(ixc^s,mom_c(id))+csound(ixc^s,id)))
      ixamax^d=ixcmax^d-kr(id,^d);
      ixamin^d=ixcmin^d-kr(id,^d);
      hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixc^s,mom_c(id))+csound(ixc^s,id)-wprim(ixa^s,mom_c(id))+csound(ixa^s,id)))
    end do
 
    do id=1,ndim
      if(id==idim) cycle
      ixamax^d=jxcmax^d+kr(id,^d);
      ixamin^d=jxcmin^d+kr(id,^d);
      hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,mom_c(id))+csound(ixa^s,id)-wprim(jxc^s,mom_c(id))+csound(jxc^s,id)))
      ixamax^d=jxcmax^d-kr(id,^d);
      ixamin^d=jxcmin^d-kr(id,^d);
      hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(jxc^s,mom_c(id))+csound(jxc^s,id)-wprim(ixa^s,mom_c(id))+csound(ixa^s,id)))
    end do
 
    ! neutrals
    ixa^l=ixo^l^ladd1;
    do id=1,ndim
      call twofl_get_csound_prim_n(wprim,x,ixi^l,ixa^l,id,tmp)
      csound(ixa^s,id)=tmp(ixa^s)
    end do
    ixcmax^d=ixomax^d;
    ixcmin^d=ixomin^d+kr(idim,^d)-1;
    jxcmax^d=ixcmax^d+kr(idim,^d);
    jxcmin^d=ixcmin^d+kr(idim,^d);
    hspeed(ixc^s,2)=0.5d0*abs(wprim(jxc^s,mom_n(idim))+csound(jxc^s,idim)-wprim(ixc^s,mom_n(idim))+csound(ixc^s,idim))
 
    do id=1,ndim
      if(id==idim) cycle
      ixamax^d=ixcmax^d+kr(id,^d);
      ixamin^d=ixcmin^d+kr(id,^d);
      hspeed(ixc^s,2)=max(hspeed(ixc^s,2),0.5d0*abs(wprim(ixa^s,mom_n(id))+csound(ixa^s,id)-wprim(ixc^s,mom_n(id))+csound(ixc^s,id)))
      ixamax^d=ixcmax^d-kr(id,^d);
      ixamin^d=ixcmin^d-kr(id,^d);
      hspeed(ixc^s,2)=max(hspeed(ixc^s,2),0.5d0*abs(wprim(ixc^s,mom_n(id))+csound(ixc^s,id)-wprim(ixa^s,mom_n(id))+csound(ixa^s,id)))
    end do
 
    do id=1,ndim
      if(id==idim) cycle
      ixamax^d=jxcmax^d+kr(id,^d);
      ixamin^d=jxcmin^d+kr(id,^d);
      hspeed(ixc^s,2)=max(hspeed(ixc^s,2),0.5d0*abs(wprim(ixa^s,mom_n(id))+csound(ixa^s,id)-wprim(jxc^s,mom_n(id))+csound(jxc^s,id)))
      ixamax^d=jxcmax^d-kr(id,^d);
      ixamin^d=jxcmin^d-kr(id,^d);
      hspeed(ixc^s,2)=max(hspeed(ixc^s,2),0.5d0*abs(wprim(jxc^s,mom_n(id))+csound(jxc^s,id)-wprim(ixa^s,mom_n(id))+csound(ixa^s,id)))
    end do
 
  end subroutine twofl_get_h_speed_species
 
  !> Estimating bounds for the minimum and maximum signal velocities
  subroutine twofl_get_cbounds_one(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
    use mod_global_parameters
    use mod_constrained_transport
 
    integer, intent(in)             :: ixI^L, ixO^L, idim
    double precision, intent(in)    :: wLC(ixI^S, nw), wRC(ixI^S, nw)
    double precision, intent(in)    :: wLp(ixI^S, nw), wRp(ixI^S, nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    double precision, intent(inout) :: cmax(ixI^S,number_species)
    double precision, intent(inout), optional :: cmin(ixI^S,number_species)
    double precision, intent(in)    :: Hspeed(ixI^S,1:number_species)
 
    double precision :: wmean(ixI^S,nw)
    double precision :: rhon(ixI^S)
    double precision :: rhoc(ixI^S)
    double precision, dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
    integer :: ix^D
 
    select case (boundspeed) 
    case (1)
      ! This implements formula (10.52) from "Riemann Solvers and Numerical
      ! Methods for Fluid Dynamics" by Toro.
      call get_rhoc_tot(wlp,x,ixi^l,ixo^l,rhoc)
      call get_rhon_tot(wlp,x,ixi^l,ixo^l,rhon)
      tmp1(ixo^s)=sqrt(abs(rhoc(ixo^s)  +rhon(ixo^s)))
 
      call get_rhoc_tot(wrp,x,ixi^l,ixo^l,rhoc)
      call get_rhon_tot(wrp,x,ixi^l,ixo^l,rhon)
      tmp2(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
 
      tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
      umean(ixo^s)=(0.5*(wlp(ixo^s,mom_n(idim))+wlp(ixo^s,mom_c(idim)))*tmp1(ixo^s) + &
                    0.5*(wrp(ixo^s,mom_n(idim))+wrp(ixo^s,mom_c(idim)))*tmp2(ixo^s))*tmp3(ixo^s)
      call twofl_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
      call twofl_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
 
      dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
       0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*(&
       0.5*(wrp(ixo^s,mom_n(idim))+wrp(ixo^s,mom_c(idim)))- & 
       0.5*(wlp(ixo^s,mom_n(idim))+wlp(ixo^s,mom_c(idim))))**2
      dmean(ixo^s)=sqrt(dmean(ixo^s))
      if(present(cmin)) then
        cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
        cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
        if(h_correction) then
          {do ix^db=ixomin^db,ixomax^db\}
            cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
            cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
          {end do\}
        end if
      else
        cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
      end if
    case (2)
    ! typeboundspeed=='cmaxmean'
      wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
      call get_rhon_tot(wmean,x,ixi^l,ixo^l,rhon)
      tmp2(ixo^s)=wmean(ixo^s,mom_n(idim))/rhon(ixo^s)
      call get_rhoc_tot(wmean,x,ixi^l,ixo^l,rhoc)
      tmp1(ixo^s)=wmean(ixo^s,mom_c(idim))/rhoc(ixo^s)
      call twofl_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
      if(present(cmin)) then
        cmax(ixo^s,1)=max(max(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) +csoundr(ixo^s),zero)
        cmin(ixo^s,1)=min(min(abs(tmp2(ixo^s)),  abs(tmp1(ixo^s)) ) -csoundr(ixo^s),zero)
        if(h_correction) then
          {do ix^db=ixomin^db,ixomax^db\}
            cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
            cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
          {end do\}
        end if
      else
        cmax(ixo^s,1)= max(abs(tmp2(ixo^s)),abs(tmp1(ixo^s)))+csoundr(ixo^s)
      end if
    case (3)
      ! Miyoshi 2005 JCP 208, 315 equation (67)
      call twofl_get_csound(wlp,x,ixi^l,ixo^l,idim,csoundl)
      call twofl_get_csound(wrp,x,ixi^l,ixo^l,idim,csoundr)
      csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
      if(present(cmin)) then
        cmin(ixo^s,1)=min(0.5*(wlp(ixo^s,mom_c(idim))+ wrp(ixo^s,mom_n(idim))),&
            0.5*(wrp(ixo^s,mom_c(idim))+ wrp(ixo^s,mom_n(idim))))-csoundl(ixo^s)
        cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,mom_c(idim))+ wrp(ixo^s,mom_n(idim))),&
            0.5*(wrp(ixo^s,mom_c(idim))+ wrp(ixo^s,mom_n(idim))))+csoundl(ixo^s)
        if(h_correction) then
          {do ix^db=ixomin^db,ixomax^db\}
            cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
            cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
          {end do\}
        end if
      else
        cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,mom_c(idim))+ wrp(ixo^s,mom_n(idim))),&
          0.5*(wrp(ixo^s,mom_c(idim))+ wrp(ixo^s,mom_n(idim))))+csoundl(ixo^s)
      end if
    end select
 
  end subroutine twofl_get_cbounds_one
 
  !> Calculate fast magnetosonic wave speed
  subroutine twofl_get_csound_prim_c(w,x,ixI^L,ixO^L,idim,csound)
    use mod_global_parameters
 
    integer, intent(in)          :: ixI^L, ixO^L, idim
    double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
    double precision, intent(out):: csound(ixI^S)
    double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
    double precision :: inv_rho(ixO^S)
    double precision :: rhoc(ixI^S)
 
    integer :: ix1,ix2
 
 
    call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
    inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
 
    if(phys_energy) then
      call twofl_get_pthermal_c_primitive(w,x,ixi^l,ixo^l,csound)
      csound(ixo^s)=twofl_gamma*csound(ixo^s)/rhoc(ixo^s)
    else
      call twofl_get_csound2_adiab_c(w,x,ixi^l,ixo^l,csound)
    endif
 
    ! store |B|^2 in v
    b2(ixo^s)        = twofl_mag_en_all(w,ixi^l,ixo^l)
    cfast2(ixo^s)   = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
    avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
         * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
         * inv_rho(ixo^s)
 
    where(avmincs2(ixo^s)<zero)
       avmincs2(ixo^s)=zero
    end where
 
    avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
 
    if (.not. twofl_hall) then
       csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
    else
       ! take the Hall velocity into account:
       ! most simple estimate, high k limit:
       ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
       kmax = dpi/min({dxlevel(^d)},bigdouble)*half
       csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
            twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
    end if
 
  end subroutine twofl_get_csound_prim_c
 
  !> Calculate fast magnetosonic wave speed
  subroutine twofl_get_csound_prim_n(w,x,ixI^L,ixO^L,idim,csound)
    use mod_global_parameters
 
    integer, intent(in)          :: ixI^L, ixO^L, idim
    double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
    double precision, intent(out):: csound(ixI^S)
    double precision :: rhon(ixI^S)
 
    if(phys_energy) then
      call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
      call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,csound)
      csound(ixo^s)=twofl_gamma*csound(ixo^s)/rhon(ixo^s)
    else
      call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound)
    endif
    csound(ixo^s) = sqrt(csound(ixo^s))
 
  end subroutine twofl_get_csound_prim_n
 
  !> Estimating bounds for the minimum and maximum signal velocities
  subroutine twofl_get_cbounds_species(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
    use mod_global_parameters
    use mod_constrained_transport
    use mod_variables
 
    integer, intent(in)             :: ixI^L, ixO^L, idim
    double precision, intent(in)    :: wLC(ixI^S, nw), wRC(ixI^S, nw)
    double precision, intent(in)    :: wLp(ixI^S, nw), wRp(ixI^S, nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    double precision, intent(inout) :: cmax(ixI^S,1:number_species)
    double precision, intent(inout), optional :: cmin(ixI^S,1:number_species)
    double precision, intent(in)    :: Hspeed(ixI^S,1:number_species)
 
    double precision :: wmean(ixI^S,nw)
    double precision :: rho(ixI^S)
    double precision, dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
    integer :: ix^D
 
    select case (boundspeed) 
    case (1)
      ! This implements formula (10.52) from "Riemann Solvers and Numerical
      ! Methods for Fluid Dynamics" by Toro.
      ! charges
      call get_rhoc_tot(wlp,x,ixi^l,ixo^l,rho)
      tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
 
      call get_rhoc_tot(wrp,x,ixi^l,ixo^l,rho)
      tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
 
      tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
      umean(ixo^s)=(wlp(ixo^s,mom_c(idim))*tmp1(ixo^s)+wrp(ixo^s,mom_c(idim))*tmp2(ixo^s))*tmp3(ixo^s)
      call twofl_get_csound_prim_c(wlp,x,ixi^l,ixo^l,idim,csoundl)
      call twofl_get_csound_prim_c(wrp,x,ixi^l,ixo^l,idim,csoundr)
 
 
      dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
       0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
       (wrp(ixo^s,mom_c(idim)) - wlp(ixo^s,mom_c(idim)))**2
      dmean(ixo^s)=sqrt(dmean(ixo^s))
      if(present(cmin)) then
        cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
        cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
        if(h_correction) then
          {do ix^db=ixomin^db,ixomax^db\}
            cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
            cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
          {end do\}
        end if
      else
        cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
      end if
    
      ! neutrals
 
      call get_rhon_tot(wlp,x,ixi^l,ixo^l,rho)
      tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
 
      call get_rhon_tot(wrp,x,ixi^l,ixo^l,rho)
      tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
 
      tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
      umean(ixo^s)=(wlp(ixo^s,mom_n(idim))*tmp1(ixo^s)+wrp(ixo^s,mom_n(idim))*tmp2(ixo^s))*tmp3(ixo^s)
      call twofl_get_csound_prim_n(wlp,x,ixi^l,ixo^l,idim,csoundl)
      call twofl_get_csound_prim_n(wrp,x,ixi^l,ixo^l,idim,csoundr)
 
 
      dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
       0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
       (wrp(ixo^s,mom_n(idim)) - wlp(ixo^s,mom_n(idim)))**2
      dmean(ixo^s)=sqrt(dmean(ixo^s))
      if(present(cmin)) then
        cmin(ixo^s,2)=umean(ixo^s)-dmean(ixo^s)
        cmax(ixo^s,2)=umean(ixo^s)+dmean(ixo^s)
        if(h_correction) then
          {do ix^db=ixomin^db,ixomax^db\}
            cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
            cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
          {end do\}
        end if
      else
        cmax(ixo^s,2)=abs(umean(ixo^s))+dmean(ixo^s)
      end if
 
    case (2)
    ! typeboundspeed=='cmaxmean'
      wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
     ! charges 
 
      call get_rhoc_tot(wmean,x,ixi^l,ixo^l,rho)
      tmp1(ixo^s)=wmean(ixo^s,mom_c(idim))/rho(ixo^s)
      call twofl_get_csound_c_idim(wmean,x,ixi^l,ixo^l,idim,csoundr)
      if(present(cmin)) then
        cmax(ixo^s,1)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
        cmin(ixo^s,1)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
        if(h_correction) then
          {do ix^db=ixomin^db,ixomax^db\}
            cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
            cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
          {end do\}
        end if
      else
        cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
      end if
      !neutrals
      
      call get_rhon_tot(wmean,x,ixi^l,ixo^l,rho)
      tmp1(ixo^s)=wmean(ixo^s,mom_n(idim))/rho(ixo^s)
      call twofl_get_csound_n(wmean,x,ixi^l,ixo^l,csoundr)
      if(present(cmin)) then
        cmax(ixo^s,2)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
        cmin(ixo^s,2)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
        if(h_correction) then
          {do ix^db=ixomin^db,ixomax^db\}
            cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
            cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
          {end do\}
        end if
      else
        cmax(ixo^s,2)= abs(tmp1(ixo^s))+csoundr(ixo^s)
      end if
    case (3)
      ! Miyoshi 2005 JCP 208, 315 equation (67)
      call twofl_get_csound_c_idim(wlp,x,ixi^l,ixo^l,idim,csoundl)
      call twofl_get_csound_c_idim(wrp,x,ixi^l,ixo^l,idim,csoundr)
      csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
      if(present(cmin)) then
        cmin(ixo^s,1)=min(wlp(ixo^s,mom_c(idim)),wrp(ixo^s,mom_c(idim)))-csoundl(ixo^s)
        cmax(ixo^s,1)=max(wlp(ixo^s,mom_c(idim)),wrp(ixo^s,mom_c(idim)))+csoundl(ixo^s)
        if(h_correction) then
          {do ix^db=ixomin^db,ixomax^db\}
            cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
            cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
          {end do\}
        end if
      else
        cmax(ixo^s,1)=max(wlp(ixo^s,mom_c(idim)),wrp(ixo^s,mom_c(idim)))+csoundl(ixo^s)
      end if
      call twofl_get_csound_n(wlp,x,ixi^l,ixo^l,csoundl)
      call twofl_get_csound_n(wrp,x,ixi^l,ixo^l,csoundr)
      csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
      if(present(cmin)) then
        cmin(ixo^s,2)=min(wlp(ixo^s,mom_n(idim)),wrp(ixo^s,mom_n(idim)))-csoundl(ixo^s)
        cmax(ixo^s,2)=max(wlp(ixo^s,mom_n(idim)),wrp(ixo^s,mom_n(idim)))+csoundl(ixo^s)
        if(h_correction) then
          {do ix^db=ixomin^db,ixomax^db\}
            cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,1)),hspeed(ix^d,2))
            cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,1)),hspeed(ix^d,2))
          {end do\}
        end if
      else
        cmax(ixo^s,2)=max(wlp(ixo^s,mom_n(idim)),wrp(ixo^s,mom_n(idim)))+csoundl(ixo^s)
      end if
 
    end select
 
  end subroutine twofl_get_cbounds_species
 
  !> prepare velocities for ct methods
  subroutine twofl_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
    use mod_global_parameters
 
    integer, intent(in)             :: ixI^L, ixO^L, idim
    double precision, intent(in)    :: wLp(ixI^S, nw), wRp(ixI^S, nw)
    double precision, intent(in)    :: cmax(ixI^S)
    double precision, intent(in), optional :: cmin(ixI^S)
    type(ct_velocity), intent(inout):: vcts
 
    integer                         :: idimE,idimN
 
    ! calculate velocities related to different UCT schemes
    select case(type_ct)
    case('average')
    case('uct_contact')
      if(.not.allocated(vcts%vnorm)) allocate(vcts%vnorm(ixi^s,1:ndim))
      ! get average normal velocity at cell faces
      vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,mom_c(idim))+wrp(ixo^s,mom_c(idim)))
    case('uct_hll')
      if(.not.allocated(vcts%vbarC)) then
        allocate(vcts%vbarC(ixi^s,1:ndir,2),vcts%vbarLC(ixi^s,1:ndir,2),vcts%vbarRC(ixi^s,1:ndir,2))
        allocate(vcts%cbarmin(ixi^s,1:ndim),vcts%cbarmax(ixi^s,1:ndim)) 
      end if
      ! Store magnitude of characteristics
      if(present(cmin)) then
        vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
        vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
      else
        vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
        vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
      end if
 
      idimn=mod(idim,ndir)+1 ! 'Next' direction
      idime=mod(idim+1,ndir)+1 ! Electric field direction
      ! Store velocities
      vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,mom_c(idimn))
      vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,mom_c(idimn))
      vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
           +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
          /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
 
      vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,mom_c(idime))
      vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,mom_c(idime))
      vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
           +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
          /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
    case default
      call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
    end select
 
  end subroutine twofl_get_ct_velocity
 
  subroutine twofl_get_csound_c_idim(w,x,ixI^L,ixO^L,idim,csound)
    use mod_global_parameters
 
    integer, intent(in)          :: ixI^L, ixO^L, idim
    double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
    double precision, intent(out):: csound(ixI^S)
    double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
    double precision :: inv_rho(ixO^S)
    double precision :: tmp(ixI^S)
#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
    double precision :: rhon(ixI^S)
#endif
    call get_rhoc_tot(w,x,ixi^l,ixo^l,tmp)
#if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
    call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
    inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+tmp(ixo^s)) 
#else
    inv_rho(ixo^s)=1.d0/tmp(ixo^s)
#endif
 
    call twofl_get_csound2_c_from_conserved(w,x,ixi^l,ixo^l,csound)
 
    ! store |B|^2 in v
    b2(ixo^s) = twofl_mag_en_all(w,ixi^l,ixo^l)
 
    cfast2(ixo^s)   = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
    avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
         * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
         * inv_rho(ixo^s)
 
    where(avmincs2(ixo^s)<zero)
       avmincs2(ixo^s)=zero
    end where
 
    avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
 
    if (.not. twofl_hall) then
       csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
    else
       ! take the Hall velocity into account:
       ! most simple estimate, high k limit:
       ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
       kmax = dpi/min({dxlevel(^d)},bigdouble)*half
       csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
            twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
    end if
 
  end subroutine twofl_get_csound_c_idim
 
  !> Calculate fast magnetosonic wave speed when cbounds_species=false
  subroutine twofl_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
    use mod_global_parameters
 
    integer, intent(in)          :: ixI^L, ixO^L, idim
    double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
    double precision, intent(out):: csound(ixI^S)
    double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
    double precision :: inv_rho(ixO^S)
    double precision :: rhoc(ixI^S)
#if (defined(A_TOT) && A_TOT == 1)
    double precision :: rhon(ixI^S)
#endif
    call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
#if  (defined(A_TOT) && A_TOT == 1)
    call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
    inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s)) 
#else
    inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
#endif
 
    call twofl_get_csound2_primitive(w,x,ixi^l,ixo^l,csound)
 
    ! store |B|^2 in v
    b2(ixo^s)        = twofl_mag_en_all(w,ixi^l,ixo^l)
    cfast2(ixo^s)   = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
    avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
         * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
         * inv_rho(ixo^s)
 
    where(avmincs2(ixo^s)<zero)
       avmincs2(ixo^s)=zero
    end where
 
    avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
 
    if (.not. twofl_hall) then
       csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
    else
       ! take the Hall velocity into account:
       ! most simple estimate, high k limit:
       ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
       kmax = dpi/min({dxlevel(^d)},bigdouble)*half
       csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
            twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
    end if
 
    contains
    !TODO copy it inside
    subroutine twofl_get_csound2_primitive(w,x,ixI^L,ixO^L,csound2)
      use mod_global_parameters
      integer, intent(in)             :: ixI^L, ixO^L
      double precision, intent(in)    :: w(ixI^S,nw)
      double precision, intent(in)    :: x(ixI^S,1:ndim)
      double precision, intent(out)   :: csound2(ixI^S)
      double precision  :: pth_c(ixI^S)
      double precision  :: pth_n(ixI^S)
  
      if(phys_energy) then
        call twofl_get_pthermal_c_primitive(w,x,ixi^l,ixo^l,pth_c)
        call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,pth_n)
        call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
      else
        call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
      endif
    end subroutine twofl_get_csound2_primitive
 
  end subroutine twofl_get_csound_prim
 
  subroutine twofl_get_csound2(w,x,ixI^L,ixO^L,csound2)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: w(ixI^S,nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    double precision, intent(out)   :: csound2(ixI^S)
    double precision  :: pth_c(ixI^S)
    double precision  :: pth_n(ixI^S)
 
    if(phys_energy) then
      call twofl_get_pthermal_c(w,x,ixi^l,ixo^l,pth_c)
      call twofl_get_pthermal_n(w,x,ixi^l,ixo^l,pth_n)
      call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
    else
      call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
    endif
  end subroutine twofl_get_csound2
 
  subroutine twofl_get_csound2_adiab(w,x,ixI^L,ixO^L,csound2)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: w(ixI^S,nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    double precision, intent(out)   :: csound2(ixI^S)
    double precision  :: rhoc(ixI^S)
    double precision  :: rhon(ixI^S)
 
    call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
    call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
    csound2(ixo^s)=twofl_gamma*twofl_adiab*&
                  max((rhoc(ixo^s)**twofl_gamma + rhon(ixo^s)**twofl_gamma)/(rhoc(ixo^s)+ rhon(ixo^s)),&
                  rhon(ixo^s)**gamma_1,rhoc(ixo^s)**gamma_1) 
  end subroutine twofl_get_csound2_adiab
 
  subroutine twofl_get_csound(w,x,ixI^L,ixO^L,idim,csound)
    use mod_global_parameters
 
    integer, intent(in)          :: ixI^L, ixO^L, idim
    double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
    double precision, intent(out):: csound(ixI^S)
    double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
    double precision :: inv_rho(ixO^S)
    double precision :: rhoc(ixI^S)
#if (defined(A_TOT) && A_TOT == 1)
    double precision :: rhon(ixI^S)
#endif
    call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
#if (defined(A_TOT) && A_TOT == 1)
    call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
    inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s)) 
#else
    inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
#endif
 
    call twofl_get_csound2(w,x,ixi^l,ixo^l,csound)
 
    ! store |B|^2 in v
    b2(ixo^s) = twofl_mag_en_all(w,ixi^l,ixo^l)
 
    cfast2(ixo^s)   = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
    avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
         * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
         * inv_rho(ixo^s)
 
    where(avmincs2(ixo^s)<zero)
       avmincs2(ixo^s)=zero
    end where
 
    avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
 
    if (.not. twofl_hall) then
       csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
    else
       ! take the Hall velocity into account:
       ! most simple estimate, high k limit:
       ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
       kmax = dpi/min({dxlevel(^d)},bigdouble)*half
       csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
            twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
    end if
 
  end subroutine twofl_get_csound
 
  subroutine twofl_get_csound2_from_pthermal(w,x,ixI^L,ixO^L,pth_c,pth_n,csound2)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: w(ixI^S,nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    double precision, intent(in)    :: pth_c(ixI^S)
    double precision, intent(in)    :: pth_n(ixI^S)
    double precision, intent(out)   :: csound2(ixI^S)
    double precision  :: csound1(ixI^S),rhon(ixI^S),rhoc(ixI^S)
 
    call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
    call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
#if !defined(C_TOT) || C_TOT == 0
    csound2(ixo^s)=twofl_gamma*max((pth_c(ixo^s) + pth_n(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s)),&
                      pth_n(ixo^s)/rhon(ixo^s), pth_c(ixo^s)/rhoc(ixo^s))
#else
    csound2(ixo^s)=twofl_gamma*(csound2(ixo^s) + csound1(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s))
 
#endif
  end subroutine twofl_get_csound2_from_pthermal
 
! end cbounds_species=false
 
  subroutine twofl_get_csound_n(w,x,ixI^L,ixO^L,csound)
    use mod_global_parameters
 
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
    double precision, intent(out):: csound(ixI^S)
    double precision :: pe_n1(ixI^S)
    call twofl_get_csound2_n_from_conserved(w,x,ixi^l,ixo^l,csound)
    csound(ixo^s) = sqrt(csound(ixo^s))
  end subroutine twofl_get_csound_n
 
  !> separate routines so that it is faster
  !> Calculate temperature=p/rho when in e_ the internal energy is stored
  subroutine twofl_get_temperature_from_eint_n(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
 
    res(ixo^s) = 1d0/rn * gamma_1 * w(ixo^s, e_n_) /w(ixo^s,rho_n_)
 
  end subroutine twofl_get_temperature_from_eint_n
 
  subroutine twofl_get_temperature_from_eint_n_with_equi(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
 
      res(ixo^s) = 1d0/rn * (gamma_1 * w(ixo^s, e_n_) + block%equi_vars(ixo^s,equi_pe_n0_,b0i)) /&
                (w(ixo^s,rho_n_) +block%equi_vars(ixo^s,equi_rho_n0_,b0i))
  end subroutine twofl_get_temperature_from_eint_n_with_equi
 
!  subroutine twofl_get_temperature_n_pert_from_tot(Te, ixI^L, ixO^L, res)
!    use mod_global_parameters
!    integer, intent(in)          :: ixI^L, ixO^L
!    double precision, intent(in) :: Te(ixI^S)
!    double precision, intent(out):: res(ixI^S)
!      res(ixO^S) = Te(ixO^S) -1d0/Rn * &
!                block%equi_vars(ixO^S,equi_pe_n0_,0)/block%equi_vars(ixO^S,equi_rho_n0_,0)
!  end subroutine twofl_get_temperature_n_pert_from_tot
 
  subroutine twofl_get_temperature_n_equi(w,x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
      res(ixo^s) = 1d0/rn * &
                block%equi_vars(ixo^s,equi_pe_n0_,b0i)/block%equi_vars(ixo^s,equi_rho_n0_,b0i)
  end subroutine twofl_get_temperature_n_equi
 
  subroutine twofl_get_rho_n_equi(w, x,ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
      res(ixo^s) = block%equi_vars(ixo^s,equi_rho_n0_,b0i)
  end subroutine twofl_get_rho_n_equi
 
  subroutine twofl_get_pe_n_equi(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
      res(ixo^s) = block%equi_vars(ixo^s,equi_pe_n0_,b0i)
  end subroutine twofl_get_pe_n_equi
 
  !> Calculate temperature=p/rho when in e_ the total energy is stored
  !> this does not check the values of twofl_energy and twofl_internal_e, 
  !>  twofl_energy = .true. and twofl_internal_e = .false.
  !> also check small_values is avoided
  subroutine twofl_get_temperature_from_etot_n(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
    res(ixo^s)=1d0/rn * (gamma_1*(w(ixo^s,e_n_)&
           - twofl_kin_en_n(w,ixi^l,ixo^l)))/w(ixo^s,rho_n_)
  end subroutine twofl_get_temperature_from_etot_n
 
  subroutine twofl_get_temperature_from_etot_n_with_equi(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
    res(ixo^s)=1d0/rn * (gamma_1*(w(ixo^s,e_n_)&
           - twofl_kin_en_n(w,ixi^l,ixo^l)) +  block%equi_vars(ixo^s,equi_pe_n0_,b0i))&
            /(w(ixo^s,rho_n_) +block%equi_vars(ixo^s,equi_rho_n0_,b0i))
            
  end subroutine twofl_get_temperature_from_etot_n_with_equi
 
  !> separate routines so that it is faster
  !> Calculate temperature=p/rho when in e_ the internal energy is stored
  subroutine twofl_get_temperature_from_eint_c(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
 
    res(ixo^s) = 1d0/rc * gamma_1 * w(ixo^s, e_c_) /w(ixo^s,rho_c_)
 
  end subroutine twofl_get_temperature_from_eint_c
 
  subroutine twofl_get_temperature_from_eint_c_with_equi(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
      res(ixo^s) = 1d0/rc * (gamma_1 * w(ixo^s, e_c_) + block%equi_vars(ixo^s,equi_pe_c0_,b0i)) /&
                (w(ixo^s,rho_c_) +block%equi_vars(ixo^s,equi_rho_c0_,b0i))
  end subroutine twofl_get_temperature_from_eint_c_with_equi
 
!  subroutine twofl_get_temperature_c_pert_from_tot(Te, ixI^L, ixO^L, res)
!    use mod_global_parameters
!    integer, intent(in)          :: ixI^L, ixO^L
!    double precision, intent(in) :: Te(ixI^S)
!    double precision, intent(out):: res(ixI^S)
!      res(ixO^S) = Te(ixO^S) -1d0/Rc * &
!                block%equi_vars(ixO^S,equi_pe_c0_,0)/block%equi_vars(ixO^S,equi_rho_c0_,0)
!  end subroutine twofl_get_temperature_c_pert_from_tot
 
  subroutine twofl_get_temperature_c_equi(w,x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
      res(ixo^s) = 1d0/rc * &
                block%equi_vars(ixo^s,equi_pe_c0_,b0i)/block%equi_vars(ixo^s,equi_rho_c0_,b0i)
  end subroutine twofl_get_temperature_c_equi
 
  subroutine twofl_get_rho_c_equi(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
      res(ixo^s) = block%equi_vars(ixo^s,equi_rho_c0_,b0i)
  end subroutine twofl_get_rho_c_equi
 
  subroutine twofl_get_pe_c_equi(w,x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
      res(ixo^s) = block%equi_vars(ixo^s,equi_pe_c0_,b0i)
  end subroutine twofl_get_pe_c_equi
 
  !> Calculate temperature=p/rho when in e_ the total energy is stored
  !> this does not check the values of twofl_energy and twofl_internal_e, 
  !>  twofl_energy = .true. and twofl_internal_e = .false.
  !> also check small_values is avoided
  subroutine twofl_get_temperature_from_etot_c(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
    res(ixo^s)=1d0/rc * (gamma_1*(w(ixo^s,e_c_)&
           - twofl_kin_en_c(w,ixi^l,ixo^l)&
           - twofl_mag_en(w,ixi^l,ixo^l)))/w(ixo^s,rho_c_)
  end subroutine twofl_get_temperature_from_etot_c
  subroutine twofl_get_temperature_from_eki_c(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
    res(ixo^s)=1d0/rc * (gamma_1*(w(ixo^s,e_c_)&
           - twofl_kin_en_c(w,ixi^l,ixo^l)))/w(ixo^s,rho_c_)
  end subroutine twofl_get_temperature_from_eki_c
 
  subroutine twofl_get_temperature_from_etot_c_with_equi(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
    res(ixo^s)=1d0/rc * (gamma_1*(w(ixo^s,e_c_)&
           - twofl_kin_en_c(w,ixi^l,ixo^l)&
           - twofl_mag_en(w,ixi^l,ixo^l)) +  block%equi_vars(ixo^s,equi_pe_c0_,b0i))&
            /(w(ixo^s,rho_c_) +block%equi_vars(ixo^s,equi_rho_c0_,b0i))
            
  end subroutine twofl_get_temperature_from_etot_c_with_equi
 
  subroutine twofl_get_temperature_from_eki_c_with_equi(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
    res(ixo^s)=1d0/rc * (gamma_1*(w(ixo^s,e_c_)&
           - twofl_kin_en_c(w,ixi^l,ixo^l)) +  block%equi_vars(ixo^s,equi_pe_c0_,b0i))&
            /(w(ixo^s,rho_c_) +block%equi_vars(ixo^s,equi_rho_c0_,b0i))
            
  end subroutine twofl_get_temperature_from_eki_c_with_equi
 
  subroutine twofl_get_csound2_adiab_n(w,x,ixI^L,ixO^L,csound2)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: w(ixI^S,nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    double precision, intent(out)   :: csound2(ixI^S)
    double precision  :: rhon(ixI^S)
 
    call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
    csound2(ixo^s)=twofl_gamma*twofl_adiab*rhon(ixo^s)**gamma_1
 
  end subroutine twofl_get_csound2_adiab_n
 
  subroutine twofl_get_csound2_n_from_conserved(w,x,ixI^L,ixO^L,csound2)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: w(ixI^S,nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    double precision, intent(out)   :: csound2(ixI^S)
    double precision                :: rhon(ixI^S)
 
    if(phys_energy) then
      call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
      call twofl_get_pthermal_n(w,x,ixi^l,ixo^l,csound2)
      csound2(ixo^s)=twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
    else
      call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
    endif
  end subroutine twofl_get_csound2_n_from_conserved
 
  !! TO DELETE
  subroutine twofl_get_csound2_n_from_primitive(w,x,ixI^L,ixO^L,csound2)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: w(ixI^S,nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    double precision, intent(out)   :: csound2(ixI^S)
    double precision                :: rhon(ixI^S)
 
    if(phys_energy) then
      call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
      call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,csound2)
      csound2(ixo^s)=twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
    else
      call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
    endif
  end subroutine twofl_get_csound2_n_from_primitive
 
  subroutine twofl_get_csound2_adiab_c(w,x,ixI^L,ixO^L,csound2)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: w(ixI^S,nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    double precision, intent(out)   :: csound2(ixI^S)
    double precision  :: rhoc(ixI^S)
 
    call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
    csound2(ixo^s)=twofl_gamma*twofl_adiab* rhoc(ixo^s)**gamma_1
 
  end subroutine twofl_get_csound2_adiab_c
 
  subroutine twofl_get_csound2_c_from_conserved(w,x,ixI^L,ixO^L,csound2)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: w(ixi^s,nw)
    double precision, intent(in)    :: x(ixi^s,1:ndim)
    double precision, intent(out)   :: csound2(ixi^s)
    double precision                :: rhoc(ixi^s)
 
    if(phys_energy) then
      call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
      call twofl_get_pthermal_c(w,x,ixi^l,ixo^l,csound2)
      csound2(ixo^s)=twofl_gamma*csound2(ixo^s)/rhoc(ixo^s)
    else
      call twofl_get_csound2_adiab_c(w,x,ixi^l,ixo^l,csound2)
    endif
  end subroutine twofl_get_csound2_c_from_conserved
 
  !> Calculate fluxes within ixO^L.
  subroutine twofl_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in)          :: ixI^L, ixO^L, idim
    ! conservative w
    double precision, intent(in) :: wC(ixI^S,nw)
    ! primitive w
    double precision, intent(in) :: w(ixI^S,nw)
    double precision, intent(in) :: x(ixI^S,1:ndim)
    double precision,intent(out) :: f(ixI^S,nwflux)
 
    double precision             :: pgas(ixO^S), ptotal(ixO^S),tmp(ixI^S)
    double precision, allocatable:: vHall(:^D&,:)
    integer                      :: idirmin, iw, idir, jdir, kdir
 
    ! value at the interfaces, idim =  block%iw0 --> b0i 
    ! reuse tmp, used afterwards
    ! value at the interface so we can't put momentum
    call get_rhoc_tot(w,x,ixi^l,ixo^l,tmp)
    ! Get flux of density
    f(ixo^s,rho_c_)=w(ixo^s,mom_c(idim))*tmp(ixo^s)
    ! pgas is time dependent only
    if(phys_energy) then
      pgas(ixo^s)=w(ixo^s,e_c_)
    else
      pgas(ixo^s)=twofl_adiab*tmp(ixo^s)**twofl_gamma
      if(has_equi_pe_c0) then
        pgas(ixo^s)=pgas(ixo^s)-block%equi_vars(ixo^s,equi_pe_c0_,b0i)
      end if
    end if
 
    if (twofl_hall) then
      allocate(vhall(ixi^s,1:ndir))
      call twofl_getv_hall(w,x,ixi^l,ixo^l,vhall)
    end if
 
    if(b0field) tmp(ixo^s)=sum(block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=ndim+1)
 
    ptotal(ixo^s) = pgas(ixo^s) + 0.5d0*sum(w(ixo^s, mag(:))**2, dim=ndim+1)
 
    ! Get flux of momentum
    ! f_i[m_k]=v_i*m_k-b_k*b_i [+ptotal if i==k]
    do idir=1,ndir
      if(idim==idir) then
        f(ixo^s,mom_c(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
        if(b0field) f(ixo^s,mom_c(idir))=f(ixo^s,mom_c(idir))+tmp(ixo^s)
      else
        f(ixo^s,mom_c(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
      end if
      if (b0field) then
        f(ixo^s,mom_c(idir))=f(ixo^s,mom_c(idir))&
             -w(ixo^s,mag(idir))*block%B0(ixo^s,idim,idim)&
             -w(ixo^s,mag(idim))*block%B0(ixo^s,idir,idim)
      end if
      f(ixo^s,mom_c(idir))=f(ixo^s,mom_c(idir))+w(ixo^s,mom_c(idim))*wc(ixo^s,mom_c(idir))
    end do
 
    ! Get flux of energy
    ! f_i[e]=v_i*e+v_i*ptotal-b_i*(b_k*v_k)
    if(phys_energy) then
      if (phys_internal_e) then
         f(ixo^s,e_c_)=w(ixo^s,mom_c(idim))*wc(ixo^s,e_c_)
      else if(twofl_eq_energy == eq_energy_ki) then
 
        f(ixo^s,e_c_)=w(ixo^s,mom_c(idim))*(wc(ixo^s,e_c_)+pgas(ixo^s))
      else
        f(ixo^s,e_c_)=w(ixo^s,mom_c(idim))*(wc(ixo^s,e_c_)+ptotal(ixo^s))&
           -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,mom_c(:)),dim=ndim+1)
 
        if (b0field) then
           f(ixo^s,e_c_) = f(ixo^s,e_c_) &
              + w(ixo^s,mom_c(idim)) * tmp(ixo^s) &
              - sum(w(ixo^s,mom_c(:))*w(ixo^s,mag(:)),dim=ndim+1) * block%B0(ixo^s,idim,idim)
        end if
 
        if (twofl_hall) then
        ! f_i[e]= f_i[e] + vHall_i*(b_k*b_k) - b_i*(vHall_k*b_k)
           if (twofl_etah>zero) then
              f(ixo^s,e_c_) = f(ixo^s,e_c_) + vhall(ixo^s,idim) * &
                 sum(w(ixo^s, mag(:))**2,dim=ndim+1) &
                 - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1)
              if (b0field) then
                 f(ixo^s,e_c_) = f(ixo^s,e_c_) &
                    + vhall(ixo^s,idim) * tmp(ixo^s) &
                    - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1) * block%B0(ixo^s,idim,idim)
              end if
           end if
        end if
      end if !total_energy
      ! add flux of equilibrium internal energy corresponding to pe_c0
      if(has_equi_pe_c0) then
#if !defined(E_RM_W0) || E_RM_W0 == 1
        f(ixo^s,e_c_)=  f(ixo^s,e_c_) &
          + w(ixo^s,mom_c(idim)) * block%equi_vars(ixo^s,equi_pe_c0_,idim) * inv_gamma_1
#else
        if(phys_internal_e) then
          f(ixo^s,e_c_)=  f(ixo^s,e_c_) &
          + w(ixo^s,mom_c(idim)) * block%equi_vars(ixo^s,equi_pe_c0_,idim) * inv_gamma_1
        else
          f(ixo^s,e_c_)=  f(ixo^s,e_c_) &
          + w(ixo^s,mom_c(idim)) * block%equi_vars(ixo^s,equi_pe_c0_,idim) * twofl_gamma * inv_gamma_1
        end if
#endif
      end if
    end if !phys_energy
 
    ! compute flux of magnetic field
    ! f_i[b_k]=v_i*b_k-v_k*b_i
    do idir=1,ndir
      if (idim==idir) then
        ! f_i[b_i] should be exactly 0, so we do not use the transport flux
        if (twofl_glm) then
           f(ixo^s,mag(idir))=w(ixo^s,psi_)
        else
           f(ixo^s,mag(idir))=zero
        end if
      else
        f(ixo^s,mag(idir))=w(ixo^s,mom_c(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,mom_c(idir))
 
        if (b0field) then
          f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
                +w(ixo^s,mom_c(idim))*block%B0(ixo^s,idir,idim)&
                -w(ixo^s,mom_c(idir))*block%B0(ixo^s,idim,idim)
        end if
 
        if (twofl_hall) then
          ! f_i[b_k] = f_i[b_k] + vHall_i*b_k - vHall_k*b_i
          if (twofl_etah>zero) then
            if (b0field) then
              f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
                   - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+block%B0(ixo^s,idim,idim)) &
                   + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+block%B0(ixo^s,idir,idim))
            else
              f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
                   - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
                   + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
            end if
          end if
        end if
 
      end if
    end do
 
    if (twofl_glm) then
      !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
      f(ixo^s,psi_)  = cmax_global**2*w(ixo^s,mag(idim))
    end if
 
    if (twofl_hall) then
      deallocate(vhall)
    end if
 
    !!neutrals
    call get_rhon_tot(w,x,ixi^l,ixo^l,tmp)
    f(ixo^s,rho_n_)=w(ixo^s,mom_n(idim))*tmp(ixo^s)
    if(phys_energy) then
      pgas(ixo^s) = w(ixo^s, e_n_)
    else
      pgas(ixo^s)=twofl_adiab*tmp(ixo^s)**twofl_gamma
      if(has_equi_pe_n0) then
        pgas(ixo^s)=pgas(ixo^s)-block%equi_vars(ixo^s,equi_pe_n0_,b0i)
      end if
    end if
    ! Momentum flux is v_i*m_i, +p in direction idim
    do idir = 1, ndir
        !if(idim==idir) then
        !  f(ixO^S,mom_c(idir)) = pgas(ixO^S)
        !else
        !  f(ixO^S,mom_c(idir)) = 0.0d0
        !end if
        !f(ixO^S,mom_c(idir))=f(ixO^S,mom_c(idir))+w(ixO^S,mom_c(idim))*wC(ixO^S,mom_c(idir))
       f(ixo^s, mom_n(idir)) = w(ixo^s,mom_n(idim)) * wc(ixo^s, mom_n(idir))
    end do
 
    f(ixo^s, mom_n(idim)) = f(ixo^s, mom_n(idim)) + pgas(ixo^s)
 
    if(phys_energy) then
      !reuse pgas for storing a in the term: div (u_n * a) and make multiplication at the end
      pgas(ixo^s) = wc(ixo^s,e_n_)
      if(.not. phys_internal_e) then
        ! add pressure perturbation
        pgas(ixo^s) = pgas(ixo^s) + w(ixo^s,e_n_)
      end if
      ! add flux of equilibrium internal energy corresponding to pe_n0
      if(has_equi_pe_n0) then
#if !defined(E_RM_W0) || E_RM_W0 == 1
        pgas(ixo^s) = pgas(ixo^s) + block%equi_vars(ixo^s,equi_pe_n0_,idim) * inv_gamma_1
#else
        pgas(ixo^s) = pgas(ixo^s) + block%equi_vars(ixo^s,equi_pe_n0_,idim) * twofl_gamma * inv_gamma_1
#endif
      end if
      ! add u_n * a in the flux
      f(ixo^s, e_n_) = w(ixo^s,mom_n(idim)) * pgas(ixo^s)
 
      ! Viscosity fluxes - viscInDiv
      !if (hd_viscosity) then
      !  call visc_get_flux_prim(w, x, ixI^L, ixO^L, idim, f, phys_energy)
      !endif
    end if
 
  end subroutine twofl_get_flux
 
  !> w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
  subroutine twofl_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
    use mod_global_parameters
    use mod_radiative_cooling, only: radiative_cooling_add_source
    use mod_viscosity, only: viscosity_add_source
    !use mod_gravity, only: gravity_add_source
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt,dtfactor
    double precision, intent(in)    :: wCT(ixI^S,1:nw),wCTprim(ixI^S,1:nw),x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    logical, intent(in)             :: qsourcesplit
    logical, intent(inout)            :: active
 
    if (.not. qsourcesplit) then
      ! Source for solving internal energy
      if(phys_internal_e) then
        active = .true.
        call internal_energy_add_source_n(qdt,ixi^l,ixo^l,wct,w,x)
        call internal_energy_add_source_c(qdt,ixi^l,ixo^l,wct,w,x,e_c_)
      else 
#if !defined(E_RM_W0) || E_RM_W0==1
        ! add -p0 div v source terms when equi are present
        if(has_equi_pe_n0) then
          active = .true.
          call add_pe_n0_divv(qdt,ixi^l,ixo^l,wct,w,x)
        endif
        if(has_equi_pe_c0) then
          active = .true.
          call add_pe_c0_divv(qdt,ixi^l,ixo^l,wct,w,x)
        endif
#endif
        if(twofl_eq_energy == eq_energy_ki) then
          active = .true.
          call add_source_lorentz_work(qdt,ixi^l,ixo^l,w,wct,x)
        endif
      endif
 
      ! Source for B0 splitting
      if (b0field) then
        active = .true.
        call add_source_b0split(qdt,ixi^l,ixo^l,wct,w,x)
      end if
 
      ! Sources for resistivity in eqs. for e, B1, B2 and B3
      if (abs(twofl_eta)>smalldouble)then
        active = .true.
        call add_source_res2(qdt,ixi^l,ixo^l,wct,w,x)
      end if
 
      if (twofl_eta_hyper>0.d0)then
        active = .true.
        call add_source_hyperres(qdt,ixi^l,ixo^l,wct,w,x)
      end if
      !it is not added in a split manner
      if(.not. use_imex_scheme .and. has_collisions()) then
        active = .true.
        call  twofl_explicit_coll_terms_update(qdt,ixi^l,ixo^l,w,wct,x)
      endif
    
      if(twofl_hyperdiffusivity) then
        active = .true.
        call  add_source_hyperdiffusive(qdt,ixi^l,ixo^l,w,wct,x)
      endif  
 
    end if
 
    {^nooned
    if(source_split_divb .eqv. qsourcesplit) then
      ! Sources related to div B
      select case (type_divb)
      case (divb_none)
        ! Do nothing
      case (divb_glm)
        active = .true.
        call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
      case (divb_powel)
        active = .true.
        call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
      case (divb_janhunen)
        active = .true.
        call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
      case (divb_linde)
        active = .true.
        call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
      case (divb_lindejanhunen)
        active = .true.
        call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
        call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
      case (divb_lindepowel)
        active = .true.
        call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
        call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
      case (divb_lindeglm)
        active = .true.
        call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
        call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
      case (divb_ct)
        continue ! Do nothing
      case (divb_multigrid)
        continue ! Do nothing
      case default
        call mpistop('Unknown divB fix')
      end select
    end if
    }
 
    if(twofl_radiative_cooling_c) then
      call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,wctprim,&
           w,x,qsourcesplit,active,rc_fl_c)
    end if
    if(twofl_radiative_cooling_n) then
      call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,wctprim,&
           w,x,qsourcesplit,active,rc_fl_n)
    end if
!
!    if(twofl_viscosity) then
!      call viscosity_add_source(qdt,ixI^L,ixO^L,wCT,&
!           w,x,phys_energy,qsourcesplit,active)
!    end if
!
    if(twofl_gravity) then
      call gravity_add_source(qdt,ixi^l,ixo^l,wct,&
           w,x,twofl_eq_energy .eq. eq_energy_ki .or. phys_total_energy,qsourcesplit,active)
    end if
 
  end subroutine twofl_add_source
 
  subroutine add_pe_n0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt
    double precision, intent(in)    :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision                :: v(ixI^S,1:ndir)
 
    call twofl_get_v_n(wct,x,ixi^l,ixi^l,v)
    call add_geom_pdivv(qdt,ixi^l,ixo^l,v,-block%equi_vars(ixi^s,equi_pe_n0_,0),w,x,e_n_)
 
  end subroutine add_pe_n0_divv
 
  subroutine add_pe_c0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt
    double precision, intent(in)    :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision                :: v(ixI^S,1:ndir)
 
    call twofl_get_v_c(wct,x,ixi^l,ixi^l,v)
    call add_geom_pdivv(qdt,ixi^l,ixo^l,v,-block%equi_vars(ixi^s,equi_pe_c0_,0),w,x,e_c_)
 
  end subroutine add_pe_c0_divv
 
  subroutine add_geom_pdivv(qdt,ixI^L,ixO^L,v,p,w,x,ind)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in)             :: ixI^L, ixO^L,ind
    double precision, intent(in)    :: qdt
    double precision, intent(in)    :: p(ixI^S), v(ixI^S,1:ndir), x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision                :: divv(ixI^S)
 
    if(slab_uniform) then
      if(nghostcells .gt. 2) then
        call divvector(v,ixi^l,ixo^l,divv,sixthorder=.true.)
      else
        call divvector(v,ixi^l,ixo^l,divv,fourthorder=.true.)
      end if
    else
     call divvector(v,ixi^l,ixo^l,divv)
    end if
    w(ixo^s,ind)=w(ixo^s,ind)+qdt*p(ixo^s)*divv(ixo^s)
  end subroutine add_geom_pdivv
 
  !> Compute the Lorentz force (JxB)
  subroutine get_lorentz(ixI^L,ixO^L,w,JxB)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: w(ixI^S,1:nw)
    double precision, intent(inout) :: JxB(ixI^S,3)
    double precision                :: a(ixI^S,3), b(ixI^S,3), tmp(ixI^S,3)
    integer                         :: idir, idirmin
    ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
    double precision :: current(ixI^S,7-2*ndir:3)
 
    b=0.0d0
    do idir = 1, ndir
      b(ixo^s, idir) = twofl_mag_i_all(w, ixi^l, ixo^l,idir)
    end do
 
    ! store J current in a
    call get_current(w,ixi^l,ixo^l,idirmin,current)
 
    a=0.0d0
    do idir=7-2*ndir,3
      a(ixo^s,idir)=current(ixo^s,idir)
    end do
 
    call cross_product(ixi^l,ixo^l,a,b,jxb)
  end subroutine get_lorentz
 
  subroutine  add_source_lorentz_work(qdt,ixI^L,ixO^L,w,wCT,x)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt
    double precision, intent(in)    :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision                :: a(ixI^S,3), b(ixI^S,1:ndir)
    
    call get_lorentz(ixi^l, ixo^l,wct,a)
    call twofl_get_v_c(wct,x,ixi^l,ixo^l,b)
    w(ixo^s,e_c_)=w(ixo^s,e_c_)+qdt*sum(a(ixo^s,1:ndir)*b(ixo^s,1:ndir),dim=ndim+1)
 
  end subroutine  add_source_lorentz_work
 
  !> Calculate v_n vector
  subroutine twofl_get_v_n(w,x,ixI^L,ixO^L,v)
    use mod_global_parameters
 
    integer, intent(in)           :: ixI^L, ixO^L
    double precision, intent(in)  :: w(ixI^S,nw), x(ixI^S,1:ndim)
    double precision, intent(out) :: v(ixI^S,ndir)
    double precision              :: rhon(ixI^S)
    integer :: idir
 
    call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
 
    do idir=1,ndir
      v(ixo^s,idir) = w(ixo^s, mom_n(idir)) / rhon(ixo^s)
    end do
 
  end subroutine twofl_get_v_n
 
  subroutine get_rhon_tot(w,x,ixI^L,ixO^L,rhon)
    use mod_global_parameters
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s,1:nw), x(ixi^s,1:ndim)
    double precision, intent(out) :: rhon(ixi^s)
    if(has_equi_rho_n0) then
      rhon(ixo^s) = w(ixo^s,rho_n_) + block%equi_vars(ixo^s,equi_rho_n0_,b0i)
    else  
      rhon(ixo^s) = w(ixo^s,rho_n_) 
    endif
 
  end subroutine get_rhon_tot
 
  subroutine twofl_get_pthermal_n(w,x,ixI^L,ixO^L,pth)
    use mod_global_parameters
    use mod_small_values, only: trace_small_values
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s,1:nw)
    double precision, intent(in)  :: x(ixi^s,1:ndim)
    double precision, intent(out) :: pth(ixi^s)
 
    integer :: ix^d, iw
 
    if(phys_energy) then
      if(phys_internal_e) then
        pth(ixo^s)=gamma_1*w(ixo^s,e_n_)
      else
        pth(ixo^s)=gamma_1*(w(ixo^s,e_n_)&
           - twofl_kin_en_n(w,ixi^l,ixo^l))
      end if
      if(has_equi_pe_n0) then
        pth(ixo^s) = pth(ixo^s) + block%equi_vars(ixo^s,equi_pe_n0_,b0i)
      endif
    else
      call get_rhon_tot(w,x,ixi^l,ixo^l,pth)
      pth(ixo^s)=twofl_adiab*pth(ixo^s)**twofl_gamma
    end if
 
    if (fix_small_values) then
      {do ix^db= ixo^lim^db\}
         if(pth(ix^d)<small_pressure) then
            pth(ix^d)=small_pressure
         end if
      {enddo^d&\}
    else if (check_small_values) then
      {do ix^db= ixo^lim^db\}
         if(pth(ix^d)<small_pressure) then
           write(*,*) "Error: small value of gas pressure",pth(ix^d),&
                " encountered when call twofl_get_pthermal_n"
           write(*,*) "Iteration: ", it, " Time: ", global_time
           write(*,*) "Location: ", x(ix^d,:)
           write(*,*) "Cell number: ", ix^d
           do iw=1,nw
             write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
           end do
           ! use erroneous arithmetic operation to crash the run
           if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
           write(*,*) "Saving status at the previous time step"
           crash=.true.
         end if
      {enddo^d&\}
    end if
 
  end subroutine twofl_get_pthermal_n
 
  subroutine twofl_get_pthermal_n_primitive(w,x,ixI^L,ixO^L,pth)
    use mod_global_parameters
    integer, intent(in)           :: ixI^L, ixO^L
    double precision, intent(in)  :: w(ixI^S,1:nw)
    double precision, intent(in)  :: x(ixI^S,1:ndim)
    double precision, intent(out) :: pth(ixI^S)
 
    if(phys_energy) then
      if(has_equi_pe_n0) then
        pth(ixo^s) = w(ixo^s,e_n_) + block%equi_vars(ixo^s,equi_pe_n0_,b0i)
      else
        pth(ixo^s) = w(ixo^s,e_n_) 
      endif
    else
      call get_rhon_tot(w,x,ixi^l,ixo^l,pth)
      pth(ixo^s)=twofl_adiab*pth(ixo^s)**twofl_gamma
    end if
  end subroutine twofl_get_pthermal_n_primitive
 
  !> Calculate v component
  subroutine twofl_get_v_n_idim(w,x,ixI^L,ixO^L,idim,v)
    use mod_global_parameters
 
    integer, intent(in)           :: ixi^l, ixo^l, idim
    double precision, intent(in)  :: w(ixi^s,nw), x(ixi^s,1:ndim)
    double precision, intent(out) :: v(ixi^s)
    double precision              :: rhon(ixi^s)
 
    call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
    v(ixo^s) = w(ixo^s, mom_n(idim)) / rhon(ixo^s)
 
  end subroutine twofl_get_v_n_idim
 
  subroutine internal_energy_add_source_n(qdt,ixI^L,ixO^L,wCT,w,x)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt
    double precision, intent(in)    :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision                :: pth(ixI^S),v(ixI^S,1:ndir),divv(ixI^S)
 
    call twofl_get_pthermal_n(wct,x,ixi^l,ixo^l,pth)
    call twofl_get_v_n(wct,x,ixi^l,ixi^l,v)
    call add_geom_pdivv(qdt,ixi^l,ixo^l,v,-pth,w,x,e_n_)
 
    if(fix_small_values .and. .not. has_equi_pe_n0) then
      call twofl_handle_small_ei_n(w,x,ixi^l,ixo^l,e_n_,'internal_energy_add_source')
    end if
  end subroutine internal_energy_add_source_n
 
  !> Calculate v_c vector
  subroutine twofl_get_v_c(w,x,ixI^L,ixO^L,v)
    use mod_global_parameters
 
    integer, intent(in)           :: ixI^L, ixO^L
    double precision, intent(in)  :: w(ixI^S,nw), x(ixI^S,1:ndim)
    double precision, intent(out) :: v(ixI^S,ndir)
    double precision              :: rhoc(ixI^S)
    integer :: idir
 
    call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
    do idir=1,ndir
      v(ixo^s,idir) = w(ixo^s, mom_c(idir)) / rhoc(ixo^s)
    end do
 
  end subroutine twofl_get_v_c
 
  subroutine get_rhoc_tot(w,x,ixI^L,ixO^L,rhoc)
    use mod_global_parameters
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s,1:nw), x(ixi^s,1:ndim)
    double precision, intent(out) :: rhoc(ixi^s)
    if(has_equi_rho_c0) then
      rhoc(ixo^s) = w(ixo^s,rho_c_) + block%equi_vars(ixo^s,equi_rho_c0_,b0i)
    else  
      rhoc(ixo^s) = w(ixo^s,rho_c_) 
    endif
 
  end subroutine get_rhoc_tot
 
  subroutine twofl_get_pthermal_c(w,x,ixI^L,ixO^L,pth)
    use mod_global_parameters
    use mod_small_values, only: trace_small_values
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s,1:nw)
    double precision, intent(in)  :: x(ixi^s,1:ndim)
    double precision, intent(out) :: pth(ixi^s)
    integer :: ix^d, iw
 
    if(phys_energy) then
      if(phys_internal_e) then
        pth(ixo^s)=gamma_1*w(ixo^s,e_c_)
      elseif(phys_total_energy) then
        pth(ixo^s)=gamma_1*(w(ixo^s,e_c_)&
           - twofl_kin_en_c(w,ixi^l,ixo^l)&
           - twofl_mag_en(w,ixi^l,ixo^l))
      else
        pth(ixo^s)=gamma_1*(w(ixo^s,e_c_)&
           - twofl_kin_en_c(w,ixi^l,ixo^l))
      end if
      if(has_equi_pe_c0) then
        pth(ixo^s) = pth(ixo^s) + block%equi_vars(ixo^s,equi_pe_c0_,b0i)
      endif
    else
      call get_rhoc_tot(w,x,ixi^l,ixo^l,pth)
      pth(ixo^s)=twofl_adiab*pth(ixo^s)**twofl_gamma
    end if
 
    if (fix_small_values) then
      {do ix^db= ixo^lim^db\}
         if(pth(ix^d)<small_pressure) then
            pth(ix^d)=small_pressure
         end if
      {enddo^d&\}
    else if (check_small_values) then
      {do ix^db= ixo^lim^db\}
         if(pth(ix^d)<small_pressure) then
           write(*,*) "Error: small value of gas pressure",pth(ix^d),&
                " encountered when call twofl_get_pe_c1"
           write(*,*) "Iteration: ", it, " Time: ", global_time
           write(*,*) "Location: ", x(ix^d,:)
           write(*,*) "Cell number: ", ix^d
           do iw=1,nw
             write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
           end do
           ! use erroneous arithmetic operation to crash the run
           if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
           write(*,*) "Saving status at the previous time step"
           crash=.true.
         end if
      {enddo^d&\}
    end if
 
  end subroutine twofl_get_pthermal_c
 
  subroutine twofl_get_pthermal_c_primitive(w,x,ixI^L,ixO^L,pth)
    use mod_global_parameters
    integer, intent(in)           :: ixI^L, ixO^L
    double precision, intent(in)  :: w(ixI^S,1:nw)
    double precision, intent(in)  :: x(ixI^S,1:ndim)
    double precision, intent(out) :: pth(ixI^S)
 
    if(phys_energy) then
      if(has_equi_pe_c0) then
        pth(ixo^s) = w(ixo^s,e_c_) + block%equi_vars(ixo^s,equi_pe_c0_,b0i)
      else
        pth(ixo^s) = w(ixo^s,e_c_) 
      endif
    else
      call get_rhoc_tot(w,x,ixi^l,ixo^l,pth)
      pth(ixo^s)=twofl_adiab*pth(ixo^s)**twofl_gamma
    end if
  end subroutine twofl_get_pthermal_c_primitive
 
  !> Calculate v_c component
  subroutine twofl_get_v_c_idim(w,x,ixI^L,ixO^L,idim,v)
    use mod_global_parameters
 
    integer, intent(in)           :: ixi^l, ixo^l, idim
    double precision, intent(in)  :: w(ixi^s,nw), x(ixi^s,1:ndim)
    double precision, intent(out) :: v(ixi^s)
    double precision              :: rhoc(ixi^s)
 
    call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
    v(ixo^s) = w(ixo^s, mom_c(idim)) / rhoc(ixo^s)
 
  end subroutine twofl_get_v_c_idim
 
  subroutine internal_energy_add_source_c(qdt,ixI^L,ixO^L,wCT,w,x,ie)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in)             :: ixI^L, ixO^L,ie
    double precision, intent(in)    :: qdt
    double precision, intent(in)    :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision                :: pth(ixI^S),v(ixI^S,1:ndir),divv(ixI^S)
 
    call twofl_get_pthermal_c(wct,x,ixi^l,ixo^l,pth)
    call twofl_get_v_c(wct,x,ixi^l,ixi^l,v)
    call add_geom_pdivv(qdt,ixi^l,ixo^l,v,-pth,w,x,ie)
    if(fix_small_values .and. .not. has_equi_pe_c0) then
      call twofl_handle_small_ei_c(w,x,ixi^l,ixo^l,ie,'internal_energy_add_source')
    end if
  end subroutine internal_energy_add_source_c
 
  !> handle small or negative internal energy
  subroutine twofl_handle_small_ei_c(w, x, ixI^L, ixO^L, ie, subname)
    use mod_global_parameters
    use mod_small_values
    integer, intent(in)             :: ixI^L,ixO^L, ie
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    character(len=*), intent(in)    :: subname
 
    integer :: idir
    logical :: flag(ixI^S,1:nw)
    double precision              :: rhoc(ixI^S)
    double precision              :: rhon(ixI^S)
 
    flag=.false.
    if(has_equi_pe_c0) then
      where(w(ixo^s,ie)+block%equi_vars(ixo^s,equi_pe_c0_,0)*inv_gamma_1<small_e)&
             flag(ixo^s,ie)=.true.
    else
      where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
    endif  
    if(any(flag(ixo^s,ie))) then
      select case (small_values_method)
      case ("replace")
        if(has_equi_pe_c0) then
          where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
                  block%equi_vars(ixo^s,equi_pe_c0_,0)*inv_gamma_1
        else
          where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
        endif
      case ("average")
        call small_values_average(ixi^l, ixo^l, w, x, flag, ie)
      case default
        ! small values error shows primitive variables
        ! to_primitive subroutine cannot be used as this error handling
        ! is also used in TC where e_to_ei is explicitly called
        w(ixo^s,e_n_)=w(ixo^s,e_n_)*gamma_1
        call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
        w(ixo^s,e_c_)=w(ixo^s,e_c_)*gamma_1
        call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
        do idir = 1, ndir
           w(ixo^s, mom_n(idir)) = w(ixo^s, mom_n(idir))/rhon(ixo^s)
           w(ixo^s, mom_c(idir)) = w(ixo^s, mom_c(idir))/rhoc(ixo^s)
        end do
        call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
      end select
    end if
 
  end subroutine twofl_handle_small_ei_c
 
  !> handle small or negative internal energy
  subroutine twofl_handle_small_ei_n(w, x, ixI^L, ixO^L, ie, subname)
    use mod_global_parameters
    use mod_small_values
    integer, intent(in)             :: ixI^L,ixO^L, ie
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    character(len=*), intent(in)    :: subname
 
    integer :: idir
    logical :: flag(ixI^S,1:nw)
    double precision              :: rhoc(ixI^S)
    double precision              :: rhon(ixI^S)
 
    flag=.false.
    if(has_equi_pe_n0) then
      where(w(ixo^s,ie)+block%equi_vars(ixo^s,equi_pe_n0_,0)*inv_gamma_1<small_e)&
                         flag(ixo^s,ie)=.true.
    else
      where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
    endif
    if(any(flag(ixo^s,ie))) then
      select case (small_values_method)
      case ("replace")
        if(has_equi_pe_n0) then
          where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
                block%equi_vars(ixo^s,equi_pe_n0_,0)*inv_gamma_1
        else
          where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
        endif
      case ("average")
        call small_values_average(ixi^l, ixo^l, w, x, flag, ie)
      case default
        ! small values error shows primitive variables
        w(ixo^s,e_n_)=w(ixo^s,e_n_)*gamma_1
        call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
        w(ixo^s,e_c_)=w(ixo^s,e_c_)*gamma_1
        call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
        do idir = 1, ndir
           w(ixo^s, mom_n(idir)) = w(ixo^s, mom_n(idir))/rhon(ixo^s)
           w(ixo^s, mom_c(idir)) = w(ixo^s, mom_c(idir))/rhoc(ixo^s)
        end do
        call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
      end select
    end if
 
  end subroutine twofl_handle_small_ei_n
 
  !> Source terms after split off time-independent magnetic field
  subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
    use mod_global_parameters
 
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
 
    double precision :: a(ixI^S,3), b(ixI^S,3), axb(ixI^S,3)
    integer :: idir
 
    a=0.d0
    b=0.d0
    ! for force-free field J0xB0 =0
    if(.not.b0field_forcefree) then
      ! store B0 magnetic field in b
      b(ixo^s,1:ndir)=block%B0(ixo^s,1:ndir,0)
 
      ! store J0 current in a
      do idir=7-2*ndir,3
        a(ixo^s,idir)=block%J0(ixo^s,idir)
      end do
      call cross_product(ixi^l,ixo^l,a,b,axb)
      axb(ixo^s,:)=axb(ixo^s,:)*qdt
      ! add J0xB0 source term in momentum equations
      w(ixo^s,mom_c(1:ndir))=w(ixo^s,mom_c(1:ndir))+axb(ixo^s,1:ndir)
    end if
 
    if(phys_total_energy) then
      a=0.d0
      ! for free-free field -(vxB0) dot J0 =0
      b(ixo^s,:)=wct(ixo^s,mag(:))
      ! store full magnetic field B0+B1 in b
      if(.not.b0field_forcefree) b(ixo^s,:)=b(ixo^s,:)+block%B0(ixo^s,:,0)
      ! store velocity in a
      do idir=1,ndir
        call twofl_get_v_c_idim(wct,x,ixi^l,ixo^l,idir,a(ixi^s,idir))
      end do
      call cross_product(ixi^l,ixo^l,a,b,axb)
      axb(ixo^s,:)=axb(ixo^s,:)*qdt
      ! add -(vxB) dot J0 source term in energy equation
      do idir=7-2*ndir,3
        w(ixo^s,e_c_)=w(ixo^s,e_c_)-axb(ixo^s,idir)*block%J0(ixo^s,idir)
      end do
    end if
 
    if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_B0')
 
  end subroutine add_source_b0split
 
  !> Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in
  !> each direction, non-conservative. If the fourthorder precompiler flag is
  !> set, uses fourth order central difference for the laplacian. Then the
  !> stencil is 5 (2 neighbours).
  subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
    use mod_global_parameters
    use mod_usr_methods
    use mod_geometry
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt
    double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
    integer :: lxO^L, kxO^L
 
    double precision :: tmp(ixI^S),tmp2(ixI^S)
 
    ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
    double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
    double precision :: gradeta(ixI^S,1:ndim), Bf(ixI^S,1:ndir)
 
    ! Calculating resistive sources involve one extra layer
    if (twofl_4th_order) then
      ixa^l=ixo^l^ladd2;
    else
      ixa^l=ixo^l^ladd1;
    end if
 
    if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
         call mpistop("Error in add_source_res1: Non-conforming input limits")
 
    ! Calculate current density and idirmin
    call get_current(wct,ixi^l,ixo^l,idirmin,current)
 
    if (twofl_eta>zero)then
       eta(ixa^s)=twofl_eta
       gradeta(ixo^s,1:ndim)=zero
    else
       call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
       ! assumes that eta is not function of current?
       do idim=1,ndim
          call gradient(eta,ixi^l,ixo^l,idim,tmp)
          gradeta(ixo^s,idim)=tmp(ixo^s)
       end do
    end if
 
    if(b0field) then
      bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+block%B0(ixi^s,1:ndir,0)
    else
      bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
    end if
 
    do idir=1,ndir
       ! Put B_idir into tmp2 and eta*Laplace B_idir into tmp
       if (twofl_4th_order) then
         tmp(ixo^s)=zero
         tmp2(ixi^s)=bf(ixi^s,idir)
         do idim=1,ndim
            lxo^l=ixo^l+2*kr(idim,^d);
            jxo^l=ixo^l+kr(idim,^d);
            hxo^l=ixo^l-kr(idim,^d);
            kxo^l=ixo^l-2*kr(idim,^d);
            tmp(ixo^s)=tmp(ixo^s)+&
                 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
                 /(12.0d0 * dxlevel(idim)**2)
         end do
       else
         tmp(ixo^s)=zero
         tmp2(ixi^s)=bf(ixi^s,idir)
         do idim=1,ndim
            jxo^l=ixo^l+kr(idim,^d);
            hxo^l=ixo^l-kr(idim,^d);
            tmp(ixo^s)=tmp(ixo^s)+&
                 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/dxlevel(idim)**2
         end do
       end if
 
       ! Multiply by eta
       tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
 
       ! Subtract grad(eta) x J = eps_ijk d_j eta J_k if eta is non-constant
       if (twofl_eta<zero)then
          do jdir=1,ndim; do kdir=idirmin,3
             if (lvc(idir,jdir,kdir)/=0)then
                if (lvc(idir,jdir,kdir)==1)then
                   tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
                else
                   tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
                end if
             end if
          end do; end do
       end if
 
       ! Add sources related to eta*laplB-grad(eta) x J to B and e
       w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
       if (phys_total_energy) then
          w(ixo^s,e_c_)=w(ixo^s,e_c_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
       end if
    end do ! idir
 
    if (phys_energy) then
       ! de/dt+=eta*J**2
      w(ixo^s,e_c_)=w(ixo^s,e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
    end if
 
    if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res1')
 
  end subroutine add_source_res1
 
  !> Add resistive source to w within ixO
  !> Uses 5 point stencil (2 neighbours) in each direction, conservative
  subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
    use mod_global_parameters
    use mod_usr_methods
    use mod_geometry
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt
    double precision, intent(in)    :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
 
    ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
    double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S),curlj(ixI^S,1:3)
    double precision :: tmpvec(ixI^S,1:3),tmp(ixO^S)
    integer :: ixA^L,idir,idirmin,idirmin1
 
    ixa^l=ixo^l^ladd2;
 
    if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
         call mpistop("Error in add_source_res2: Non-conforming input limits")
 
    ixa^l=ixo^l^ladd1;
    ! Calculate current density within ixL: J=curl B, thus J_i=eps_ijk*d_j B_k
    ! Determine exact value of idirmin while doing the loop.
    call get_current(wct,ixi^l,ixa^l,idirmin,current)
 
    if (twofl_eta>zero)then
       eta(ixa^s)=twofl_eta
    else
       call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
    end if
 
    ! dB/dt= -curl(J*eta), thus B_i=B_i-eps_ijk d_j Jeta_k
    tmpvec(ixa^s,1:ndir)=zero
    do idir=idirmin,3
       tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
    end do
    curlj=0.d0
    call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
    if(stagger_grid.and.ndim==2.and.ndir==3) then
      ! if 2.5D
      w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
    else
      w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
    end if
 
    if(phys_energy) then
      if(phys_total_energy) then
        ! de/dt= +div(B x Jeta) = eta J^2 - B dot curl(eta J)
        ! de1/dt= eta J^2 - B1 dot curl(eta J)
        w(ixo^s,e_c_)=w(ixo^s,e_c_)+qdt*(eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)-&
          sum(wct(ixo^s,mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1))
      else
        ! add eta*J**2 source term in the internal energy equation
        w(ixo^s,e_c_)=w(ixo^s,e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
      end if
 
    end if
 
    if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res2')
  end subroutine add_source_res2
 
  !> Add Hyper-resistive source to w within ixO
  !> Uses 9 point stencil (4 neighbours) in each direction.
  subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt
    double precision, intent(in)    :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    !.. local ..
    double precision                :: current(ixI^S,7-2*ndir:3)
    double precision                :: tmpvec(ixI^S,1:3),tmpvec2(ixI^S,1:3),tmp(ixI^S),ehyper(ixI^S,1:3)
    integer                         :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
 
    ixa^l=ixo^l^ladd3;
    if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
         call mpistop("Error in add_source_hyperres: Non-conforming input limits")
 
    call get_current(wct,ixi^l,ixa^l,idirmin,current)
    tmpvec(ixa^s,1:ndir)=zero
    do jdir=idirmin,3
       tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
    end do
 
    ixa^l=ixo^l^ladd2;
    call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
 
    ixa^l=ixo^l^ladd1;
    tmpvec(ixa^s,1:ndir)=zero
    call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
    ehyper(ixa^s,1:ndir) = - tmpvec(ixa^s,1:ndir)*twofl_eta_hyper
 
    ixa^l=ixo^l;
    tmpvec2(ixa^s,1:ndir)=zero
    call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
 
    do idir=1,ndir
      w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
    end do
 
    if (phys_energy) then
      ! de/dt= +div(B x Ehyper)
      ixa^l=ixo^l^ladd1;
      tmpvec2(ixa^s,1:ndir)=zero
      do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
        tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
        + lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
      end do; end do; end do
      tmp(ixo^s)=zero
      call divvector(tmpvec2,ixi^l,ixo^l,tmp)
      w(ixo^s,e_c_)=w(ixo^s,e_c_)+tmp(ixo^s)*qdt
    end if
 
    if (fix_small_values)  call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_hyperres')
 
  end subroutine add_source_hyperres
 
  subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
    ! Add divB related sources to w within ixO
    ! corresponding to Dedner JCP 2002, 175, 645 _equation 24_
    ! giving the EGLM-MHD scheme
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision:: divb(ixI^S)
    integer          :: idim,idir
    double precision :: gradPsi(ixI^S)
 
    ! We calculate now div B
    call get_divb(wct,ixi^l,ixo^l,divb, twofl_divb_4thorder)
 
    ! dPsi/dt =  - Ch^2/Cp^2 Psi
    if (twofl_glm_alpha < zero) then
      w(ixo^s,psi_) = abs(twofl_glm_alpha)*wct(ixo^s,psi_)
    else
      ! implicit update of Psi variable
      ! equation (27) in Mignone 2010 J. Com. Phys. 229, 2117
      if(slab_uniform) then
        w(ixo^s,psi_) = dexp(-qdt*cmax_global*twofl_glm_alpha/minval(dxlevel(:)))*w(ixo^s,psi_)
      else
        w(ixo^s,psi_) = dexp(-qdt*cmax_global*twofl_glm_alpha/minval(block%ds(ixo^s,:),dim=ndim+1))*w(ixo^s,psi_)
      end if
    end if
 
    ! gradient of Psi
    do idim=1,ndim
       select case(typegrad)
       case("central")
          call gradient(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
       case("limited")
          call gradients(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
       end select
       if (phys_total_energy) then
       ! e  = e  -qdt (b . grad(Psi))
         w(ixo^s,e_c_) = w(ixo^s,e_c_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
       end if
    end do
 
    ! m = m - qdt b div b
    do idir=1,ndir
      w(ixo^s,mom_c(idir))=w(ixo^s,mom_c(idir))-qdt*twofl_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
    end do
 
    if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_glm')
 
  end subroutine add_source_glm
 
  !> Add divB related sources to w within ixO corresponding to Powel
  subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
    use mod_global_parameters
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt,   wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision                :: divb(ixI^S),v(ixI^S,1:ndir)
    integer                         :: idir
 
    ! We calculate now div B
    call get_divb(wct,ixi^l,ixo^l,divb, twofl_divb_4thorder)
 
    ! calculate velocity
    call twofl_get_v_c(wct,x,ixi^l,ixo^l,v)
 
    if (phys_total_energy) then
      ! e = e - qdt (v . b) * div b
      w(ixo^s,e_c_)=w(ixo^s,e_c_)-&
           qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=ndim+1)*divb(ixo^s)
    end if
 
    ! b = b - qdt v * div b
    do idir=1,ndir
      w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
    end do
 
    ! m = m - qdt b div b
    do idir=1,ndir
      w(ixo^s,mom_c(idir))=w(ixo^s,mom_c(idir))-qdt*twofl_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
    end do
 
    if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_powel')
 
  end subroutine add_source_powel
 
  subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
    ! Add divB related sources to w within ixO
    ! corresponding to Janhunen, just the term in the induction equation.
    use mod_global_parameters
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt,   wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision                :: divb(ixI^S),vel(ixI^S)
    integer                         :: idir
 
    ! We calculate now div B
    call get_divb(wct,ixi^l,ixo^l,divb, twofl_divb_4thorder)
 
    ! b = b - qdt v * div b
    do idir=1,ndir
      call twofl_get_v_c_idim(wct,x,ixi^l,ixo^l,idir,vel)
      w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s)*divb(ixo^s)
    end do
 
    if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_janhunen')
 
  end subroutine add_source_janhunen
 
  subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
    ! Add Linde's divB related sources to wnew within ixO
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    integer :: idim, idir, ixp^L, i^D, iside
    double precision :: divb(ixI^S),graddivb(ixI^S)
    logical, dimension(-1:1^D&) :: leveljump
 
    ! Calculate div B
    ixp^l=ixo^l^ladd1;
    call get_divb(wct,ixi^l,ixp^l,divb, twofl_divb_4thorder)
 
    ! for AMR stability, retreat one cell layer from the boarders of level jump
    {do i^db=-1,1\}
      if(i^d==0|.and.) cycle
      if(neighbor_type(i^d,block%igrid)==2 .or. neighbor_type(i^d,block%igrid)==4) then
        leveljump(i^d)=.true.
      else
        leveljump(i^d)=.false.
      end if
    {end do\}
 
    ixp^l=ixo^l;
    do idim=1,ndim
      select case(idim)
       {case(^d)
          do iside=1,2
            i^dd=kr(^dd,^d)*(2*iside-3);
            if (leveljump(i^dd)) then
              if (iside==1) then
                ixpmin^d=ixomin^d-i^d
              else
                ixpmax^d=ixomax^d-i^d
              end if
            end if
          end do
       \}
      end select
    end do
 
    ! Add Linde's diffusive terms
    do idim=1,ndim
       ! Calculate grad_idim(divb)
       select case(typegrad)
       case("central")
         call gradient(divb,ixi^l,ixp^l,idim,graddivb)
       case("limited")
         call gradients(divb,ixi^l,ixp^l,idim,graddivb)
       end select
 
       ! Multiply by Linde's eta*dt = divbdiff*(c_max*dx)*dt = divbdiff*dx**2
       if (slab_uniform) then
          graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
       else
          graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
                          /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
       end if
 
       w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
 
       if (typedivbdiff=='all' .and. phys_total_energy) then
         ! e += B_idim*eta*grad_idim(divb)
         w(ixp^s,e_c_)=w(ixp^s,e_c_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
       end if
    end do
 
    if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_linde')
 
  end subroutine add_source_linde
 
 
  !> get dimensionless div B = |divB| * volume / area / |B|
  subroutine get_normalized_divb(w,ixI^L,ixO^L,divb)
 
    use mod_global_parameters
 
    integer, intent(in)                :: ixi^l, ixo^l
    double precision, intent(in)       :: w(ixi^s,1:nw)
    double precision                   :: divb(ixi^s), dsurface(ixi^s)
 
    double precision :: invb(ixo^s)
    integer :: ixa^l,idims
 
    call get_divb(w,ixi^l,ixo^l,divb)
    invb(ixo^s)=sqrt(twofl_mag_en_all(w,ixi^l,ixo^l))
    where(invb(ixo^s)/=0.d0)
      invb(ixo^s)=1.d0/invb(ixo^s)
    end where
    if(slab_uniform) then
      divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/dxlevel(:))
    else
      ixamin^d=ixomin^d-1;
      ixamax^d=ixomax^d-1;
      dsurface(ixo^s)= sum(block%surfaceC(ixo^s,:),dim=ndim+1)
      do idims=1,ndim
        ixa^l=ixo^l-kr(idims,^d);
        dsurface(ixo^s)=dsurface(ixo^s)+block%surfaceC(ixa^s,idims)
      end do
      divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
      block%dvolume(ixo^s)/dsurface(ixo^s)
    end if
 
  end subroutine get_normalized_divb
 
  !> Calculate idirmin and the idirmin:3 components of the common current array
  !> make sure that dxlevel(^D) is set correctly.
  subroutine get_current(w,ixI^L,ixO^L,idirmin,current)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in)  :: ixo^l, ixi^l
    double precision, intent(in) :: w(ixi^s,1:nw)
    integer, intent(out) :: idirmin
    integer :: idir, idirmin0
 
    ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
    double precision :: current(ixi^s,7-2*ndir:3),bvec(ixi^s,1:ndir)
 
    idirmin0 = 7-2*ndir
 
    bvec(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))
 
    call curlvector(bvec,ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
 
    if(b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
        block%J0(ixo^s,idirmin0:3)
 
  end subroutine get_current
 
  ! copied from gravity
  !> w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
  subroutine gravity_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
       energy,qsourcesplit,active)
    use mod_global_parameters
    use mod_usr_methods
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt, x(ixI^S,1:ndim)
    double precision, intent(in)    :: wCT(ixI^S,1:nw)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    logical, intent(in) :: energy,qsourcesplit
    logical, intent(inout) :: active
    double precision       :: vel(ixI^S)
    integer                         :: idim
 
    double precision :: gravity_field(ixI^S,ndim)
 
    if(qsourcesplit .eqv. grav_split) then
      active = .true.
 
      if (.not. associated(usr_gravity)) then
        write(*,*) "mod_usr.t: please point usr_gravity to a subroutine"
        write(*,*) "like the phys_gravity in mod_usr_methods.t"
        call mpistop("gravity_add_source: usr_gravity not defined")
      else
        call usr_gravity(ixi^l,ixo^l,wct,x,gravity_field)
      end if
  
      do idim = 1, ndim
        w(ixo^s,mom_n(idim)) = w(ixo^s,mom_n(idim)) &
              + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,rho_n_)
        w(ixo^s,mom_c(idim)) = w(ixo^s,mom_c(idim)) &
              + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,rho_c_)
        if(energy) then
#if !defined(E_RM_W0) || E_RM_W0 == 1
          call twofl_get_v_n_idim(wct,x,ixi^l,ixo^l,idim,vel)
          w(ixo^s,e_n_)=w(ixo^s,e_n_) &
              + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,rho_n_)
          call twofl_get_v_c_idim(wct,x,ixi^l,ixo^l,idim,vel)
          w(ixo^s,e_c_)=w(ixo^s,e_c_) &
              + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,rho_c_)
#else
          w(ixo^s,e_n_)=w(ixo^s,e_n_) &
              + qdt * gravity_field(ixo^s,idim) *  wct(ixo^s,mom_n(idim))
          w(ixo^s,e_c_)=w(ixo^s,e_c_) &
              + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,mom_c(idim))
#endif
 
 
        end if
      end do
    end if
 
  end subroutine gravity_add_source
 
  subroutine gravity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
    use mod_global_parameters
    use mod_usr_methods
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: dx^D, x(ixI^S,1:ndim), w(ixI^S,1:nw)
    double precision, intent(inout) :: dtnew
 
    double precision                :: dxinv(1:ndim), max_grav
    integer                         :: idim
 
    double precision :: gravity_field(ixI^S,ndim)
 
    ^d&dxinv(^d)=one/dx^d;
 
    if(.not. associated(usr_gravity)) then
      write(*,*) "mod_usr.t: please point usr_gravity to a subroutine"
      write(*,*) "like the phys_gravity in mod_usr_methods.t"
      call mpistop("gravity_get_dt: usr_gravity not defined")
    else
      call usr_gravity(ixi^l,ixo^l,w,x,gravity_field)
    end if
 
    do idim = 1, ndim
      max_grav = maxval(abs(gravity_field(ixo^s,idim)))
      max_grav = max(max_grav, epsilon(1.0d0))
      dtnew = min(dtnew, 1.0d0 / sqrt(max_grav * dxinv(idim)))
    end do
 
  end subroutine gravity_get_dt
 
  !> If resistivity is not zero, check diffusion time limit for dt
  subroutine twofl_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
    use mod_global_parameters
    use mod_usr_methods
    use mod_radiative_cooling, only: cooling_get_dt
    !use mod_viscosity, only: viscosity_get_dt
    !use mod_gravity, only: gravity_get_dt
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(inout) :: dtnew
    double precision, intent(in)    :: dx^D
    double precision, intent(in)    :: w(ixI^S,1:nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
 
    integer                       :: idirmin,idim
    double precision              :: dxarr(ndim)
    double precision              :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
 
    dtnew = bigdouble
 
    ^d&dxarr(^d)=dx^d;
    if (twofl_eta>zero)then
       dtnew=dtdiffpar*minval(dxarr(1:ndim))**2/twofl_eta
    else if (twofl_eta<zero)then
       call get_current(w,ixi^l,ixo^l,idirmin,current)
       call usr_special_resistivity(w,ixi^l,ixo^l,idirmin,x,current,eta)
       dtnew=bigdouble
       do idim=1,ndim
         if(slab_uniform) then
           dtnew=min(dtnew,&
                dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
         else
           dtnew=min(dtnew,&
                dtdiffpar/(smalldouble+maxval(eta(ixo^s)/block%ds(ixo^s,idim)**2)))
         end if
       end do
    end if
 
    if(twofl_eta_hyper>zero) then
      if(slab_uniform) then
        dtnew=min(dtdiffpar*minval(dxarr(1:ndim))**4/twofl_eta_hyper,dtnew)
      else
        dtnew=min(dtdiffpar*minval(block%ds(ixo^s,1:ndim))**4/twofl_eta_hyper,dtnew)
      end if
    end if
 
    ! the timestep related to coll terms: 1/(rho_n rho_c alpha)
    if(dtcollpar>0d0 .and. has_collisions()) then
        call coll_get_dt(w,x,ixi^l,ixo^l,dtnew)
    endif
 
    if(twofl_radiative_cooling_c) then
      call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x,rc_fl_c)
    end if
    if(twofl_radiative_cooling_n) then
      call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x,rc_fl_n)
    end if
!
!    if(twofl_viscosity) then
!      call viscosity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
!    end if
!
    if(twofl_gravity) then
      call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
    end if
    if(twofl_hyperdiffusivity) then
      call hyperdiffusivity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
    end if
 
 
  end subroutine twofl_get_dt
 
  pure function has_collisions() result(res)
    logical :: res
    res = .not. twofl_alpha_coll_constant .or. twofl_alpha_coll >0d0
  end function has_collisions
 
  subroutine coll_get_dt(w,x,ixI^L,ixO^L,dtnew)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: w(ixi^s,1:nw)
    double precision, intent(in)    :: x(ixi^s,1:ndim)
    double precision, intent(inout) :: dtnew
 
    double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
    double precision, allocatable :: gamma_rec(:^d&), gamma_ion(:^D&)
    double precision :: max_coll_rate
 
    call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
    call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
 
    call get_alpha_coll(ixi^l, ixo^l, w, x, alpha)
    max_coll_rate = maxval(alpha(ixo^s) * max(rhon(ixo^s), rhoc(ixo^s)))
 
    if(twofl_coll_inc_ionrec) then
       allocate(gamma_ion(ixi^s), gamma_rec(ixi^s)) 
       call get_gamma_ion_rec(ixi^l, ixo^l, w, x, gamma_rec, gamma_ion)
       max_coll_rate=max(max_coll_rate, maxval(gamma_ion(ixo^s)), maxval(gamma_rec(ixo^s))) 
       deallocate(gamma_ion, gamma_rec) 
    endif
    dtnew = min(dtcollpar/max_coll_rate, dtnew)
 
  end subroutine coll_get_dt
 
  ! Add geometrical source terms to w
  subroutine twofl_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt, dtfactor,x(ixI^S,1:ndim)
    double precision, intent(inout) :: wCT(ixI^S,1:nw), wprim(ixI^S,1:nw), w(ixI^S,1:nw)
 
    integer          :: iw,idir, h1x^L{^NOONED, h2x^L}
    double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),rho(ixI^S)
 
    integer :: mr_,mphi_ ! Polar var. names
    integer :: br_,bphi_
 
    ! charges
 
    mr_=mom_c(1); mphi_=mom_c(1)-1+phi_  ! Polar var. names
    br_=mag(1); bphi_=mag(1)-1+phi_
    call get_rhoc_tot(wct,x,ixi^l,ixo^l,rho)
 
    select case (coordinate)
    case (cylindrical)
      call twofl_get_p_c_total(wct,x,ixi^l,ixo^l,tmp)
 
      if(phi_>0) then
        w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
                  wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/rho(ixo^s))
        w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
                 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/rho(ixo^s) &
                 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
        if(.not.stagger_grid) then
          w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
                   (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
                   -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
                   /rho(ixo^s)
        end if
      else
        w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
      end if
      if(twofl_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,psi_)/x(ixo^s,1)
    case (spherical)
       h1x^l=ixo^l-kr(1,^d); {^nooned h2x^l=ixo^l-kr(2,^d);}
       call twofl_get_p_c_total(wct,x,ixi^l,ixo^l,tmp1)
       tmp(ixo^s)=tmp1(ixo^s)
       if(b0field) then
         tmp2(ixo^s)=sum(block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=ndim+1)
         tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
       end if
       ! m1
       tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
                  *(block%surfaceC(ixo^s,1)-block%surfaceC(h1x^s,1))/block%dvolume(ixo^s)
       if(ndir>1) then
         do idir=2,ndir
           tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom_c(idir))**2/rho(ixo^s)-wct(ixo^s,mag(idir))**2
           if(b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
         end do
       end if
       w(ixo^s,mom_c(1))=w(ixo^s,mom_c(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
       ! b1
       if(twofl_glm) then
         w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,psi_)
       end if
 
       {^nooned
       ! m2
       tmp(ixo^s)=tmp1(ixo^s)
       if(b0field) then
         tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
       end if
       ! This will make hydrostatic p=const an exact solution
       w(ixo^s,mom_c(2))=w(ixo^s,mom_c(2))+qdt*tmp(ixo^s) &
            *(block%surfaceC(ixo^s,2)-block%surfaceC(h2x^s,2)) &
            /block%dvolume(ixo^s)
       tmp(ixo^s)=-(wct(ixo^s,mom_c(1))*wct(ixo^s,mom_c(2))/rho(ixo^s) &
            -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
       if (b0field) then
          tmp(ixo^s)=tmp(ixo^s)+block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
               +wct(ixo^s,mag(1))*block%B0(ixo^s,2,0)
       end if
       if(ndir==3) then
         tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom_c(3))**2/rho(ixo^s) &
              -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
         if (b0field) then
            tmp(ixo^s)=tmp(ixo^s)-2.0d0*block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
                 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
         end if
       end if
       w(ixo^s,mom_c(2))=w(ixo^s,mom_c(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
       ! b2
       if(.not.stagger_grid) then
         tmp(ixo^s)=(wct(ixo^s,mom_c(1))*wct(ixo^s,mag(2)) &
              -wct(ixo^s,mom_c(2))*wct(ixo^s,mag(1)))/rho(ixo^s)
         if(b0field) then
           tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom_c(1))*block%B0(ixo^s,2,0) &
                -wct(ixo^s,mom_c(2))*block%B0(ixo^s,1,0))/rho(ixo^s)
         end if
         if(twofl_glm) then
           tmp(ixo^s)=tmp(ixo^s) &
                + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,psi_)
         end if
         w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
       end if
       }
 
       if(ndir==3) then
         ! m3
         tmp(ixo^s)=-(wct(ixo^s,mom_c(3))*wct(ixo^s,mom_c(1))/rho(ixo^s) &
              -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
              -(wct(ixo^s,mom_c(2))*wct(ixo^s,mom_c(3))/rho(ixo^s) &
              -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
              *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
         if (b0field) then
            tmp(ixo^s)=tmp(ixo^s)+block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
                 +wct(ixo^s,mag(1))*block%B0(ixo^s,3,0) {^nooned &
                 +(block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
                 +wct(ixo^s,mag(2))*block%B0(ixo^s,3,0)) &
                 *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
         end if
         w(ixo^s,mom_c(3))=w(ixo^s,mom_c(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
         ! b3
         if(.not.stagger_grid) then
           tmp(ixo^s)=(wct(ixo^s,mom_c(1))*wct(ixo^s,mag(3)) &
                -wct(ixo^s,mom_c(3))*wct(ixo^s,mag(1)))/rho(ixo^s) {^nooned &
                -(wct(ixo^s,mom_c(3))*wct(ixo^s,mag(2)) &
                -wct(ixo^s,mom_c(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
                /(rho(ixo^s)*dsin(x(ixo^s,2))) }
           if (b0field) then
              tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom_c(1))*block%B0(ixo^s,3,0) &
                   -wct(ixo^s,mom_c(3))*block%B0(ixo^s,1,0))/rho(ixo^s){^nooned &
                   -(wct(ixo^s,mom_c(3))*block%B0(ixo^s,2,0) &
                   -wct(ixo^s,mom_c(2))*block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
                   /(rho(ixo^s)*dsin(x(ixo^s,2))) }
           end if
           w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
         end if
       end if
    end select
 
    ! neutrals 
    !TODO no dust: see and implement them from hd/mod_hd_phys !
    !uncomment cartesian expansion
    call get_rhon_tot(wct,x,ixi^l,ixo^l,rho)
    call twofl_get_pthermal_n(wct, x, ixi^l, ixo^l, tmp1)
 
    select case (coordinate)
!    case(Cartesian_expansion)
!      !the user provides the functions of exp_factor and del_exp_factor
!      if(associated(usr_set_surface)) call usr_set_surface(ixI^L,x,block%dx,exp_factor,del_exp_factor,exp_factor_primitive)
!      tmp(ixO^S) = tmp1(ixO^S)*del_exp_factor(ixO^S)/exp_factor(ixO^S)
!      w(ixO^S,mom(1)) = w(ixO^S,mom(1)) + qdt*tmp(ixO^S)
 
    case (cylindrical)
      mr_   = mom_n(r_)
      if (phi_ > 0) then
         where (rho(ixo^s) > 0d0)
            tmp(ixo^s) = tmp1(ixo^s) + wct(ixo^s, mphi_)**2 / rho(ixo^s)
            w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, r_)
         end where
         ! s[mphi]=(-mphi*mr/rho)/radius
         where (rho(ixo^s) > 0d0)
             tmp(ixo^s) = -wct(ixo^s, mphi_) * wct(ixo^s, mr_) / rho(ixo^s)
             w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * tmp(ixo^s) / x(ixo^s, r_)
         end where
      else
         ! s[mr]=2pthermal/radius
         w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp1(ixo^s) / x(ixo^s, r_)
      end if
    case (spherical)
       if(phi_>0) mphi_ = mom_n(phi_)
       h1x^l=ixo^l-kr(1,^d); {^nooned h2x^l=ixo^l-kr(2,^d);}
       ! s[mr]=((mtheta**2+mphi**2)/rho+2*p)/r
       tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
            *(block%surfaceC(ixo^s, 1) - block%surfaceC(h1x^s, 1)) &
            /block%dvolume(ixo^s)
       if (ndir > 1) then
         do idir = 2, ndir
           tmp(ixo^s) = tmp(ixo^s) + wct(ixo^s, mom_n(idir))**2 / rho(ixo^s)
         end do
       end if
       w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, 1)
 
       {^nooned
       ! s[mtheta]=-(mr*mtheta/rho)/r+cot(theta)*(mphi**2/rho+p)/r
       tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
            * (block%surfaceC(ixo^s, 2) - block%surfaceC(h2x^s, 2)) &
            / block%dvolume(ixo^s)
       if (ndir == 3) then
          tmp(ixo^s) = tmp(ixo^s) + (wct(ixo^s, mom_n(3))**2 / rho(ixo^s)) / tan(x(ixo^s, 2))
       end if
       tmp(ixo^s) = tmp(ixo^s) - (wct(ixo^s, mom_n(2)) * wct(ixo^s, mr_)) / rho(ixo^s)
       w(ixo^s, mom_n(2)) = w(ixo^s, mom_n(2)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
 
       if (ndir == 3) then
         ! s[mphi]=-(mphi*mr/rho)/r-cot(theta)*(mtheta*mphi/rho)/r
         tmp(ixo^s) = -(wct(ixo^s, mom_n(3)) * wct(ixo^s, mr_)) / rho(ixo^s)&
                     - (wct(ixo^s, mom_n(2)) * wct(ixo^s, mom_n(3))) / rho(ixo^s) / tan(x(ixo^s, 2))
         w(ixo^s, mom_n(3)) = w(ixo^s, mom_n(3)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
       end if
       }
    end select
 
    contains
      subroutine twofl_get_p_c_total(w,x,ixI^L,ixO^L,p)
        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)   :: p(ixI^S)
    
        call twofl_get_pthermal_c(w,x,ixi^l,ixo^l,p)
    
        p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
    
      end subroutine twofl_get_p_c_total
 
  end subroutine twofl_add_source_geom
 
  subroutine twofl_get_temp_c_pert_from_etot(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
 
    ! store pe1 in res
    res(ixo^s)=(gamma_1*(w(ixo^s,e_c_)&
           - twofl_kin_en_c(w,ixi^l,ixo^l)&
           - twofl_mag_en(w,ixi^l,ixo^l)))
    if(has_equi_pe_c0) then
      res(ixo^s) = res(ixo^s) + block%equi_vars(ixo^s,equi_pe_c0_,b0i)
      if(has_equi_rho_c0) then
        res(ixo^s) = res(ixo^s)/(rc * (w(ixo^s,rho_c_)+ block%equi_vars(ixo^s,equi_rho_c0_,b0i))) - &
                      block%equi_vars(ixo^s,equi_pe_c0_,b0i)/(rc * block%equi_vars(ixo^s,equi_rho_c0_,b0i))
      else
        ! infinite equi temperature with p0 and 0 density
        res(ixo^s) = 0d0
      endif
    else
      res(ixo^s) = res(ixo^s)/(rc * w(ixo^s,rho_c_))
    endif
 
  end subroutine twofl_get_temp_c_pert_from_etot
 
  !> Compute 2 times total magnetic energy
  function twofl_mag_en_all(w, ixI^L, ixO^L) result(mge)
    use mod_global_parameters
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s, nw)
    double precision              :: mge(ixo^s)
 
    if (b0field) then
      mge(ixo^s) = sum((w(ixo^s, mag(:))+block%B0(ixo^s,:,b0i))**2, dim=ndim+1)
    else
      mge(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=ndim+1)
    end if
  end function twofl_mag_en_all
 
  !> Compute full magnetic field by direction
  function twofl_mag_i_all(w, ixI^L, ixO^L,idir) result(mgf)
    use mod_global_parameters
    integer, intent(in)           :: ixi^l, ixo^l, idir
    double precision, intent(in)  :: w(ixi^s, nw)
    double precision              :: mgf(ixo^s)
 
    if (b0field) then
      mgf(ixo^s) = w(ixo^s, mag(idir))+block%B0(ixo^s,idir,b0i)
    else
      mgf(ixo^s) = w(ixo^s, mag(idir))
    end if
  end function twofl_mag_i_all
 
  !> Compute evolving magnetic energy
  function twofl_mag_en(w, ixI^L, ixO^L) result(mge)
    use mod_global_parameters, only: nw, ndim
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s, nw)
    double precision              :: mge(ixo^s)
 
    mge(ixo^s) = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
  end function twofl_mag_en
 
  !> compute kinetic energy of neutrals
  function twofl_kin_en_n(w, ixI^L, ixO^L) result(ke)
    use mod_global_parameters, only: nw, ndim,block
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s, nw)
    double precision              :: ke(ixo^s)
 
    if(has_equi_rho_n0) then
      ke(ixo^s) = 0.5d0 * sum(w(ixo^s, mom_n(:))**2, dim=ndim+1) / (w(ixo^s, rho_n_) + block%equi_vars(ixo^s,equi_rho_n0_,0))
    else
      ke(ixo^s) = 0.5d0 * sum(w(ixo^s, mom_n(:))**2, dim=ndim+1) / w(ixo^s, rho_n_)
    endif
 
  end function twofl_kin_en_n
 
  subroutine twofl_get_temp_n_pert_from_etot(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out):: res(ixI^S)
 
    ! store pe1 in res
    res(ixo^s)=(gamma_1*(w(ixo^s,e_c_)- twofl_kin_en_c(w,ixi^l,ixo^l)))
    if(has_equi_pe_n0) then
      res(ixo^s) = res(ixo^s) + block%equi_vars(ixo^s,equi_pe_n0_,b0i)
      if(has_equi_rho_n0) then
        res(ixo^s) = res(ixo^s)/(rn * (w(ixo^s,rho_n_)+ block%equi_vars(ixo^s,equi_rho_n0_,b0i))) - &
                      block%equi_vars(ixo^s,equi_pe_n0_,b0i)/(rn * block%equi_vars(ixo^s,equi_rho_n0_,b0i))
      else
        ! infinite equi temperature with p0 and 0 density
        res(ixo^s) = 0d0
      endif
    else
      res(ixo^s) = res(ixo^s)/(rn * w(ixo^s,rho_n_))
    endif
 
  end subroutine twofl_get_temp_n_pert_from_etot
 
  !> compute kinetic energy of charges
  !> w are conserved variables
  function twofl_kin_en_c(w, ixI^L, ixO^L) result(ke)
    use mod_global_parameters, only: nw, ndim,block
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s, nw)
    double precision              :: ke(ixo^s)
 
    if(has_equi_rho_c0) then
      ke(ixo^s) = 0.5d0 * sum(w(ixo^s, mom_c(:))**2, dim=ndim+1) / (w(ixo^s, rho_c_) + block%equi_vars(ixo^s,equi_rho_c0_,0))
    else
      ke(ixo^s) = 0.5d0 * sum(w(ixo^s, mom_c(:))**2, dim=ndim+1) / w(ixo^s, rho_c_)
    endif
  end function twofl_kin_en_c
 
  subroutine twofl_getv_hall(w,x,ixI^L,ixO^L,vHall)
    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(inout) :: vHall(ixI^S,1:3)
 
    integer          :: idir, idirmin
    double precision :: current(ixI^S,7-2*ndir:3)
    double precision :: rho(ixI^S)
 
    call get_rhoc_tot(w,x,ixi^l,ixo^l,rho)
    ! Calculate current density and idirmin
    call get_current(w,ixi^l,ixo^l,idirmin,current)
    vhall(ixo^s,1:3) = zero
    vhall(ixo^s,idirmin:3) = - twofl_etah*current(ixo^s,idirmin:3)
    do idir = idirmin, 3
       vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
    end do
 
  end subroutine twofl_getv_hall
 
! the following not used
!  subroutine twofl_getdt_Hall(w,x,ixI^L,ixO^L,dx^D,dthall)
!    use mod_global_parameters
!
!    integer, intent(in) :: ixI^L, ixO^L
!    double precision, intent(in)    :: dx^D
!    double precision, intent(in)    :: w(ixI^S,1:nw)
!    double precision, intent(in)    :: x(ixI^S,1:ndim)
!    double precision, intent(out)   :: dthall
!    !.. local ..
!    double precision :: dxarr(ndim)
!    double precision :: bmag(ixI^S)
!
!    dthall=bigdouble
!
!    ! because we have that in cmax now:
!    return
!
!    ^D&dxarr(^D)=dx^D;
!
!    if (.not. B0field) then
!       bmag(ixO^S)=sqrt(sum(w(ixO^S,mag(:))**2, dim=ndim+1))
!       bmag(ixO^S)=sqrt(sum((w(ixO^S,mag(:)) + block%B0(ixO^S,1:ndir,b0i))**2))
!    end if
!
!    if(slab_uniform) then
!      dthall=dtdiffpar*minval(dxarr(1:ndim))**2.0d0/(twofl_etah*maxval(bmag(ixO^S)/w(ixO^S,rho_c_)))
!    else
!      dthall=dtdiffpar*minval(block%ds(ixO^S,1:ndim))**2.0d0/(twofl_etah*maxval(bmag(ixO^S)/w(ixO^S,rho_c_)))
!    end if
!
!  end subroutine twofl_getdt_Hall
 
  subroutine twofl_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
    use mod_global_parameters
    use mod_usr_methods
    integer, intent(in)             :: ixI^L, ixO^L, idir
    double precision, intent(in)    :: qt
    double precision, intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
    double precision, intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
    type(state)                     :: s
    double precision                :: dB(ixI^S), dPsi(ixI^S)
 
    if(stagger_grid) then
      wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
      wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
      wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
      wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
    else
      ! Solve the Riemann problem for the linear 2x2 system for normal
      ! B-field and GLM_Psi according to Dedner 2002:
      ! This implements eq. (42) in Dedner et al. 2002 JcP 175
      ! Gives the Riemann solution on the interface
      ! for the normal B component and Psi in the GLM-MHD system.
      ! 23/04/2013 Oliver Porth
      db(ixo^s)   = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
      dpsi(ixo^s) = wrp(ixo^s,psi_) - wlp(ixo^s,psi_)
 
      wlp(ixo^s,mag(idir))   = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
           - 0.5d0/cmax_global * dpsi(ixo^s)
      wlp(ixo^s,psi_)       = 0.5d0 * (wrp(ixo^s,psi_) + wlp(ixo^s,psi_)) &
           - 0.5d0*cmax_global * db(ixo^s)
 
      wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
      wrp(ixo^s,psi_) = wlp(ixo^s,psi_)
 
      if(phys_total_energy) then
        wrc(ixo^s,e_c_)=wrc(ixo^s,e_c_)-half*wrc(ixo^s,mag(idir))**2
        wlc(ixo^s,e_c_)=wlc(ixo^s,e_c_)-half*wlc(ixo^s,mag(idir))**2
      end if
      wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
      wrc(ixo^s,psi_) = wlp(ixo^s,psi_)
      wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
      wlc(ixo^s,psi_) = wlp(ixo^s,psi_)
      ! modify total energy according to the change of magnetic field
      if(phys_total_energy) then
        wrc(ixo^s,e_c_)=wrc(ixo^s,e_c_)+half*wrc(ixo^s,mag(idir))**2
        wlc(ixo^s,e_c_)=wlc(ixo^s,e_c_)+half*wlc(ixo^s,mag(idir))**2
      end if
    end if
 
    if(associated(usr_set_wlr)) call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
 
  end subroutine twofl_modify_wlr
 
  subroutine twofl_boundary_adjust(igrid,psb)
    use mod_global_parameters
    integer, intent(in) :: igrid
    type(state), target :: psb(max_blocks)
 
    integer :: iB, idims, iside, ixO^L, i^D
 
    block=>ps(igrid)
    ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
    do idims=1,ndim
       ! to avoid using as yet unknown corner info in more than 1D, we
       ! fill only interior mesh ranges of the ghost cell ranges at first,
       ! and progressively enlarge the ranges to include corners later
       do iside=1,2
          i^d=kr(^d,idims)*(2*iside-3);
          if (neighbor_type(i^d,igrid)/=1) cycle
          ib=(idims-1)*2+iside
          if(.not.boundary_divbfix(ib)) cycle
          if(any(typeboundary(:,ib)==bc_special)) then
            ! MF nonlinear force-free B field extrapolation and data driven
            ! require normal B of the first ghost cell layer to be untouched by
            ! fixdivB=0 process, set boundary_divbfix_skip(iB)=1 in par file
            select case (idims)
            {case (^d)
               if (iside==2) then
                  ! maximal boundary
                  ixomin^dd=ixghi^d+1-nghostcells+boundary_divbfix_skip(2*^d)^d%ixOmin^dd=ixglo^dd;
                  ixomax^dd=ixghi^dd;
               else
                  ! minimal boundary
                  ixomin^dd=ixglo^dd;
                  ixomax^dd=ixglo^d-1+nghostcells-boundary_divbfix_skip(2*^d-1)^d%ixOmax^dd=ixghi^dd;
               end if \}
            end select
            call fixdivb_boundary(ixg^ll,ixo^l,psb(igrid)%w,psb(igrid)%x,ib)
          end if
       end do
    end do
 
  end subroutine twofl_boundary_adjust
 
  subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
    use mod_global_parameters
 
    integer, intent(in) :: ixG^L,ixO^L,iB
    double precision, intent(inout) :: w(ixG^S,1:nw)
    double precision, intent(in) :: x(ixG^S,1:ndim)
 
    double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
    integer :: ix^D,ixF^L
 
    select case(ib)
     case(1)
       ! 2nd order CD for divB=0 to set normal B component better
       {^iftwod
       ixfmin1=ixomin1+1
       ixfmax1=ixomax1+1
       ixfmin2=ixomin2+1
       ixfmax2=ixomax2-1
       if(slab_uniform) then
         dx1x2=dxlevel(1)/dxlevel(2)
         do ix1=ixfmax1,ixfmin1,-1
           w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
            +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
                    w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
         enddo
       else
         do ix1=ixfmax1,ixfmin1,-1
           w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
             w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
           +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
           -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
             block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
            /block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
         end do
       end if
       }
       {^ifthreed
       ixfmin1=ixomin1+1
       ixfmax1=ixomax1+1
       ixfmin2=ixomin2+1
       ixfmax2=ixomax2-1
       ixfmin3=ixomin3+1
       ixfmax3=ixomax3-1
       if(slab_uniform) then
         dx1x2=dxlevel(1)/dxlevel(2)
         dx1x3=dxlevel(1)/dxlevel(3)
         do ix1=ixfmax1,ixfmin1,-1
           w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
                     w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
             +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
                     w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
             +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
                     w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
         end do
       else
         do ix1=ixfmax1,ixfmin1,-1
           w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
          ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
           +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
           -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
             w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
             block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
           +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
           -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
            /block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
         end do
       end if
       }
     case(2)
       {^iftwod
       ixfmin1=ixomin1-1
       ixfmax1=ixomax1-1
       ixfmin2=ixomin2+1
       ixfmax2=ixomax2-1
       if(slab_uniform) then
         dx1x2=dxlevel(1)/dxlevel(2)
         do ix1=ixfmin1,ixfmax1
           w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
            -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
                    w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
         enddo
       else
         do ix1=ixfmin1,ixfmax1
           w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
             w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
           -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
           +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
             block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
            /block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
         end do
       end if
       }
       {^ifthreed
       ixfmin1=ixomin1-1
       ixfmax1=ixomax1-1
       ixfmin2=ixomin2+1
       ixfmax2=ixomax2-1
       ixfmin3=ixomin3+1
       ixfmax3=ixomax3-1
       if(slab_uniform) then
         dx1x2=dxlevel(1)/dxlevel(2)
         dx1x3=dxlevel(1)/dxlevel(3)
         do ix1=ixfmin1,ixfmax1
           w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
                     w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
             -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
                     w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
             -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
                     w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
         end do
       else
         do ix1=ixfmin1,ixfmax1
           w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
          ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
             block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
           -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
           +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
             w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
             block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
           -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
           +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
            /block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
         end do
       end if
       }
     case(3)
       {^iftwod
       ixfmin1=ixomin1+1
       ixfmax1=ixomax1-1
       ixfmin2=ixomin2+1
       ixfmax2=ixomax2+1
       if(slab_uniform) then
         dx2x1=