mod_rmhd_phys.t

Go to the documentation of this file.

!> Radiation-magneto-hydrodynamics module
module mod_rmhd_phys
 
#include "amrvac.h"
 
  use mod_global_parameters, only: std_len, const_c
  use mod_thermal_conduction, only: tc_fluid
  use mod_thermal_emission, only: te_fluid
  use mod_physics
  use mod_comm_lib, only: mpistop
  use mod_functions_bfield, only: get_divb, mag
 
  implicit none
  private
 
  !> The adiabatic index
  double precision, public                :: rmhd_gamma = 5.d0/3.0d0
  !> The adiabatic constant
  double precision, public                :: rmhd_adiab = 1.0d0
  !> The MHD resistivity
  double precision, public                :: rmhd_eta = 0.0d0
  !> The MHD hyper-resistivity
  double precision, public                :: rmhd_eta_hyper = 0.0d0
  !> Hall resistivity
  double precision, public                :: rmhd_etah = 0.0d0
  !> The small_est allowed energy
  double precision, protected             :: small_e
  !> The smallest allowed radiation energy
  double precision, public, protected     :: small_r_e = 0.d0
  !> Height of the mask used in the TRAC method
  double precision, public, protected     :: rmhd_trac_mask = 0.d0
  !> GLM-MHD parameter: ratio of the diffusive and advective time scales for div b
  !> taking values within [0, 1]
  double precision, public                :: rmhd_glm_alpha = 0.5d0
  !> The thermal conductivity kappa in hyperbolic thermal conduction
  double precision, public                :: hypertc_kappa
  !> Coefficient of diffusive divB cleaning
  double precision                        :: divbdiff     = 0.8d0
  !> Helium abundance over Hydrogen
  double precision, public, protected  :: he_abundance=0.1d0
  !> Ionization fraction of H
  !> H_ion_fr = H+/(H+ + H)
  double precision, public, protected  :: h_ion_fr=1d0
  !> Ionization fraction of He
  !> He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
  double precision, public, protected  :: he_ion_fr=1d0
  !> Ratio of number He2+ / number He+ + He2+
  !> He_ion_fr2 = He2+/(He2+ + He+)
  double precision, public, protected  :: he_ion_fr2=1d0
  ! used for eq of state when it is not defined by units,
  ! the units do not contain terms related to ionization fraction
  ! and it is p = RR * rho * T
  double precision, public, protected  :: rr=1d0
  !> gamma minus one and its inverse
  double precision :: gamma_1, inv_gamma_1
  !> inverse of squared speed of light c0 and reduced speed of light c
  double precision :: inv_squared_c0, inv_squared_c
  !> equi vars indices in the state%equi_vars array
  integer, public :: equi_rho0_ = -1
  integer, public :: equi_pe0_ = -1
  !> Number of tracer species
  integer, public, protected              :: rmhd_n_tracer = 0
  !> Index of the density (in the w array)
  integer, public, protected              :: rho_
  !> Indices of the momentum density
  integer, allocatable, public, protected :: mom(:)
  !> Indices of the momentum density for the form of better vectorization
  integer, public, protected              :: ^c&m^C_
  !> Index of the energy density (-1 if not present)
  integer, public, protected              :: e_
  !> Index of the radiation energy
  integer, public, protected              :: r_e
  !> Indices of the magnetic field for the form of better vectorization
  integer, public, protected              :: ^c&b^C_
  !> Index of the gas pressure (-1 if not present) should equal e_
  integer, public, protected              :: p_
  !> Index of the heat flux q
  integer, public, protected :: q_
  !> Indices of the GLM psi
  integer, public, protected :: psi_
  !> Indices of temperature
  integer, public, protected :: te_
  !> Index of the cutoff temperature for the TRAC method
  integer, public, protected              :: tcoff_
  integer, public, protected              :: tweight_
  !> Indices of the tracers
  integer, allocatable, public, protected :: tracer(:)
  !> The number of waves
  integer :: nwwave=8
  !> Method type in a integer for good performance
  integer :: type_divb
  !> To skip * layer of ghost cells during divB=0 fix for boundary
  integer, public, protected :: boundary_divbfix_skip(2*^nd)=0
  ! 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
  !> Whether an energy equation is used
  logical, public, protected              :: rmhd_energy = .true.
  !> Whether thermal conduction is used
  logical, public, protected              :: rmhd_thermal_conduction = .false.
  !> Whether thermal conduction is used
  logical, public, protected              :: rmhd_hyperbolic_thermal_conduction = .false.
  !> Whether viscosity is added
  logical, public, protected              :: rmhd_viscosity = .false.
  !> Whether gravity is added
  logical, public, protected              :: rmhd_gravity = .false.
  !> Whether particles module is added
  logical, public, protected              :: rmhd_particles = .false.
  !> Whether GLM-MHD is used to control div B
  logical, public, protected              :: rmhd_glm = .false.
  !> Whether extended GLM-MHD is used with additional sources
  logical, public, protected              :: rmhd_glm_extended = .true.
  !> Whether TRAC method is used
  logical, public, protected              :: rmhd_trac = .false.
  !> Which TRAC method is used
  integer, public, protected              :: rmhd_trac_type=1
  !> Distance between two adjacent traced magnetic field lines (in finest cell size)
  integer, public, protected              :: rmhd_trac_finegrid=4
  !> Whether divB cleaning sources are added splitting from fluid solver
  logical, public, protected              :: source_split_divb = .false.
  !> Whether plasma is partially ionized
  logical, public, protected              :: rmhd_partial_ionization = .false.
  !> Whether CAK radiation line force is activated
  logical, public, protected              :: rmhd_cak_force = .false.
  !> MHD fourth order
  logical, public, protected              :: rmhd_4th_order = .false.
  !> whether split off equilibrium density
  logical, public :: has_equi_rho0 = .false.
  !> whether split off equilibrium thermal pressure
  logical, public :: has_equi_pe0 = .false.
  logical, public :: rmhd_equi_thermal = .false.
  !> whether dump full variables (when splitting is used) in a separate dat file
  logical, public, protected :: rmhd_dump_full_vars = .false.
  !> Whether divB is computed with a fourth order approximation
  integer, public, protected :: rmhd_divb_nth = 1
  !> Use a compact way to add resistivity
  logical :: compactres = .false.
  !> Add divB wave in Roe solver
  logical, public :: divbwave = .true.
  !> clean initial divB
  logical, public :: clean_initial_divb = .false.
  !> Formalism to treat radiation
  character(len=8), public :: rmhd_radiation_formalism = 'fld'
  !> In the case of no rmhd_energy, how to compute pressure
  character(len=8), public :: rmhd_pressure = 'Trad'
  !> Treat radiation fld_Rad_force
  logical, public, protected :: rmhd_radiation_force = .true.
  !> Treat radiation-gas energy interaction
  logical, public, protected :: rmhd_energy_interact = .true.
  !> Treat radiation energy diffusion
  logical, public, protected :: rmhd_radiation_diffusion = .true.
  !> Treat radiation advection
  logical, public, protected :: rmhd_radiation_advection = .true.
  !> Do a running mean over the radiation pressure when determining dt
  logical, protected :: radio_acoustic_filter = .false.
  integer, protected :: size_ra_filter = 1
  !> kb/(m_p mu)* 1/a_rad**4,
  double precision, public :: kbmpmua4
  !> Use the speed of light to calculate the timestep, usefull for debugging
  logical :: dt_c = .false.
  ! remove the below flag  and assume default value = .false.
  ! when eq state properly implemented everywhere
  ! and not anymore through units
  logical, public, protected :: eq_state_units = .true.
  !> To control divB=0 fix for boundary
  logical, public, protected :: boundary_divbfix(2*^nd)=.true.
  !> B0 field is force-free
  logical, public, protected :: b0field_forcefree=.true.
  !> Whether an total energy equation is used
  logical :: total_energy = .true.
  !> Whether an internal or hydrodynamic energy equation is used
  logical, public :: partial_energy = .false.
  !> Whether gravity work is included in energy equation
  logical :: gravity_energy
  !> gravity work is calculated use density times velocity or conservative momentum
  logical :: gravity_rhov = .false.
  !> 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'
  !> Update all equations due to divB cleaning
  character(len=std_len) ::    typedivbdiff = 'all'
  !> type of fluid for thermal conduction
  type(tc_fluid), public, allocatable     :: tc_fl
  !> type of fluid for thermal emission synthesis
  type(te_fluid), public, allocatable     :: te_fl_rmhd
 
  procedure(sub_convert), pointer      :: rmhd_to_primitive  => null()
  procedure(sub_convert), pointer      :: rmhd_to_conserved  => null()
  procedure(sub_small_values), pointer :: rmhd_handle_small_values => null()
  procedure(sub_get_pthermal), pointer :: rmhd_get_pthermal  => null()
  procedure(sub_get_pthermal), pointer :: rmhd_get_rfactor   => null()
  procedure(sub_get_pthermal), pointer :: rmhd_get_temperature=> null()
  ! Public methods
  public :: rmhd_phys_init
  public :: rmhd_get_pthermal
  public :: rmhd_get_temperature
  public :: rmhd_get_v
  public :: rmhd_get_rho
  public :: rmhd_to_conserved
  public :: rmhd_to_primitive
  public :: rmhd_e_to_ei
  public :: rmhd_ei_to_e
  public :: rmhd_face_to_center
  public :: get_divb
  public :: get_current
  public :: get_normalized_divb
  public :: b_from_vector_potential
  public :: rmhd_mag_en_all
  {^nooned
  public :: rmhd_clean_divb_multigrid
  }
  public :: rmhd_get_pradiation
  public :: rmhd_get_pthermal_plus_pradiation
  public :: rmhd_get_tgas
  public :: rmhd_get_trad
  public :: rmhd_set_mg_bounds
 
contains
 
  !> Read this module"s parameters from a file
  subroutine rmhd_read_params(files)
    use mod_global_parameters
    use mod_particles, only: particles_eta, particles_etah
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /rmhd_list/ rmhd_energy, rmhd_n_tracer, rmhd_gamma, rmhd_adiab,&
      rmhd_eta, rmhd_eta_hyper, rmhd_etah, rmhd_glm_alpha, rmhd_glm_extended,&
      rmhd_thermal_conduction, rmhd_gravity,&
      rmhd_viscosity, rmhd_4th_order, typedivbfix, source_split_divb, divbdiff,&
      typedivbdiff, type_ct, compactres, divbwave, he_abundance,&
      h_ion_fr, he_ion_fr, he_ion_fr2, eq_state_units, si_unit, b0field ,rmhd_dump_full_vars,&
      b0field_forcefree, bdip, bquad, boct, busr, rmhd_particles, rmhd_partial_ionization,&
      particles_eta, particles_etah,has_equi_rho0, has_equi_pe0,rmhd_equi_thermal,&
      boundary_divbfix, boundary_divbfix_skip, rmhd_divb_nth, clean_initial_divb,&
      rmhd_trac, rmhd_trac_type, rmhd_trac_mask, rmhd_trac_finegrid, rmhd_cak_force,&
      rmhd_hyperbolic_thermal_conduction, &
      rmhd_pressure, rmhd_radiation_formalism,&
      rmhd_radiation_force, rmhd_energy_interact, rmhd_radiation_diffusion,&
      rmhd_radiation_advection, radio_acoustic_filter, size_ra_filter, dt_c
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status="old")
       read(unitpar, rmhd_list, end=111)
111    close(unitpar)
    end do
 
  end subroutine rmhd_read_params
 
  !> Write this module's parameters to a snapsoht
  subroutine rmhd_write_info(fh)
    use mod_global_parameters
    integer, intent(in)                 :: fh
    integer                             :: er
    integer, parameter                  :: n_par = 1
    double precision                    :: values(n_par)
    integer, dimension(MPI_STATUS_SIZE) :: st
    character(len=name_len)             :: names(n_par)
 
    call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
 
    names(1) = "gamma"
    values(1) = rmhd_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 rmhd_write_info
 
  subroutine rmhd_phys_init()
    use mod_global_parameters
    use mod_thermal_conduction
    use mod_viscosity, only: viscosity_init
    use mod_gravity, only: gravity_init
    use mod_particles, only: particles_init, particles_eta, particles_etah
    use mod_fld
    use mod_afld
    use mod_supertimestepping, only: sts_init, add_sts_method,&
            set_conversion_methods_to_head, set_error_handling_to_head
    use mod_cak_force, only: cak_init
    use mod_ionization_degree
    use mod_usr_methods, only: usr_rfactor
    {^nooned
    use mod_multigrid_coupling
    }
 
    integer :: itr, idir
 
    call rmhd_read_params(par_files)
 
    if(.not.eq_state_units) then
      if(rmhd_partial_ionization) then
        rmhd_partial_ionization=.false.
        if(mype==0) write(*,*) 'WARNING: set rmhd_partial_ionization=F when eq_state_units=F'
      end if
    end if
 
    if(rmhd_hyperbolic_thermal_conduction) then
      rmhd_thermal_conduction=.false.
      if(mype==0) write(*,*) 'WARNING: turn off parabolic TC when using hyperbolic TC'
    end if
 
    physics_type = "rmhd"
    phys_energy=rmhd_energy
    phys_trac=rmhd_trac
    phys_trac_type=rmhd_trac_type
    phys_partial_ionization=rmhd_partial_ionization
 
    phys_gamma = rmhd_gamma
    phys_trac_finegrid=rmhd_trac_finegrid
 
    if(rmhd_energy) then
      partial_energy=.false.
      total_energy=.true.
    else
      total_energy=.false.
    end if
    phys_total_energy=total_energy
    if(rmhd_energy) then
      gravity_energy=.true.
      if(has_equi_rho0) then
        gravity_rhov=.true.
      end if
    else
      gravity_energy=.false.
    end if
 
    {^ifoned
    if(rmhd_trac .and. rmhd_trac_type .gt. 2) then
      rmhd_trac_type=1
      if(mype==0) write(*,*) 'WARNING: reset rmhd_trac_type=1 for 1D simulation'
    end if
    }
    if(rmhd_trac .and. rmhd_trac_type .le. 4) then
      rmhd_trac_mask=bigdouble
      if(mype==0) write(*,*) 'WARNING: set rmhd_trac_mask==bigdouble for global TRAC method'
    end if
    phys_trac_mask=rmhd_trac_mask
 
    ! set default gamma for polytropic/isothermal process
    use_particles=rmhd_particles
    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 => rmhd_clean_divb_multigrid
    }
    case ('glm')
      rmhd_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')
      rmhd_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),stop_indices(number_species))
    ! set the index of the first flux variable for species 1
    start_indices(1)=1
    ! Determine flux variables
    rho_ = var_set_rho()
 
    allocate(mom(ndir))
    mom(:) = var_set_momentum(ndir)
    m^c_=mom(^c);
 
    ! Set index of energy variable
    if (rmhd_energy) then
      nwwave = 8
      e_     = var_set_energy() ! energy density
      p_     = e_               ! gas pressure
    else
      nwwave = 7
      e_     = -1
      p_     = -1
    end if
 
    allocate(mag(ndir))
    mag(:) = var_set_bfield(ndir)
    b^c_=mag(^c);
 
    if (rmhd_glm) then
      psi_ = var_set_fluxvar('psi', 'psi', need_bc=.false.)
    else
      psi_ = -1
    end if
 
    !> set radiation energy
    r_e = var_set_radiation_energy()
 
    if(rmhd_hyperbolic_thermal_conduction) then
      ! hyperbolic thermal conduction flux q
      q_ = var_set_q()
      need_global_cmax=.true.
    else
      q_=-1
    end if
 
    allocate(tracer(rmhd_n_tracer))
    ! Set starting index of tracers
    do itr = 1, rmhd_n_tracer
      tracer(itr) = var_set_fluxvar("trc", "trp", itr, need_bc=.false.)
    end do
 
    !if(rmhd_hyperbolic_thermal_conduction) then
    !  ! hyperbolic thermal conduction flux q
    !  q_ = var_set_auxvar('q','q')
    !  need_global_cmax=.true.
    !else
    !  q_=-1
    !end if
 
    !  set temperature as an auxiliary variable to get ionization degree
    if(rmhd_partial_ionization) then
      te_ = var_set_auxvar('Te','Te')
    else
      te_ = -1
    end if
 
    ! set number of variables which need update ghostcells
    nwgc=nwflux+nwaux
 
    ! set the index of the last flux variable for species 1
    stop_indices(1)=nwflux
 
    ! set cutoff temperature when using the TRAC method, as well as an auxiliary weight
    tweight_ = -1
    if(rmhd_trac) then
      tcoff_ = var_set_wextra()
      iw_tcoff=tcoff_
      if(rmhd_trac_type .ge. 3) then
        tweight_ = var_set_wextra()
      endif
    else
      tcoff_ = -1
    end if
 
    ! set indices of equi vars and update number_equi_vars
    number_equi_vars = 0
    if(has_equi_rho0) then
      number_equi_vars = number_equi_vars + 1
      equi_rho0_ = number_equi_vars
      iw_equi_rho = equi_rho0_
    endif
    if(has_equi_pe0) then
      number_equi_vars = number_equi_vars + 1
      equi_pe0_ = number_equi_vars
      iw_equi_p = equi_pe0_
    endif
    ! determine number of stagger variables
    nws=ndim
 
    nvector      = 2 ! No. vector vars
    allocate(iw_vector(nvector))
    iw_vector(1) = mom(1) - 1   ! TODO: why like this?
    iw_vector(2) = mag(1) - 1   ! TODO: why like this?
 
    ! Check whether custom flux types have been defined
    if (.not. allocated(flux_type)) then
       allocate(flux_type(ndir, nwflux))
       flux_type = flux_default
    else if (any(shape(flux_type) /= [ndir, nwflux])) then
       call mpistop("phys_check error: flux_type has wrong shape")
    end if
 
    if(nwflux>mag(ndir)) then
      ! for flux of tracers, using hll flux
      flux_type(:,mag(ndir)+1:nwflux)=flux_hll
    end if
 
    if(ndim>1) then
      if(rmhd_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_rho             => rmhd_get_rho
    phys_get_dt              => rmhd_get_dt
    phys_get_cmax            => rmhd_get_cmax_origin
    phys_get_a2max           => rmhd_get_a2max
    phys_get_tcutoff         => rmhd_get_tcutoff
    phys_get_h_speed         => rmhd_get_h_speed
    if(has_equi_rho0) then
      phys_get_cbounds       => rmhd_get_cbounds_split_rho
    else
      phys_get_cbounds       => rmhd_get_cbounds
    end if
    if(has_equi_rho0) then
      phys_to_primitive      => rmhd_to_primitive_split_rho
      rmhd_to_primitive      => rmhd_to_primitive_split_rho
      phys_to_conserved      => rmhd_to_conserved_split_rho
      rmhd_to_conserved      => rmhd_to_conserved_split_rho
    else
      phys_to_primitive      => rmhd_to_primitive_origin
      rmhd_to_primitive      => rmhd_to_primitive_origin
      phys_to_conserved      => rmhd_to_conserved_origin
      rmhd_to_conserved      => rmhd_to_conserved_origin
    end if
    if(b0field.or.has_equi_rho0.or.has_equi_pe0) then
      phys_get_flux            => rmhd_get_flux_split
    else
      phys_get_flux            => rmhd_get_flux
    end if
    phys_get_v                 => rmhd_get_v
    if(b0field.or.has_equi_rho0) then
      phys_add_source_geom     => rmhd_add_source_geom_split
    else
      phys_add_source_geom     => rmhd_add_source_geom
    end if
    phys_add_source          => rmhd_add_source
    phys_check_params        => rmhd_check_params
    phys_write_info          => rmhd_write_info
 
    phys_handle_small_values => rmhd_handle_small_values_origin
    rmhd_handle_small_values => rmhd_handle_small_values_origin
    phys_check_w             => rmhd_check_w_origin
 
    phys_set_mg_bounds       => rmhd_set_mg_bounds
    phys_get_trad            => rmhd_get_trad
    phys_get_tgas            => rmhd_get_tgas
 
    phys_get_pthermal        => rmhd_get_pthermal_origin
    rmhd_get_pthermal         => rmhd_get_pthermal_origin
 
    if(number_equi_vars>0) then
      phys_set_equi_vars => set_equi_vars_grid
    endif
 
    if(type_divb==divb_glm) then
      phys_modify_wlr => rmhd_modify_wlr
    end if
 
    ! choose Rfactor in ideal gas law
    if(rmhd_partial_ionization) then
      rmhd_get_rfactor=>rfactor_from_temperature_ionization
      phys_update_temperature => rmhd_update_temperature
    else if(associated(usr_rfactor)) then
      rmhd_get_rfactor=>usr_rfactor
    else
      rmhd_get_rfactor=>rfactor_from_constant_ionization
    end if
 
    if(rmhd_partial_ionization) then
      rmhd_get_temperature => rmhd_get_temperature_from_te
    else
      if(has_equi_pe0 .and. has_equi_rho0) then
        rmhd_get_temperature => rmhd_get_temperature_from_etot_with_equi
      else
        rmhd_get_temperature => rmhd_get_temperature_from_etot
      end if
    end if
 
    ! if using ct stagger grid, boundary divb=0 is not done here
    if(stagger_grid) then
      phys_get_ct_velocity => rmhd_get_ct_velocity
      phys_update_faces => rmhd_update_faces
      phys_face_to_center => rmhd_face_to_center
      phys_modify_wlr => rmhd_modify_wlr
    else if(ndim>1) then
      phys_boundary_adjust => rmhd_boundary_adjust
    end if
 
    {^nooned
    ! clean initial divb
    if(clean_initial_divb) phys_clean_divb => rmhd_clean_divb_multigrid
    }
 
    ! derive units from basic units
    call rmhd_physical_units()
 
    !> Initiate radiation-closure module
    select case(rmhd_radiation_formalism)
    case('fld')
      call fld_init(he_abundance, rmhd_radiation_diffusion, rmhd_energy_interact, rmhd_gamma)
    case('afld')
      call afld_init(he_abundance, rmhd_radiation_diffusion, rmhd_gamma)
    case default
      call mpistop('Radiation formalism unknown')
    end select
 
    if(rmhd_hyperbolic_thermal_conduction) then
      hypertc_kappa=8.d-7*unit_temperature**3.5d0/unit_length/unit_density/unit_velocity**3
    end if
    if(.not. rmhd_energy .and. rmhd_thermal_conduction) then
      call mpistop("thermal conduction needs rmhd_energy=T")
    end if
    if(.not. rmhd_energy .and. rmhd_hyperbolic_thermal_conduction) then
      call mpistop("hyperbolic thermal conduction needs rmhd_energy=T")
    end if
 
    ! initialize thermal conduction module
    if (rmhd_thermal_conduction) then
      call sts_init()
      call tc_init_params(rmhd_gamma)
 
      allocate(tc_fl)
      call tc_get_mhd_params(tc_fl,tc_params_read_rmhd)
      call add_sts_method(rmhd_get_tc_dt_rmhd,rmhd_sts_set_source_tc_rmhd,e_,1,e_,1,.false.)
      if(has_equi_pe0 .and. has_equi_rho0) then
        tc_fl%get_temperature_from_conserved => rmhd_get_temperature_from_etot_with_equi
      else
        tc_fl%get_temperature_from_conserved => rmhd_get_temperature_from_etot
      end if
      if(has_equi_pe0 .and. has_equi_rho0) then
        tc_fl%get_temperature_from_eint => rmhd_get_temperature_from_eint_with_equi
        if(rmhd_equi_thermal) then
          tc_fl%has_equi = .true.
          tc_fl%get_temperature_equi => rmhd_get_temperature_equi
          tc_fl%get_rho_equi => rmhd_get_rho_equi
        else
          tc_fl%has_equi = .false.
        end if
      else
        tc_fl%get_temperature_from_eint => rmhd_get_temperature_from_eint
      end if
      call set_conversion_methods_to_head(rmhd_e_to_ei, rmhd_ei_to_e)
      call set_error_handling_to_head(rmhd_tc_handle_small_e)
      tc_fl%get_rho => rmhd_get_rho
      tc_fl%e_ = e_
      tc_fl%Tcoff_ = tcoff_
    end if
 
    allocate(te_fl_rmhd)
    te_fl_rmhd%get_rho=> rmhd_get_rho
    te_fl_rmhd%get_pthermal=> rmhd_get_pthermal
    te_fl_rmhd%get_var_Rfactor => rmhd_get_rfactor
{^ifthreed
    phys_te_images => rmhd_te_images
}
    ! Initialize viscosity module
    if (rmhd_viscosity) call viscosity_init(phys_wider_stencil)
 
    ! Initialize gravity module
    if(rmhd_gravity) then
      call gravity_init()
    end if
 
    ! Initialize particles module
    if(rmhd_particles) then
      call particles_init()
      if (particles_eta  < zero) particles_eta = rmhd_eta
      if (particles_etah < zero) particles_eta = rmhd_etah
      if(mype==0) then
         write(*,*) '*****Using particles:        with rmhd_eta, rmhd_etah :', rmhd_eta, rmhd_etah
         write(*,*) '*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
      end if
    end if
 
    ! initialize ionization degree table
    if(rmhd_partial_ionization) call ionization_degree_init()
 
    ! Initialize CAK radiation force module
    if (rmhd_cak_force) call cak_init(rmhd_gamma)
  end subroutine rmhd_phys_init
 
{^ifthreed
  subroutine rmhd_te_images
    use mod_global_parameters
    use mod_thermal_emission
 
    select case(convert_type)
      case('EIvtiCCmpi','EIvtuCCmpi')
        call get_euv_image(unitconvert,te_fl_rmhd)
      case('ESvtiCCmpi','ESvtuCCmpi')
        call get_euv_spectrum(unitconvert,te_fl_rmhd)
      case('SIvtiCCmpi','SIvtuCCmpi')
        call get_sxr_image(unitconvert,te_fl_rmhd)
      case('WIvtiCCmpi','WIvtuCCmpi')
        call get_whitelight_image(unitconvert,te_fl_rmhd)
      case default
        call mpistop("Error in synthesize emission: Unknown convert_type")
      end select
  end subroutine rmhd_te_images
}
 
!!start th cond
  ! wrappers for STS functions in thermal_conductivity module
  ! which take as argument the tc_fluid (defined in the physics module)
  subroutine  rmhd_sts_set_source_tc_rmhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
    use mod_global_parameters
    use mod_fix_conserve
    use mod_thermal_conduction, only: sts_set_source_tc_mhd
    integer, intent(in) :: ixi^l, ixo^l, igrid, nflux
    double precision, intent(in) ::  x(ixi^s,1:ndim)
    double precision, intent(inout) ::  wres(ixi^s,1:nw), w(ixi^s,1:nw)
    double precision, intent(in) :: my_dt
    logical, intent(in) :: fix_conserve_at_step
    call sts_set_source_tc_mhd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl)
  end subroutine rmhd_sts_set_source_tc_rmhd
 
  function rmhd_get_tc_dt_rmhd(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
    !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
    !where                      tc_k_para_i=tc_k_para*B_i**2/B**2
    !and                        T=p/rho
    use mod_global_parameters
    use mod_thermal_conduction, only: get_tc_dt_mhd
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
    double precision, intent(in) :: w(ixi^s,1:nw)
    double precision :: dtnew
 
    dtnew=get_tc_dt_mhd(w,ixi^l,ixo^l,dx^d,x,tc_fl)
  end function rmhd_get_tc_dt_rmhd
 
  subroutine rmhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l,ixo^l
    double precision, intent(inout) :: w(ixi^s,1:nw)
    double precision, intent(in)    :: x(ixi^s,1:ndim)
    integer, intent(in)    :: step
    character(len=140) :: error_msg
 
    write(error_msg,"(a,i3)") "Thermal conduction step ", step
    call rmhd_handle_small_ei(w,x,ixi^l,ixo^l,e_,error_msg)
  end subroutine rmhd_tc_handle_small_e
 
  ! fill in tc_fluid fields from namelist
  subroutine tc_params_read_rmhd(fl)
    use mod_global_parameters, only: unitpar,par_files
    type(tc_fluid), intent(inout) :: fl
    double precision :: tc_k_para=0d0
    double precision :: tc_k_perp=0d0
    integer                      :: n
    ! list parameters
    logical :: tc_perpendicular=.false.
    logical :: tc_saturate=.false.
    character(len=std_len)  :: tc_slope_limiter="MC"
 
    namelist /tc_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_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_params_read_rmhd
!!end th cond
 
  !> sets the equilibrium variables
  subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
    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, add splitted variables back to wnew
  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)
 
    if(has_equi_rho0) then
      wnew(ixo^s,rho_)=w(ixo^s,rho_)+block%equi_vars(ixo^s,equi_rho0_,0)
    else
      wnew(ixo^s,rho_)=w(ixo^s,rho_)
    endif
    wnew(ixo^s,mom(:))=w(ixo^s,mom(:))
 
    if (b0field) then
      ! add background magnetic field B0 to B
      wnew(ixo^s,mag(1:ndir))=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,0)
    else
      wnew(ixo^s,mag(1:ndir))=w(ixo^s,mag(1:ndir))
    end if
 
    if(rmhd_energy) then
      wnew(ixo^s,e_)=w(ixo^s,e_)
      if(has_equi_pe0) then
        wnew(ixo^s,e_)=wnew(ixo^s,e_)+block%equi_vars(ixo^s,equi_pe0_,0)*inv_gamma_1
      end if
      if(b0field .and. total_energy) then
        wnew(ixo^s,e_)=wnew(ixo^s,e_)+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)
      end if
    end if
  end function convert_vars_splitting
 
  subroutine rmhd_check_params
    use mod_global_parameters
    use mod_usr_methods
    use mod_convert, only: add_convert_method
 
    ! after user parameter setting
    gamma_1=rmhd_gamma-1.d0
    if (.not. rmhd_energy) then
       if (rmhd_gamma <= 0.0d0) call mpistop ("Error: rmhd_gamma <= 0")
       if (rmhd_adiab < 0.0d0) call mpistop ("Error: rmhd_adiab < 0")
       small_pressure = rmhd_adiab*small_density**rmhd_gamma
    else
       if (rmhd_gamma <= 0.0d0 .or. rmhd_gamma == 1.0d0) &
            call mpistop ("Error: rmhd_gamma <= 0 or rmhd_gamma == 1")
       inv_gamma_1=1.d0/gamma_1
       small_e = small_pressure * inv_gamma_1
    end if
 
    small_r_e = small_pressure/(rmhd_gamma - 1.0d0)
 
    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(rmhd_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
      endif
    endif
 
    if (use_multigrid) call rmhd_set_mg_bounds()
  end subroutine rmhd_check_params
 
  !> Set the boundaries for the diffusion of E
  subroutine rmhd_set_mg_bounds
    use mod_global_parameters
    use mod_multigrid_coupling
    use mod_usr_methods
    integer :: ib
 
    ! Set boundary conditions for the multigrid solver
    do ib = 1, 2*ndim
       select case (typeboundary(r_e, ib))
       case (bc_symm)
       ! d/dx u = 0
           mg%bc(ib, mg_iphi)%bc_type = mg_bc_neumann
           mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
       case (bc_asymm)
           ! u = 0
           mg%bc(ib, mg_iphi)%bc_type = mg_bc_dirichlet
           mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
       case (bc_cont)
           ! d/dx u = 0
           ! mg%bc(iB, mg_iphi)%bc_type = mg_bc_continuous
           mg%bc(ib, mg_iphi)%bc_type = mg_bc_neumann
           mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
       case (bc_periodic)
           ! Nothing to do here
       case (bc_noinflow)
           call usr_special_mg_bc(ib)
       case (bc_special)
           call usr_special_mg_bc(ib)
       case default
           call mpistop("divE_multigrid warning: unknown b.c. ")
       end select
    end do
  end subroutine rmhd_set_mg_bounds
 
  subroutine rmhd_physical_units()
    use mod_global_parameters
    double precision :: mp,kb,miu0,c_lightspeed
    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 ! G^2 cm^2 dyne^-1
      c_lightspeed=const_c
    end if
    if(eq_state_units) then
      a=1d0+4d0*he_abundance
      if(rmhd_partial_ionization) then
        b=1d0+h_ion_fr+he_abundance*(he_ion_fr*(he_ion_fr2+1d0)+1d0)
      else
        b=2d0+3d0*he_abundance
      end if
      rr=1d0
    else
      a=1d0
      b=1d0
      rr=(1d0+h_ion_fr+he_abundance*(he_ion_fr*(he_ion_fr2+1d0)+1d0))/(1d0+4d0*he_abundance)
    end if
    if(unit_density/=1.d0 .or. unit_numberdensity/=1.d0) then
      if(unit_density/=1.d0) then
        unit_numberdensity=unit_density/(a*mp)
      else if(unit_numberdensity/=1.d0) then
        unit_density=a*mp*unit_numberdensity
      end if
      if(unit_temperature/=1.d0) then
        unit_pressure=b*unit_numberdensity*kb*unit_temperature
        unit_velocity=sqrt(unit_pressure/unit_density)
        unit_magneticfield=sqrt(miu0*unit_pressure)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_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)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_pressure/=1.d0) then
        unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
        unit_velocity=sqrt(unit_pressure/unit_density)
        unit_magneticfield=sqrt(miu0*unit_pressure)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_velocity/=1.d0) then
        unit_pressure=unit_density*unit_velocity**2
        unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
        unit_magneticfield=sqrt(miu0*unit_pressure)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_time/=1.d0) then
        unit_velocity=unit_length/unit_time
        unit_pressure=unit_density*unit_velocity**2
        unit_magneticfield=sqrt(miu0*unit_pressure)
        unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
      end if
    else if(unit_temperature/=1.d0) then
      ! units of temperature and velocity are dependent
      if(unit_magneticfield/=1.d0) then
        unit_pressure=unit_magneticfield**2/miu0
        unit_numberdensity=unit_pressure/(b*unit_temperature*kb)
        unit_density=a*mp*unit_numberdensity
        unit_velocity=sqrt(unit_pressure/unit_density)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_pressure/=1.d0) then
        unit_magneticfield=sqrt(miu0*unit_pressure)
        unit_numberdensity=unit_pressure/(b*unit_temperature*kb)
        unit_density=a*mp*unit_numberdensity
        unit_velocity=sqrt(unit_pressure/unit_density)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      end if
    else if(unit_magneticfield/=1.d0) then
      ! units of magnetic field and pressure are dependent
      if(unit_velocity/=1.d0) then
        unit_pressure=unit_magneticfield**2/miu0
        unit_numberdensity=unit_pressure/(b*unit_temperature*kb)
        unit_density=a*mp*unit_numberdensity
        unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_time/=0.d0) then
        unit_pressure=unit_magneticfield**2/miu0
        unit_velocity=unit_length/unit_time
        unit_density=unit_pressure/unit_velocity**2
        unit_numberdensity=unit_density/(a*mp)
        unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
      end if
    else if(unit_pressure/=1.d0) then
      if(unit_velocity/=1.d0) then
        unit_magneticfield=sqrt(miu0*unit_pressure)
        unit_density=unit_pressure/unit_velocity**2
        unit_numberdensity=unit_density/(a*mp)
        unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_time/=0.d0) then
        unit_magneticfield=sqrt(miu0*unit_pressure)
        unit_velocity=unit_length/unit_time
        unit_density=unit_pressure/unit_velocity**2
        unit_numberdensity=unit_density/(a*mp)
        unit_temperature=unit_pressure/(b*unit_numberdensity*kb)
      end if
    end if
    ! 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
 
    !> Units for radiative flux and opacity
    unit_radflux = unit_velocity*unit_pressure
    unit_opacity = one/(unit_density*unit_length)
  end subroutine rmhd_physical_units
 
  subroutine rmhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
    use mod_global_parameters
    logical, intent(in) :: primitive
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s,nw)
    logical, intent(inout) :: flag(ixi^s,1:nw)
    double precision :: tmp
    integer :: ix^d
 
    flag=.false.
    {do ix^db=ixomin^db,ixomax^db\}
      if(has_equi_rho0) then
        tmp=w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,0)
      else
        tmp=w(ix^d,rho_)
      end if
      if(tmp<small_density) flag(ix^d,rho_) = .true.
      if(primitive) then
        if(has_equi_pe0) then
          if(w(ix^d,p_)+block%equi_vars(ix^d,equi_pe0_,0)<small_pressure) flag(ix^d,e_) = .true.
        else
          if(w(ix^d,p_)<small_pressure) flag(ix^d,e_) = .true.
        end if
      else
        tmp=w(ix^d,e_)-half*((^c&w(ix^d,m^c_)**2+)/tmp+(^c&w(ix^d,b^c_)**2+))
        if(has_equi_pe0) then
          if(tmp+block%equi_vars(ix^d,equi_pe0_,0)*inv_gamma_1<small_e) flag(ix^d,e_) = .true.
        else
          if(tmp<small_e) flag(ix^d,e_) = .true.
        end if
      end if
      if(w(ix^d,r_e)<small_r_e) flag(ix^d,r_e) = .true.
    {end do\}
  end subroutine rmhd_check_w_origin
 
  !> Transform primitive variables into conservative ones
  subroutine rmhd_to_conserved_origin(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
    integer :: ix^d
 
    {do ix^db=ixomin^db,ixomax^db\}
      ! Calculate total energy from pressure, kinetic and magnetic energy
      w(ix^d,e_)=w(ix^d,p_)*inv_gamma_1&
                 +half*((^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)&
                 +(^c&w(ix^d,b^c_)**2+))
      ! Convert velocity to momentum
      ^c&w(ix^d,m^c_)=w(ix^d,rho_)*w(ix^d,m^c_)\
    {end do\}
  end subroutine rmhd_to_conserved_origin
 
  !> Transform primitive variables into conservative ones
  subroutine rmhd_to_conserved_split_rho(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)
    double precision :: rho
    integer :: ix^d
 
    {do ix^db=ixomin^db,ixomax^db\}
      rho=w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i)
      ! Calculate total energy from pressure, kinetic and magnetic energy
      w(ix^d,e_)=w(ix^d,p_)*inv_gamma_1&
                 +half*((^c&w(ix^d,m^c_)**2+)*rho&
                       +(^c&w(ix^d,b^c_)**2+))
      ! Convert velocity to momentum
      ^c&w(ix^d,m^c_)=rho*w(ix^d,m^c_)\
    {end do\}
  end subroutine rmhd_to_conserved_split_rho
 
  !> Transform conservative variables into primitive ones
  subroutine rmhd_to_primitive_origin(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
    double precision                :: inv_rho
    integer :: ix^d
 
    if (fix_small_values) then
      ! fix small values preventing NaN numbers in the following converting
      call rmhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'rmhd_to_primitive_origin')
    end if
 
    {do ix^db=ixomin^db,ixomax^db\}
      inv_rho = 1.d0/w(ix^d,rho_)
      ! Convert momentum to velocity
      ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
      ! Calculate pressure = (gamma-1) * (e-ek-eb)
      w(ix^d,p_)=gamma_1*(w(ix^d,e_)&
                -half*(w(ix^d,rho_)*(^c&w(ix^d,m^c_)**2+)&
                  +(^c&w(ix^d,b^c_)**2+)))
    {end do\}
  end subroutine rmhd_to_primitive_origin
 
  !> Transform conservative variables into primitive ones
  subroutine rmhd_to_primitive_split_rho(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)
    double precision :: inv_rho
    integer :: ix^d
 
    if (fix_small_values) then
      ! fix small values preventing NaN numbers in the following converting
      call rmhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'rmhd_to_primitive_split_rho')
    end if
 
    {do ix^db=ixomin^db,ixomax^db\}
      inv_rho=1.d0/(w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i))
      ! Convert momentum to velocity
      ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
      ! Calculate pressure = (gamma-1) * (e-ek-eb)
      w(ix^d,p_)=gamma_1*(w(ix^d,e_)&
                  -half*((w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i))*&
                  (^c&w(ix^d,m^c_)**2+)+(^c&w(ix^d,b^c_)**2+)))
    {end do\}
  end subroutine rmhd_to_primitive_split_rho
 
  !> Transform internal energy to total energy
  subroutine rmhd_ei_to_e(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
    integer :: ix^d
 
    if(has_equi_rho0) then
      {do ix^db=ixomin^db,ixomax^db\}
        ! Calculate e = ei + ek + eb
        w(ix^d,e_)=w(ix^d,e_)&
                  +half*((^c&w(ix^d,m^c_)**2+)/&
         (w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,0))&
                  +(^c&w(ix^d,b^c_)**2+))
      {end do\}
    else
      {do ix^db=ixomin^db,ixomax^db\}
        ! Calculate e = ei + ek + eb
        w(ix^d,e_)=w(ix^d,e_)&
                  +half*((^c&w(ix^d,m^c_)**2+)/w(ix^d,rho_)&
                        +(^c&w(ix^d,b^c_)**2+))
      {end do\}
    end if
  end subroutine rmhd_ei_to_e
 
  !> Transform total energy to internal energy
  subroutine rmhd_e_to_ei(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
    integer :: ix^d
 
    if(has_equi_rho0) then
      {do ix^db=ixomin^db,ixomax^db\}
        ! Calculate ei = e - ek - eb
        w(ix^d,e_)=w(ix^d,e_)&
                  -half*((^c&w(ix^d,m^c_)**2+)/&
         (w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,0))&
                        +(^c&w(ix^d,b^c_)**2+))
      {end do\}
    else
      {do ix^db=ixomin^db,ixomax^db\}
        ! Calculate ei = e - ek - eb
        w(ix^d,e_)=w(ix^d,e_)&
                  -half*((^c&w(ix^d,m^c_)**2+)/w(ix^d,rho_)&
                        +(^c&w(ix^d,b^c_)**2+))
      {end do\}
    end if
 
    if(fix_small_values) then
      call rmhd_handle_small_ei(w,x,ixi^l,ixi^l,e_,'rmhd_e_to_ei')
    end if
  end subroutine rmhd_e_to_ei
 
  subroutine rmhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
    use mod_global_parameters
    use mod_small_values
    logical, intent(in)             :: primitive
    integer, intent(in)             :: ixi^l,ixo^l
    double precision, intent(inout) :: w(ixi^s,1:nw)
    double precision, intent(in)    :: x(ixi^s,1:ndim)
    character(len=*), intent(in)    :: subname
    double precision :: rho
    integer :: idir, ix^d
    logical :: flag(ixi^s,1:nw)
 
    call phys_check_w(primitive, ixi^l, ixi^l, w, flag)
    if(any(flag)) then
      select case (small_values_method)
      case ("replace")
       {do ix^db=ixomin^db,ixomax^db\}
          if(has_equi_rho0) then
            rho=w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,0)
            if(flag(ix^d,rho_)) w(ix^d,rho_)=small_density-block%equi_vars(ix^d,equi_rho0_,0)
          else
            rho=w(ix^d,rho_)
            if(flag(ix^d,rho_)) w(ix^d,rho_)=small_density
          end if
         {
          if(small_values_fix_iw(m^c_)) then
            if(flag({ix^d},rho_)) w({ix^d},m^c_)=0.0d0
          end if
          \}
          if(primitive) then
            if(has_equi_pe0) then
              if(flag(ix^d,e_)) w(ix^d,p_)=small_pressure-block%equi_vars(ix^d,equi_pe0_,0)
            else
              if(flag(ix^d,e_)) w(ix^d,p_)=small_pressure
            end if
          else
            if(has_equi_pe0) then
              if(flag(ix^d,e_)) &
                w(ix^d,e_)=small_e+half*((^c&w(ix^d,m^c_)**2+)/rho+(^c&w(ix^d,b^c_)**2+))&
                -block%equi_vars(ix^d,equi_pe0_,0)*inv_gamma_1
            else
              if(flag(ix^d,e_)) &
                w(ix^d,e_)=small_e+half*((^c&w(ix^d,m^c_)**2+)/rho+(^c&w(ix^d,b^c_)**2+))
            end if
          end if
          if(flag(ix^d,r_e)) w(ix^d,r_e)=small_r_e
       {end do\}
      case ("average")
        ! do averaging of density
        call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
        if(primitive)then
          call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
        else
          ! do averaging of internal energy
          {do ix^db=iximin^db,iximax^db\}
            if(has_equi_rho0) then
              rho=w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,0)
            else
              rho=w(ix^d,rho_)
            end if
            w(ix^d,e_)=w(ix^d,e_)&
                -half*((^c&w(ix^d,m^c_)**2+)/rho+(^c&w(ix^d,b^c_)**2+))
          {end do\}
          call small_values_average(ixi^l, ixo^l, w, x, flag, e_)
          ! convert back
          {do ix^db=iximin^db,iximax^db\}
            if(has_equi_rho0) then
              rho=w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,0)
            else
              rho=w(ix^d,rho_)
            end if
            w(ix^d,e_)=w(ix^d,e_)&
                +half*((^c&w(ix^d,m^c_)**2+)/rho+(^c&w(ix^d,b^c_)**2+))
          {end do\}
        end if
        call small_values_average(ixi^l, ixo^l, w, x, flag, r_e)
      case default
        if(.not.primitive) then
          !convert w to primitive
          ! do averaging of internal energy
         {do ix^db=iximin^db,iximax^db\}
            if(has_equi_rho0) then
              rho=w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,0)
            else
              rho=w(ix^d,rho_)
            end if
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)/rho\
            w(ix^d,p_)=gamma_1*(w(ix^d,e_)&
                -half*((^c&w(ix^d,m^c_)**2+)/rho+(^c&w(ix^d,b^c_)**2+)))
         {end do\}
        end if
        call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
      end select
    end if
  end subroutine rmhd_handle_small_values_origin
 
  !> Calculate v vector
  subroutine rmhd_get_v(w,x,ixI^L,ixO^L,v)
    use mod_global_parameters
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s,nw), x(ixi^s,1:ndim)
    double precision, intent(out) :: v(ixi^s,ndir)
    double precision :: rho(ixi^s)
    integer :: idir
 
    call rmhd_get_rho(w,x,ixi^l,ixo^l,rho)
    rho(ixo^s)=1.d0/rho(ixo^s)
    ! Convert momentum to velocity
    do idir = 1, ndir
       v(ixo^s, idir) = w(ixo^s, mom(idir))*rho(ixo^s)
    end do
  end subroutine rmhd_get_v
 
  !> Calculate cmax_idim=csound+abs(v_idim) within ixO^L
  subroutine rmhd_get_cmax_origin(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 :: rho, inv_rho, cfast2, avmincs2, b2, kmax
    integer :: ix^d
 
    if(b0field) then
     {do ix^db=ixomin^db,ixomax^db \}
        if(has_equi_rho0) then
          rho=(w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i))
        else
          rho=w(ix^d,rho_)
        end if
        inv_rho=1.d0/rho
        ! sound speed**2 
        cmax(ix^d)=rmhd_gamma*w(ix^d,p_)*inv_rho
        ! store |B|^2 in v
        b2=(^c&(w(ix^d,b^c_)+block%B0(ix^d,^c,b0i))**2+)
        cfast2=b2*inv_rho+cmax(ix^d)
        avmincs2=cfast2**2-4.0d0*cmax(ix^d)*(w(ix^d,mag(idim))+block%B0(ix^d,idim,b0i))**2*inv_rho
        if(avmincs2<zero) avmincs2=zero
        cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
        cmax(ix^d)=abs(w(ix^d,mom(idim)))+cmax(ix^d)
     {end do\}
    else
     {do ix^db=ixomin^db,ixomax^db \}
        if(has_equi_rho0) then
          rho=(w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i))
        else
          rho=w(ix^d,rho_)
        end if
        inv_rho=1.d0/rho
        ! sound speed**2 
        cmax(ix^d)=rmhd_gamma*w(ix^d,p_)*inv_rho
        ! store |B|^2 in v
        b2=(^c&w(ix^d,b^c_)**2+)
        cfast2=b2*inv_rho+cmax(ix^d)
        avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
        if(avmincs2<zero) avmincs2=zero
        cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
        cmax(ix^d)=abs(w(ix^d,mom(idim)))+cmax(ix^d)
     {end do\}
    end if
  end subroutine rmhd_get_cmax_origin
 
  subroutine rmhd_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 rmhd_get_a2max
 
  !> get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
  subroutine rmhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in) :: ixi^l,ixo^l
    double precision, intent(in) :: x(ixi^s,1:ndim)
    double precision, intent(out) :: tco_local,tmax_local
    ! in primitive form
    double precision, intent(inout) :: w(ixi^s,1:nw)
    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 :: ltrc,ltrp,altr(ixi^s)
    integer :: idims,jxo^l,hxo^l,ixa^d,ixb^d,ix^d
    integer :: jxp^l,hxp^l,ixp^l,ixq^l
    logical :: lrlt(ixi^s)
 
    if(rmhd_partial_ionization) then
      call rmhd_get_temperature_from_te(w,x,ixi^l,ixi^l,te)
    else
      call rmhd_get_rfactor(w,x,ixi^l,ixi^l,te)
      te(ixi^s)=w(ixi^s,p_)/(te(ixi^s)*w(ixi^s,rho_))
    end if
    tco_local=zero
    tmax_local=maxval(te(ixo^s))
 
    {^ifoned
    select case(rmhd_trac_type)
    case(0)
      !> test case, fixed cutoff temperature
      block%wextra(ixi^s,tcoff_)=2.5d5/unit_temperature
    case(1)
      hxo^l=ixo^l-1;
      jxo^l=ixo^l+1;
      lts(ixo^s)=0.5d0*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=4.d0
      ixp^l=ixo^l^ladd1;
      hxo^l=ixo^l-1;
      jxo^l=ixo^l+1;
      hxp^l=ixp^l-1;
      jxp^l=ixp^l+1;
      lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
      lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
      lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
      block%wextra(ixo^s,tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
    case default
      call mpistop("rmhd_trac_type not allowed for 1D simulation")
    end select
    }
    {^nooned
    select case(rmhd_trac_type)
    case(0)
      !> test case, fixed cutoff temperature
      block%wextra(ixi^s,tcoff_)=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(rmhd_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^ladd2;
      ! temperature gradient at cell centers
      do idims=1,ndim
        ixq^l=ixp^l;
        hxp^l=ixp^l;
        jxp^l=ixp^l;
        select case(idims)
          {case(^d)
           ixqmin^d=ixqmin^d+1
           ixqmax^d=ixqmax^d-1
           hxpmax^d=ixpmin^d
           jxpmin^d=ixpmax^d
          \}
        end select
        call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
        call gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
        call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
      end do
      ! B vector
     {do ix^db=ixpmin^db,ixpmax^db\}
        if(b0field) then
          ^c&bunitvec(ix^d,^c)=w(ix^d,iw_mag(^c))+block%B0(ix^d,^c,0)\
        else
          ^c&bunitvec(ix^d,^c)=w(ix^d,iw_mag(^c))\
        end if
        tmp1(ix^d)=1.d0/(dsqrt(^c&bunitvec(ix^d,^c)**2+)+smalldouble)
        ! b unit vector: magnetic field direction vector
        ^d&bunitvec({ix^d},^d)=bunitvec({ix^d},^d)*tmp1({ix^d})\
        ! temperature length scale inversed
        lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec({ix^d},^d)+)/te(ix^d)
        ! fraction of cells size to temperature length scale
        if(slab_uniform) then
          lts(ix^d)=min(^d&dxlevel(^d))*lts(ix^d)
        else
          lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
        end if
        lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
     {end do\}
      ! need one ghost layer for thermal conductivity
      ixp^l=ixo^l^ladd1;
      do idims=1,ndim
        hxo^l=ixp^l-kr(idims,^d);
        jxo^l=ixp^l+kr(idims,^d);
        if(idims==1) then
          altr(ixp^s)=0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
        else
          altr(ixp^s)=altr(ixp^s)+0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
        end if
      end do
      block%wextra(ixp^s,tcoff_)=te(ixp^s)*altr(ixp^s)**0.4d0
    case(3,5)
      !> do nothing here
    case default
      call mpistop("unknown rmhd_trac_type")
    end select
    }
  end subroutine rmhd_get_tcutoff
 
  !> get H speed for H-correction to fix the carbuncle problem at grid-aligned shock front
  subroutine rmhd_get_h_speed(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)
    double precision, allocatable :: tmp(:^d&)
    integer :: jxc^l, ixc^l, ixa^l, id, ix^d
 
    hspeed=0.d0
    ixa^l=ixo^l^ladd1;
    allocate(tmp(ixa^s))
    do id=1,ndim
      call rmhd_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(wprim(jxc^s,mom(idim))+csound(jxc^s,idim)-wprim(ixc^s,mom(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(id))+csound(ixa^s,id)-wprim(ixc^s,mom(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(id))+csound(ixc^s,id)-wprim(ixa^s,mom(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(id))+csound(ixa^s,id)-wprim(jxc^s,mom(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(id))+csound(jxc^s,id)-wprim(ixa^s,mom(id))+csound(ixa^s,id)))
    end do
    deallocate(tmp)
 
  end subroutine rmhd_get_h_speed
 
  !> Estimating bounds for the minimum and maximum signal velocities without split
  subroutine rmhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l, idim
    double precision, intent(in)    :: wlc(ixi^s, nw), wrc(ixi^s, nw)
    double precision, intent(in)    :: wlp(ixi^s, nw), wrp(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s,1:ndim)
    double precision, intent(inout) :: cmax(ixi^s,1:number_species)
    double precision, intent(inout), optional :: cmin(ixi^s,1:number_species)
    double precision, intent(in)    :: hspeed(ixi^s,1:number_species)
 
    double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
    double precision :: umean, dmean, 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 rmhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
      call rmhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
      if(present(cmin)) then
       {do ix^db=ixomin^db,ixomax^db\}
          tmp1=sqrt(wlp(ix^d,rho_))
          tmp2=sqrt(wrp(ix^d,rho_))
          tmp3=1.d0/(tmp1+tmp2)
          umean=(wlp(ix^d,mom(idim))*tmp1+wrp(ix^d,mom(idim))*tmp2)*tmp3
          dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
           half*tmp1*tmp2*tmp3**2*(wrp(ix^d,mom(idim))-wlp(ix^d,mom(idim)))**2)
          cmin(ix^d,1)=umean-dmean
          cmax(ix^d,1)=umean+dmean
       {end do\}
        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
       {do ix^db=ixomin^db,ixomax^db\}
          tmp1=sqrt(wlp(ix^d,rho_))
          tmp2=sqrt(wrp(ix^d,rho_))
          tmp3=1.d0/(tmp1+tmp2)
          umean=(wlp(ix^d,mom(idim))*tmp1+wrp(ix^d,mom(idim))*tmp2)*tmp3
          dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
           half*tmp1*tmp2*tmp3**2*(wrp(ix^d,mom(idim))-wlp(ix^d,mom(idim)))**2)
          cmax(ix^d,1)=abs(umean)+dmean
       {end do\}
      end if
    case (2)
      wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
      call rmhd_get_csound_prim(wmean,x,ixi^l,ixo^l,idim,csoundr)
      if(present(cmin)) then
       {do ix^db=ixomin^db,ixomax^db\}
          cmax(ix^d,1)=max(wmean(ix^d,mom(idim))+csoundr(ix^d),zero)
          cmin(ix^d,1)=min(wmean(ix^d,mom(idim))-csoundr(ix^d),zero)
       {end do\}
        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(wmean(ixo^s,mom(idim)))+csoundr(ixo^s)
      end if
    case (3)
      ! Miyoshi 2005 JCP 208, 315 equation (67)
      call rmhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
      call rmhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
      if(present(cmin)) then
       {do ix^db=ixomin^db,ixomax^db\}
          csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
          cmin(ix^d,1)=min(wlp(ix^d,mom(idim)),wrp(ix^d,mom(idim)))-csoundl(ix^d)
          cmax(ix^d,1)=max(wlp(ix^d,mom(idim)),wrp(ix^d,mom(idim)))+csoundl(ix^d)
       {end do\}
        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
       {do ix^db=ixomin^db,ixomax^db\}
          csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
          cmax(ix^d,1)=max(wlp(ix^d,mom(idim)),wrp(ix^d,mom(idim)))+csoundl(ix^d)
       {end do\}
      end if
    end select
  end subroutine rmhd_get_cbounds
 
  !> Estimating bounds for the minimum and maximum signal velocities with rho split
  subroutine rmhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l, idim
    double precision, intent(in)    :: wlc(ixi^s, nw), wrc(ixi^s, nw)
    double precision, intent(in)    :: wlp(ixi^s, nw), wrp(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s,1:ndim)
    double precision, intent(inout) :: cmax(ixi^s,1:number_species)
    double precision, intent(inout), optional :: cmin(ixi^s,1:number_species)
    double precision, intent(in)    :: hspeed(ixi^s,1:number_species)
    double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
    double precision :: umean, dmean, 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 rmhd_get_csound_prim_split(wlp,x,ixi^l,ixo^l,idim,csoundl)
      call rmhd_get_csound_prim_split(wrp,x,ixi^l,ixo^l,idim,csoundr)
      if(present(cmin)) then
       {do ix^db=ixomin^db,ixomax^db\}
          tmp1=sqrt(wlp(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i))
          tmp2=sqrt(wrp(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i))
          tmp3=1.d0/(tmp1+tmp2)
          umean=(wlp(ix^d,mom(idim))*tmp1+wrp(ix^d,mom(idim))*tmp2)*tmp3
          dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
           half*tmp1*tmp2*tmp3**2*(wrp(ix^d,mom(idim))-wlp(ix^d,mom(idim)))**2)
          cmin(ix^d,1)=umean-dmean
          cmax(ix^d,1)=umean+dmean
       {end do\}
        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
       {do ix^db=ixomin^db,ixomax^db\}
          tmp1=sqrt(wlp(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i))
          tmp2=sqrt(wrp(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i))
          tmp3=1.d0/(tmp1+tmp2)
          umean=(wlp(ix^d,mom(idim))*tmp1+wrp(ix^d,mom(idim))*tmp2)*tmp3
          dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
           half*tmp1*tmp2*tmp3**2*(wrp(ix^d,mom(idim))-wlp(ix^d,mom(idim)))**2)
          cmax(ix^d,1)=abs(umean)+dmean
       {end do\}
      end if
    case (2)
      wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
      call rmhd_get_csound_prim_split(wmean,x,ixi^l,ixo^l,idim,csoundr)
      if(present(cmin)) then
       {do ix^db=ixomin^db,ixomax^db\}
          cmax(ix^d,1)=max(wmean(ix^d,mom(idim))+csoundr(ix^d),zero)
          cmin(ix^d,1)=min(wmean(ix^d,mom(idim))-csoundr(ix^d),zero)
       {end do\}
        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(wmean(ixo^s,mom(idim)))+csoundr(ixo^s)
      end if
    case (3)
      ! Miyoshi 2005 JCP 208, 315 equation (67)
      call rmhd_get_csound_prim_split(wlp,x,ixi^l,ixo^l,idim,csoundl)
      call rmhd_get_csound_prim_split(wrp,x,ixi^l,ixo^l,idim,csoundr)
      if(present(cmin)) then
       {do ix^db=ixomin^db,ixomax^db\}
          csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
          cmin(ix^d,1)=min(wlp(ix^d,mom(idim)),wrp(ix^d,mom(idim)))-csoundl(ix^d)
          cmax(ix^d,1)=max(wlp(ix^d,mom(idim)),wrp(ix^d,mom(idim)))+csoundl(ix^d)
       {end do\}
        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
       {do ix^db=ixomin^db,ixomax^db\}
          csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
          cmax(ix^d,1)=max(wlp(ix^d,mom(idim)),wrp(ix^d,mom(idim)))+csoundl(ix^d)
       {end do\}
      end if
    end select
  end subroutine rmhd_get_cbounds_split_rho
 
  !> prepare velocities for ct methods
  subroutine rmhd_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(idim))+wrp(ixo^s,mom(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(idimn))
      vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,mom(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(idime))
      vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,mom(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 rmhd_get_ct_velocity
 
  !> Calculate fast magnetosonic wave speed
  subroutine rmhd_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(ixo^s)
    double precision :: inv_rho, cfast2, avmincs2, b2, kmax
    double precision :: prad_tensor(ixo^s, 1:ndim, 1:ndim)
    double precision :: prad_max(ixo^s)
    integer :: ix^d
 
    call rmhd_get_pradiation(w, x, ixi^l, ixo^l, prad_tensor, nghostcells-1)
    !> filter cmax
    if(radio_acoustic_filter) then
      call rmhd_radio_acoustic_filter(x, ixi^l, ixo^l, prad_max)
    endif
    ! store |B|^2 in v
    if(b0field) then
      {do ix^db=ixomin^db,ixomax^db \}
        inv_rho=1.d0/w(ix^d,rho_)
        prad_max(ix^d) = maxval(prad_tensor(ix^d,:,:))
        if(rmhd_energy) then
          csound(ix^d)=max(rmhd_gamma,4.d0/3.d0)*(w(ix^d,p_)+prad_max(ix^d))*inv_rho
        else
          csound(ix^d)=rmhd_gamma*rmhd_adiab*w(ix^d,rho_)**gamma_1
        end if
        b2=(^c&(w(ix^d,b^c_)+block%B0(ix^d,^c,b0i))**2+)
        cfast2=b2*inv_rho+csound(ix^d)
        avmincs2=cfast2**2-4.0d0*csound(ix^d)*(w(ix^d,mag(idim))+&
         block%B0(ix^d,idim,b0i))**2*inv_rho
        if(avmincs2<zero) avmincs2=zero
        csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
      {end do\}
    else
      {do ix^db=ixomin^db,ixomax^db \}
        inv_rho=1.d0/w(ix^d,rho_)
        prad_max(ix^d)=maxval(prad_tensor(ix^d,:,:))
        if(rmhd_energy) then
          csound(ix^d)=max(rmhd_gamma,4.d0/3.d0)*(w(ix^d,p_)+prad_max(ix^d))*inv_rho
        else
          csound(ix^d)=rmhd_gamma*rmhd_adiab*w(ix^d,rho_)**gamma_1
        end if
        b2=(^c&w(ix^d,b^c_)**2+)
        cfast2=b2*inv_rho+csound(ix^d)
        avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
        if(avmincs2<zero) avmincs2=zero
        csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
      {end do\}
    end if
  end subroutine rmhd_get_csound_prim
 
  !> Calculate fast magnetosonic wave speed
  subroutine rmhd_get_csound_prim_split(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(ixo^s)
    double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
    double precision :: prad_tensor(ixo^s, 1:ndim, 1:ndim)
    double precision :: prad_max(ixo^s)
    integer :: ix^d
 
    call rmhd_get_pradiation(w, x, ixi^l, ixo^l, prad_tensor, nghostcells-1)
    !> filter cmax
    if (radio_acoustic_filter) then
      call rmhd_radio_acoustic_filter(x, ixi^l, ixo^l, prad_max)
    endif
 
    ! store |B|^2 in v
    if(b0field) then
      {do ix^db=ixomin^db,ixomax^db \}
        rho=(w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i))
        inv_rho=1.d0/rho
        prad_max(ix^d) = maxval(prad_tensor(ix^d,:,:))
        if(has_equi_pe0) then
          csound(ix^d)=max(rmhd_gamma,4.d0/3.d0)*(w(ix^d,p_)+block%equi_vars(ix^d,equi_pe0_,b0i)+prad_max(ix^d))*inv_rho
        end if
        b2=(^c&(w(ix^d,b^c_)+block%B0(ix^d,^c,b0i))**2+)
        cfast2=b2*inv_rho+csound(ix^d)
        avmincs2=cfast2**2-4.0d0*csound(ix^d)*(w(ix^d,mag(idim))+&
         block%B0(ix^d,idim,b0i))**2*inv_rho
        if(avmincs2<zero) avmincs2=zero
        csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
      {end do\}
    else
      {do ix^db=ixomin^db,ixomax^db \}
        rho=(w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i))
        inv_rho=1.d0/rho
        prad_max(ix^d) = maxval(prad_tensor(ix^d,:,:))
        if(has_equi_pe0) then
          csound(ix^d)=max(rmhd_gamma,4.d0/3.d0)*(w(ix^d,p_)+block%equi_vars(ix^d,equi_pe0_,b0i)+prad_max(ix^d))*inv_rho
        end if
        b2=(^c&w(ix^d,b^c_)**2+)
        cfast2=b2*inv_rho+csound(ix^d)
        avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
        if(avmincs2<zero) avmincs2=zero
        csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
      {end do\}
    end if
  end subroutine rmhd_get_csound_prim_split
 
  !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L
  subroutine rmhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
    use mod_global_parameters
    use mod_small_values, only: trace_small_values
 
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s,nw)
    double precision, intent(in) :: x(ixi^s,1:ndim)
    double precision, intent(out):: pth(ixi^s)
 
    integer :: iw, ix^d
 
   {do ix^db=ixomin^db,ixomax^db\}
      if(has_equi_rho0) then
        pth(ix^d)=gamma_1*(w(ix^d,e_)-half*((^c&w(ix^d,m^c_)**2+)/(w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,0))&
             +(^c&w(ix^d,b^c_)**2+)))+block%equi_vars(ix^d,equi_pe0_,0)
      else
        pth(ix^d)=gamma_1*(w(ix^d,e_)-half*((^c&w(ix^d,m^c_)**2+)/w(ix^d,rho_)&
             +(^c&w(ix^d,b^c_)**2+)))
      end if
      if(fix_small_values.and.pth(ix^d)<small_pressure) pth(ix^d)=small_pressure
   {end do\}
 
    if(check_small_values.and..not.fix_small_values) then
     {do ix^db=ixomin^db,ixomax^db\}
        if(pth(ix^d)<small_pressure) then
          write(*,*) "Error: small value of gas pressure",pth(ix^d),&
               " encountered when call rmhd_get_pthermal"
          write(*,*) "Iteration: ", it, " Time: ", global_time
          write(*,*) "Location: ", x(ix^d,:)
          write(*,*) "Cell number: ", ix^d
          do iw=1,nw
            write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
          end do
          ! use erroneous arithmetic operation to crash the run
          if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
          write(*,*) "Saving status at the previous time step"
          crash=.true.
        end if
     {end do\}
    end if
  end subroutine rmhd_get_pthermal_origin
 
  !> copy temperature from stored Te variable
  subroutine rmhd_get_temperature_from_te(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out):: res(ixi^s)
    res(ixo^s) = w(ixo^s, te_)
  end subroutine rmhd_get_temperature_from_te
 
  !> Calculate temperature=p/rho when in e_ the internal energy is stored
  subroutine rmhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out):: res(ixi^s)
    double precision :: r(ixi^s)
 
    call rmhd_get_rfactor(w,x,ixi^l,ixo^l,r)
    res(ixo^s) = gamma_1 * w(ixo^s, e_)/(w(ixo^s,rho_)*r(ixo^s))
  end subroutine rmhd_get_temperature_from_eint
 
  !> Calculate temperature=p/rho when in e_ the total energy is stored
  subroutine rmhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
    use mod_global_parameters
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out):: res(ixi^s)
    double precision :: r(ixi^s)
 
    call rmhd_get_rfactor(w,x,ixi^l,ixo^l,r)
    call rmhd_get_pthermal(w,x,ixi^l,ixo^l,res)
    res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,rho_))
  end subroutine rmhd_get_temperature_from_etot
 
  subroutine rmhd_get_temperature_from_etot_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)
    double precision :: r(ixi^s)
 
    call rmhd_get_rfactor(w,x,ixi^l,ixo^l,r)
    call rmhd_get_pthermal(w,x,ixi^l,ixo^l,res)
    res(ixo^s)=res(ixo^s)/(r(ixo^s)*(w(ixo^s,rho_)+block%equi_vars(ixo^s,equi_rho0_,b0i)))
  end subroutine rmhd_get_temperature_from_etot_with_equi
 
  subroutine rmhd_get_temperature_from_eint_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)
    double precision :: r(ixi^s)
 
    call rmhd_get_rfactor(w,x,ixi^l,ixo^l,r)
    res(ixo^s) = (gamma_1 * w(ixo^s, e_) + block%equi_vars(ixo^s,equi_pe0_,b0i)) /&
                ((w(ixo^s,rho_) +block%equi_vars(ixo^s,equi_rho0_,b0i))*r(ixo^s))
  end subroutine rmhd_get_temperature_from_eint_with_equi
 
  subroutine rmhd_get_temperature_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)
    double precision :: r(ixi^s)
 
    call rmhd_get_rfactor(w,x,ixi^l,ixo^l,r)
    res(ixo^s)= block%equi_vars(ixo^s,equi_pe0_,b0i)/(block%equi_vars(ixo^s,equi_rho0_,b0i)*r(ixo^s))
  end subroutine rmhd_get_temperature_equi
 
  subroutine rmhd_get_rho_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_rho0_,b0i)
  end subroutine rmhd_get_rho_equi
 
  subroutine rmhd_get_pe_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_pe0_,b0i)
  end subroutine rmhd_get_pe_equi
 
  !> Calculate total pressure within ixO^L including magnetic pressure
  subroutine rmhd_get_p_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 rmhd_get_pthermal(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 rmhd_get_p_total
 
  !> Calculate radiation pressure within ixO^L
  subroutine rmhd_get_pradiation(w, x, ixI^L, ixO^L, prad, nth)
    use mod_global_parameters
    use mod_fld
    use mod_afld
    integer, intent(in)          :: ixi^l, ixo^l, nth
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out):: prad(ixo^s, 1:ndim, 1:ndim)
 
    select case (rmhd_radiation_formalism)
    case('fld')
      call fld_get_radpress(w, x, ixi^l, ixo^l, prad, nth)
    case('afld')
      call afld_get_radpress(w, x, ixi^l, ixo^l, prad, nth)
    case default
      call mpistop('Radiation formalism unknown')
    end select
  end subroutine rmhd_get_pradiation
 
  !> Calculates the sum of the gas pressure and the max Prad tensor element
  subroutine rmhd_get_pthermal_plus_pradiation(w, x, ixI^L, ixO^L, pth_plus_prad)
    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              :: pth(ixi^s)
    double precision              :: prad_tensor(ixo^s, 1:ndim, 1:ndim)
    double precision              :: prad_max(ixo^s)
    double precision, intent(out) :: pth_plus_prad(ixi^s)
    integer :: ix^d
 
    call rmhd_get_pthermal(w, x, ixi^l, ixo^l, pth)
    call rmhd_get_pradiation(w, x, ixi^l, ixo^l, prad_tensor, nghostcells)
    {do ix^d = ixomin^d,ixomax^d\}
      prad_max(ix^d) = maxval(prad_tensor(ix^d,:,:))
    {enddo\}
    !> filter cmax
    if (radio_acoustic_filter) then
      call rmhd_radio_acoustic_filter(x, ixi^l, ixo^l, prad_max)
    endif
    pth_plus_prad(ixo^s) = pth(ixo^s) + prad_max(ixo^s)
  end subroutine rmhd_get_pthermal_plus_pradiation
 
  !> Filter peaks in cmax due to radiation energy density, used for debugging
  subroutine rmhd_radio_acoustic_filter(x, ixI^L, ixO^L, prad_max)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
    double precision, intent(inout) :: prad_max(ixo^s)
    double precision :: tmp_prad(ixi^s)
    integer :: ix^d, filter, idim
 
    if (size_ra_filter .lt. 1) call mpistop("ra filter of size < 1 makes no sense")
    if (size_ra_filter .gt. nghostcells) call mpistop("ra filter of size < nghostcells makes no sense")
 
    tmp_prad(ixi^s) = zero
    tmp_prad(ixo^s) = prad_max(ixo^s)
    do filter = 1,size_ra_filter
      do idim = 1,ndim
        ! {do ix^D = ixOmin^D+filter,ixOmax^D-filter\}
        {do ix^d = ixomin^d,ixomax^d\}
            prad_max(ix^d) = min(tmp_prad(ix^d),tmp_prad(ix^d+filter*kr(idim,^d)))
            prad_max(ix^d) = min(tmp_prad(ix^d),tmp_prad(ix^d-filter*kr(idim,^d)))
        {enddo\}
      enddo
    enddo
  end subroutine rmhd_radio_acoustic_filter
 
  !> Calculates gas temperature
  subroutine rmhd_get_tgas(w, x, ixI^L, ixO^L, tgas)
    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             :: pth(ixi^s)
    double precision, intent(out):: tgas(ixi^s)
 
    call rmhd_get_pthermal(w, x, ixi^l, ixo^l, pth)
        tgas(ixi^s) = pth(ixi^s)/w(ixi^s,rho_)
  end subroutine rmhd_get_tgas
 
  !> Calculates radiation temperature
  subroutine rmhd_get_trad(w, x, ixI^L, ixO^L, trad)
    use mod_global_parameters
    use mod_constants
 
    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):: trad(ixi^s)
 
    trad(ixi^s) = (w(ixi^s,r_e)*unit_pressure&
    /const_rad_a)**(1.d0/4.d0)/unit_temperature
  end subroutine rmhd_get_trad
 
  !> Calculate fluxes within ixO^L without any splitting
  subroutine rmhd_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             :: vhall(ixi^s,1:ndir)
    double precision             :: ptotal
    integer                      :: iw, ix^d
 
    {do ix^db=ixomin^db,ixomax^db\}
       ! Get flux of density
       f(ix^d,rho_)=w(ix^d,mom(idim))*w(ix^d,rho_)
       ! f_i[m_k]=v_i*m_k-b_k*b_i
       ^c&f(ix^d,m^c_)=wc(ix^d,mom(idim))*w(ix^d,m^c_)-w(ix^d,mag(idim))*w(ix^d,b^c_)\
       ptotal=w(ix^d,p_)+half*(^c&w(ix^d,b^c_)**2+)
       ! normal one includes total pressure
       f(ix^d,mom(idim))=f(ix^d,mom(idim))+ptotal
       ! Get flux of total energy
       ! f_i[e]=v_i*e+v_i*ptotal-b_i*(b_k*v_k)
       f(ix^d,e_)=w(ix^d,mom(idim))*(wc(ix^d,e_)+ptotal)&
          -w(ix^d,mag(idim))*(^c&w(ix^d,b^c_)*w(ix^d,m^c_)+)
       ! f_i[b_k]=v_i*b_k-v_k*b_i
       ^c&f(ix^d,b^c_)=w(ix^d,mom(idim))*w(ix^d,b^c_)-w(ix^d,mag(idim))*w(ix^d,m^c_)\
    {end do\}
    if(rmhd_glm) then
     {do ix^db=ixomin^db,ixomax^db\}
        f(ix^d,mag(idim))=w(ix^d,psi_)
        !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
        f(ix^d,psi_)=cmax_global**2*w(ix^d,mag(idim))
     {end do\}
    end if
    if(rmhd_radiation_advection) then
      {do ix^db=ixomin^db,ixomax^db\}
        f(ix^d,r_e)=w(ix^d,mom(idim))*wc(ix^d,r_e)
      {end do\}
    else
      f(ixo^s,r_e)=zero
    endif
    ! Get flux of tracer
    do iw=1,rmhd_n_tracer
     {do ix^db=ixomin^db,ixomax^db\}
        f(ix^d,tracer(iw))=w(ix^d,mom(idim))*w(ix^d,tracer(iw))
     {end do\}
    end do
    if(rmhd_hyperbolic_thermal_conduction) then
     {do ix^db=ixomin^db,ixomax^db\}
        f(ix^d,e_)=f(ix^d,e_)+w(ix^d,q_)*w(ix^d,mag(idim))/(dsqrt(^d&w({ix^d},b^d_)**2+)+smalldouble)
        f(ix^d,q_)=zero
     {end do\}
    end if
  end subroutine rmhd_get_flux
 
  !> Calculate fluxes within ixO^L with possible splitting
  subroutine rmhd_get_flux_split(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             :: vhall(ixi^s,1:ndir)
    double precision             :: ptotal, btotal(ixo^s,1:ndir)
    integer                      :: iw, ix^d
 
    {do ix^db=ixomin^db,ixomax^db\}
      ! Get flux of density
      if(has_equi_rho0) then
        f(ix^d,rho_)=w(ix^d,mom(idim))*(w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i))
      else
        f(ix^d,rho_)=w(ix^d,mom(idim))*w(ix^d,rho_)
      endif
      if(rmhd_energy) then
        ptotal=w(ix^d,p_)+half*(^c&w(ix^d,b^c_)**2+)
      else
        ptotal=rmhd_adiab*w(ix^d,rho_)**rmhd_gamma+half*(^c&w(ix^d,b^c_)**2+)
        if(has_equi_pe0) then
          ptotal=ptotal-block%equi_vars(ix^d,equi_pe0_,b0i)
        end if
      end if
      if(b0field) then
        ^c&btotal(ix^d,^c)=w(ix^d,b^c_)+block%B0(ix^d,^c,idim)\
        ptotal=ptotal+(^c&w(ix^d,b^c_)*block%B0(ix^d,^c,idim)+)
        ! Get flux of momentum and magnetic field
        ! f_i[m_k]=v_i*m_k-b_k*b_i
        ^c&f(ix^d,m^c_)=wc(ix^d,mom(idim))*w(ix^d,m^c_)-&
                          btotal(ix^d,idim)*w(ix^d,b^c_)-w(ix^d,mag(idim))*block%B0(ix^d,^c,idim)\
        f(ix^d,mom(idim))=f(ix^d,mom(idim))+ptotal
      else
        ^c&btotal(ix^d,^c)=w(ix^d,b^c_)\
        ! Get flux of momentum and magnetic field
        ! f_i[m_k]=v_i*m_k-b_k*b_i
        ^c&f(ix^d,m^c_)=wc(ix^d,mom(idim))*w(ix^d,m^c_)-w(ix^d,mag(idim))*w(ix^d,b^c_)\
        f(ix^d,mom(idim))=f(ix^d,mom(idim))+ptotal
      end if
      ! f_i[b_k]=v_i*b_k-v_k*b_i
      ^c&f(ix^d,b^c_)=w(ix^d,mom(idim))*btotal(ix^d,^c)-btotal(ix^d,idim)*w(ix^d,m^c_)\
      ! Get flux of energy
      ! f_i[e]=v_i*e+v_i*ptotal-b_i*(b_k*v_k)
      if(rmhd_energy) then
        f(ix^d,e_)=w(ix^d,mom(idim))*(wc(ix^d,e_)+ptotal)&
           -btotal(ix^d,idim)*(^c&w(ix^d,b^c_)*w(ix^d,m^c_)+)
      end if
    {end do\}
    if(rmhd_glm) then
      {do ix^db=ixomin^db,ixomax^db\}
        f(ix^d,mag(idim))=w(ix^d,psi_)
        !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
        f(ix^d,psi_) = cmax_global**2*w(ix^d,mag(idim))
      {end do\}
    end if
    if(rmhd_radiation_advection) then
      {do ix^db=ixomin^db,ixomax^db\}
        f(ix^d,r_e)=w(ix^d,mom(idim))*wc(ix^d,r_e)
      {end do\}
    else
      f(ixo^s,r_e)=zero
    endif
    ! Get flux of tracer
    do iw=1,rmhd_n_tracer
     {do ix^db=ixomin^db,ixomax^db\}
        f(ix^d,tracer(iw))=w(ix^d,mom(idim))*w(ix^d,tracer(iw))
     {end do\}
    end do
    if(rmhd_hyperbolic_thermal_conduction) then
     {do ix^db=ixomin^db,ixomax^db\}
        f(ix^d,e_)=f(ix^d,e_)+w(ix^d,q_)*btotal(ix^d,idim)/(dsqrt(^c&btotal(ix^d,^c)**2+)+smalldouble)
        f(ix^d,q_)=zero
     {end do\}
    end if
  end subroutine rmhd_get_flux_split
 
  !> use cell-center flux to get cell-face flux
  !> and get the source term as the divergence of the flux
  subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
    use mod_global_parameters
 
    integer, intent(in) :: ixi^l, ixo^l
    double precision, dimension(:^D&,:), intent(inout) :: ff
    double precision, intent(out) :: src(ixi^s)
 
    double precision :: ffc(ixi^s,1:ndim)
    double precision :: dxinv(ndim)
    integer :: idims, ix^d, ixa^l, ixb^l, ixc^l
 
    ixa^l=ixo^l^ladd1;
    dxinv=1.d0/dxlevel
    ! cell corner flux in ffc
    ffc=0.d0
    ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
    {do ix^db=0,1\}
      ixbmin^d=ixcmin^d+ix^d;
      ixbmax^d=ixcmax^d+ix^d;
      ffc(ixc^s,1:ndim)=ffc(ixc^s,1:ndim)+ff(ixb^s,1:ndim)
    {end do\}
    ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
    ! flux at cell face
    ff(ixi^s,1:ndim)=0.d0
    do idims=1,ndim
      ixb^l=ixo^l-kr(idims,^d);
      ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
      {do ix^db=0,1 \}
         if({ ix^d==0 .and. ^d==idims | .or.}) then
           ixbmin^d=ixcmin^d-ix^d;
           ixbmax^d=ixcmax^d-ix^d;
           ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
         end if
      {end do\}
      ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
    end do
    src=0.d0
    if(slab_uniform) then
      do idims=1,ndim
        ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
        ixb^l=ixo^l-kr(idims,^d);
        src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
      end do
    else
      do idims=1,ndim
        ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
        ixb^l=ixo^l-kr(idims,^d);
        src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
      end do
      src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
    end if
  end subroutine get_flux_on_cell_face
 
  !> w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
  subroutine rmhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
    use mod_global_parameters
    use mod_viscosity, only: viscosity_add_source
    use mod_gravity, only: gravity_add_source
    use mod_cak_force, only: cak_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
 
    ! TODO local_timestep support is only added for splitting
    ! but not for other nonideal terms such gravity, RC, viscosity,..
    ! it will also only work for divbfix  'linde', which does not require
    ! modification as it does not use dt in the update
    if (.not. qsourcesplit) then
      if(has_equi_pe0) then
        active = .true.
        call add_pe0_divv(qdt,dtfactor,ixi^l,ixo^l,wctprim,w,x)
      end if
      if(rmhd_hyperbolic_thermal_conduction) then
        call add_hypertc_source(qdt,ixi^l,ixo^l,wct,w,x,wctprim)
      end if
      ! Source for B0 splitting
      if (b0field) then
        active = .true.
        call add_source_b0split(qdt,dtfactor,ixi^l,ixo^l,wctprim,w,x)
      end if
      ! Sources for resistivity in eqs. for e, B1, B2 and B3
      if (abs(rmhd_eta)>smalldouble)then
        active = .true.
        call add_source_res2(qdt,ixi^l,ixo^l,wct,w,x)
      end if
      if (rmhd_eta_hyper>0.d0)then
        active = .true.
        call add_source_hyperres(qdt,ixi^l,ixo^l,wct,w,x)
      end if
    end if
    {^nooned
    if(source_split_divb .eqv. qsourcesplit) then
      ! Sources related to div B
      select case (type_divb)
      case (divb_ct)
        continue ! Do nothing
      case (divb_linde)
        active = .true.
        call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
      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,wctprim,w,x)
      case (divb_janhunen)
        active = .true.
        call add_source_janhunen(qdt,ixi^l,ixo^l,wctprim,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,wctprim,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,wctprim,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_multigrid)
        continue ! Do nothing
      case (divb_none)
        ! Do nothing
      case default
        call mpistop('Unknown divB fix')
      end select
    end if
    }
    if(rmhd_viscosity) then
      call viscosity_add_source(qdt,ixi^l,ixo^l,wct,&
           w,x,rmhd_energy,qsourcesplit,active)
    end if
    if(rmhd_gravity) then
      call gravity_add_source(qdt,ixi^l,ixo^l,wct,wctprim,&
           w,x,gravity_energy,gravity_rhov,qsourcesplit,active)
    end if
    if (rmhd_cak_force) then
      call cak_add_source(qdt,ixi^l,ixo^l,wct,w,x,rmhd_energy,qsourcesplit,active)
    end if
    !> This is where the radiation force and heating/cooling are added
    call rmhd_add_radiation_source(qdt,ixi^l,ixo^l,wct,w,x,qsourcesplit,active)
    ! update temperature from new pressure, density, and old ionization degree
    if(rmhd_partial_ionization) then
      if(.not.qsourcesplit) then
        active = .true.
        call rmhd_update_temperature(ixi^l,ixo^l,wct,w,x)
      end if
    end if
  end subroutine rmhd_add_source
 
  subroutine rmhd_add_radiation_source(qdt,ixI^L,ixO^L,wCT,w,x,qsourcesplit,active)
    use mod_constants
    use mod_global_parameters
    use mod_usr_methods
    use mod_fld
    use mod_afld
    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) :: qsourcesplit
    logical, intent(inout) :: active
    double precision :: cmax(ixi^s)
 
    select case(rmhd_radiation_formalism)
    case('fld')
      if(fld_diff_scheme .eq. 'mg') call fld_get_diffcoef_central(w, wct, x, ixi^l, ixo^l)
      !> radiation force
      if(rmhd_radiation_force) call get_fld_rad_force(qdt,ixi^l,ixo^l,wct,w,x,rmhd_energy,qsourcesplit,active)
      call rmhd_handle_small_values(.true., w, x, ixi^l, ixo^l, 'fld_e_interact')
    case('afld')
      if(fld_diff_scheme .eq. 'mg') call afld_get_diffcoef_central(w, wct, x, ixi^l, ixo^l)
      !> radiation force
      if(rmhd_radiation_force) call get_afld_rad_force(qdt,ixi^l,ixo^l,wct,w,x,rmhd_energy,qsourcesplit,active)
      call rmhd_handle_small_values(.true., w, x, ixi^l, ixo^l, 'fld_e_interact')
      !> photon tiring, heating and cooling
      if(rmhd_energy) then
        if (rmhd_energy_interact) call get_afld_energy_interact(qdt,ixi^l,ixo^l,wct,w,x,rmhd_energy,qsourcesplit,active)
      endif
    case default
      call mpistop('Radiation formalism unknown')
    end select
  end subroutine rmhd_add_radiation_source
 
  subroutine add_pe0_divv(qdt,dtfactor,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,dtfactor
    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                :: divv(ixi^s)
 
    if(slab_uniform) then
      if(nghostcells .gt. 2) then
        call divvector(wct(ixi^s,mom(1:ndir)),ixi^l,ixo^l,divv,3)
      else
        call divvector(wct(ixi^s,mom(1:ndir)),ixi^l,ixo^l,divv,2)
      end if
    else
     call divvector(wct(ixi^s,mom(1:ndir)),ixi^l,ixo^l,divv)
    end if
    if(local_timestep) then
      w(ixo^s,e_)=w(ixo^s,e_)-dtfactor*block%dt(ixo^s)*block%equi_vars(ixo^s,equi_pe0_,0)*divv(ixo^s)
    else
      w(ixo^s,e_)=w(ixo^s,e_)-qdt*block%equi_vars(ixo^s,equi_pe0_,0)*divv(ixo^s)
    end if
  end subroutine add_pe0_divv
 
  subroutine get_tau(ixI^L,ixO^L,w,Te,tau,sigT5)
    use mod_global_parameters
    integer, intent(in) :: ixi^l, ixo^l
    double precision, dimension(ixI^S,1:nw), intent(in) :: w
    double precision, dimension(ixI^S), intent(in) :: te
    double precision, dimension(ixI^S), intent(out) :: tau,sigt5
    double precision :: dxmin,taumin
    double precision, dimension(ixI^S) :: sigt7,eint
    integer :: ix^d
 
    taumin=4.d0
    !> w supposed to be wCTprim here
    if(rmhd_trac) then
      where(te(ixo^s) .lt. block%wextra(ixo^s,tcoff_))
        sigt5(ixo^s)=hypertc_kappa*sqrt(block%wextra(ixo^s,tcoff_)**5)
        sigt7(ixo^s)=sigt5(ixo^s)*block%wextra(ixo^s,tcoff_)
      else where
        sigt5(ixo^s)=hypertc_kappa*sqrt(te(ixo^s)**5)
        sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
      end where
    else
      sigt5(ixo^s)=hypertc_kappa*sqrt(te(ixo^s)**5)
      sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
    end if
    eint(ixo^s)=w(ixo^s,p_)/(rmhd_gamma-one)
    tau(ixo^s)=max(taumin*dt,sigt7(ixo^s)/eint(ixo^s)/cmax_global**2)
  end subroutine get_tau
 
  subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
    use mod_global_parameters
    integer, intent(in) :: ixi^l,ixo^l
    double precision, intent(in) :: qdt
    double precision, dimension(ixI^S,1:ndim), intent(in) :: x
    double precision, dimension(ixI^S,1:nw), intent(in) :: wct,wctprim
    double precision, dimension(ixI^S,1:nw), intent(inout) :: w
    double precision :: invdx
    double precision, dimension(ixI^S) :: te,tau,sigt,htc_qsrc,tface,r
    double precision, dimension(ixI^S) :: htc_esrc,bsum,bunit
    double precision, dimension(ixI^S,1:ndim) :: btot
    integer :: idims
    integer :: hxc^l,hxo^l,ixc^l,jxc^l,jxo^l,kxc^l
 
    call rmhd_get_rfactor(wctprim,x,ixi^l,ixi^l,r)
    !Te(ixI^S)=wCTprim(ixI^S,p_)/wCT(ixI^S,rho_)
    te(ixi^s)=wctprim(ixi^s,p_)/(r(ixi^s)*w(ixi^s,rho_))
    call get_tau(ixi^l,ixo^l,wctprim,te,tau,sigt)
    htc_qsrc=zero
    do idims=1,ndim
      if(b0field) then
        btot(ixo^s,idims)=wct(ixo^s,mag(idims))+block%B0(ixo^s,idims,0)
      else
        btot(ixo^s,idims)=wct(ixo^s,mag(idims))
      endif
    enddo
    bsum(ixo^s)=sqrt(sum(btot(ixo^s,:)**2,dim=ndim+1))+smalldouble
    do idims=1,ndim
      invdx=1.d0/dxlevel(idims)
      ixc^l=ixo^l;
      ixcmin^d=ixomin^d-kr(idims,^d);ixcmax^d=ixomax^d;
      jxc^l=ixc^l+kr(idims,^d);
      kxc^l=jxc^l+kr(idims,^d);
      hxc^l=ixc^l-kr(idims,^d);
      hxo^l=ixo^l-kr(idims,^d);
      tface(ixc^s)=(7.d0*(te(ixc^s)+te(jxc^s))-(te(hxc^s)+te(kxc^s)))/12.d0
      bunit(ixo^s)=btot(ixo^s,idims)/bsum(ixo^s)
      htc_qsrc(ixo^s)=htc_qsrc(ixo^s)+sigt(ixo^s)*bunit(ixo^s)*(tface(ixo^s)-tface(hxo^s))*invdx
    end do
    htc_qsrc(ixo^s)=(htc_qsrc(ixo^s)+wct(ixo^s,q_))/tau(ixo^s)
    w(ixo^s,q_)=w(ixo^s,q_)-qdt*htc_qsrc(ixo^s)
  end subroutine add_hypertc_source
 
  !> Compute the Lorentz force (JxB)
  subroutine get_lorentz_force(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)
    ! 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)
    integer                         :: idir, idirmin
 
    b=0.0d0
    if(b0field) then
      do idir = 1, ndir
        b(ixo^s, idir) = w(ixo^s,mag(idir))+block%B0(ixo^s,idir,0)
      end do
    else
      do idir = 1, ndir
        b(ixo^s, idir) = w(ixo^s,mag(idir))
      end do
    end if
    ! 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_force
 
  !> Compute 1/(1+v_A^2/c^2) for semirelativistic MHD, where v_A is the Alfven
  !> velocity
  subroutine rmhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
    use mod_global_parameters
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s, nw)
    double precision, intent(out) :: gamma_a2(ixo^s)
    double precision              :: rho(ixi^s)
 
    ! rmhd_get_rho cannot be used as x is not a param
    if(has_equi_rho0) then
      rho(ixo^s) = w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,b0i)
    else
      rho(ixo^s) = w(ixo^s,rho_)
    endif
    ! Compute the inverse of 1 + B^2/(rho * c^2)
    gamma_a2(ixo^s) = 1.0d0/(1.0d0+rmhd_mag_en_all(w, ixi^l, ixo^l)/rho(ixo^s)*inv_squared_c)
  end subroutine rmhd_gamma2_alfven
 
  !> Compute 1/sqrt(1+v_A^2/c^2) for semirelativisitic MHD, where v_A is the
  !> Alfven velocity
  function rmhd_gamma_alfven(w, ixI^L, ixO^L) result(gamma_A)
    use mod_global_parameters
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s, nw)
    double precision              :: gamma_a(ixo^s)
 
    call rmhd_gamma2_alfven(ixi^l, ixo^l, w, gamma_a)
    gamma_a = sqrt(gamma_a)
  end function rmhd_gamma_alfven
 
  subroutine rmhd_get_rho(w,x,ixI^L,ixO^L,rho)
    use mod_global_parameters
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s,1:nw),x(ixi^s,1:ndim)
    double precision, intent(out) :: rho(ixi^s)
 
    if(has_equi_rho0) then
      rho(ixo^s) = w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,b0i)
    else
      rho(ixo^s) = w(ixo^s,rho_)
    endif
  end subroutine rmhd_get_rho
 
  !> handle small or negative internal energy
  subroutine rmhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
    use mod_global_parameters
    use mod_small_values
    integer, intent(in)             :: ixi^l,ixo^l, ie
    double precision, intent(inout) :: w(ixi^s,1:nw)
    double precision, intent(in)    :: x(ixi^s,1:ndim)
    character(len=*), intent(in)    :: subname
    double precision :: rho(ixi^s)
    integer :: idir
    logical :: flag(ixi^s,1:nw)
 
    flag=.false.
    if(has_equi_pe0) then
      where(w(ixo^s,ie)+block%equi_vars(ixo^s,equi_pe0_,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_pe0) then
          where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
                  block%equi_vars(ixo^s,equi_pe0_,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_)=w(ixo^s,e_)*gamma_1
        call rmhd_get_rho(w,x,ixi^l,ixo^l,rho)
        do idir = 1, ndir
           w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))/rho(ixo^s)
        end do
        call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
      end select
    end if
  end subroutine rmhd_handle_small_ei
 
  subroutine rmhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
    use mod_global_parameters
    use mod_ionization_degree
 
    integer, intent(in)             :: ixi^l, ixo^l
    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 :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
 
    call ionization_degree_from_temperature(ixi^l,ixo^l,wct(ixi^s,te_),iz_h,iz_he)
 
    call rmhd_get_pthermal(w,x,ixi^l,ixo^l,pth)
 
    w(ixo^s,te_)=(2.d0+3.d0*he_abundance)*pth(ixo^s)/(w(ixo^s,rho_)*(1.d0+iz_h(ixo^s)+&
     he_abundance*(iz_he(ixo^s)*(iz_he(ixo^s)+1.d0)+1.d0)))
  end subroutine rmhd_update_temperature
 
  !> Source terms after split off time-independent magnetic field
  subroutine add_source_b0split(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
    use mod_global_parameters
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: qdt, dtfactor,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)
      if(local_timestep) then
        do idir=1,3
          axb(ixo^s,idir)=axb(ixo^s,idir)*block%dt(ixo^s)*dtfactor
        enddo
      else
        axb(ixo^s,:)=axb(ixo^s,:)*qdt
      endif
      ! add J0xB0 source term in momentum equations
      w(ixo^s,mom(1:ndir))=w(ixo^s,mom(1:ndir))+axb(ixo^s,1:ndir)
    end if
    if(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
      a(ixi^s,1:ndir)=wct(ixi^s,mom(1:ndir))
      ! -E = a x b
      call cross_product(ixi^l,ixo^l,a,b,axb)
      if(local_timestep) then
        do idir=1,3
          axb(ixo^s,idir)=axb(ixo^s,idir)*block%dt(ixo^s)*dtfactor
        enddo
      else
        axb(ixo^s,:)=axb(ixo^s,:)*qdt
      endif
      ! add -(vxB) dot J0 source term in energy equation
      do idir=7-2*ndir,3
        w(ixo^s,e_)=w(ixo^s,e_)-axb(ixo^s,idir)*block%J0(ixo^s,idir)
      end do
    end if
    if (fix_small_values) call rmhd_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 (rmhd_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 (rmhd_eta>zero)then
       eta(ixa^s)=rmhd_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 (rmhd_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 (rmhd_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(total_energy) then
          w(ixo^s,e_)=w(ixo^s,e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
       end if
    end do ! idir
    if(rmhd_energy) then
      ! de/dt+=eta*J**2
      w(ixo^s,e_)=w(ixo^s,e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
    end if
    if (fix_small_values) call rmhd_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)
    tmpvec=zero
    if(rmhd_eta>zero)then
      do idir=idirmin,3
        tmpvec(ixa^s,idir)=current(ixa^s,idir)*rmhd_eta
      end do
    else
      call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
      do idir=idirmin,3
        tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
      end do
    end if
    ! dB/dt= -curl(J*eta), thus B_i=B_i-eps_ijk d_j Jeta_k
    call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
    if(stagger_grid) then
      if(ndim==2.and.ndir==3) then
        ! if 2.5D
        w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
      end if
    else
      w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
    end if
    if(rmhd_energy) then
      if(rmhd_eta>zero)then
        tmp(ixo^s)=qdt*rmhd_eta*sum(current(ixo^s,:)**2,dim=ndim+1)
      else
        tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
      end if
      if(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_)=w(ixo^s,e_)+tmp(ixo^s)-&
        qdt*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_)=w(ixo^s,e_)+tmp(ixo^s)
      end if
    end if
    if (fix_small_values) call rmhd_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)*rmhd_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(total_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_)=w(ixo^s,e_)+tmp(ixo^s)*qdt
    end if
    if (fix_small_values)  call rmhd_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 or GLM-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), gradpsi(ixi^s), ba(ixo^s,1:ndir)
    integer          :: idir
 
    ! dPsi/dt =  - Ch^2/Cp^2 Psi
    if (rmhd_glm_alpha < zero) then
      w(ixo^s,psi_) = abs(rmhd_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*rmhd_glm_alpha/minval(dxlevel(:)))*w(ixo^s,psi_)
      else
        w(ixo^s,psi_) = dexp(-qdt*cmax_global*rmhd_glm_alpha/minval(block%ds(ixo^s,:),dim=ndim+1))*w(ixo^s,psi_)
      end if
    end if
    if(rmhd_glm_extended) then
      if(b0field) then
        ba(ixo^s,1:ndir)=wct(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,0)
      else
        ba(ixo^s,1:ndir)=wct(ixo^s,mag(1:ndir))
      end if
      ! gradient of Psi
      if(total_energy) then
        do idir=1,ndim
          select case(typegrad)
          case("central")
            call gradient(wct(ixi^s,psi_),ixi^l,ixo^l,idir,gradpsi)
          case("limited")
            call gradientl(wct(ixi^s,psi_),ixi^l,ixo^l,idir,gradpsi)
          end select
          ! e  = e  -qdt (b . grad(Psi))
          w(ixo^s,e_) = w(ixo^s,e_)-qdt*ba(ixo^s,idir)*gradpsi(ixo^s)
        end do
      end if
      ! We calculate now div B
      call get_divb(wct,ixi^l,ixo^l,divb,rmhd_divb_nth)
      ! m = m - qdt b div b
      do idir=1,ndir
        w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-qdt*ba(ixo^s,idir)*divb(ixo^s)
      end do
    end if
    if (fix_small_values) call rmhd_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), ba(1:ndir)
    integer                         :: idir, ix^d
 
    ! calculate div B
    call get_divb(wct,ixi^l,ixo^l,divb,rmhd_divb_nth)
    if(b0field) then
     {do ix^db=ixomin^db,ixomax^db\}
        ! b = b - qdt v * div b
        ^c&w(ix^d,b^c_)=w(ix^d,b^c_)-qdt*wct(ix^d,m^c_)*divb(ix^d)\
        ! m = m - qdt b div b
        ^c&w(ix^d,m^c_)=w(ix^d,m^c_)-qdt*(wct(ix^d,b^c_)+block%B0(ix^d,^c,0))*divb(ix^d)\
        if (total_energy) then
          ! e = e - qdt (v . b) * div b
          w(ix^d,e_)=w(ix^d,e_)-qdt*(^c&wct(ix^d,m^c_)*(wct(ix^d,b^c_)+block%B0(ix^d,^c,0))+)*divb(ix^d)
        end if
     {end do\}
    else
     {do ix^db=ixomin^db,ixomax^db\}
        ! b = b - qdt v * div b
        ^c&w(ix^d,b^c_)=w(ix^d,b^c_)-qdt*wct(ix^d,m^c_)*divb(ix^d)\
        ! m = m - qdt b div b
        ^c&w(ix^d,m^c_)=w(ix^d,m^c_)-qdt*wct(ix^d,b^c_)*divb(ix^d)\
        if (total_energy) then
          ! e = e - qdt (v . b) * div b
          w(ix^d,e_)=w(ix^d,e_)-qdt*(^c&wct(ix^d,m^c_)*wct(ix^d,b^c_)+)*divb(ix^d)
        end if
     {end do\}
    end if
    if (fix_small_values) call rmhd_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)
    integer                         :: idir, ix^d
 
    ! calculate div B
    call get_divb(wct,ixi^l,ixo^l,divb,rmhd_divb_nth)
    {do ix^db=ixomin^db,ixomax^db\}
      ! b = b - qdt v * div b
      ^c&w(ix^d,b^c_)=w(ix^d,b^c_)-qdt*wct(ix^d,m^c_)*divb(ix^d)\
    {end do\}
    if (fix_small_values) call rmhd_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)
    double precision :: divb(ixi^s),graddivb(ixi^s)
    integer :: idim, idir, ixp^l, i^d, iside
    logical, dimension(-1:1^D&) :: leveljump
 
    ! Calculate div B
    ixp^l=ixo^l^ladd1;
    call get_divb(wct,ixi^l,ixp^l,divb,rmhd_divb_nth)
    ! 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 gradientl(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. total_energy) then
         ! e += B_idim*eta*grad_idim(divb)
         w(ixp^s,e_)=w(ixp^s,e_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
       end if
    end do
    if (fix_small_values) call rmhd_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(rmhd_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
    ! 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)
    integer :: idir, idirmin0
 
    idirmin0 = 7-2*ndir
    call curlvector(w(ixi^s,mag(1:ndir)),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
 
  !> If resistivity is not zero, check diffusion time limit for dt
  subroutine rmhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
    use mod_global_parameters
    use mod_usr_methods
    use mod_viscosity, only: viscosity_get_dt
    use mod_gravity, only: gravity_get_dt
    use mod_cak_force, only: cak_get_dt
    use mod_fld, only: fld_radforce_get_dt
    use mod_afld, only: afld_radforce_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)
    double precision              :: dxarr(ndim)
    double precision              :: current(ixi^s,7-2*ndir:3),eta(ixi^s)
    integer                       :: idirmin,idim
 
    dtnew = bigdouble
 
    if (.not. dt_c) then
      ^d&dxarr(^d)=dx^d;
      if (rmhd_eta>zero)then
         dtnew=dtdiffpar*minval(dxarr(1:ndim))**2/rmhd_eta
      else if (rmhd_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(rmhd_eta_hyper>zero) then
        if(slab_uniform) then
          dtnew=min(dtdiffpar*minval(dxarr(1:ndim))**4/rmhd_eta_hyper,dtnew)
        else
          dtnew=min(dtdiffpar*minval(block%ds(ixo^s,1:ndim))**4/rmhd_eta_hyper,dtnew)
        end if
      end if
      if(rmhd_radiation_force) then
        select case(rmhd_radiation_formalism)
        case('fld')
          call fld_radforce_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
        case('afld')
          call afld_radforce_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
        case default
          call mpistop('Radiation formalism unknown')
        end select
      endif
      if(rmhd_viscosity) then
        call viscosity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
      end if
      if(rmhd_gravity) then
        call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
      end if
      if (rmhd_cak_force) then
        call cak_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
      end if
    else
      {^ifoned dtnew = dx1*unit_velocity/const_c}
      {^nooned dtnew = min(dx^d*unit_velocity/const_c)}
    endif
  end subroutine rmhd_get_dt
 
  ! Add geometrical source terms to w
  subroutine rmhd_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)
    double precision :: tmp,tmp1,invr,cot
    integer          :: ix^d
    integer :: mr_,mphi_ ! Polar var. names
    integer :: br_,bphi_
 
    mr_=mom(1); mphi_=mom(1)-1+phi_  ! Polar var. names
    br_=mag(1); bphi_=mag(1)-1+phi_
    select case (coordinate)
    case (cylindrical)
     {do ix^db=ixomin^db,ixomax^db\}
        ! include dt in invr, invr is always used with qdt
        if(local_timestep) then
          invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
        else
          invr=qdt/x(ix^d,1)
        end if
        if(rmhd_energy) then
          tmp=wprim(ix^d,p_)+half*(^c&wprim(ix^d,b^c_)**2+)
        else
          tmp=rmhd_adiab*wprim(ix^d,rho_)**rmhd_gamma+half*(^c&wprim(ix^d,b^c_)**2+)
        end if
        if(phi_>0) then
          w(ix^d,mr_)=w(ix^d,mr_)+invr*(tmp-&
                    wprim(ix^d,bphi_)**2+wprim(ix^d,mphi_)*wct(ix^d,mphi_))
          w(ix^d,mphi_)=w(ix^d,mphi_)+invr*(&
                   -wct(ix^d,mphi_)*wprim(ix^d,mr_) &
                   +wprim(ix^d,bphi_)*wprim(ix^d,br_))
          if(.not.stagger_grid) then
            w(ix^d,bphi_)=w(ix^d,bphi_)+invr*&
                     (wprim(ix^d,bphi_)*wprim(ix^d,mr_) &
                     -wprim(ix^d,br_)*wprim(ix^d,mphi_))
          end if
        else
          w(ix^d,mr_)=w(ix^d,mr_)+invr*tmp
        end if
        if(rmhd_glm) w(ix^d,br_)=w(ix^d,br_)+wprim(ix^d,psi_)*invr
     {end do\}
    case (spherical)
     {do ix^db=ixomin^db,ixomax^db\}
        ! include dt in invr, invr is always used with qdt
        if(local_timestep) then
          invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
        else
          invr=qdt/x(ix^d,1)
        end if
        if(rmhd_energy) then
          tmp1=wprim(ix^d,p_)+half*(^c&wprim(ix^d,b^c_)**2+)
        else
          tmp1=rmhd_adiab*wprim(ix^d,rho_)**rmhd_gamma+half*(^c&wprim(ix^d,b^c_)**2+)
        end if
        ! m1
        {^ifonec
        w(ix^d,mom(1))=w(ix^d,mom(1))+two*tmp1*invr
        }
        {^noonec
        w(ix^d,mom(1))=w(ix^d,mom(1))+invr*&
         (two*tmp1+(^ce&wprim(ix^d,m^ce_)*wct(ix^d,m^ce_)-wprim(ix^d,b^ce_)**2+))
        }
        ! b1
        if(rmhd_glm) then
          w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,psi_)
        end if
        {^ifoned
        cot=0.d0
        }
        {^nooned
        cot=1.d0/tan(x(ix^d,2))
        }
        {^iftwoc
        ! m2
        w(ix^d,mom(2))=w(ix^d,mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
          +wprim(ix^d,b1_)*wprim(ix^d,b2_))
        ! b2
        if(.not.stagger_grid) then
          tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
          if(rmhd_glm) then
            tmp=tmp+wprim(ix^d,psi_)*cot
          end if
          w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
        end if
        }
        {^ifthreec
        ! m2
        w(ix^d,mom(2))=w(ix^d,mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
          +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
          +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
        ! b2
        if(.not.stagger_grid) then
          tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
          if(rmhd_glm) then
            tmp=tmp+wprim(ix^d,psi_)*cot
          end if
          w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
        end if
        ! m3
        w(ix^d,mom(3))=w(ix^d,mom(3))-invr*&
             (wprim(ix^d,m3_)*wct(ix^d,m1_) &
             -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
            +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
             -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
        ! b3
        if(.not.stagger_grid) then
          w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
             (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
             -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
            -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
             -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
        end if
        }
     {end do\}
    end select
  end subroutine rmhd_add_source_geom
 
  ! Add geometrical source terms to w
  subroutine rmhd_add_source_geom_split(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)
    double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s),invrho(ixo^s),invr(ixo^s)
    integer          :: iw,idir, h1x^l{^nooned, h2x^l}
    integer :: mr_,mphi_ ! Polar var. names
    integer :: br_,bphi_
 
    mr_=mom(1); mphi_=mom(1)-1+phi_  ! Polar var. names
    br_=mag(1); bphi_=mag(1)-1+phi_
    if(has_equi_rho0) then
      invrho(ixo^s) = 1d0/(wct(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,b0i))
    else
      invrho(ixo^s) = 1d0/wct(ixo^s,rho_)
    end if
    ! include dt in invr, invr is always used with qdt
    if(local_timestep) then
      invr(ixo^s) = block%dt(ixo^s) * dtfactor/x(ixo^s,1)
    else
      invr(ixo^s) = qdt/x(ixo^s,1)
    end if
 
    select case (coordinate)
    case (cylindrical)
      call rmhd_get_p_total(wct,x,ixi^l,ixo^l,tmp)
      if(phi_>0) then
        w(ixo^s,mr_)=w(ixo^s,mr_)+invr(ixo^s)*(tmp(ixo^s)-&
                  wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2*invrho(ixo^s))
        w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt*invr(ixo^s)*(&
                 -wct(ixo^s,mphi_)*wct(ixo^s,mr_)*invrho(ixo^s) &
                 +wct(ixo^s,bphi_)*wct(ixo^s,br_))
        if(.not.stagger_grid) then
          w(ixo^s,bphi_)=w(ixo^s,bphi_)+invr(ixo^s)*&
                   (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
                   -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
                   *invrho(ixo^s)
        end if
      else
        w(ixo^s,mr_)=w(ixo^s,mr_)+invr(ixo^s)*tmp(ixo^s)
      end if
      if(rmhd_glm) w(ixo^s,br_)=w(ixo^s,br_)+wct(ixo^s,psi_)*invr(ixo^s)
    case (spherical)
       h1x^l=ixo^l-kr(1,^d); {^nooned h2x^l=ixo^l-kr(2,^d);}
       call rmhd_get_p_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(idir))**2*invrho(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(1))=w(ixo^s,mom(1))+tmp(ixo^s)*invr(ixo^s)
       ! b1
       if(rmhd_glm) then
         w(ixo^s,mag(1))=w(ixo^s,mag(1))+invr(ixo^s)*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
       if(local_timestep) then
         tmp1(ixo^s) = block%dt(ixo^s) * tmp(ixo^s)
       else
         tmp1(ixo^s) = qdt * tmp(ixo^s)
       endif  
       ! This will make hydrostatic p=const an exact solution
       w(ixo^s,mom(2))=w(ixo^s,mom(2))+tmp1(ixo^s) &
            *(block%surfaceC(ixo^s,2)-block%surfaceC(h2x^s,2)) &
            /block%dvolume(ixo^s)
       tmp(ixo^s)=-(wct(ixo^s,mom(1))*wct(ixo^s,mom(2))*invrho(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(3))**2*invrho(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(2))=w(ixo^s,mom(2))+tmp(ixo^s)*invr(ixo^s)
       ! b2
       if(.not.stagger_grid) then
         tmp(ixo^s)=(wct(ixo^s,mom(1))*wct(ixo^s,mag(2)) &
              -wct(ixo^s,mom(2))*wct(ixo^s,mag(1)))*invrho(ixo^s)
         if(b0field) then
           tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom(1))*block%B0(ixo^s,2,0) &
                -wct(ixo^s,mom(2))*block%B0(ixo^s,1,0))*invrho(ixo^s)
         end if
         if(rmhd_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))+tmp(ixo^s)*invr(ixo^s)
       end if
       }
       if(ndir==3) then
         ! m3
         tmp(ixo^s)=-(wct(ixo^s,mom(3))*wct(ixo^s,mom(1))*invrho(ixo^s) &
              -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
              -(wct(ixo^s,mom(2))*wct(ixo^s,mom(3))*invrho(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(3))=w(ixo^s,mom(3))+tmp(ixo^s)*invr(ixo^s)
         ! b3
         if(.not.stagger_grid) then
           tmp(ixo^s)=(wct(ixo^s,mom(1))*wct(ixo^s,mag(3)) &
                -wct(ixo^s,mom(3))*wct(ixo^s,mag(1)))*invrho(ixo^s) {^nooned &
                -(wct(ixo^s,mom(3))*wct(ixo^s,mag(2)) &
                -wct(ixo^s,mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
                *invrho(ixo^s)/dsin(x(ixo^s,2)) }
           if (b0field) then
              tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom(1))*block%B0(ixo^s,3,0) &
                   -wct(ixo^s,mom(3))*block%B0(ixo^s,1,0))*invrho(ixo^s){^nooned &
                   -(wct(ixo^s,mom(3))*block%B0(ixo^s,2,0) &
                   -wct(ixo^s,mom(2))*block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
                   *invrho(ixo^s)/dsin(x(ixo^s,2)) }
           end if
           w(ixo^s,mag(3))=w(ixo^s,mag(3))+tmp(ixo^s)*invr(ixo^s)
         end if
       end if
    end select
  end subroutine rmhd_add_source_geom_split
 
  !> Compute 2 times total magnetic energy
  function rmhd_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 = sum((w(ixo^s, mag(:))+block%B0(ixo^s,:,b0i))**2, dim=ndim+1)
    else
      mge = sum(w(ixo^s, mag(:))**2, dim=ndim+1)
    end if
  end function rmhd_mag_en_all
 
  subroutine rmhd_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(ixo^s), dpsi(ixo^s)
    integer :: ix^d
 
    if(stagger_grid) then
     {do ix^db=ixomin^db,ixomax^db\}
        wlc(ix^d,mag(idir))=s%ws(ix^d,idir)
        wrc(ix^d,mag(idir))=s%ws(ix^d,idir)
        wlp(ix^d,mag(idir))=s%ws(ix^d,idir)
        wrp(ix^d,mag(idir))=s%ws(ix^d,idir)
     {end do\}
    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
     {do ix^db=ixomin^db,ixomax^db\}
        db(ix^d)=wrp(ix^d,mag(idir))-wlp(ix^d,mag(idir))
        dpsi(ix^d)=wrp(ix^d,psi_)-wlp(ix^d,psi_)
        wlp(ix^d,mag(idir))=half*(wrp(ix^d,mag(idir))+wlp(ix^d,mag(idir))-dpsi(ix^d)/cmax_global)
        wlp(ix^d,psi_)=half*(wrp(ix^d,psi_)+wlp(ix^d,psi_)-db(ix^d)*cmax_global)
        wrp(ix^d,mag(idir))=wlp(ix^d,mag(idir))
        wrp(ix^d,psi_)=wlp(ix^d,psi_)
        if(total_energy) then
          wrc(ix^d,e_)=wrc(ix^d,e_)-half*wrc(ix^d,mag(idir))**2
          wlc(ix^d,e_)=wlc(ix^d,e_)-half*wlc(ix^d,mag(idir))**2
        end if
        wrc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
        wrc(ix^d,psi_)=wlp(ix^d,psi_)
        wlc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
        wlc(ix^d,psi_)=wlp(ix^d,psi_)
        ! modify total energy according to the change of magnetic field
        if(total_energy) then
          wrc(ix^d,e_)=wrc(ix^d,e_)+half*wrc(ix^d,mag(idir))**2
          wlc(ix^d,e_)=wlc(ix^d,e_)+half*wlc(ix^d,mag(idir))**2
        end if
     {end do\}
    end if
    if(associated(usr_set_wlr)) call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
  end subroutine rmhd_modify_wlr
 
  subroutine rmhd_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 rmhd_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
       if(total_energy) call rmhd_to_primitive(ixg^l,ixo^l,w,x)
       {^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
       }
       if(total_energy) call rmhd_to_conserved(ixg^l,ixo^l,w,x)
     case(2)
       if(total_energy) call rmhd_to_primitive(ixg^l,ixo^l,w,x)
       {^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
       }
       if(total_energy) call rmhd_to_conserved(ixg^l,ixo^l,w,x)
     case(3)
       if(total_energy) call rmhd_to_primitive(ixg^l,ixo^l,w,x)
       {^iftwod
       ixfmin1=ixomin1+1
       ixfmax1=ixomax1-1
       ixfmin2=ixomin2+1
       ixfmax2=ixomax2+1
       if(slab_uniform) then
         dx2x1=dxlevel(2)/dxlevel(1)
         do ix2=ixfmax2,ixfmin2,-1
           w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
            +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
                    w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
         enddo
       else
         do ix2=ixfmax2,ixfmin2,-1
           w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
             w(ixfmin1:ixfmax1,ix2,mag(2)))*block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
           +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
             block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
           -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
             block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
            /block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
         end do
       end if
       }
       {^ifthreed
       ixfmin1=ixomin1+1
       ixfmax1=ixomax1-1
       ixfmin3=ixomin3+1
       ixfmax3=ixomax3-1
       ixfmin2=ixomin2+1
       ixfmax2=ixomax2+1
       if(slab_uniform) then
         dx2x1=dxlevel(2)/dxlevel(1)
         dx2x3=dxlevel(2)/dxlevel(3)
         do ix2=ixfmax2,ixfmin2,-1
           w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
             ix2+1,ixfmin3:ixfmax3,mag(2)) &
             +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
                     w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
             +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
                     w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
         end do
       else
         do ix2=ixfmax2,ixfmin2,-1
           w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
          ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
             w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
             block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
           +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
             w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
             block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
           -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
             w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
             block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
           +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
             w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
             block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
           -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
             w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
             block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
            /block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
             w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
         end do
       end if
       }
       if(total_energy) call rmhd_to_conserved(ixg^l,ixo^l,w,x)
     case(4)
       if(total_energy) call rmhd_to_primitive(ixg^l,ixo^l,w,x)
       {^iftwod
       ixfmin1=ixomin1+1
       ixfmax1=ixomax1-1
       ixfmin2=ixomin2-1
       ixfmax2=ixomax2-1
       if(slab_uniform) then
         dx2x1=dxlevel(2)/dxlevel(1)
         do ix2=ixfmin2,ixfmax2
           w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
            -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
                    w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
         end do
       else
         do ix2=ixfmin2,ixfmax2
           w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
             w(ixfmin1:ixfmax1,ix2,mag(2)))*block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
           -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
             block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
           +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
             block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
            /block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
         end do
       end if
       }
       {^ifthreed
       ixfmin1=ixomin1+1
       ixfmax1=ixomax1-1
       ixfmin3=ixomin3+1
       ixfmax3=ixomax3-1
       ixfmin2=ixomin2-1
       ixfmax2=ixomax2-1
       if(slab_uniform) then
         dx2x1=dxlevel(2)/dxlevel(1)
         dx2x3=dxlevel(2)/dxlevel(3)
         do ix2=ixfmin2,ixfmax2
           w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
             ix2-1,ixfmin3:ixfmax3,mag(2)) &
             -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
                     w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
             -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
                     w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
         end do
       else
         do ix2=ixfmin2,ixfmax2
           w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
          ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
             w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
             block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
           -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
             w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
             block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
           +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
             w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
             block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
           -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
             w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
             block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
           +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
             w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
             block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
            /block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
             w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
         end do
       end if
       }
       if(total_energy) call rmhd_to_conserved(ixg^l,ixo^l,w,x)
     {^ifthreed
     case(5)
       if(total_energy) call rmhd_to_primitive(ixg^l,ixo^l,w,x)
       ixfmin1=ixomin1+1
       ixfmax1=ixomax1-1
       ixfmin2=ixomin2+1
       ixfmax2=ixomax2-1
       ixfmin3=ixomin3+1
       ixfmax3=ixomax3+1
       if(slab_uniform) then
         dx3x1=dxlevel(3)/dxlevel(1)
         dx3x2=dxlevel(3)/dxlevel(2)
         do ix3=ixfmax3,ixfmin3,-1
           w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
             ixfmin2:ixfmax2,ix3+1,mag(3)) &
             +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
                     w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
             +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
                     w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
         end do
       else
         do ix3=ixfmax3,ixfmin3,-1
           w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
          ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
             w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
             block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
           +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
             w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
             block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
           -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
             w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
             block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
           +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
             w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
             block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
           -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
             w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
             block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
            /block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
             w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
         end do
       end if
       if(total_energy) call rmhd_to_conserved(ixg^l,ixo^l,w,x)
     case(6)
       if(total_energy) call rmhd_to_primitive(ixg^l,ixo^l,w,x)
       ixfmin1=ixomin1+1
       ixfmax1=ixomax1-1
       ixfmin2=ixomin2+1
       ixfmax2=ixomax2-1
       ixfmin3=ixomin3-1
       ixfmax3=ixomax3-1
       if(slab_uniform) then
         dx3x1=dxlevel(3)/dxlevel(1)
         dx3x2=dxlevel(3)/dxlevel(2)
         do ix3=ixfmin3,ixfmax3
           w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
             ixfmin2:ixfmax2,ix3-1,mag(3)) &
             -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
                     w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
             -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
                     w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
         end do
       else
         do ix3=ixfmin3,ixfmax3
           w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
          ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
             w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
             block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
           -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
             w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
             block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
           +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
             w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
             block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
           -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
             w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
             block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
           +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
             w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
             block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
            /block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
             w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
         end do
       end if
       if(total_energy) call rmhd_to_conserved(ixg^l,ixo^l,w,x)
     }
     case default
       call mpistop("Special boundary is not defined for this region")
    end select
  end subroutine fixdivb_boundary
 
  {^nooned
  subroutine rmhd_clean_divb_multigrid(qdt, qt, active)
    use mod_forest
    use mod_global_parameters
    use mod_multigrid_coupling
    use mod_geometry
    double precision, intent(in) :: qdt    !< Current time step
    double precision, intent(in) :: qt     !< Current time
    logical, intent(inout)       :: active !< Output if the source is active
    integer                      :: id
    integer, parameter           :: max_its      = 50
    double precision             :: residual_it(max_its), max_divb
    double precision             :: tmp(ixg^t), grad(ixg^t, ndim)
    double precision             :: res
    double precision, parameter  :: max_residual = 1d-3
    double precision, parameter  :: residual_reduction = 1d-10
    integer                      :: iigrid, igrid
    integer                      :: n, nc, lvl, ix^l, ixc^l, idim
    type(tree_node), pointer     :: pnode
 
    mg%operator_type = mg_laplacian
    ! Set boundary conditions
    do n = 1, 2*ndim
       idim = (n+1)/2
       select case (typeboundary(mag(idim), n))
       case (bc_symm)
          ! d/dx B = 0, take phi = 0
          mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
          mg%bc(n, mg_iphi)%bc_value = 0.0_dp
       case (bc_asymm)
          ! B = 0, so grad(phi) = 0
          mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
          mg%bc(n, mg_iphi)%bc_value = 0.0_dp
       case (bc_cont)
          mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
          mg%bc(n, mg_iphi)%bc_value = 0.0_dp
       case (bc_special)
          ! Assume Dirichlet boundary conditions, derivative zero
          mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
          mg%bc(n, mg_iphi)%bc_value = 0.0_dp
       case (bc_periodic)
          ! Nothing to do here
       case default
          write(*,*) "rmhd_clean_divb_multigrid warning: unknown boundary type"
          mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
          mg%bc(n, mg_iphi)%bc_value = 0.0_dp
       end select
    end do
 
    ix^l=ixm^ll^ladd1;
    max_divb = 0.0d0
    ! Store divergence of B as right-hand side
    do iigrid = 1, igridstail
       igrid =  igrids(iigrid);
       pnode => igrid_to_node(igrid, mype)%node
       id    =  pnode%id
       lvl   =  mg%boxes(id)%lvl
       nc    =  mg%box_size_lvl(lvl)
 
       ! Geometry subroutines expect this to be set
       block => ps(igrid)
       ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
 
       call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^ll, ixm^ll, tmp, &
            rmhd_divb_nth)
       mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(ixm^t)
       max_divb = max(max_divb, maxval(abs(tmp(ixm^t))))
    end do
 
    ! Solve laplacian(phi) = divB
    if(stagger_grid) then
      call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
           mpi_max, icomm, ierrmpi)
 
      if (mype == 0) print *, "Performing multigrid divB cleaning"
      if (mype == 0) print *, "iteration vs residual"
      ! Solve laplacian(phi) = divB
      do n = 1, max_its
         call mg_fas_fmg(mg, n>1, max_res=residual_it(n))
         if (mype == 0) write(*, "(I4,E11.3)") n, residual_it(n)
         if (residual_it(n) < residual_reduction * max_divb) exit
      end do
      if (mype == 0 .and. n > max_its) then
         print *, "divb_multigrid warning: not fully converged"
         print *, "current amplitude of divb: ", residual_it(max_its)
         print *, "multigrid smallest grid: ", &
              mg%domain_size_lvl(:, mg%lowest_lvl)
         print *, "note: smallest grid ideally has <= 8 cells"
         print *, "multigrid dx/dy/dz ratio: ", mg%dr(:, 1)/mg%dr(1, 1)
         print *, "note: dx/dy/dz should be similar"
      end if
    else
      do n = 1, max_its
         call mg_fas_vcycle(mg, max_res=res)
         if (res < max_residual) exit
      end do
      if (res > max_residual) call mpistop("divb_multigrid: no convergence")
    end if
 
    ! Correct the magnetic field
    do iigrid = 1, igridstail
       igrid = igrids(iigrid);
       pnode => igrid_to_node(igrid, mype)%node
       id    =  pnode%id
       ! Geometry subroutines expect this to be set
       block => ps(igrid)
       ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
       ! Compute the gradient of phi
       tmp(ix^s) = mg%boxes(id)%cc({:,}, mg_iphi)
       if(stagger_grid) then
         do idim =1, ndim
           ixcmin^d=ixmlo^d-kr(idim,^d);
           ixcmax^d=ixmhi^d;
           call gradientf(tmp,ps(igrid)%x,ixg^ll,ixc^l,idim,grad(ixg^t,idim))
           ! Apply the correction B* = B - gradient(phi)
           ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
         end do
         ! store cell-center magnetic energy
         tmp(ixm^t) = sum(ps(igrid)%w(ixm^t, mag(1:ndim))**2, dim=ndim+1)
         ! change cell-center magnetic field
         call rmhd_face_to_center(ixm^ll,ps(igrid))
       else
         do idim = 1, ndim
            call gradient(tmp,ixg^ll,ixm^ll,idim,grad(ixg^t, idim))
         end do
         ! store cell-center magnetic energy
         tmp(ixm^t) = sum(ps(igrid)%w(ixm^t, mag(1:ndim))**2, dim=ndim+1)
         ! Apply the correction B* = B - gradient(phi)
         ps(igrid)%w(ixm^t, mag(1:ndim)) = &
              ps(igrid)%w(ixm^t, mag(1:ndim)) - grad(ixm^t, :)
       end if
       if(total_energy) then
         ! Determine magnetic energy difference
         tmp(ixm^t) = 0.5_dp * (sum(ps(igrid)%w(ixm^t, &
              mag(1:ndim))**2, dim=ndim+1) - tmp(ixm^t))
         ! Keep thermal pressure the same
         ps(igrid)%w(ixm^t, e_) = ps(igrid)%w(ixm^t, e_) + tmp(ixm^t)
       end if
    end do
    active = .true.
  end subroutine rmhd_clean_divb_multigrid
  }
 
  subroutine rmhd_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s,vcts)
    use mod_global_parameters
    integer, intent(in)                :: ixi^l, ixo^l
    double precision, intent(in)       :: qt,qdt
    ! cell-center primitive variables
    double precision, intent(in)       :: wprim(ixi^s,1:nw)
    type(state)                        :: sct, s
    type(ct_velocity)                  :: vcts
    double precision, intent(in)       :: fc(ixi^s,1:nwflux,1:ndim)
    double precision, intent(inout)    :: fe(ixi^s,sdim:3)
 
    select case(type_ct)
    case('average')
      call update_faces_average(ixi^l,ixo^l,qt,qdt,fc,fe,sct,s)
    case('uct_contact')
      call update_faces_contact(ixi^l,ixo^l,qt,qdt,wprim,fc,fe,sct,s,vcts)
    case('uct_hll')
      call update_faces_hll(ixi^l,ixo^l,qt,qdt,fe,sct,s,vcts)
    case default
      call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
    end select
  end subroutine rmhd_update_faces
 
  !> get electric field though averaging neighors to update faces in CT
  subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
    use mod_global_parameters
    use mod_usr_methods
    integer, intent(in)                :: ixi^l, ixo^l
    double precision, intent(in)       :: qt, qdt
    type(state)                        :: sct, s
    double precision, intent(in)       :: fc(ixi^s,1:nwflux,1:ndim)
    double precision, intent(inout)    :: fe(ixi^s,sdim:3)
    double precision                   :: circ(ixi^s,1:ndim)
    ! non-ideal electric field on cell edges
    double precision, dimension(ixI^S,sdim:3) :: e_resi
    integer                            :: ix^d,ixc^l,ixa^l,i1kr^d,i2kr^d
    integer                            :: idim1,idim2,idir,iwdim1,iwdim2
 
    associate(bfaces=>s%ws,x=>s%x)
    ! Calculate contribution to FEM of each edge,
    ! that is, estimate value of line integral of
    ! electric field in the positive idir direction.
    ! if there is resistivity, get eta J
    if(rmhd_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,e_resi)
    do idim1=1,ndim
      iwdim1 = mag(idim1)
      i1kr^d=kr(idim1,^d);
      do idim2=1,ndim
        iwdim2 = mag(idim2)
        i2kr^d=kr(idim2,^d);
        do idir=sdim,3! Direction of line integral
          ! Allow only even permutations
          if (lvc(idim1,idim2,idir)==1) then
            ixcmax^d=ixomax^d;
            ixcmin^d=ixomin^d+kr(idir,^d)-1;
            ! average cell-face electric field to cell edges
           {do ix^db=ixcmin^db,ixcmax^db\}
              fe(ix^d,idir)=quarter*&
                (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
                -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
              ! add resistive electric field at cell edges E=-vxB+eta J
              if(rmhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
              ! times time step and edge length
              fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
           {end do\}
          end if
        end do
      end do
    end do
    ! allow user to change inductive electric field, especially for boundary driven applications
    if(associated(usr_set_electric_field)) &
      call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
    circ(ixi^s,1:ndim)=zero
    ! Calculate circulation on each face
    do idim1=1,ndim ! Coordinate perpendicular to face
      ixcmax^d=ixomax^d;
      ixcmin^d=ixomin^d-kr(idim1,^d);
      do idim2=1,ndim
        ixa^l=ixc^l-kr(idim2,^d);
        do idir=sdim,3 ! Direction of line integral
          ! Assemble indices
          if(lvc(idim1,idim2,idir)==1) then
            ! Add line integrals in direction idir
            circ(ixc^s,idim1)=circ(ixc^s,idim1)&
                             +(fe(ixc^s,idir)&
                              -fe(ixa^s,idir))
          else if(lvc(idim1,idim2,idir)==-1) then
            ! Add line integrals in direction idir
            circ(ixc^s,idim1)=circ(ixc^s,idim1)&
                             -(fe(ixc^s,idir)&
                              -fe(ixa^s,idir))
          end if
        end do
      end do
      ! Divide by the area of the face to get dB/dt
      where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
        circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
      elsewhere
        circ(ixc^s,idim1)=zero
      end where
      ! Time update cell-face magnetic field component
      bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
    end do
    end associate
  end subroutine update_faces_average
 
  !> update faces using UCT contact mode by Gardiner and Stone 2005 JCP 205, 509
  subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
    use mod_global_parameters
    use mod_usr_methods
    use mod_geometry
    integer, intent(in)                :: ixi^l, ixo^l
    double precision, intent(in)       :: qt, qdt
    ! cell-center primitive variables
    double precision, intent(in)       :: wp(ixi^s,1:nw)
    type(state)                        :: sct, s
    type(ct_velocity)                  :: vcts
    double precision, intent(in)       :: fc(ixi^s,1:nwflux,1:ndim)
    double precision, intent(inout)    :: fe(ixi^s,sdim:3)
    double precision                   :: circ(ixi^s,1:ndim)
    ! electric field at cell centers
    double precision                   :: ecc(ixi^s,sdim:3)
    double precision                   :: ein(ixi^s,sdim:3)
    ! gradient of E at left and right side of a cell face
    double precision                   :: el(ixi^s),er(ixi^s)
    ! gradient of E at left and right side of a cell corner
    double precision                   :: elc,erc
    ! non-ideal electric field on cell edges
    double precision, dimension(ixI^S,sdim:3) :: e_resi
    ! current on cell edges
    double precision :: jce(ixi^s,sdim:3)
    ! location at cell faces
    double precision :: xs(ixgs^t,1:ndim)
    double precision :: gradi(ixgs^t)
    integer :: ixc^l,ixa^l
    integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^d,i1kr^d,i2kr^d
 
    associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
    ! if there is resistivity, get eta J
    if(rmhd_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,e_resi)
    if(b0field) then
     {do ix^db=iximin^db,iximax^db\}
        ! Calculate electric field at cell centers
       {^ifthreed
        ecc(ix^d,1)=(wp(ix^d,b2_)+block%B0(ix^d,2,0))*wp(ix^d,m3_)-(wp(ix^d,b3_)+block%B0(ix^d,3,0))*wp(ix^d,m2_)
        ecc(ix^d,2)=(wp(ix^d,b3_)+block%B0(ix^d,3,0))*wp(ix^d,m1_)-(wp(ix^d,b1_)+block%B0(ix^d,1,0))*wp(ix^d,m3_)
        ecc(ix^d,3)=(wp(ix^d,b1_)+block%B0(ix^d,1,0))*wp(ix^d,m2_)-(wp(ix^d,b2_)+block%B0(ix^d,2,0))*wp(ix^d,m1_)
       }
       {^iftwod
        ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
       }
       {^ifoned
        ecc(ix^d,3)=0.d0
       }
     {end do\}
    else
     {do ix^db=iximin^db,iximax^db\}
        ! Calculate electric field at cell centers
       {^ifthreed
        ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
        ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
        ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
       }
       {^iftwod
        ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
       }
       {^ifoned
        ecc(ix^d,3)=0.d0
       }
     {end do\}
    end if
 
    ! Calculate contribution to FEM of each edge,
    ! that is, estimate value of line integral of
    ! electric field in the positive idir direction.
    ! evaluate electric field along cell edges according to equation (41)
    do idim1=1,ndim
      iwdim1 = mag(idim1)
      i1kr^d=kr(idim1,^d);
      do idim2=1,ndim
        iwdim2 = mag(idim2)
        i2kr^d=kr(idim2,^d);
        do idir=sdim,3 ! Direction of line integral
          ! Allow only even permutations
          if (lvc(idim1,idim2,idir)==1) then
            ixcmax^d=ixomax^d;
            ixcmin^d=ixomin^d+kr(idir,^d)-1;
            ! Assemble indices
            ! average cell-face electric field to cell edges
           {do ix^db=ixcmin^db,ixcmax^db\}
              fe(ix^d,idir)=quarter*&
              (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
              -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
              if(partial_energy) ein(ix^d,idir)=fe(ix^d,idir)
           {end do\}
            ! add slope in idim2 direction from equation (50)
            ixamin^d=ixcmin^d;
            ixamax^d=ixcmax^d+i1kr^d;
           {do ix^db=ixamin^db,ixamax^db\}
              el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
              er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
           {end do\}
           {!dir$ ivdep
            do ix^db=ixcmin^db,ixcmax^db\}
              if(vnorm(ix^d,idim1)>0.d0) then
                elc=el(ix^d)
              else if(vnorm(ix^d,idim1)<0.d0) then
                elc=el({ix^d+i1kr^d})
              else
                elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
              end if
              if(vnorm({ix^d+i2kr^d},idim1)>0.d0) then
                erc=er(ix^d)
              else if(vnorm({ix^d+i2kr^d},idim1)<0.d0) then
                erc=er({ix^d+i1kr^d})
              else
                erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
              end if
              fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
           {end do\}
            ! add slope in idim1 direction from equation (50)
            ixamin^d=ixcmin^d;
            ixamax^d=ixcmax^d+i2kr^d;
           {do ix^db=ixamin^db,ixamax^db\}
              el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
              er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
           {end do\}
           {!dir$ ivdep
            do ix^db=ixcmin^db,ixcmax^db\}
              if(vnorm(ix^d,idim2)>0.d0) then
                elc=el(ix^d)
              else if(vnorm(ix^d,idim2)<0.d0) then
                elc=el({ix^d+i2kr^d})
              else
                elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
              end if
              if(vnorm({ix^d+i1kr^d},idim2)>0.d0) then
                erc=er(ix^d)
              else if(vnorm({ix^d+i1kr^d},idim2)<0.d0) then
                erc=er({ix^d+i2kr^d})
              else
                erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
              end if
              fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
              ! difference between average and upwind interpolated E
              if(partial_energy) ein(ix^d,idir)=fe(ix^d,idir)-ein(ix^d,idir)
              ! add resistive electric field at cell edges E=-vxB+eta J
              if(rmhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
              ! times time step and edge length
              fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
           {end do\}
          end if
        end do
      end do
    end do
    if(partial_energy) then
      ! add upwind diffused magnetic energy back to energy
      ! calculate current density at cell edges
      jce=0.d0
      do idim1=1,ndim
        do idim2=1,ndim
          do idir=sdim,3
            if (lvc(idim1,idim2,idir)==0) cycle
            ixcmax^d=ixomax^d;
            ixcmin^d=ixomin^d+kr(idir,^d)-1;
            ixamax^d=ixcmax^d-kr(idir,^d)+1;
            ixamin^d=ixcmin^d;
            ! current at transverse faces
            xs(ixa^s,:)=x(ixa^s,:)
            xs(ixa^s,idim2)=x(ixa^s,idim2)+half*s%dx(ixa^s,idim2)
            call gradientf(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi)
            if (lvc(idim1,idim2,idir)==1) then
              jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
            else
              jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
            end if
          end do
        end do
      end do
      do idir=sdim,3
        ixcmax^d=ixomax^d;
        ixcmin^d=ixomin^d+kr(idir,^d)-1;
        ! E dot J on cell edges
        ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
        ! average from cell edge to cell center
       {^ifthreed
        if(idir==1) then
         {do ix^db=ixomin^db,ixomax^db\}
            jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1,ix2-1,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
                          +ein(ix1,ix2-1,ix3-1,idir))
            if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
            w(ix^d,e_)=w(ix^d,e_)+qdt*jce(ix^d,idir)
         {end do\}
        else if(idir==2) then
         {do ix^db=ixomin^db,ixomax^db\}
            jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
                          +ein(ix1-1,ix2,ix3-1,idir))
            if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
            w(ix^d,e_)=w(ix^d,e_)+qdt*jce(ix^d,idir)
         {end do\}
        else
         {do ix^db=ixomin^db,ixomax^db\}
            jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2-1,ix3,idir)&
                          +ein(ix1-1,ix2-1,ix3,idir))
            if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
            w(ix^d,e_)=w(ix^d,e_)+qdt*jce(ix^d,idir)
         {end do\}
        end if
       }
       {^iftwod
        !idir=3
       {do ix^db=ixomin^db,ixomax^db\}
          jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,idir)+ein(ix1,ix2-1,idir)&
                        +ein(ix1-1,ix2-1,idir))
          if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
          w(ix^d,e_)=w(ix^d,e_)+qdt*jce(ix^d,idir)
       {end do\}
       }
        ! save additional numerical resistive heating to an extra variable
        if(nwextra>0) then
          block%w(ixo^s,nw)=block%w(ixo^s,nw)+jce(ixo^s,idir)
        end if
      end do
    end if
    ! allow user to change inductive electric field, especially for boundary driven applications
    if(associated(usr_set_electric_field)) &
      call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
    circ(ixi^s,1:ndim)=zero
    ! Calculate circulation on each face
    do idim1=1,ndim ! Coordinate perpendicular to face
      ixcmax^d=ixomax^d;
      ixcmin^d=ixomin^d-kr(idim1,^d);
      do idim2=1,ndim
        ixa^l=ixc^l-kr(idim2,^d);
        do idir=sdim,3 ! Direction of line integral
          ! Assemble indices
          if(lvc(idim1,idim2,idir)==1) then
            ! Add line integrals in direction idir
            circ(ixc^s,idim1)=circ(ixc^s,idim1)&
                             +(fe(ixc^s,idir)&
                              -fe(ixa^s,idir))
          else if(lvc(idim1,idim2,idir)==-1) then
            ! Add line integrals in direction idir
            circ(ixc^s,idim1)=circ(ixc^s,idim1)&
                             -(fe(ixc^s,idir)&
                              -fe(ixa^s,idir))
          end if
        end do
      end do
      ! Divide by the area of the face to get dB/dt
      where(s%surfaceC(ixc^s,idim1) > smalldouble)
        circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
      elsewhere
        circ(ixc^s,idim1)=zero
      end where
      ! Time update cell-face magnetic field component
      bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
    end do
    end associate
  end subroutine update_faces_contact
 
  !> update faces
  subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s,vcts)
    use mod_global_parameters
    use mod_constrained_transport
    use mod_usr_methods
    integer, intent(in)                :: ixi^l, ixo^l
    double precision, intent(in)       :: qt, qdt
    double precision, intent(inout)    :: fe(ixi^s,sdim:3)
    type(state)                        :: sct, s
    type(ct_velocity)                  :: vcts
    double precision                   :: vtill(ixi^s,2)
    double precision                   :: vtilr(ixi^s,2)
    double precision                   :: bfacetot(ixi^s,ndim)
    double precision                   :: btill(ixi^s,ndim)
    double precision                   :: btilr(ixi^s,ndim)
    double precision                   :: cp(ixi^s,2)
    double precision                   :: cm(ixi^s,2)
    double precision                   :: circ(ixi^s,1:ndim)
    ! non-ideal electric field on cell edges
    double precision, dimension(ixI^S,sdim:3) :: e_resi
    integer                            :: hxc^l,ixc^l,ixcp^l,jxc^l,ixcm^l
    integer                            :: idim1,idim2,idir
 
    associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
      cbarmax=>vcts%cbarmax)
    ! Calculate contribution to FEM of each edge,
    ! that is, estimate value of line integral of
    ! electric field in the positive idir direction.
 
    ! Loop over components of electric field
 
    ! idir: electric field component we need to calculate
    ! idim1: directions in which we already performed the reconstruction
    ! idim2: directions in which we perform the reconstruction
 
    ! if there is resistivity, get eta J
    if(rmhd_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,e_resi)
 
    do idir=sdim,3
      ! Indices
      ! idir: electric field component
      ! idim1: one surface
      ! idim2: the other surface
      ! cyclic permutation: idim1,idim2,idir=1,2,3
      ! Velocity components on the surface
      ! follow cyclic premutations:
      ! Sx(1),Sx(2)=y,z ; Sy(1),Sy(2)=z,x ; Sz(1),Sz(2)=x,y
      ixcmax^d=ixomax^d;
      ixcmin^d=ixomin^d-1+kr(idir,^d);
      ! Set indices and directions
      idim1=mod(idir,3)+1
      idim2=mod(idir+1,3)+1
      jxc^l=ixc^l+kr(idim1,^d);
      ixcp^l=ixc^l+kr(idim2,^d);
      ! Reconstruct transverse transport velocities
      call reconstruct(ixi^l,ixc^l,idim2,vbarc(ixi^s,idim1,1),&
               vtill(ixi^s,2),vtilr(ixi^s,2))
      call reconstruct(ixi^l,ixc^l,idim1,vbarc(ixi^s,idim2,2),&
               vtill(ixi^s,1),vtilr(ixi^s,1))
      ! Reconstruct magnetic fields
      ! Eventhough the arrays are larger, reconstruct works with
      ! the limits ixG.
      if(b0field) then
        bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+block%B0(ixi^s,idim1,idim1)
        bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+block%B0(ixi^s,idim2,idim2)
      else
        bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
        bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
      end if
      call reconstruct(ixi^l,ixc^l,idim2,bfacetot(ixi^s,idim1),&
               btill(ixi^s,idim1),btilr(ixi^s,idim1))
      call reconstruct(ixi^l,ixc^l,idim1,bfacetot(ixi^s,idim2),&
               btill(ixi^s,idim2),btilr(ixi^s,idim2))
      ! Take the maximum characteristic
      cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
      cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
      cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
      cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
      ! Calculate eletric field
      fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
                     + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
                     - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
                     /(cp(ixc^s,1)+cm(ixc^s,1)) &
                     +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
                     + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
                     - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
                     /(cp(ixc^s,2)+cm(ixc^s,2))
      ! add resistive electric field at cell edges E=-vxB+eta J
      if(rmhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
      fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
      if (.not.slab) then
        where(abs(x(ixc^s,r_)+half*dxlevel(r_)).lt.1.0d-9)
          fe(ixc^s,idir)=zero
        end where
      end if
    end do
    ! allow user to change inductive electric field, especially for boundary driven applications
    if(associated(usr_set_electric_field)) &
      call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
    circ(ixi^s,1:ndim)=zero
    ! Calculate circulation on each face: interal(fE dot dl)
    do idim1=1,ndim ! Coordinate perpendicular to face
      ixcmax^d=ixomax^d;
      ixcmin^d=ixomin^d-kr(idim1,^d);
      do idim2=1,ndim
        do idir=sdim,3 ! Direction of line integral
          ! Assemble indices
          if(lvc(idim1,idim2,idir)/=0) then
            hxc^l=ixc^l-kr(idim2,^d);
            ! Add line integrals in direction idir
            circ(ixc^s,idim1)=circ(ixc^s,idim1)&
                             +lvc(idim1,idim2,idir)&
                             *(fe(ixc^s,idir)&
                              -fe(hxc^s,idir))
          end if
        end do
      end do
      ! Divide by the area of the face to get dB/dt
      where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
        circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
      elsewhere
        circ(ixc^s,idim1)=zero
      end where
      ! Time update cell-face magnetic field component
      bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
    end do
    end associate
  end subroutine update_faces_hll
 
  !> calculate eta J at cell edges
  subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
    use mod_global_parameters
    use mod_usr_methods
    use mod_geometry
    integer, intent(in)                :: ixi^l, ixo^l
    type(state), intent(in)            :: sct, s
    ! current on cell edges
    double precision :: jce(ixi^s,sdim:3)
    ! current on cell centers
    double precision :: jcc(ixi^s,7-2*ndir:3)
    ! location at cell faces
    double precision :: xs(ixgs^t,1:ndim)
    ! resistivity
    double precision :: eta(ixi^s)
    double precision :: gradi(ixgs^t)
    integer :: ix^d,ixc^l,ixa^l,ixb^l,idir,idirmin,idim1,idim2
 
    associate(x=>s%x,dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
    ! calculate current density at cell edges
    jce=0.d0
    do idim1=1,ndim
      do idim2=1,ndim
        do idir=sdim,3
          if (lvc(idim1,idim2,idir)==0) cycle
          ixcmax^d=ixomax^d;
          ixcmin^d=ixomin^d+kr(idir,^d)-1;
          ixbmax^d=ixcmax^d-kr(idir,^d)+1;
          ixbmin^d=ixcmin^d;
          ! current at transverse faces
          xs(ixb^s,:)=x(ixb^s,:)
          xs(ixb^s,idim2)=x(ixb^s,idim2)+half*dx(ixb^s,idim2)
          call gradientf(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi,2)
          if (lvc(idim1,idim2,idir)==1) then
            jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
          else
            jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
          end if
        end do
      end do
    end do
    ! get resistivity
    if(rmhd_eta>zero)then
      jce(ixi^s,:)=jce(ixi^s,:)*rmhd_eta
    else
      ixa^l=ixo^l^ladd1;
      call get_current(wct,ixi^l,ixa^l,idirmin,jcc)
      call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,jcc,eta)
      ! calcuate eta on cell edges
      do idir=sdim,3
        ixcmax^d=ixomax^d;
        ixcmin^d=ixomin^d+kr(idir,^d)-1;
        jcc(ixc^s,idir)=0.d0
       {do ix^db=0,1\}
          if({ ix^d==1 .and. ^d==idir | .or.}) cycle
          ixamin^d=ixcmin^d+ix^d;
          ixamax^d=ixcmax^d+ix^d;
          jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
       {end do\}
        jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
        jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
      end do
    end if
    end associate
  end subroutine get_resistive_electric_field
 
  !> calculate cell-center values from face-center values
  subroutine rmhd_face_to_center(ixO^L,s)
    use mod_global_parameters
    ! Non-staggered interpolation range
    integer, intent(in) :: ixo^l
    type(state) :: s
    integer :: ix^d
 
    ! calculate cell-center values from face-center values in 2nd order
    ! because the staggered arrays have an additional place to the left.
    ! Interpolate to cell barycentre using arithmetic average
    ! This might be done better later, to make the method less diffusive.
    {!dir$ ivdep
    do ix^db=ixomin^db,ixomax^db\}
      {^ifthreed
      s%w(ix^d,b1_)=half/s%surface(ix^d,1)*(s%ws(ix^d,1)*s%surfaceC(ix^d,1)&
        +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
      s%w(ix^d,b2_)=half/s%surface(ix^d,2)*(s%ws(ix^d,2)*s%surfaceC(ix^d,2)&
        +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
      s%w(ix^d,b3_)=half/s%surface(ix^d,3)*(s%ws(ix^d,3)*s%surfaceC(ix^d,3)&
        +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
      }
      {^iftwod
      s%w(ix^d,b1_)=half/s%surface(ix^d,1)*(s%ws(ix^d,1)*s%surfaceC(ix^d,1)&
        +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
      s%w(ix^d,b2_)=half/s%surface(ix^d,2)*(s%ws(ix^d,2)*s%surfaceC(ix^d,2)&
        +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
      }
   {end do\}
    ! calculate cell-center values from face-center values in 4th order
    !do idim=1,ndim
    !  gxO^L=ixO^L-2*kr(idim,^D);
    !  hxO^L=ixO^L-kr(idim,^D);
    !  jxO^L=ixO^L+kr(idim,^D);
 
    !  ! Interpolate to cell barycentre using fourth order central formula
    !  w(ixO^S,mag(idim))=(0.0625d0/s%surface(ixO^S,idim))*&
    !         ( -ws(gxO^S,idim)*s%surfaceC(gxO^S,idim) &
    !     +9.0d0*ws(hxO^S,idim)*s%surfaceC(hxO^S,idim) &
    !     +9.0d0*ws(ixO^S,idim)*s%surfaceC(ixO^S,idim) &
    !           -ws(jxO^S,idim)*s%surfaceC(jxO^S,idim) )
    !end do
 
    ! calculate cell-center values from face-center values in 6th order
    !do idim=1,ndim
    !  fxO^L=ixO^L-3*kr(idim,^D);
    !  gxO^L=ixO^L-2*kr(idim,^D);
    !  hxO^L=ixO^L-kr(idim,^D);
    !  jxO^L=ixO^L+kr(idim,^D);
    !  kxO^L=ixO^L+2*kr(idim,^D);
 
    !  ! Interpolate to cell barycentre using sixth order central formula
    !  w(ixO^S,mag(idim))=(0.00390625d0/s%surface(ixO^S,idim))* &
    !     (  +3.0d0*ws(fxO^S,idim)*s%surfaceC(fxO^S,idim) &
    !       -25.0d0*ws(gxO^S,idim)*s%surfaceC(gxO^S,idim) &
    !      +150.0d0*ws(hxO^S,idim)*s%surfaceC(hxO^S,idim) &
    !      +150.0d0*ws(ixO^S,idim)*s%surfaceC(ixO^S,idim) &
    !       -25.0d0*ws(jxO^S,idim)*s%surfaceC(jxO^S,idim) &
    !        +3.0d0*ws(kxO^S,idim)*s%surfaceC(kxO^S,idim) )
    !end do
  end subroutine rmhd_face_to_center
 
  !> calculate magnetic field from vector potential
  subroutine b_from_vector_potential(ixIs^L, ixI^L, ixO^L, ws, x)
    use mod_global_parameters
    use mod_constrained_transport
    integer, intent(in)                :: ixis^l, ixi^l, ixo^l
    double precision, intent(inout)    :: ws(ixis^s,1:nws)
    double precision, intent(in)       :: x(ixi^s,1:ndim)
    double precision                   :: adummy(ixis^s,1:3)
 
    call b_from_vector_potentiala(ixis^l, ixi^l, ixo^l, ws, x, adummy)
  end subroutine b_from_vector_potential
 
  subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
    use mod_global_parameters
    use mod_ionization_degree
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s,1:nw)
    double precision, intent(in) :: x(ixi^s,1:ndim)
    double precision, intent(out):: rfactor(ixi^s)
    double precision :: iz_h(ixo^s),iz_he(ixo^s)
 
    call ionization_degree_from_temperature(ixi^l,ixo^l,w(ixi^s,te_),iz_h,iz_he)
    ! assume the first and second ionization of Helium have the same degree
    rfactor(ixo^s)=(1.d0+iz_h(ixo^s)+0.1d0*(1.d0+iz_he(ixo^s)*(1.d0+iz_he(ixo^s))))/(2.d0+3.d0*he_abundance)
  end subroutine rfactor_from_temperature_ionization
 
  subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
    use mod_global_parameters
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s,1:nw)
    double precision, intent(in) :: x(ixi^s,1:ndim)
    double precision, intent(out):: rfactor(ixi^s)
 
    rfactor(ixo^s)=rr
  end subroutine rfactor_from_constant_ionization
end module mod_rmhd_phys