mod_mf_phys.t

Go to the documentation of this file.

!> Magnetofriction module
module mod_mf_phys
  use mod_global_parameters, only: std_len
  use mod_physics
  use mod_functions_bfield, only: get_divb,mag
  use mod_comm_lib, only: mpistop
 
 
  implicit none
  private
 
  !> viscosity coefficient in s cm^-2 for solar corona (Cheung 2012 ApJ)
  double precision, public                :: mf_nu = 1.d-15
 
  !> maximal limit of magnetofrictional velocity in cm s^-1 (Pomoell 2019 A&A)
  double precision, public                :: mf_vmax = 3.d6
 
  !> decay scale of frictional velocity near boundaries
  double precision, public                :: mf_decay_scale(2*^nd)=0.d0
 
  !> GLM-MHD parameter: ratio of the diffusive and advective time scales for div b
  !> taking values within [0, 1]
  double precision, public                :: mf_glm_alpha = 0.5d0
 
  !> The resistivity
  double precision, public                :: mf_eta = 0.0d0
 
  !> The hyper-resistivity
  double precision, public                :: mf_eta_hyper = 0.0d0
 
  !> Coefficient of diffusive divB cleaning
  double precision :: divbdiff     = 0.8d0
 
  !> Helium abundance over Hydrogen
  double precision, public, protected  :: he_abundance=0.1d0
 
  !> Method type in a integer for good performance
  integer :: type_divb
 
  !> 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_
 
  !> Indices of the magnetic field for the form of better vectorization
  integer, public, protected              :: ^c&b^C_
 
  !> Indices of the GLM psi
  integer, public, protected :: psi_
 
  !> 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 particles module is added
  logical, public, protected              :: mf_particles = .false.
 
  !> Whether GLM-MHD is used
  logical, public, protected              :: mf_glm = .false.
 
  !> Whether divB cleaning sources are added splitting from fluid solver
  logical, public, protected              :: source_split_divb = .false.
 
  !> set to true if need to record electric field on cell edges
  logical, public, protected              :: mf_record_electric_field = .false.
 
  !> Whether divB is computed with a fourth order approximation
  integer, public, protected :: mf_divb_nth = 1
 
  !> To control divB=0 fix for boundary
  logical, public, protected :: boundary_divbfix(2*^nd)=.true.
 
  !> Use a compact way to add resistivity
  logical :: compactres   = .false.
 
  !> Add divB wave in Roe solver
  logical, public :: divbwave     = .true.
 
  !> clean divb in the initial condition
  logical, public, protected :: clean_initial_divb=.false.
 
  !> Method type to clean divergence of B
  character(len=std_len), public, protected :: typedivbfix  = 'ct'
 
  !> Method type of constrained transport
  character(len=std_len), public, protected :: type_ct  = 'average'
 
  !> Update all equations due to divB cleaning
  character(len=std_len) ::    typedivbdiff = 'all'
 
 
  ! Public methods
  public :: mf_phys_init
  public :: mf_to_conserved
  public :: mf_to_primitive
  public :: mf_face_to_center
  public :: get_divb
  public :: get_current
  public :: get_normalized_divb
  public :: b_from_vector_potential
  public :: record_force_free_metrics
  {^nooned
  public :: mf_clean_divb_multigrid
  }
 
contains
 
  !> Read this module's parameters from a file
  subroutine mf_read_params(files)
    use mod_global_parameters
    use mod_particles, only: particles_eta
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /mf_list/ mf_nu, mf_vmax, mf_decay_scale, &
      mf_eta, mf_eta_hyper, mf_glm_alpha, mf_particles,&
      particles_eta, mf_record_electric_field,&
      typedivbfix, source_split_divb, divbdiff,&
      typedivbdiff, type_ct, compactres, divbwave, he_abundance, si_unit, &
      bdip, bquad, boct, busr, clean_initial_divb, &
      boundary_divbfix, boundary_divbfix_skip, mf_divb_nth
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status="old")
       read(unitpar, mf_list, end=111)
111    close(unitpar)
    end do
 
  end subroutine mf_read_params
 
  !> Write this module's parameters to a snapsoht
  subroutine mf_write_info(fh)
    use mod_global_parameters
    integer, intent(in)                 :: fh
    integer, parameter                  :: n_par = 1
    double precision                    :: values(n_par)
    character(len=name_len)             :: names(n_par)
    integer, dimension(MPI_STATUS_SIZE) :: st
    integer                             :: er
 
    call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
 
    names(1) = "nu"
    values(1) = mf_nu
    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 mf_write_info
 
  subroutine mf_phys_init()
    use mod_global_parameters
    {^nooned
    use mod_multigrid_coupling
    }
 
    integer :: itr, idir
 
    call mf_read_params(par_files)
 
    physics_type = "mf"
    phys_energy = .false.
    use_particles=mf_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 => mf_clean_divb_multigrid
    }
    case ('glm')
      mf_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')
      mf_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
 
    ! set velocity field as flux variables
    allocate(mom(ndir))
    mom(:) = var_set_momentum(ndir)
    m^c_=mom(^c);
 
    ! set magnetic field as flux variables
    allocate(mag(ndir))
    mag(:) = var_set_bfield(ndir)
    b^c_=mag(^c);
 
    ! start with magnetic field and skip velocity when update ghostcells
    iwstart=mag(1)
    allocate(start_indices(number_species),stop_indices(number_species))
    ! set the index of the first flux variable for species 1
    start_indices(1)=mag(1)
 
    if (mf_glm) then
      psi_ = var_set_fluxvar('psi', 'psi', need_bc=.false.)
    else
      psi_ = -1
    end if
 
    ! set number of variables which need update ghostcells
    nwgc=nwflux-ndir
 
    ! set the index of the last flux variable for species 1
    stop_indices(1)=nwflux
 
    ! 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, nw))
       flux_type = flux_default
    else if (any(shape(flux_type) /= [ndir, nw])) then
       call mpistop("phys_check error: flux_type has wrong shape")
    end if
    do idir=1,ndir
      if(ndim>1) flux_type(idir,mag(idir))=flux_tvdlf
    end do
    if(mf_glm .and. ndim>1) flux_type(:,psi_)=flux_tvdlf
 
    phys_get_dt              => mf_get_dt
    phys_get_cmax            => mf_get_cmax
    phys_get_cbounds         => mf_get_cbounds
    phys_get_flux            => mf_get_flux
    phys_add_source_geom     => mf_add_source_geom
    phys_add_source          => mf_add_source
    phys_to_conserved        => mf_to_conserved
    phys_to_primitive        => mf_to_primitive
    phys_check_params        => mf_check_params
    phys_write_info          => mf_write_info
    phys_special_advance     => mf_velocity_update
 
    if(type_divb==divb_glm) then
      phys_modify_wlr => mf_modify_wlr
    end if
 
    ! pass to global variable to record electric field
    record_electric_field=mf_record_electric_field
 
 
    ! if using ct stagger grid, boundary divb=0 is not done here
    if(stagger_grid) then
      phys_face_to_center => mf_face_to_center
      select case(type_ct)
      case('average')
        transverse_ghost_cells = 1
        phys_get_ct_velocity => mf_get_ct_velocity_average
        phys_update_faces => mf_update_faces_average
      case('uct_contact')
        transverse_ghost_cells = 1
        phys_get_ct_velocity => mf_get_ct_velocity_contact
        phys_update_faces => mf_update_faces_contact
      case('uct_hll')
        transverse_ghost_cells = 2
        phys_get_ct_velocity => mf_get_ct_velocity_hll
        phys_update_faces => mf_update_faces_hll
      case default
        call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
      end select
    else if(ndim>1) then
      phys_boundary_adjust => mf_boundary_adjust
    end if
 
    {^nooned
    ! clean initial divb
    if(clean_initial_divb) phys_clean_divb => mf_clean_divb_multigrid
    }
 
    ! derive units from basic units
    call mf_physical_units()
 
  end subroutine mf_phys_init
 
  subroutine mf_check_params
    use mod_global_parameters
    use mod_geometry, only: coordinate
    use mod_particles, only: particles_init, particles_eta, particles_etah
    use mod_particles, only: npayload,nusrpayload, &
                ngridvars,num_particles,physics_type_particles
 
    ! Initialize particles module
    if(mf_particles) then
      call particles_init()
      ! never allow Hall effects in particles when doing magnetofrictional
      particles_etah=0.0d0
    end if
 
    if(mype==0)then 
           write(*,*)'====MF run with settings===================='
           write(*,*)'Using mod_mf_phys with settings:'
           write(*,*)'SI_unit=',si_unit
           write(*,*)'Dimensionality   :',ndim
           write(*,*)'vector components:',ndir
           write(*,*)'coordinate set to type,slab:',coordinate,slab
           write(*,*)'number of variables          nw=',nw
           write(*,*)'    start index         iwstart=',iwstart
           write(*,*)'number of      vector variables=',nvector
           write(*,*)'number of stagger variables nws=',nws
           write(*,*)'number of    variables with BCs=',nwgc
           write(*,*)'number of      vars with fluxes=',nwflux
           write(*,*)'number of   vars with flux + BC=',nwfluxbc
           write(*,*)'number of   auxiliary variables=',nwaux
           write(*,*)'number of extra vars without flux=',nwextra
           write(*,*)'number of extra vars   for wextra=',nw_extra
           write(*,*)'number of auxiliary I/O variables=',nwauxio
           write(*,*)'    mf_eta=',mf_eta,' nonzero implies resistivity'
           write(*,*)'    mf_eta_hyper=',mf_eta_hyper
           write(*,*)'    mf_particles=',mf_particles
           if(mf_particles) then
              write(*,*) '*****Using particles:     with mf_eta, mf_eta_hyper :', mf_eta, mf_eta_hyper
              write(*,*) '*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
              write(*,*) '*****Using particles: npayload,ngridvars :', npayload,ngridvars
              write(*,*) '*****Using particles:        nusrpayload :', nusrpayload
              write(*,*) '*****Using particles:      num_particles :', num_particles
              write(*,*) '*****Using particles: physics_type_particles=',physics_type_particles
           end if
           write(*,*)'number of             ghostcells=',nghostcells
           write(*,*)'number due to phys_wider_stencil=',phys_wider_stencil
           write(*,*)'==========================================='
    endif
 
 
  end subroutine mf_check_params
 
  subroutine mf_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
      c_lightspeed=const_c
    end if
    a=1.d0+4.d0*he_abundance
    b=2.d0+3.d0*he_abundance
    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
 
    ! get dimensionless mf nu
    mf_nu=mf_nu/unit_time*unit_length**2
 
    ! get dimensionless maximal mf velocity limit
    mf_vmax=mf_vmax/unit_velocity
 
  end subroutine mf_physical_units
 
  !> Transform primitive variables into conservative ones
  subroutine mf_to_conserved(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
    ! nothing to do for mf
  end subroutine mf_to_conserved
 
  !> Transform conservative variables into primitive ones
  subroutine mf_to_primitive(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
    ! nothing to do for mf
  end subroutine mf_to_primitive
 
  !> Calculate cmax_idim=csound+abs(v_idim) within ixO^L
  subroutine mf_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
    use mod_global_parameters
 
    integer, intent(in)          :: ixi^l, ixo^l, idim
    double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
    double precision, intent(inout) :: cmax(ixi^s)
 
    cmax(ixo^s)=abs(w(ixo^s,mom(idim)))+one
 
  end subroutine mf_get_cmax
 
  !> Estimating bounds for the minimum and maximum signal velocities
  subroutine mf_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
    use mod_global_parameters
    use mod_variables
 
    integer, intent(in)             :: ixi^l, ixo^l, idim
    double precision, intent(in)    :: wlc(ixi^s, nw), wrc(ixi^s, nw)
    double precision, intent(in)    :: wlp(ixi^s, nw), wrp(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s,1:ndim)
    double precision, intent(inout) :: cmax(ixi^s,1:number_species)
    double precision, intent(inout), optional :: cmin(ixi^s,1:number_species)
    double precision, intent(in)    :: hspeed(ixi^s,1:number_species)
 
    double precision, dimension(ixI^S) :: tmp1
 
    tmp1(ixo^s)=0.5d0*(wlc(ixo^s,mom(idim))+wrc(ixo^s,mom(idim)))
    if(present(cmin)) then
      cmax(ixo^s,1)=max(tmp1(ixo^s)+one,zero)
      cmin(ixo^s,1)=min(tmp1(ixo^s)-one,zero)
    else
      cmax(ixo^s,1)=abs(tmp1(ixo^s))+one
    end if
 
  end subroutine mf_get_cbounds
 
  !> prepare velocities for ct methods
  subroutine mf_get_ct_velocity_average(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
 
  end subroutine mf_get_ct_velocity_average
 
  subroutine mf_get_ct_velocity_contact(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
 
    ! calculate velocities related to different UCT schemes
    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)))
 
  end subroutine mf_get_ct_velocity_contact
 
  subroutine mf_get_ct_velocity_hll(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
 
    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))
 
  end subroutine mf_get_ct_velocity_hll
 
  !> Calculate fluxes within ixO^L.
  subroutine mf_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
    use mod_global_parameters
    use mod_usr_methods
    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)
 
    integer :: ix^d
 
   {do ix^db=ixomin^db,ixomax^db\}
      ! compute flux of magnetic field
      ! 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(mf_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
 
  end subroutine mf_get_flux
 
  !> Add global source terms to update frictional velocity on complete domain
  subroutine mf_velocity_update(qt,psa)
    use mod_global_parameters
    use mod_ghostcells_update
    double precision, intent(in) :: qt     !< Current time
    type(state), target :: psa(max_blocks) !< Compute based on this state
 
    integer :: iigrid,igrid
    logical :: stagger_flag
    logical :: firstmf=.true.
 
    if(firstmf) then
      ! point bc mpi datatype to partial type for velocity field
      type_send_srl=>type_send_srl_p1
      type_recv_srl=>type_recv_srl_p1
      type_send_r=>type_send_r_p1
      type_recv_r=>type_recv_r_p1
      type_send_p=>type_send_p_p1
      type_recv_p=>type_recv_p_p1
      call create_bc_mpi_datatype(mom(1),ndir)
      firstmf=.false.
    end if
 
    !$OMP PARALLEL DO PRIVATE(igrid)
    do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
      block=>psa(igrid)
      call frictional_velocity(psa(igrid)%w,psa(igrid)%x,ixg^ll,ixm^ll)
    end do
    !$OMP END PARALLEL DO
 
    ! only update mf velocity in ghost cells
    stagger_flag=stagger_grid
    type_send_srl=>type_send_srl_p1
    type_recv_srl=>type_recv_srl_p1
    type_send_r=>type_send_r_p1
    type_recv_r=>type_recv_r_p1
    type_send_p=>type_send_p_p1
    type_recv_p=>type_recv_p_p1
    stagger_grid=.false.
    bcphys=.false.
    call getbc(qt,0.d0,psa,mom(1),ndir)
    bcphys=.true.
    type_send_srl=>type_send_srl_f
    type_recv_srl=>type_recv_srl_f
    type_send_r=>type_send_r_f
    type_recv_r=>type_recv_r_f
    type_send_p=>type_send_p_f
    type_recv_p=>type_recv_p_f
    stagger_grid=stagger_flag
 
  end subroutine mf_velocity_update
 
  !> w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
  subroutine mf_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
    use mod_global_parameters
 
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: qdt, dtfactor
    double precision, intent(in)    :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:ndim)
    double precision, intent(inout) :: w(ixi^s,1:nw)
    logical, intent(in)             :: qsourcesplit
    logical, intent(inout)            :: active
 
    if (.not. qsourcesplit) then
 
      ! Sources for resistivity in eqs. for e, B1, B2 and B3
      if (abs(mf_eta)>smalldouble)then
        active = .true.
        call add_source_res2(qdt,ixi^l,ixo^l,wct,w,x)
      end if
 
      if (mf_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_none)
        ! Do nothing
      case (divb_glm)
        active = .true.
        call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
      case (divb_powel)
        active = .true.
        call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
      case (divb_janhunen)
        active = .true.
        call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
      case (divb_linde)
        active = .true.
        call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
      case (divb_lindejanhunen)
        active = .true.
        call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
        call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
      case (divb_lindepowel)
        active = .true.
        call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
        call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
      case (divb_lindeglm)
        active = .true.
        call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
        call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
      case (divb_ct)
        continue ! Do nothing
      case (divb_multigrid)
        continue ! Do nothing
      case default
        call mpistop('Unknown divB fix')
      end select
    end if
    }
 
  end subroutine mf_add_source
 
  subroutine frictional_velocity(w,x,ixI^L,ixO^L)
    use mod_global_parameters
 
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: x(ixi^s,1:ndim)
    double precision, intent(inout) :: w(ixi^s,1:nw)
 
    double precision :: xmin(ndim),xmax(ndim)
    double precision :: decay(ixo^s)
    double precision :: current(ixi^s,7-2*ndir:3),tmp(ixo^s)
    integer :: ix^d, idirmin,idir,jdir,kdir
 
    call get_current(w,ixi^l,ixo^l,idirmin,current)
    ! extrapolate current for the outmost layer, 
    ! because extrapolation of B at boundaries introduces artificial current
{
    if(block%is_physical_boundary(2*^d-1)) then
      if(slab) then
        ! exclude boundary 5 in 3D Cartesian
        if(2*^d-1==5.and.ndim==3) then
        else
          current(ixomin^d^d%ixO^s,:)=current(ixomin^d+1^d%ixO^s,:)
        end if
      else
        ! exclude boundary 1 spherical/cylindrical
        if(2*^d-1>1) then
          current(ixomin^d^d%ixO^s,:)=current(ixomin^d+1^d%ixO^s,:)
        end if
      end if
    end if
    if(block%is_physical_boundary(2*^d)) then
      current(ixomax^d^d%ixO^s,:)=current(ixomax^d-1^d%ixO^s,:)
    end if
\}
    w(ixo^s,mom(:))=0.d0
    ! calculate Lorentz force
    do idir=1,ndir; do jdir=idirmin,3; do kdir=1,ndir
      if(lvc(idir,jdir,kdir)/=0) then
        if(lvc(idir,jdir,kdir)==1) then
          w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+current(ixo^s,jdir)*w(ixo^s,mag(kdir))
        else
          w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-current(ixo^s,jdir)*w(ixo^s,mag(kdir))
        end if
      end if
    end do; end do; end do
 
    tmp(ixo^s)=sum(w(ixo^s,mag(:))**2,dim=ndim+1)
    ! frictional coefficient
    where(tmp(ixo^s)/=0.d0)
      tmp(ixo^s)=1.d0/(tmp(ixo^s)*mf_nu)
    endwhere
 
    do idir=1,ndir
      w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*tmp(ixo^s)
    end do
 
    ! decay frictional velocity near selected boundaries
    ^d&xmin(^d)=xprobmin^d\
    ^d&xmax(^d)=xprobmax^d\
    decay(ixo^s)=1.d0
    do idir=1,ndim
      if(mf_decay_scale(2*idir-1)>0.d0) decay(ixo^s)=decay(ixo^s)*&
         (1.d0-exp(-(x(ixo^s,idir)-xmin(idir))/mf_decay_scale(2*idir-1)))
      if(mf_decay_scale(2*idir)>0.d0) decay(ixo^s)=decay(ixo^s)*&
         (1.d0-exp((x(ixo^s,idir)-xmax(idir))/mf_decay_scale(2*idir)))
    end do
 
    ! saturate mf velocity at mf_vmax
    tmp(ixo^s)=sqrt(sum(w(ixo^s,mom(:))**2,dim=ndim+1))/mf_vmax+1.d-12
    tmp(ixo^s)=dtanh(tmp(ixo^s))/tmp(ixo^s)
    do idir=1,ndir
      w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*tmp(ixo^s)*decay(ixo^s)
    end do
 
  end subroutine frictional_velocity
 
  !> 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)
 
    ! 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)
    double precision :: tmp(ixi^s),tmp2(ixi^s)
    integer :: ixa^l,idir,jdir,kdir,idirmin,idim,jxo^l,hxo^l,ix
    integer :: lxo^l, kxo^l
 
    ! Calculating resistive sources involve one extra layer
      ixa^l=ixo^l^ladd2;
 
    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 (mf_eta>zero)then
       eta(ixa^s)=mf_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
 
    bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
 
    do idir=1,ndir
       ! Put B_idir into tmp2 and eta*Laplace B_idir into tmp
         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
 
       ! 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 (mf_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)
 
    end do ! idir
 
  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)
    integer :: ixa^l,idir,idirmin,idirmin1
 
    ! resistive term already added as an electric field component
    if(stagger_grid .and. ndim==ndir) return
 
    ixa^l=ixo^l^ladd2;
 
    if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
         call mpistop("Error in add_source_res2: Non-conforming input limits")
 
    ixa^l=ixo^l^ladd1;
    ! Calculate current density within ixL: J=curl B, thus J_i=eps_ijk*d_j B_k
    ! Determine exact value of idirmin while doing the loop.
    call get_current(wct,ixi^l,ixa^l,idirmin,current)
 
    if (mf_eta>zero)then
       eta(ixa^s)=mf_eta
    else
       call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
    end if
 
    ! dB/dt= -curl(J*eta), thus B_i=B_i-eps_ijk d_j Jeta_k
    tmpvec(ixa^s,1:ndir)=zero
    do idir=idirmin,3
       tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
    end do
    curlj=0.d0
    call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
    if(stagger_grid) 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
 
  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)*mf_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
 
  end subroutine add_source_hyperres
 
  subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
    ! Add divB related sources to w within ixO
    ! corresponding to Dedner JCP 2002, 175, 645 _equation 24_
    ! giving the EGLM-MHD scheme
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:ndim)
    double precision, intent(inout) :: w(ixi^s,1:nw)
    double precision:: divb(ixi^s)
    double precision :: gradpsi(ixi^s)
    integer          :: idim,idir
 
    ! We calculate now div B
    call get_divb(wct,ixi^l,ixo^l,divb,mf_divb_nth)
 
    ! dPsi/dt =  - Ch^2/Cp^2 Psi
    if (mf_glm_alpha < zero) then
      w(ixo^s,psi_) = abs(mf_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*mf_glm_alpha/minval(dxlevel(:)))*w(ixo^s,psi_)
      else
        w(ixo^s,psi_) = dexp(-qdt*cmax_global*mf_glm_alpha/minval(block%ds(ixo^s,:),dim=ndim+1))*w(ixo^s,psi_)
      end if
    end if
 
    ! gradient of Psi
    do idim=1,ndim
       select case(typegrad)
       case("central")
          call gradient(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
       case("limited")
          call gradientl(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
       end select
    end do
  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)
    integer                         :: idir
 
    ! We calculate now div B
    call get_divb(wct,ixi^l,ixo^l,divb,mf_divb_nth)
 
    ! b = b - qdt v * div b
    do idir=1,ndir
      w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,mom(idir))*divb(ixo^s)
    end do
 
  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
 
    ! We calculate now div B
    call get_divb(wct,ixi^l,ixo^l,divb,mf_divb_nth)
 
    ! b = b - qdt v * div b
    do idir=1,ndir
      w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,mom(idir))*divb(ixo^s)
    end do
 
  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,mf_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)
    end do
 
  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(sum(w(ixo^s, mag(:))**2, dim=ndim+1))
    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,fourthorder)
    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
    logical, intent(in), optional :: fourthorder
 
    ! 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
 
    if(present(fourthorder)) then
      call curlvector(w(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,idirmin0,ndir,fourthorder)
    else
      call curlvector(w(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
    end if
 
  end subroutine get_current
 
  !> If resistivity is not zero, check diffusion time limit for dt
  subroutine mf_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
    use mod_global_parameters
    use mod_usr_methods
 
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(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
 
    ^d&dxarr(^d)=dx^d;
    if (mf_eta>zero)then
       dtnew=dtdiffpar*minval(dxarr(1:ndim))**2/mf_eta
    else if (mf_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(mf_eta_hyper>zero) then
      if(slab_uniform) then
        dtnew=min(dtdiffpar*minval(dxarr(1:ndim))**4/mf_eta_hyper,dtnew)
      else
        dtnew=min(dtdiffpar*minval(block%ds(ixo^s,1:ndim))**4/mf_eta_hyper,dtnew)
      end if
    end if
  end subroutine mf_get_dt
 
  ! Add geometrical source terms to w
  subroutine mf_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,invr,cs2
    integer          :: iw,idir,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)
      if(phi_>0) then
        if(.not.stagger_grid) then
         {do ix^db=ixomin^db,ixomax^db\}
            w(ix^d,bphi_)=w(ix^d,bphi_)+qdt/x(ix^d,1)*&
                     (wct(ix^d,bphi_)*wct(ix^d,mr_) &
                     -wct(ix^d,br_)*wct(ix^d,mphi_))
         {end do\}
        end if
      end if
      if(mf_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,psi_)/x(ixo^s,1)
    case (spherical)
       {do ix^db=ixomin^db,ixomax^db\}
          invr=qdt/x(ix^d,1)
         ! b1
         if(mf_glm) then
           w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wct(ix^d,psi_)
         end if
 
         {^nooned
         ! b2
         if(.not.stagger_grid) then
           tmp=wct(ix^d,mom(1))*wct(ix^d,mag(2))-wct(ix^d,mom(2))*wct(ix^d,mag(1))
           cs2=1.d0/tan(x(ix^d,2))
           if(mf_glm) then
             tmp=tmp+cs2*wct(ix^d,psi_)
           end if
           w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
         end if
         }
 
         if(ndir==3) then
          {^ifoned
          ! b3
          if(.not.stagger_grid) then
            w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
               (wct(ix^d,mom(1))*wct(ix^d,mag(3)) &
               -wct(ix^d,mom(3))*wct(ix^d,mag(1)))
          end if
          }
          {^nooned
          ! b3
          if(.not.stagger_grid) then
            w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
               (wct(ix^d,mom(1))*wct(ix^d,mag(3)) &
               -wct(ix^d,mom(3))*wct(ix^d,mag(1)) &
              -(wct(ix^d,mom(3))*wct(ix^d,mag(2)) &
               -wct(ix^d,mom(2))*wct(ix^d,mag(3)))*cs2)
          end if
          }
         end if
       {end do\}
    end select
 
  end subroutine mf_add_source_geom
 
  subroutine mf_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
    use mod_global_parameters
    use mod_usr_methods
    integer, intent(in)             :: ixi^l, ixo^l, idir
    double precision, intent(in)    :: qt
    double precision, intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
    double precision, intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
    type(state)                     :: s
    double precision                :: db(ixi^s), dpsi(ixi^s)
 
    ! Solve the Riemann problem for the linear 2x2 system for normal
    ! B-field and GLM_Psi according to Dedner 2002:
    ! This implements eq. (42) in Dedner et al. 2002 JcP 175
    ! Gives the Riemann solution on the interface
    ! for the normal B component and Psi in the GLM-MHD system.
    ! 23/04/2013 Oliver Porth
    db(ixo^s)   = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
    dpsi(ixo^s) = wrp(ixo^s,psi_) - wlp(ixo^s,psi_)
 
    wlp(ixo^s,mag(idir))   = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
         - 0.5d0/cmax_global * dpsi(ixo^s)
    wlp(ixo^s,psi_)       = 0.5d0 * (wrp(ixo^s,psi_) + wlp(ixo^s,psi_)) &
         - 0.5d0*cmax_global * db(ixo^s)
 
    wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
    wrp(ixo^s,psi_) = wlp(ixo^s,psi_)
 
  end subroutine mf_modify_wlr
 
  subroutine mf_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 mf_boundary_adjust
 
  subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
    use mod_global_parameters
 
    integer, intent(in) :: ixg^l,ixo^l,ib
    double precision, intent(inout) :: w(ixg^s,1:nw)
    double precision, intent(in) :: x(ixg^s,1:ndim)
 
    double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
    integer :: ix^d,ixf^l
 
    select case(ib)
     case(1)
       ! 2nd order CD for divB=0 to set normal B component better
       {^iftwod
       ixfmin1=ixomin1+1
       ixfmax1=ixomax1+1
       ixfmin2=ixomin2+1
       ixfmax2=ixomax2-1
       if(slab_uniform) then
         dx1x2=dxlevel(1)/dxlevel(2)
         do ix1=ixfmax1,ixfmin1,-1
           w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
            +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
                    w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
         enddo
       else
         do ix1=ixfmax1,ixfmin1,-1
           w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
             w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
           +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
           -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
             block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
            /block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
         end do
       end if
       }
       {^ifthreed
       ixfmin1=ixomin1+1
       ixfmax1=ixomax1+1
       ixfmin2=ixomin2+1
       ixfmax2=ixomax2-1
       ixfmin3=ixomin3+1
       ixfmax3=ixomax3-1
       if(slab_uniform) then
         dx1x2=dxlevel(1)/dxlevel(2)
         dx1x3=dxlevel(1)/dxlevel(3)
         do ix1=ixfmax1,ixfmin1,-1
           w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
                     w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
             +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
                     w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
             +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
                     w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
         end do
       else
         do ix1=ixfmax1,ixfmin1,-1
           w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
          ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
           +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
           -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
             w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
             block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
           +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
           -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
            /block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
         end do
       end if
       }
     case(2)
       {^iftwod
       ixfmin1=ixomin1-1
       ixfmax1=ixomax1-1
       ixfmin2=ixomin2+1
       ixfmax2=ixomax2-1
       if(slab_uniform) then
         dx1x2=dxlevel(1)/dxlevel(2)
         do ix1=ixfmin1,ixfmax1
           w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
            -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
                    w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
         enddo
       else
         do ix1=ixfmin1,ixfmax1
           w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
             w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
           -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
           +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
             block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
            /block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
         end do
       end if
       }
       {^ifthreed
       ixfmin1=ixomin1-1
       ixfmax1=ixomax1-1
       ixfmin2=ixomin2+1
       ixfmax2=ixomax2-1
       ixfmin3=ixomin3+1
       ixfmax3=ixomax3-1
       if(slab_uniform) then
         dx1x2=dxlevel(1)/dxlevel(2)
         dx1x3=dxlevel(1)/dxlevel(3)
         do ix1=ixfmin1,ixfmax1
           w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
                     w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
             -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
                     w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
             -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
                     w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
         end do
       else
         do ix1=ixfmin1,ixfmax1
           w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
          ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
             block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
           -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
           +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
             w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
             block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
           -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
           +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
             w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
             block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
            /block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
             w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
         end do
       end if
       }
     case(3)
       {^iftwod
       ixfmin1=ixomin1+1
       ixfmax1=ixomax1-1
       ixfmin2=ixomin2+1
       ixfmax2=ixomax2+1
       if(slab_uniform) then
         dx2x1=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
       }
     case(4)
       {^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
       }
     {^ifthreed
     case(5)
       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
     case(6)
       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
     }
     case default
       call mpistop("Special boundary is not defined for this region")
    end select
 
  end subroutine fixdivb_boundary
 
  {^nooned
  subroutine mf_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
 
    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                      :: id
    integer, parameter           :: max_its      = 50
    double precision             :: residual_it(max_its), max_divb
    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_neumann
          mg%bc(n, mg_iphi)%bc_value = 0.0_dp
       case (bc_periodic)
          ! Nothing to do here
       case default
          write(*,*) "mf_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, &
            mf_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))
           ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
         end do
         call mf_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
         ! 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
 
    end do
 
    active = .true.
 
  end subroutine mf_clean_divb_multigrid
  }
 
  !> get electric field though averaging neighors to update faces in CT
  subroutine mf_update_faces_average(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
    use mod_global_parameters
    use mod_usr_methods
 
    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)
    double precision                   :: e_resi(ixi^s,sdim:3)
    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(mf_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(mf_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
              if(record_electric_field) s%we(ix^d,idir)=fe(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
      circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
      ! Time update
      bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
    end do
 
    end associate
 
  end subroutine mf_update_faces_average
 
  !> update faces using UCT contact mode by Gardiner and Stone 2005 JCP 205, 509
  subroutine mf_update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
    use mod_global_parameters
    use mod_usr_methods
 
    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)
    ! 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(ixi^s),erc(ixi^s)
    ! current on cell edges
    double precision                   :: jce(ixi^s,sdim:3)
    integer                            :: hxc^l,ixc^l,jxc^l,ixa^l,ixb^l
    integer                            :: idim1,idim2,idir,iwdim1,iwdim2
 
    associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
 
    ecc=0.d0
    ! Calculate electric field at cell centers
    do idim1=1,ndim; do idim2=1,ndim; do idir=sdim,3
      if(lvc(idim1,idim2,idir)==1)then
         ecc(ixi^s,idir)=ecc(ixi^s,idir)+wp(ixi^s,mag(idim1))*wp(ixi^s,mom(idim2))
      else if(lvc(idim1,idim2,idir)==-1) then
         ecc(ixi^s,idir)=ecc(ixi^s,idir)-wp(ixi^s,mag(idim1))*wp(ixi^s,mom(idim2))
      endif
    enddo; enddo; enddo
 
    ! if there is resistivity, get eta J
    if(mf_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,jce)
 
    ! 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)
      do idim2=1,ndim
        iwdim2 = mag(idim2)
        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
            jxc^l=ixc^l+kr(idim1,^d);
            hxc^l=ixc^l+kr(idim2,^d);
            ! average cell-face electric field to cell edges
            fe(ixc^s,idir)=quarter*&
            (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
            -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
 
            ! add slope in idim2 direction from equation (50)
            ixamin^d=ixcmin^d;
            ixamax^d=ixcmax^d+kr(idim1,^d);
            el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
            hxc^l=ixa^l+kr(idim2,^d);
            er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
            where(vnorm(ixc^s,idim1)>0.d0)
              elc(ixc^s)=el(ixc^s)
            else where(vnorm(ixc^s,idim1)<0.d0)
              elc(ixc^s)=el(jxc^s)
            else where
              elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
            end where
            hxc^l=ixc^l+kr(idim2,^d);
            where(vnorm(hxc^s,idim1)>0.d0)
              erc(ixc^s)=er(ixc^s)
            else where(vnorm(hxc^s,idim1)<0.d0)
              erc(ixc^s)=er(jxc^s)
            else where
              erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
            end where
            fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
 
            ! add slope in idim1 direction from equation (50)
            jxc^l=ixc^l+kr(idim2,^d);
            ixamin^d=ixcmin^d;
            ixamax^d=ixcmax^d+kr(idim2,^d);
            el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
            hxc^l=ixa^l+kr(idim1,^d);
            er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
            where(vnorm(ixc^s,idim2)>0.d0)
              elc(ixc^s)=el(ixc^s)
            else where(vnorm(ixc^s,idim2)<0.d0)
              elc(ixc^s)=el(jxc^s)
            else where
              elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
            end where
            hxc^l=ixc^l+kr(idim1,^d);
            where(vnorm(hxc^s,idim2)>0.d0)
              erc(ixc^s)=er(ixc^s)
            else where(vnorm(hxc^s,idim2)<0.d0)
              erc(ixc^s)=er(jxc^s)
            else where
              erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
            end where
            fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
 
            ! add current component of electric field at cell edges E=-vxB+eta J
            if(mf_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
 
            if(record_electric_field) s%we(ixc^s,idir)=fe(ixc^s,idir)
            ! times time step and edge length 
            fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
            if (.not.slab) then
              where(abs(x(ixc^s,r_)+half*dxlevel(r_))<1.0d-9)
                fe(ixc^s,idir)=zero
              end where
            end if
          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
        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 mf_update_faces_contact
 
  !> update faces
  subroutine mf_update_faces_hll(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
    use mod_global_parameters
    use mod_usr_methods
    use mod_constrained_transport
 
    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                   :: vtill(ixi^s,2)
    double precision                   :: vtilr(ixi^s,2)
    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)
    ! current on cell edges
    double precision                   :: jce(ixi^s,sdim:3)
    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(mf_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,jce)
 
    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.
      call reconstruct(ixi^l,ixc^l,idim2,bfacesct(ixi^s,idim1),&
               btill(ixi^s,idim1),btilr(ixi^s,idim1))
 
      call reconstruct(ixi^l,ixc^l,idim1,bfacesct(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 current component of electric field at cell edges E=-vxB+eta J
      if(mf_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
 
      if(record_electric_field) s%we(ixc^s,idir)=fe(ixc^s,idir)
      ! times time step and edge length 
      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
      bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
    end do
 
    end associate
  end subroutine mf_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+1;
          ixcmin^d=ixomin^d+kr(idir,^d)-2;
          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)
          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(mf_eta>zero)then
      jce(ixi^s,:)=jce(ixi^s,:)*mf_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)
      enddo
    end if
 
    end associate
  end subroutine get_resistive_electric_field
 
  !> calculate cell-center values from face-center values
  subroutine mf_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
   {!dir$ ivdep
    do ix^db=ixomin^db,ixomax^db\}
      {^ifthreed
      s%w(ix^d,mag(1))=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,mag(2))=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,mag(3))=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,mag(1))=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,mag(2))=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 mf_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 record_force_free_metrics()
    use mod_global_parameters
 
    double precision :: sum_jbb_ipe, sum_j_ipe, sum_l_ipe, f_i_ipe, volume_pe
    double precision :: sum_jbb, sum_j, sum_l, f_i, volume, cw_sin_theta
    integer :: iigrid, igrid, ix^d
    integer :: amode, istatus(mpi_status_size)
    integer, save :: fhmf
    logical :: patchwi(ixg^t)
    logical, save :: logmfopened=.false.
    character(len=800) :: filename,filehead
    character(len=800) :: line,datastr
 
    sum_jbb_ipe = 0.d0
    sum_j_ipe = 0.d0
    sum_l_ipe = 0.d0
    f_i_ipe = 0.d0
    volume_pe=0.d0
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      block=>ps(igrid)
      ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
      call mask_inner(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x,patchwi,volume_pe)
      sum_jbb_ipe = sum_jbb_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
        ps(igrid)%x,1,patchwi)
      sum_j_ipe   = sum_j_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
        ps(igrid)%x,2,patchwi)
      f_i_ipe=f_i_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
        ps(igrid)%x,3,patchwi)
      sum_l_ipe   = sum_l_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
        ps(igrid)%x,4,patchwi)
    end do
    call mpi_reduce(sum_jbb_ipe,sum_jbb,1,mpi_double_precision,&
       mpi_sum,0,icomm,ierrmpi)
    call mpi_reduce(sum_j_ipe,sum_j,1,mpi_double_precision,mpi_sum,0,&
       icomm,ierrmpi)
    call mpi_reduce(f_i_ipe,f_i,1,mpi_double_precision,mpi_sum,0,&
       icomm,ierrmpi)
    call mpi_reduce(sum_l_ipe,sum_l,1,mpi_double_precision,mpi_sum,0,&
       icomm,ierrmpi)
    call mpi_reduce(volume_pe,volume,1,mpi_double_precision,mpi_sum,0,&
       icomm,ierrmpi)
 
    if(mype==0) then
      ! current- and volume-weighted average of the sine of the angle
      ! between the magnetic field B and the current density J
      cw_sin_theta = sum_jbb/sum_j
      ! volume-weighted average of the absolute value of the fractional
      ! magnetic flux change
      f_i = f_i/volume
      sum_j=sum_j/volume
      sum_l=sum_l/volume
      if(.not.logmfopened) then
        ! generate filename
        write(filename,"(a,a)") trim(base_filename), "_mf_metrics.csv"
 
        ! Delete the log when not doing a restart run
        if (restart_from_file == undefined) then
           open(unit=20,file=trim(filename),status='replace')
           close(20, status='delete')
        end if
 
        amode=ior(mpi_mode_create,mpi_mode_wronly)
        amode=ior(amode,mpi_mode_append)
        call mpi_file_open(mpi_comm_self,filename,amode,mpi_info_null,fhmf,ierrmpi)
        logmfopened=.true.
        filehead="    it, time, CW_sin_theta, current, Lorenz_force, f_i"
        if (it == 0) then
          call mpi_file_write(fhmf,filehead,len_trim(filehead), &
                              mpi_character,istatus,ierrmpi)
          call mpi_file_write(fhmf,achar(10),1,mpi_character,istatus,ierrmpi)
        end if
      end if
      line=''
      write(datastr,'(i8,a)') it,','
      line=trim(line)//trim(datastr)
      write(datastr,'(es13.6,a)') global_time,','
      line=trim(line)//trim(datastr)
      write(datastr,'(es13.6,a)') cw_sin_theta,','
      line=trim(line)//trim(datastr)
      write(datastr,'(es13.6,a)') sum_j,','
      line=trim(line)//trim(datastr)
      write(datastr,'(es13.6,a)') sum_l,','
      line=trim(line)//trim(datastr)
      write(datastr,'(es13.6)') f_i
      line=trim(line)//trim(datastr)//new_line('A')
      call mpi_file_write(fhmf,line,len_trim(line),mpi_character,istatus,ierrmpi)
      if(.not.time_advance) then
        call mpi_file_close(fhmf,ierrmpi)
        logmfopened=.false.
      end if
    end if
 
  end subroutine record_force_free_metrics
 
  subroutine mask_inner(ixI^L,ixO^L,w,x,patchwi,volume_pe)
    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) :: volume_pe
    logical, intent(out)        :: patchwi(ixi^s)
 
    double precision            :: xo^l
    integer                     :: ix^d
 
    {xomin^d = xprobmin^d + 0.05d0*(xprobmax^d-xprobmin^d)\}
    {xomax^d = xprobmax^d - 0.05d0*(xprobmax^d-xprobmin^d)\}
    if(slab) then
      xomin^nd = xprobmin^nd
    else
      xomin1 = xprobmin1
    end if
 
    {do ix^db=ixomin^db,ixomax^db\}
        if({ x(ix^dd,^d) > xomin^d .and. x(ix^dd,^d) < xomax^d | .and. }) then
          patchwi(ix^d)=.true.
          volume_pe=volume_pe+block%dvolume(ix^d)
        else
          patchwi(ix^d)=.false.
        endif
    {end do\}
 
  end subroutine mask_inner
 
  function integral_grid_mf(ixI^L,ixO^L,w,x,iw,patchwi)
    use mod_global_parameters
 
    integer, intent(in)                :: ixi^l,ixo^l,iw
    double precision, intent(in)       :: x(ixi^s,1:ndim)
    double precision                   :: w(ixi^s,nw+nwauxio)
    logical, intent(in) :: patchwi(ixi^s)
 
    double precision, dimension(ixI^S,7-2*ndir:3) :: current
    double precision, dimension(ixI^S,1:ndir) :: bvec,qvec
    double precision :: integral_grid_mf,tmp(ixi^s),bm2
    integer :: ix^d,idirmin0,idirmin,idir,jdir,kdir
 
    integral_grid_mf=0.d0
    select case(iw)
     case(1)
      ! Sum(|JxB|/|B|*dvolume)
      bvec(ixo^s,:)=w(ixo^s,mag(:))
      call get_current(w,ixi^l,ixo^l,idirmin,current)
      ! calculate Lorentz force
      qvec(ixo^s,1:ndir)=zero
      do idir=1,ndir; do jdir=idirmin,3; do kdir=1,ndir
         if(lvc(idir,jdir,kdir)/=0)then
           if(lvc(idir,jdir,kdir)==1)then
             qvec(ixo^s,idir)=qvec(ixo^s,idir)+current(ixo^s,jdir)*bvec(ixo^s,kdir)
           else
             qvec(ixo^s,idir)=qvec(ixo^s,idir)-current(ixo^s,jdir)*bvec(ixo^s,kdir)
           end if
         end if
      end do; end do; end do
 
      {do ix^db=ixomin^db,ixomax^db\}
         if(patchwi(ix^d)) then
           bm2=sum(bvec(ix^d,:)**2)
           if(bm2/=0.d0) bm2=1.d0/bm2
           integral_grid_mf=integral_grid_mf+sqrt(sum(qvec(ix^d,:)**2)*&
                           bm2)*block%dvolume(ix^d)
         end if
      {end do\}
     case(2)
      ! Sum(|J|*dvolume)
      call get_current(w,ixi^l,ixo^l,idirmin,current)
      {do ix^db=ixomin^db,ixomax^db\}
         if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+sqrt(sum(current(ix^d,:)**2))*&
                           block%dvolume(ix^d)
      {end do\}
     case(3)
      ! f_i solenoidal property of B: (dvolume |div B|)/(dsurface |B|)
      ! Sum(f_i*dvolume)
      call get_normalized_divb(w,ixi^l,ixo^l,tmp)
      {do ix^db=ixomin^db,ixomax^db\}
         if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+tmp(ix^d)*block%dvolume(ix^d)
      {end do\}
     case(4)
      ! Sum(|JxB|*dvolume)
      bvec(ixo^s,:)=w(ixo^s,mag(:))
      call get_current(w,ixi^l,ixo^l,idirmin,current)
      ! calculate Lorentz force
      qvec(ixo^s,1:ndir)=zero
      do idir=1,ndir; do jdir=idirmin,3; do kdir=1,ndir
         if(lvc(idir,jdir,kdir)/=0)then
           if(lvc(idir,jdir,kdir)==1)then
             qvec(ixo^s,idir)=qvec(ixo^s,idir)+current(ixo^s,jdir)*bvec(ixo^s,kdir)
           else
             qvec(ixo^s,idir)=qvec(ixo^s,idir)-current(ixo^s,jdir)*bvec(ixo^s,kdir)
           end if
         end if
      end do; end do; end do
 
      {do ix^db=ixomin^db,ixomax^db\}
         if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+&
                            sqrt(sum(qvec(ix^d,:)**2))*block%dvolume(ix^d)
      {end do\}
    end select
    return
  end function integral_grid_mf
 
end module mod_mf_phys