mod_magnetofriction.t

Go to the documentation of this file.

!> module mod_magnetofriction.t
!> Purpose: use magnetofrictional method to relax 3D magnetic field to
!>          force-free field
!> 01.04.2016 developed by Chun Xia and Yang Guo
!> 04.10.2017 modulized by Chun Xia
!> Usage:
!> in amrvac.par:
!>   &methodlist
!>    time_stepper='onestep' ! time marching scheme, or 'twostep','threestep'
!>    flux_method=13*'cd4' ! or 'tvdlf', 'fd'
!>    limiter= 13*'koren' ! or 'vanleer','cada3','mp5' so on
!>   /
!>   &meshlist
!>    ditregrid=20 ! set iteration interval for adjusting AMR
!>   /
!>   &mhd_list
!>    mhd_magnetofriction=.true.
!>   /
!>   &mf_list
!>    mf_it_max=60000  ! set the maximum iteration number
!>    mf_ditsave=20000 ! set iteration interval for data output
!>    mf_cc=0.3     ! stability coefficient controls numerical stability
!>    mf_cy=0.2     ! frictional velocity coefficient
!>    mf_cdivb=0.01 ! divb cleaning coefficient controls diffusion speed of divb
!>   /
module mod_magnetofriction
  implicit none
 
  !> stability coefficient controls numerical stability
  double precision :: mf_cc
  !> frictional velocity coefficient
  double precision :: mf_cy, mf_cy_max
  !> divb cleaning coefficient controls diffusion speed of divb
  double precision :: mf_cdivb, mf_cdivb_max
  !> TVDLF dissipation coefficient controls the dissipation term
  double precision :: mf_tvdlfeps, mf_tvdlfeps_min
  !> time in magnetofriction process
  double precision :: tmf
  !> maximal speed for fd scheme
  double precision :: cmax_mype
  !> maximal speed for fd scheme
  double precision :: cmax_global
  !> maximal limit of magnetofrictional velocity in cm s^-1 (Pomoell 2019 A&A)
  double precision, public :: mf_vmax = 3.d6
 
  !> Index of the density (in the w array)
  integer, private, protected              :: rho_
 
  !> Indices of the momentum density
  integer, allocatable, private, protected :: mom(:)
 
  !> Indices of the magnetic field
  integer, allocatable, private, protected :: mag(:)
 
  integer :: mf_ditsave
  integer :: mf_it_max
  integer :: mf_it
  !> How to open the MF diagnostics CSV: auto, append, or replace.
  character(len=16) :: mf_log_mode
  !> Optional diagnostics filename; empty uses <base_filename>_mflog.csv.
  character(len=256) :: mf_log_filename
  !> Write the legacy, method-specific MF diagnostics history.
  logical :: mf_write_detailed_history
  logical :: mf_advance
  logical :: fix_conserve_at_step = .true.
 
  !> Whether this run continues an existing MF diagnostics series.
  logical :: mf_continue_run = .false.
 
contains
  !> Read this module"s parameters from a file
  subroutine mf_params_read(files)
    use mod_global_parameters, only: unitpar
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /mf_list/ mf_ditsave, mf_it_max, mf_it, mf_cc, mf_cy, mf_cy_max, &
                       mf_cdivb, mf_cdivb_max, mf_tvdlfeps, mf_tvdlfeps_min, &
                       mf_log_mode, mf_log_filename, mf_write_detailed_history
 
    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_params_read
 
  !> Initialize the module
  subroutine magnetofriction_init()
    use mod_global_parameters
    use mod_usr_methods, only: usr_before_main_loop
 
    rho_=iw_rho
    allocate(mom(ndir))
    mom=iw_mom
    allocate(mag(ndir))
    mag=iw_mag
 
    mf_it=0          ! set the initial iteration number
    mf_it_max=60000  ! set the maximum iteration number
    mf_ditsave=20000 ! set iteration interval for data output
    mf_cc=0.3d0      ! stability coefficient controls numerical stability
    mf_cy=0.2d0      ! frictional velocity coefficient. The default value is mf_cy itself.
    mf_cy_max=mf_cy  ! maximum of the frictional velocity coefficient
    mf_cdivb=0.01d0  ! divb cleaning coefficient controls diffusion speed of divb
    mf_cdivb_max=mf_cdivb ! maximum of the divb cleaning coefficient
    mf_tvdlfeps=1.d0 ! coefficient to control the TVDLF dissipation
    mf_tvdlfeps_min = mf_tvdlfeps ! minimum of the TVDLF dissipation coefficient
    mf_log_mode='auto'
    mf_log_filename=''
    mf_write_detailed_history=.false.
    ! get dimensionless maximal mf velocity limit
    mf_vmax=mf_vmax/unit_velocity
 
    call mf_params_read(par_files)
 
    if(mf_cy_max < mf_cy) mf_cy_max = mf_cy
    if(mf_cdivb_max < mf_cdivb) mf_cdivb_max = mf_cdivb
    if(mf_tvdlfeps_min > mf_tvdlfeps) mf_tvdlfeps_min = mf_tvdlfeps
 
    usr_before_main_loop => magnetofriction
 
  end subroutine magnetofriction_init
 
  subroutine magnetofriction
    use mod_global_parameters
    use mod_physics
    use mod_ghostcells_update
    use mod_input_output
    use mod_amr_grid, only: resettree
    use mod_comm_lib, only: mpistop
    {^ifthreed
    use mod_nlfff_diagnostics, only: nlfff_physical_metrics,&
       evaluate_nlfff_metrics_amrvac,write_nlfff_metrics_header,&
       write_nlfff_metrics_row
    }
 
    double precision :: dvolume(ixG^T),dsurface(ixG^T),dvone
    double precision :: dtfff,dtfff_pe,dtnew,dx^D
    double precision :: cwsin_theta_new,cwsin_theta_old
    double precision :: sum_jbb,sum_jbb_ipe,sum_j,sum_j_ipe,sum_l_ipe,sum_l
    double precision :: f_i_ipe,f_i,volumepe,volume,tmpt,time_in
    double precision, external :: integral_grid
    {^ifthreed
    type(nlfff_physical_metrics) :: physical_metrics
    }
    integer :: i,iigrid, igrid, idims,ix^D,hxM^LL,fhmf,tmpit,i^D
    integer :: common_metrics_unit
    logical :: patchwi(ixG^T), stagger_flag,common_metrics_exists
 
    ! not do fix conserve and getbc for staggered values if stagger is used
    stagger_flag=stagger_grid
    stagger_grid=.false.
 
    time_in=mpi_wtime()
    if(mype==0) write(*,*) 'Evolving to force-free field using magnetofricitonal method...'
    if(prolongprimitive) call mpistop('use prolongprimitive=.false. in MF module')
    mf_advance=.false.
    dtfff=1.d-2
    tmpt=global_time
    tmpit=it
    tmf=global_time
    i=mf_it
    common_metrics_unit=-1
    ! MF checkpoints store the MF iteration in the standard snapshot header.
    ! Recover it automatically unless the user explicitly supplied mf_it.
    if(i==0 .and. it>0) i=it
    select case(trim(mf_log_mode))
    case('auto')
      ! Potential-field snapshots normally have it=0, whereas MF checkpoints
      ! store their positive MF iteration in the standard snapshot header.
      mf_continue_run=(i>0)
    case('append')
      mf_continue_run=.true.
    case('replace')
      mf_continue_run=.false.
    case default
      call mpistop("mf_log_mode must be 'auto', 'append', or 'replace'")
    end select
    if(snapshotini==0 .and. i==0) then
      call saveamrfile(1)
      call saveamrfile(2)
    end if
    mf_advance=.true.
    ! point bc mpi datatype to partial type for magnetic 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
    ! create bc mpi datatype for ghostcells update
    call create_bc_mpi_datatype(mag(1),ndir)
    ! point bc mpi datatype to partial type for velocity field
    type_send_srl=>type_send_srl_p2
    type_recv_srl=>type_recv_srl_p2
    type_send_r=>type_send_r_p2
    type_recv_r=>type_recv_r_p2
    type_send_p=>type_send_p_p2
    type_recv_p=>type_recv_p_p2
    ! create bc mpi datatype for ghostcells update
    call create_bc_mpi_datatype(mom(1),ndir)
    ! convert conservative variables to primitive ones which are used during MF
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       call phys_to_primitive(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
    end do
    ! calculate magnetofrictional velocity
    call mf_velocity_update(dtfff)
    ! update velocity in ghost cells
    call getbc(tmf,0.d0,ps,mom(1),ndir)
    ! calculate initial values of metrics
    if(i==0) then
      call metrics
      call printlog_mf
    end if
    {^ifthreed
    common_metrics_exists=.false.
    if(mype==0) then
      inquire(file=trim(base_filename)//'_nlfff_metrics.csv',&
         exist=common_metrics_exists)
    end if
    call mpi_bcast(common_metrics_exists,1,mpi_logical,0,icomm,ierrmpi)
    if(mype==0) then
      if(mf_continue_run .and. common_metrics_exists) then
        open(newunit=common_metrics_unit,file=trim(base_filename)//&
           '_nlfff_metrics.csv',status='old',position='append',action='write')
      else
        open(newunit=common_metrics_unit,file=trim(base_filename)//&
           '_nlfff_metrics.csv',status='replace',action='write')
        call write_nlfff_metrics_header(common_metrics_unit)
      end if
    end if
    if(.not.mf_continue_run .or. .not.common_metrics_exists) then
      call evaluate_nlfff_metrics_amrvac(mag,physical_metrics)
      if(mype==0) call write_nlfff_metrics_row(common_metrics_unit,i,&
         physical_metrics)
    end if
    }
    ! magnetofrictional loops
    do
      ! calculate time step based on Cmax= Alfven speed + abs(frictional speed)
      dtfff_pe=bigdouble
      cmax_mype=zero
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        block=>ps(igrid)
        ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
        call getdtfff_courant(ps(igrid)%w,ps(igrid)%x,ixg^ll,ixm^ll,dtnew)
        dtfff_pe=min(dtfff_pe,dtnew)
      end do
      call mpi_allreduce(dtfff_pe,dtfff,1,mpi_double_precision,mpi_min, &
                         icomm,ierrmpi)
      call mpi_allreduce(cmax_mype,cmax_global,1,mpi_double_precision,&
                 mpi_max,icomm,ierrmpi)
 
      ! =======
      ! evolve
      ! =======
      call advectmf(1,ndim,tmf,dtfff)
 
      if(i>=10000)  then
        mf_tvdlfeps = 0.9998d0 * mf_tvdlfeps
        mf_tvdlfeps = max(mf_tvdlfeps_min,mf_tvdlfeps)
      end if
      if(i<=60000) then
        mf_cy=1.0001d0*mf_cy
        mf_cy = min(mf_cy_max,mf_cy)
      end if
      if(i>=100000) then
        mf_cdivb=1.0000018d0*mf_cdivb
        mf_cdivb=min(mf_cdivb_max,mf_cdivb)
      end if
 
      i=i+1
      tmf=tmf+dtfff
      if(mod(i,10)==0) then
        ! calculate metrics
        call metrics
        call printlog_mf
        {^ifthreed
        call evaluate_nlfff_metrics_amrvac(mag,physical_metrics)
        if(mype==0) call write_nlfff_metrics_row(common_metrics_unit,i,&
           physical_metrics)
        }
      end if
      if(mod(i,mf_ditsave)==0) then
        it=i
        global_time=tmf
        do iigrid=1,igridstail; igrid=igrids(iigrid);
          call phys_to_conserved(ixg^ll,ixg^ll,ps(igrid)%w,ps(igrid)%x)
        end do
        mf_advance=.false.
        call saveamrfile(1)
        call saveamrfile(2)
        do iigrid=1,igridstail; igrid=igrids(iigrid);
           call phys_to_primitive(ixg^ll,ixg^ll,ps(igrid)%w,ps(igrid)%x)
        end do
        mf_advance=.true.
        if(mype==0) then
          write(*,*) "itmf=",i
          write(*,*) '<CW sin theta>:',cwsin_theta_new
          write(*,*) '<f_i>:',f_i
          write(*,*) '----------------------------------------------------------'
        end if
      end if
      ! reconstruct AMR grid every 10 step
      if(mod(i,ditregrid)==0 .and. refine_max_level>1) call resettree
      if (i>=mf_it_max) then
        if(mod(i,10)/=0) then
          ! calculate metrics
          call metrics
          call printlog_mf
          {^ifthreed
          call evaluate_nlfff_metrics_amrvac(mag,physical_metrics)
          if(mype==0) call write_nlfff_metrics_row(common_metrics_unit,i,&
             physical_metrics)
          }
        end if
        if(mype==0) then
          write (*,*) 'Reach maximum iteration step!'
          write (*,*) 'The total iteration step is:', i
        end if
        exit
      end if
    enddo
    ! point bc mpi data type back to full type for MHD
    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
    bcphys=.true.
    ! set velocity back to zero and convert primitive variables back to conservative ones
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       ps(igrid)%w(ixg^t,mom(:))=zero
       call phys_to_conserved(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
    end do
    global_time=tmpt
    it=tmpit
    if (mype==0) then
      if(mf_write_detailed_history) call mpi_file_close(fhmf,ierrmpi)
      {^ifthreed close(common_metrics_unit)}
    end if
    mf_advance=.false.
    ! restore stagger_grid value
    stagger_grid=stagger_flag
    if(mype==0) write(*,*) 'Magnetofriction phase took : ',mpi_wtime()-time_in,' sec'
    contains
 
      subroutine metrics
 
        sum_jbb_ipe = 0.d0
        sum_j_ipe = 0.d0
        sum_l_ipe = 0.d0
        f_i_ipe = 0.d0
        volumepe=0.d0
        dsurface=zero
        do iigrid=1,igridstail; igrid=igrids(iigrid);
          block=>ps(igrid)
          dvolume(ixm^t)=block%dvolume(ixm^t)
          do idims=1,ndim
            hxm^ll=ixm^ll-kr(idims,^d);
            dsurface(ixm^t)=dsurface(ixm^t)+block%surfaceC(hxm^t,idims)
          end do
          ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
          call mask_inner(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
          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_allreduce(sum_jbb_ipe,sum_jbb,1,mpi_double_precision,&
           mpi_sum,icomm,ierrmpi)
        call mpi_allreduce(sum_j_ipe,sum_j,1,mpi_double_precision,mpi_sum,&
           icomm,ierrmpi)
        call mpi_allreduce(f_i_ipe,f_i,1,mpi_double_precision,mpi_sum,&
           icomm,ierrmpi)
        call mpi_allreduce(sum_l_ipe,sum_l,1,mpi_double_precision,mpi_sum,&
           icomm,ierrmpi)
        call mpi_allreduce(volumepe,volume,1,mpi_double_precision,mpi_sum,&
           icomm,ierrmpi)
        ! current- and volume-weighted average of the sine of the angle
        ! between the magnetic field B and the current density J
        cwsin_theta_new = 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
      end subroutine metrics
 
      subroutine mask_inner(ixI^L,ixO^L,w,x)
 
        integer, intent(in)         :: ixI^L,ixO^L
        double precision, intent(in):: w(ixI^S,nw),x(ixI^S,1:ndim)
        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.
              volumepe=volumepe+dvolume(ix^d)
            else
              patchwi(ix^d)=.false.
            endif
        {end do\}
 
      end subroutine mask_inner
 
      subroutine printlog_mf
        integer :: amode, status(MPI_STATUS_SIZE)
        integer :: truncate_unit, ios
        character(len=800) :: filename,filehead
        character(len=2048) :: line,datastr
        logical, save :: logmfopened=.false.
        logical :: logfile_exists
 
        if(.not.mf_write_detailed_history) return
 
        if(mype==0) then
          if(.not.logmfopened) then
            ! generate filename
            if(len_trim(mf_log_filename)>0) then
              filename=trim(mf_log_filename)
            else
              write(filename,"(a,a)") trim(base_filename), "_mflog.csv"
            end if
 
            inquire(file=trim(filename),exist=logfile_exists)
            if(.not.mf_continue_run) then
              ! A potential-field (or otherwise unrelated) restart begins a
              ! new MF relaxation, so discard diagnostics from an older run.
              open(newunit=truncate_unit,file=trim(filename),status='replace', &
                   action='write',iostat=ios)
              if(ios/=0) call mpistop('Unable to replace magnetofriction log file')
              close(truncate_unit)
              logfile_exists=.false.
            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.
            ! Preserve a single header when continuing an existing MF run.
            if(.not.logfile_exists) then
              filehead="  itmf,  dt,  <f_i>,  <CW sin theta>,  <Current>,  <Lorenz force>"
              call mpi_file_write(fhmf,filehead,len_trim(filehead), &
                                  mpi_character,status,ierrmpi)
              call mpi_file_write(fhmf,achar(10),1,mpi_character,status,ierrmpi)
            end if
          end if
          line=''
          write(datastr,'(i6,a)') i,','
          line=trim(line)//trim(datastr)
          write(datastr,'(es13.6,a)') dtfff,','
          line=trim(line)//trim(datastr)
          write(datastr,'(es13.6,a)') f_i,','
          line=trim(line)//trim(datastr)
          write(datastr,'(es13.6,a)') cwsin_theta_new,','
          line=trim(line)//trim(datastr)
          write(datastr,'(es13.6,a)') sum_j,','
          line=trim(line)//trim(datastr)
          write(datastr,'(es13.6)') sum_l
          line=trim(line)//trim(datastr)//new_line('A')
          call mpi_file_write(fhmf,line,len_trim(line),mpi_character,status,ierrmpi)
        end if
 
      end subroutine printlog_mf
 
      function integral_grid_mf(ixI^L,ixO^L,w,x,iw,patchwi)
        use mod_geometry
 
        integer, intent(in)                :: ixI^L,ixO^L,iw
        double precision, intent(in)       :: x(ixI^S,1:ndim)
        double precision, intent(in)       :: w(ixI^S,nw+nwauxio)
        logical, intent(in) :: patchwi(ixI^S)
 
        double precision, dimension(ixI^S,1:ndir) :: bvec,qvec,current
        double precision :: integral_grid_mf,tmp(ixI^S),b_mag(ixI^S)
        integer :: ix^D,i,idirmin,idir,jdir,kdir
 
        integral_grid_mf=0.d0
        select case(iw)
         case(1)
          ! Sum(dvolume*|JxB|/|B|)
          if(b0field) then
            bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
          else
            bvec(ixi^s,:)=w(ixi^s,mag(:))
          endif
          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=1,ndir; do kdir=idirmin,3
             if(lvc(idir,jdir,kdir)/=0)then
                tmp(ixo^s)=current(ixo^s,jdir)*bvec(ixo^s,kdir)
                if(lvc(idir,jdir,kdir)==1)then
                   qvec(ixo^s,idir)=qvec(ixo^s,idir)+tmp(ixo^s)
                else
                   qvec(ixo^s,idir)=qvec(ixo^s,idir)-tmp(ixo^s)
                endif
             endif
          enddo; enddo; enddo
 
          {do ix^db=ixomin^db,ixomax^db\}
             if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+sqrt(sum(qvec(ix^d,:)**2)/&
                               sum(bvec(ix^d,:)**2))*dvolume(ix^d)
          {end do\}
         case(2)
          ! Sum(dvolume*|J|)
          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))*&
                               dvolume(ix^d)
          {end do\}
         case(3)
          ! f_i solenoidal property of B: (dvolume |div B|)/(dsurface |B|)
          ! Sum(dvolume*f_i)
          if(b0field) then
            bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
          else
            bvec(ixi^s,:)=w(ixi^s,mag(:))
          endif
          call divvector(bvec,ixi^l,ixo^l,tmp)
          {do ix^db=ixomin^db,ixomax^db\}
             if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+abs(tmp(ix^d))*&
                   dvolume(ix^d)**2/sqrt(sum(bvec(ix^d,:)**2))/dsurface(ix^d)
          {end do\}
         case(4)
          ! Sum(|JxB|)
          if(b0field) then
            bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
          else
            bvec(ixi^s,:)=w(ixi^s,mag(:))
          endif
          call curlvector(bvec,ixi^l,ixo^l,current,idirmin,1,ndir)
          ! calculate Lorentz force
          qvec(ixo^s,1:ndir)=zero
          do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
             if(lvc(idir,jdir,kdir)/=0)then
                tmp(ixo^s)=current(ixo^s,jdir)*bvec(ixo^s,kdir)
                if(lvc(idir,jdir,kdir)==1)then
                   qvec(ixo^s,idir)=qvec(ixo^s,idir)+tmp(ixo^s)
                else
                   qvec(ixo^s,idir)=qvec(ixo^s,idir)-tmp(ixo^s)
                endif
             endif
          enddo; enddo; enddo
 
          {do ix^db=ixomin^db,ixomax^db\}
             if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+sqrt(sum(qvec(ix^d,:)**2))*dvolume(ix^d)
          {end do\}
        end select
        return
      end function integral_grid_mf
 
  end subroutine magnetofriction
 
  subroutine mf_velocity_update(dtfff)
    use mod_global_parameters
 
    double precision, intent(in) :: dtfff
    double precision :: vhatmax,vhatmax_pe,vhatmaxgrid
    integer :: i,iigrid, igrid
 
    vhatmax_pe=smalldouble
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      block=>ps(igrid)
      ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
      call vhat(ps(igrid)%w,ps(igrid)%x,ixg^ll,ixm^ll,vhatmaxgrid)
      vhatmax_pe=max(vhatmax_pe,vhatmaxgrid)
    end do
    call mpi_allreduce(vhatmax_pe,vhatmax,1,mpi_double_precision,mpi_max, &
                           icomm,ierrmpi)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      block=>ps(igrid)
      ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
      ! calculate frictional velocity
      call frictional_velocity(ps(igrid)%w,ps(igrid)%x,ixg^ll,ixm^ll,vhatmax,dtfff)
    end do
 
  end subroutine mf_velocity_update
 
  subroutine vhat(w,x,ixI^L,ixO^L,vhatmaxgrid)
    ! Calculate v_hat
    use mod_global_parameters
 
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(inout)  :: w(ixI^S,nw)
    double precision, intent(in)  :: x(ixI^S,1:ndim)
    double precision, intent(out) :: vhatmaxgrid
 
    double precision :: current(ixI^S,7-2*ndir:3),tmp(ixI^S),dxhm
    double precision :: dxhms(ixO^S)
    integer :: idirmin,idir,jdir,kdir
 
    call get_current(w,ixi^l,ixo^l,idirmin,current)
    w(ixi^s,mom(:))=0.d0
    ! calculate Lorentz force
    do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
       if(lvc(idir,jdir,kdir)/=0)then
          if(b0field) then
            tmp(ixo^s)=current(ixo^s,jdir)*(w(ixo^s,mag(kdir))+block%b0(ixo^s,kdir,0))
          else
            tmp(ixo^s)=current(ixo^s,jdir)*w(ixo^s,mag(kdir))
          endif
          if(lvc(idir,jdir,kdir)==1)then
             w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+tmp(ixo^s)
          else
             w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-tmp(ixo^s)
          endif
       endif
    enddo; enddo; enddo
 
    ! 1/B**2
    if(b0field) then
      tmp(ixo^s)=1.d0/(sum((w(ixo^s,mag(:))+block%b0(ixo^s,:,0))**2,dim=ndim+1)+smalldouble)
    else
      tmp(ixo^s)=1.d0/(sum(w(ixo^s,mag(:))**2,dim=ndim+1)+smalldouble)
    endif
 
    if(slab_uniform) then
      dxhm=dble(ndim)/(^d&1.0d0/dxlevel(^d)+)
      do idir=1,ndir
        w(ixo^s,mom(idir))=dxhm*w(ixo^s,mom(idir))*tmp(ixo^s)
      end do
    else
      dxhms(ixo^s)=dble(ndim)/sum(1.d0/block%dx(ixo^s,:),dim=ndim+1)
      do idir=1,ndir
        w(ixo^s,mom(idir))=dxhms(ixo^s)*w(ixo^s,mom(idir))*tmp(ixo^s)
      end do
    end if
    vhatmaxgrid=maxval(sqrt(sum(w(ixo^s,mom(:))**2,dim=ndim+1)))
 
  end subroutine vhat
 
  subroutine frictional_velocity(w,x,ixI^L,ixO^L,qvmax,qdt)
    use mod_global_parameters
 
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: x(ixI^S,1:ndim),qdt,qvmax
    double precision, intent(inout) :: w(ixI^S,1:nw)
 
    double precision :: dxhm,disbd(6),bfzone^D
    double precision :: dxhms(ixO^S)
    integer :: ix^D, idir
    logical :: buffer
 
    if(slab_uniform) then
      dxhm=dble(ndim)/(^d&1.0d0/dxlevel(^d)+)
      dxhm=mf_cc*mf_cy/qvmax*dxhm/qdt
      ! dxhm=mf_cc*mf_cy/qvmax
      w(ixo^s,mom(:))=w(ixo^s,mom(:))*dxhm
    else
      dxhms(ixo^s)=dble(ndim)/sum(1.d0/block%dx(ixo^s,:),dim=ndim+1)
      dxhms(ixo^s)=mf_cc*mf_cy/qvmax*dxhms(ixo^s)/qdt
      ! dxhm=mf_cc*mf_cy/qvmax
      do idir=1,ndir
        w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*dxhms(ixo^s)
      end do
    end if
 
{^ifthreed
    bfzone1=0.05d0*(xprobmax1-xprobmin1)
    bfzone2=0.05d0*(xprobmax2-xprobmin2)
    bfzone3=0.05d0*(xprobmax3-xprobmin3)
    {do ix^db=ixomin^db,ixomax^db\}
       disbd(1)=x(ix^d,1)-xprobmin1
       disbd(2)=xprobmax1-x(ix^d,1)
       disbd(3)=x(ix^d,2)-xprobmin2
       disbd(4)=xprobmax2-x(ix^d,2)
       disbd(5)=x(ix^d,3)-xprobmin1
       disbd(6)=xprobmax3-x(ix^d,3)
 
       if(slab) then
         if(disbd(1)<bfzone1) then
           w(ix^d,mom(:))=(1.d0-((bfzone1-disbd(1))/bfzone1)**2)*w(ix^d,mom(:))
         endif
       else
         if(disbd(5)<bfzone3) then
           w(ix^d,mom(:))=(1.d0-((bfzone3-disbd(5))/bfzone3)**2)*w(ix^d,mom(:))
         endif
       end if
       if(disbd(2)<bfzone1) then
         w(ix^d,mom(:))=(1.d0-((bfzone1-disbd(2))/bfzone1)**2)*w(ix^d,mom(:))
       endif
       if(disbd(3)<bfzone2) then
         w(ix^d,mom(:))=(1.d0-((bfzone2-disbd(3))/bfzone2)**2)*w(ix^d,mom(:))
       endif
       if(disbd(4)<bfzone2) then
         w(ix^d,mom(:))=(1.d0-((bfzone2-disbd(4))/bfzone2)**2)*w(ix^d,mom(:))
       endif
       if(disbd(6)<bfzone3) then
         w(ix^d,mom(:))=(1.d0-((bfzone3-disbd(6))/bfzone3)**2)*w(ix^d,mom(:))
       endif
    {end do\}
}
 
    ! saturate mf velocity at mf_vmax
    dxhms(ixo^s)=sqrt(sum(w(ixo^s,mom(:))**2,dim=ndim+1))/mf_vmax+1.d-12
    dxhms(ixo^s)=dtanh(dxhms(ixo^s))/dxhms(ixo^s)
    do idir=1,ndir
      w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*dxhms(ixo^s)
    end do
  end subroutine frictional_velocity
 
  subroutine advectmf(idim^LIM,qt,qdt)
    !  integrate all grids by one step of its delta(global_time)
    ! This subroutine is in VAC terminology equivalent to
    ! `advect' (with the difference that it will `advect' all grids)
    use mod_global_parameters
    use mod_fix_conserve
    use mod_comm_lib, only: mpistop
 
    integer, intent(in) :: idim^LIM
    double precision, intent(in) :: qt, qdt
 
    integer :: iigrid, igrid
 
    call init_comm_fix_conserve(idim^lim,ndir)
    fix_conserve_at_step = mf_advance .and. levmax>levmin
 
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      ps1(igrid)%w=ps(igrid)%w
    end do
 
    istep=0
 
    select case (t_stepper)
     case (onestep)
       call advect1mf(flux_method,qdt,one,idim^lim,qt,ps1,qt,ps)
     case (twostep)
       ! predictor step
       fix_conserve_at_step = .false.
       call advect1mf(typepred1,qdt,half,idim^lim,qt,ps,qt,ps1)
       ! corrector step
       fix_conserve_at_step = mf_advance .and. levmax>levmin
       call advect1mf(flux_method,qdt,one,idim^lim,qt+half*qdt,ps1,qt,ps)
     case (threestep)
       ! three step Runge-Kutta in accordance with Gottlieb & Shu 1998
       call advect1mf(flux_method,qdt,one,idim^lim,qt,ps,qt,ps1)
 
       do iigrid=1,igridstail; igrid=igrids(iigrid);
          ps2(igrid)%w(ixg^t,1:nwflux)=0.75d0*ps(igrid)%w(ixg^t,1:nwflux)+0.25d0*&
               ps1(igrid)%w(ixg^t,1:nwflux)
          if (nw>nwflux) ps2(igrid)%w(ixg^t,nwflux+1:nw) = &
               ps(igrid)%w(ixg^t,nwflux+1:nw)
       end do
 
       call advect1mf(flux_method,qdt,0.25d0,idim^lim,qt+qdt,ps1,qt+dt*0.25d0,ps2)
 
       do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
          ps(igrid)%w(ixg^t,1:nwflux)=1.0d0/3.0d0*ps(igrid)%w(ixg^t,1:nwflux)+&
               2.0d0/3.0d0*ps2(igrid)%w(ixg^t,1:nwflux)
       end do
 
       call advect1mf(flux_method,qdt,2.0d0/3.0d0,idim^lim,qt+qdt/2.0d0,ps2,qt+qdt/3.0d0,ps)
     case default
       call mpistop("unkown time_stepper in advectmf")
    end select
 
  end subroutine advectmf
 
  subroutine advect1mf(method,dtin,dtfactor,idim^LIM,qtC,psa,qt,psb)
    ! Integrate all grids by one partial step
    ! This subroutine is equivalent to VAC's `advect1', but does
    ! the advection for all grids
    use mod_global_parameters
    use mod_ghostcells_update
    use mod_fix_conserve
 
    integer, intent(in) :: idim^LIM
    type(state) :: psa(max_blocks)! Compute fluxes based on this state
    type(state) :: psb(max_blocks) ! update on this state
    double precision, intent(in) :: dtin,dtfactor, qtC, qt
    integer, intent(in) :: method(nlevelshi)
 
    double precision :: qdt
    integer :: iigrid, igrid, level, i^D
    logical :: setigrid
 
    istep=istep+1
 
    ! loop over all grids to arrive at equivalent
    qdt=dtfactor*dtin
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       block=>ps(igrid)
       level=node(plevel_,igrid)
 
       call process1_gridmf(method(level),igrid,qdt,ixg^ll,idim^lim,qtc,&
                       psa(igrid)%w,qt,psb(igrid)%w)
    end do
 
    ! opedit: Send flux for all grids, expects sends for all
    ! nsend_fc(^D), set in connectivity.t.
 
    if (fix_conserve_at_step) then
      call recvflux(idim^lim)
      call sendflux(idim^lim)
      call fix_conserve(psb,idim^lim,mag(1),ndir)
    end if
    ! point bc mpi datatype to partial type for magnetic 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
    ! update B in ghost cells
    call getbc(qt+qdt,qdt,psb,mag(1),ndir)
    ! calculate magnetofrictional velocity
    call mf_velocity_update(qdt)
    ! point bc mpi datatype to partial type for velocity field
    type_send_srl=>type_send_srl_p2
    type_recv_srl=>type_recv_srl_p2
    type_send_r=>type_send_r_p2
    type_recv_r=>type_recv_r_p2
    type_send_p=>type_send_p_p2
    type_recv_p=>type_recv_p_p2
    ! update magnetofrictional velocity in ghost cells
    call getbc(qt+qdt,qdt,psb,mom(1),ndir)
 
  end subroutine advect1mf
 
  subroutine process1_gridmf(method,igrid,qdt,ixG^L,idim^LIM,qtC,wCT,qt,w)
    ! This subroutine is equivalent to VAC's `advect1' for one grid
    use mod_global_parameters
    use mod_fix_conserve
    use mod_comm_lib, only: mpistop
 
    integer, intent(in) :: method
    integer, intent(in) :: igrid, ixG^L, idim^LIM
    double precision, intent(in) :: qdt, qtC, qt
    double precision :: wCT(ixG^S,1:nw), w(ixG^S,1:nw)
    double precision :: dx^D, fC(ixG^S,1:ndir,1:ndim)
    integer :: ixO^L
 
    dx^d=rnode(rpdx^d_,igrid);
    ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
    fc=0.d0
 
    ixo^l=ixg^l^lsubnghostcells;
    select case (method)
     case (fs_cd4)
       !================================
       ! 4th order central difference
       !================================
       call centdiff4mf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,fc,dx^d,ps(igrid)%x)
     case (fs_tvdlf)
       !================================
       ! TVDLF
       !================================
       call tvdlfmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,fc,dx^d,ps(igrid)%x)
     case (fs_hancock)
       ! hancock predict (first) step for twostep tvdlf and tvdmu scheme
       call hancockmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,dx^d,ps(igrid)%x)
     case (fs_fd)
       !================================
       ! finite difference
       !================================
       call fdmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,fc,dx^d,ps(igrid)%x)
    case default
       call mpistop("unknown flux scheme in advect1_gridmf")
    end select
 
    if (fix_conserve_at_step) then
      call store_flux(igrid,fc,idim^lim,ndir)
    end if
 
  end subroutine process1_gridmf
 
  subroutine upwindlrmf(ixI^L,ixL^L,ixR^L,idim,w,wCT,wLC,wRC,x)
    ! Determine the upwinded wLC(ixL) and wRC(ixR) from w.
    ! the wCT is only used when PPM is exploited.
    use mod_global_parameters
    use mod_limiter
 
    integer, intent(in) :: ixI^L, ixL^L, ixR^L, idim
    double precision, dimension(ixI^S,1:nw) :: w, wCT
    double precision, dimension(ixI^S,1:nw) :: wLC, wRC
    double precision, dimension(ixI^S,1:ndim) :: x
 
    double precision   :: ldw(ixI^S), rdw(ixI^S), dwC(ixI^S)
    integer :: jxR^L, ixC^L, jxC^L, iw
 
    if (type_limiter(block%level) == limiter_mp5) then
       call mp5limiter(ixi^l,ixl^l,idim,w,wlc,wrc)
    else if (type_limiter(block%level) == limiter_ppm) then
       call ppmlimiter(ixi^l,ixm^ll,idim,w,wct,wlc,wrc)
    else
       jxr^l=ixr^l+kr(idim,^d);
       ixcmax^d=jxrmax^d; ixcmin^d=ixlmin^d-kr(idim,^d);
       jxc^l=ixc^l+kr(idim,^d);
 
       do iw=1,nwflux
          if (loglimit(iw)) then
             w(ixcmin^d:jxcmax^d,iw)=dlog10(w(ixcmin^d:jxcmax^d,iw))
             wlc(ixl^s,iw)=dlog10(wlc(ixl^s,iw))
             wrc(ixr^s,iw)=dlog10(wrc(ixr^s,iw))
          end if
 
          dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
 
          ! limit flux from left and/or right
          call dwlimiter2(dwc,ixi^l,ixc^l,idim,type_limiter(block%level),ldw,rdw)
          wlc(ixl^s,iw)=wlc(ixl^s,iw)+half*ldw(ixl^s)
          wrc(ixr^s,iw)=wrc(ixr^s,iw)-half*rdw(jxr^s)
 
          if (loglimit(iw)) then
             w(ixcmin^d:jxcmax^d,iw)=10.0d0**w(ixcmin^d:jxcmax^d,iw)
             wlc(ixl^s,iw)=10.0d0**wlc(ixl^s,iw)
             wrc(ixr^s,iw)=10.0d0**wrc(ixr^s,iw)
          end if
       end do
 
    endif
 
  end subroutine upwindlrmf
 
  subroutine getfluxmf(w,x,ixI^L,ixO^L,idir,idim,f)
    ! Calculate lux f_idim[idir] within ixO^L.
    use mod_global_parameters
 
    integer, intent(in)             :: ixI^L, ixO^L, idir, idim
    double precision, intent(in)    :: w(ixI^S,nw)
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    double precision,intent(out)    :: f(ixI^S)
 
    ! compute flux of magnetic field
    ! f_i[b_k]=v_i*b_k-v_k*b_i
    if (idim==idir) then
      f(ixo^s)=zero
    else
      f(ixo^s)=w(ixo^s,mom(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,mom(idir))
      if (b0field) then
        f(ixo^s)=f(ixo^s)&
              +w(ixo^s,mom(idim))*block%B0(ixo^s,idir,idim)&
              -w(ixo^s,mom(idir))*block%B0(ixo^s,idim,idim)
      end if
    end if
 
  end subroutine getfluxmf
 
  subroutine tvdlfmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,fC,dx^D,x)
    ! method=='tvdlf'  --> 2nd order TVD-Lax-Friedrich scheme.
    ! method=='tvdlf1' --> 1st order TVD-Lax-Friedrich scheme.
    use mod_global_parameters
    use mod_comm_lib, only: mpistop
 
    double precision, intent(in)                         :: qdt, qtC, qt, dx^D
    integer, intent(in)                                  :: ixI^L, ixO^L, idim^LIM
    double precision, dimension(ixI^S,1:ndim), intent(in) ::  x
    double precision, dimension(ixI^S,1:nw)               :: wCT, wnew
    double precision, dimension(ixI^S,1:ndir,1:ndim)        :: fC
 
    double precision, dimension(ixI^S,1:nw) :: wLC, wRC, wmean
    double precision, dimension(ixI^S)      :: fLC, fRC
    double precision, dimension(ixI^S)      :: cmaxC
    double precision :: dxinv(1:ndim), inv_volume(ixO^S)
    integer :: idims, idir, ix^L, hxO^L, ixC^L, ixCR^L, jxC^L, kxC^L, kxR^L
 
    ! The flux calculation contracts by one in the idim direction it is applied.
    ! The limiter contracts the same directions by one more, so expand ixO by 2.
    ix^l=ixo^l;
    do idims= idim^lim
       ix^l=ix^l^ladd2*kr(idims,^d);
    end do
    if (ixi^l^ltix^l|.or.|.or.) &
       call mpistop("Error in tvdlfmf: Nonconforming input limits")
 
    ^d&dxinv(^d)=-qdt/dx^d;
    fc=0.d0
    do idims= idim^lim
       b0i=idims
 
       hxo^l=ixo^l-kr(idims,^d);
       ! ixC is centered index in the idim direction from ixOmin-1/2 to ixOmax+1/2
       ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
       ! Calculate wRC=uR_{j+1/2} and wLC=uL_j+1/2
       jxc^l=ixc^l+kr(idims,^d);
       kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
       kxr^l=kxc^l+kr(idims,^d);
       ixcr^l=ixc^l;
 
       wrc(kxc^s,1:nwflux)=wct(kxr^s,1:nwflux)
       wlc(kxc^s,1:nwflux)=wct(kxc^s,1:nwflux)
 
       call upwindlrmf(ixi^l,ixcr^l,ixcr^l,idims,wct,wct,wlc,wrc,x)
 
       ! For the high order Lax-Friedrich TVDLF scheme the limiter is based on
       ! the maximum eigenvalue, it is calculated in advance.
       ! determine mean state and store in wLC
       wmean=0.5d0*(wlc+wrc)
       call getcmaxfff(wmean,ixg^ll,ixc^l,idims,cmaxc)
 
       ! Calculate fLC=f(uL_j+1/2) and fRC=f(uR_j+1/2) for each idir
       do idir=1,ndir
          call getfluxmf(wlc,x,ixg^ll,ixc^l,idir,idims,flc)
          call getfluxmf(wrc,x,ixg^ll,ixc^l,idir,idims,frc)
          ! To save memory we use fLC to store (F_L+F_R)/2=half*(fLC+fRC)
          flc(ixc^s)=half*(flc(ixc^s)+frc(ixc^s))
 
          ! Add TVDLF dissipation to the flux
          if (idir==idims) then
            flc(ixc^s)=flc(ixc^s)-mf_tvdlfeps*tvdlfeps*cmaxc(ixc^s)*half*(wrc(ixc^s,mag(idir))-wlc(ixc^s,mag(idir)))
          end if
          if (slab_uniform) then
            fc(ixc^s,idir,idims)=flc(ixc^s)
          else
            fc(ixc^s,idir,idims)=block%surfaceC(ixc^s,idims)*flc(ixc^s)
          end if
 
       end do ! Next idir
    end do ! Next idims
    b0i=0
 
    !Now update the state:
    do idims= idim^lim
       hxo^l=ixo^l-kr(idims,^d);
       ! Multiply the fluxes by -dt/dx since Flux fixing expects this
       if (slab_uniform) then
          fc(ixi^s,:,idims)=dxinv(idims)*fc(ixi^s,:,idims)
          wnew(ixo^s,mag(:))=wnew(ixo^s,mag(:)) &
               + (fc(ixo^s,:,idims)-fc(hxo^s,:,idims))
       else
          inv_volume = 1.0d0/block%dvolume(ixo^s)
          fc(ixi^s,:,idims)=-qdt*fc(ixi^s,:,idims)
 
          do idir = 1, ndir
             wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir)) + (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims)) * &
                  inv_volume
          end do
       end if
 
    end do ! Next idims
 
    if (.not.slab) call addgeometrymf(qdt,ixi^l,ixo^l,wct,wnew,x)
    call divbclean(qdt,ixi^l,ixo^l,wct,wnew,x)
 
  end subroutine tvdlfmf
 
  subroutine hancockmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,dx^D,x)
    ! The non-conservative Hancock predictor for TVDLFmf
    ! on entry:
    ! input available on ixI^L=ixG^L asks for output on ixO^L=ixG^L^LSUBnghostcells
    ! one entry: (predictor): wCT -- w_n        wnew -- w_n   qdt=dt/2
    ! on exit :  (predictor): wCT -- w_n        wnew -- w_n+1/2
    use mod_global_parameters
    use mod_comm_lib, only: mpistop
 
    integer, intent(in) :: ixI^L, ixO^L, idim^LIM
    double precision, intent(in) :: qdt, qtC, qt, dx^D, x(ixI^S,1:ndim)
    double precision, intent(inout) :: wCT(ixI^S,1:nw), wnew(ixI^S,1:nw)
 
    double precision, dimension(ixI^S,1:nw) :: wLC, wRC
    double precision, dimension(ixI^S) :: fLC, fRC
    double precision :: dxinv(1:ndim)
    integer :: idims, idir, ix^L, hxO^L, ixtest^L
 
    ! Expand limits in each idims direction in which fluxes are added
    ix^l=ixo^l;
    do idims= idim^lim
       ix^l=ix^l^laddkr(idims,^d);
    end do
    if (ixi^l^ltix^l|.or.|.or.) &
       call mpistop("Error in Hancockmf: Nonconforming input limits")
 
    ^d&dxinv(^d)=-qdt/dx^d;
    do idims= idim^lim
       b0i=idims
       ! Calculate w_j+g_j/2 and w_j-g_j/2
       ! First copy all variables, then upwind wLC and wRC.
       ! wLC is to the left of ixO, wRC is to the right of wCT.
       hxo^l=ixo^l-kr(idims,^d);
 
       wrc(hxo^s,1:nwflux)=wct(ixo^s,1:nwflux)
       wlc(ixo^s,1:nwflux)=wct(ixo^s,1:nwflux)
 
       call upwindlrmf(ixi^l,ixo^l,hxo^l,idims,wct,wct,wlc,wrc,x)
 
       ! Advect mag(idir)
       do idir=1,ndir
          ! Calculate the fLC and fRC fluxes
          call getfluxmf(wrc,x,ixi^l,hxo^l,idir,idims,frc)
          call getfluxmf(wlc,x,ixi^l,ixo^l,idir,idims,flc)
 
          if (slab_uniform) then
             wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+dxinv(idims)* &
                              (flc(ixo^s)-frc(hxo^s))
          else
             wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))-qdt/block%dvolume(ixo^s) &
                   *(block%surfaceC(ixo^s,idims)*flc(ixo^s) &
                    -block%surfaceC(hxo^s,idims)*frc(hxo^s))
          end if
       end do
    end do ! next idims
    b0i=0
 
    if (.not.slab) call addgeometrymf(qdt,ixi^l,ixo^l,wct,wnew,x)
 
  end subroutine hancockmf
 
  subroutine fdmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,fC,dx^D,x)
    use mod_global_parameters
    double precision, intent(in)                                     :: qdt, qtC, qt, dx^D
    integer, intent(in)                                              :: ixI^L, ixO^L, idim^LIM
    double precision, dimension(ixI^S,1:ndim), intent(in)            :: x
 
    double precision, dimension(ixI^S,1:nw), intent(inout)           :: wCT, wnew
    double precision, dimension(ixI^S,1:ndir,1:ndim), intent(out)  :: fC
 
    double precision, dimension(ixI^S)                               :: fCT
    double precision, dimension(ixI^S,1:nw)                          :: fm, fp, fmR, fpL
    double precision, dimension(ixI^S)                               :: v
    double precision                                                 :: dxinv(1:ndim)
    integer                                                          :: idims, idir, ixC^L, ix^L, hxO^L, ixCR^L
 
    ^d&dxinv(^d)=-qdt/dx^d;
    do idims= idim^lim
 
       ! Get fluxes for the whole grid (mesh+nghostcells)
       {^d& ixcmin^d = ixomin^d - nghostcells * kr(idims,^d)\}
       {^d& ixcmax^d = ixomax^d + nghostcells * kr(idims,^d)\}
 
       hxo^l=ixo^l-kr(idims,^d);
       ! ix is centered index in the idim direction from ixOmin-1/2 to ixOmax+1/2
       ixmax^d=ixomax^d; ixmin^d=hxomin^d;
       ixcr^l=ixc^l;
 
       do idir=1,ndir
          call getfluxmf(wct,x,ixg^ll,ixcr^l,idir,idims,fct)
          ! Lax-Friedrich splitting:
          fp(ixcr^s,mag(idir)) = half * (fct(ixcr^s) + mf_tvdlfeps * tvdlfeps * cmax_global * wct(ixcr^s,mag(idir)))
          fm(ixcr^s,mag(idir)) = half * (fct(ixcr^s) - mf_tvdlfeps * tvdlfeps * cmax_global * wct(ixcr^s,mag(idir)))
       end do ! iw loop
 
       ! now do the reconstruction of fp and fm:
       call reconstructlmf(ixi^l,ix^l,idims,fp,fpl)
       call reconstructrmf(ixi^l,ix^l,idims,fm,fmr)
 
       do idir=1,ndir
          if (slab_uniform) then
             fc(ix^s,idir,idims) = dxinv(idims) * (fpl(ix^s,mag(idir)) + fmr(ix^s,mag(idir)))
             wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+ &
                  (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))
          else
             fc(ix^s,idir,idims)=-qdt*block%surfaceC(ix^s,idims) * (fpl(ix^s,mag(idir)) + fmr(ix^s,mag(idir)))
             wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+ &
               (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))/block%dvolume(ixo^s)
          end if
       end do ! iw loop
 
    end do !idims loop
 
    if (.not.slab) call addgeometrymf(qdt,ixi^l,ixo^l,wct,wnew,x)
    call divbclean(qdt,ixi^l,ixo^l,wct,wnew,x)
 
  end subroutine fdmf
 
  subroutine reconstructlmf(ixI^L,iL^L,idims,w,wLC)
    use mod_global_parameters
    use mod_limiter
 
    integer, intent(in)             :: ixI^L, iL^L, idims
    double precision, intent(in)    :: w(ixI^S,1:nw)
 
    double precision, intent(out)   :: wLC(ixI^S,1:nw)
 
    double precision                :: ldw(ixI^S), dwC(ixI^S)
    integer                         :: jxR^L, ixC^L, jxC^L, kxC^L, iw
 
    select case (type_limiter(block%level))
    case (limiter_mp5)
       call mp5limiterl(ixi^l,il^l,idims,w,wlc)
    case (limiter_weno5)
       call weno5limiterl(ixi^l,il^l,idims,w,wlc,1)
    case (limiter_wenoz5)
       call weno5limiterl(ixi^l,il^l,idims,w,wlc,2)
    case default
 
       kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
 
       wlc(kxc^s,1:nwflux) = w(kxc^s,1:nwflux)
 
       jxr^l=il^l+kr(idims,^d);
 
       ixcmax^d=jxrmax^d; ixcmin^d=ilmin^d-kr(idims,^d);
       jxc^l=ixc^l+kr(idims,^d);
 
       do iw=1,nwflux
          dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
 
          call dwlimiter2(dwc,ixi^l,ixc^l,idims,type_limiter(block%level),ldw=ldw)
 
          wlc(il^s,iw)=wlc(il^s,iw)+half*ldw(il^s)
       end do
    end select
 
  end subroutine reconstructlmf
 
  subroutine reconstructrmf(ixI^L,iL^L,idims,w,wRC)
    use mod_global_parameters
    use mod_limiter
 
    integer, intent(in)             :: ixI^L, iL^L, idims
    double precision, intent(in)    :: w(ixI^S,1:nw)
 
    double precision, intent(out)   :: wRC(ixI^S,1:nw)
 
    double precision                :: rdw(ixI^S), dwC(ixI^S)
    integer                         :: jxR^L, ixC^L, jxC^L, kxC^L, kxR^L, iw
 
    select case (type_limiter(block%level))
    case (limiter_mp5)
       call mp5limiterr(ixi^l,il^l,idims,w,wrc)
    case (limiter_weno5)
       call weno5limiterr(ixi^l,il^l,idims,w,wrc,1)
    case (limiter_wenoz5)
       call weno5limiterr(ixi^l,il^l,idims,w,wrc,2)
    case default
 
       kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
       kxr^l=kxc^l+kr(idims,^d);
 
       wrc(kxc^s,1:nwflux)=w(kxr^s,1:nwflux)
 
       jxr^l=il^l+kr(idims,^d);
       ixcmax^d=jxrmax^d; ixcmin^d=ilmin^d-kr(idims,^d);
       jxc^l=ixc^l+kr(idims,^d);
 
       do iw=1,nwflux
          dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
          call dwlimiter2(dwc,ixi^l,ixc^l,idims,type_limiter(block%level),rdw=rdw)
 
          wrc(il^s,iw)=wrc(il^s,iw)-half*rdw(jxr^s)
       end do
    end select
 
  end subroutine reconstructrmf
 
  subroutine centdiff4mf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,w,fC,dx^D,x)
    ! Advance the flow variables from global_time to global_time+qdt within ixO^L by
    ! fourth order centered differencing in space
    ! for the dw/dt+dF_i(w)/dx_i=S type equation.
    ! wCT contains the time centered variables at time qtC for flux and source.
    ! w is the old value at qt on input and the new value at qt+qdt on output.
    use mod_global_parameters
    use mod_comm_lib, only: mpistop
 
    integer, intent(in) :: ixI^L, ixO^L, idim^LIM
    double precision, intent(in) :: qdt, qtC, qt, dx^D
    double precision :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
    double precision, intent(in) :: x(ixI^S,1:ndim)
    double precision :: fC(ixI^S,1:ndir,1:ndim)
 
    double precision :: v(ixI^S,ndim), f(ixI^S)
    double precision, dimension(ixI^S,1:nw) :: wLC, wRC
    double precision, dimension(ixI^S)      :: vLC, vRC,cmaxLC,cmaxRC
    double precision :: dxinv(1:ndim)
    integer :: idims, idir, idirmin,ix^D
    integer :: ix^L, hxO^L, ixC^L, jxC^L, hxC^L, kxC^L, kkxC^L, kkxR^L
 
    ! two extra layers are needed in each direction for which fluxes are added.
    ix^l=ixo^l;
    do idims= idim^lim
       ix^l=ix^l^ladd2*kr(idims,^d);
    end do
 
    if (ixi^l^ltix^l|.or.|.or.) then
       call mpistop("Error in evolve_CentDiff4: Non-conforming input limits")
    end if
    ^d&dxinv(^d)=-qdt/dx^d;
 
    ! Add fluxes to w
    do idims= idim^lim
       b0i=idims
       ix^l=ixo^l^ladd2*kr(idims,^d);
       hxo^l=ixo^l-kr(idims,^d);
 
       ixcmin^d=hxomin^d; ixcmax^d=ixomax^d;
       hxc^l=ixc^l-kr(idims,^d);
       jxc^l=ixc^l+kr(idims,^d);
       kxc^l=ixc^l+2*kr(idims,^d);
 
       kkxcmin^d=iximin^d; kkxcmax^d=iximax^d-kr(idims,^d);
       kkxr^l=kkxc^l+kr(idims,^d);
       wrc(kkxc^s,1:nwflux)=wct(kkxr^s,1:nwflux)
       wlc(kkxc^s,1:nwflux)=wct(kkxc^s,1:nwflux)
 
       call upwindlrmf(ixi^l,ixc^l,ixc^l,idims,wct,wct,wlc,wrc,x)
 
       ! Calculate velocities from upwinded values
       call getcmaxfff(wlc,ixg^ll,ixc^l,idims,cmaxlc)
       call getcmaxfff(wrc,ixg^ll,ixc^l,idims,cmaxrc)
       ! now take the maximum of left and right states
       vlc(ixc^s)=max(cmaxrc(ixc^s),cmaxlc(ixc^s))
 
       do idir=1,ndir
          ! Get non-transported flux
          call getfluxmf(wct,x,ixi^l,ix^l,idir,idims,f)
          ! Center flux to interface
          ! f_i+1/2= (-f_(i+2) +7 f_(i+1) + 7 f_i - f_(i-1))/12
          fc(ixc^s,idir,idims)=(-f(kxc^s)+7.0d0*(f(jxc^s)+f(ixc^s))-f(hxc^s))/12.0d0
          ! add rempel dissipative flux, only second order version for now
          ! one could gradually reduce the dissipative flux to improve solutions
          ! for computing steady states (Keppens et al. 2012, JCP)
          fc(ixc^s,idir,idims)=fc(ixc^s,idir,idims)-mf_tvdlfeps*tvdlfeps*half*vlc(ixc^s) &
                                         *(wrc(ixc^s,mag(idir))-wlc(ixc^s,mag(idir)))
 
          if (slab_uniform) then
             fc(ixc^s,idir,idims)=dxinv(idims)*fc(ixc^s,idir,idims)
             ! result: f_(i+1/2)-f_(i-1/2) = [-f_(i+2)+8(f_(i+1)+f_(i-1))-f_(i-2)]/12
             w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+(fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))
          else
             fc(ixc^s,idir,idims)=-qdt*block%surfaceC(ixc^s,idims)*fc(ixc^s,idir,idims)
             w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+ &
                  (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))/block%dvolume(ixo^s)
          end if
       end do    !next idir
    end do       !next idims
    b0i=0
 
    if (.not.slab) call addgeometrymf(qdt,ixi^l,ixo^l,wct,w,x)
    call divbclean(qdt,ixi^l,ixo^l,wct,w,x)
 
  end subroutine centdiff4mf
 
  subroutine getdtfff_courant(w,x,ixI^L,ixO^L,dtnew)
    ! compute CFL limited dt (for variable time stepping)
    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), dtnew
 
    double precision :: courantmax, dxinv(1:ndim)
    double precision :: cmax(ixI^S),tmp(ixI^S),alfven(ixI^S)
    integer :: idims
 
    dtnew=bigdouble
    courantmax=0.d0
    ^d&dxinv(^d)=1.d0/dxlevel(^d);
 
    do idims=1,ndim
       call getcmaxfff(w,ixi^l,ixo^l,idims,cmax)
       cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
       if (.not.slab_uniform) then
          tmp(ixo^s)=cmax(ixo^s)/block%dx(ixo^s,idims)
          courantmax=max(courantmax,maxval(tmp(ixo^s)))
       else
          tmp(ixo^s)=cmax(ixo^s)*dxinv(idims)
          courantmax=max(courantmax,maxval(tmp(ixo^s)))
       end if
    end do
    ! courantmax='max( c/dx)'
    if (courantmax>smalldouble)  dtnew=min(dtnew,mf_cc/courantmax)
 
  end subroutine getdtfff_courant
 
  subroutine getcmaxfff(w,ixI^L,ixO^L,idims,cmax)
    use mod_global_parameters
 
    logical :: new_cmax,needcmin
    integer, intent(in) :: ixI^L, ixO^L, idims
    double precision, intent(in)    :: w(ixI^S,1:nw)
    double precision, intent(out) :: cmax(ixI^S)
 
    ! calculate alfven speed
    if(b0field) then
      cmax(ixo^s)=sqrt(sum((w(ixo^s,mag(:))+block%b0(ixo^s,:,0))**2,dim=ndim+1)/w(ixo^s,rho_))
    else
      cmax(ixo^s)=sqrt(sum(w(ixo^s,mag(:))**2,dim=ndim+1)/w(ixo^s,rho_))
    endif
    cmax(ixo^s)=cmax(ixo^s)+abs(w(ixo^s,mom(idims)))
 
  end subroutine getcmaxfff
 
  !> Clean divergence of magnetic field by Janhunen's and Linde's source terms
  subroutine divbclean(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)    :: x(ixI^S,1:ndim),wCT(ixI^S,1:nw),qdt
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision :: divb(ixI^S),graddivb(ixI^S),bdivb(ixI^S,1:ndir)
    integer :: idims, ix^L, ixp^L, i^D, iside
 
    ! Calculate div B
    ix^l=ixo^l^ladd1;
    call get_divb(wct,ixi^l,ix^l,divb)
 
    ixp^l=ixo^l;
 
    ! Add Linde's diffusive terms
    do idims=1,ndim
       ! Calculate grad_idim(divb)
       call gradient(divb,ixi^l,ixp^l,idims,graddivb)
 
       ! Multiply by Linde's eta*dt = divbdiff*(c_max*dx)*dt = divbdiff*dx**2
       if (slab_uniform) then
          graddivb(ixp^s)=graddivb(ixp^s)*mf_cdivb/(^d&1.0d0/dxlevel(^d)**2+)
       else
          graddivb(ixp^s)=graddivb(ixp^s)*mf_cdivb &
                          /(^d&1.0d0/block%dx(ixp^s,^d)**2+)
       end if
       ! B_idim += eta*grad_idim(divb)
       ! with Janhunen's term
       !w(ixp^S,mag(idims))=w(ixp^S,mag(idims))+&
       !      graddivb(ixp^S)-qdt*w(ixp^S,mom(idims))*divb(ixp^S)
       ! without Janjunen's term
       w(ixp^s,mag(idims))=w(ixp^s,mag(idims))+&
             graddivb(ixp^s)
    end do
 
  end subroutine divbclean
 
  subroutine addgeometrymf(qdt,ixI^L,ixO^L,wCT,w,x)
    ! Add geometrical source terms to w
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt, x(ixI^S,1:ndim)
    double precision, intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
    !.. local ..
    double precision :: tmp(ixI^S)
    integer          :: iw
    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
        ! s[Bphi]=(Bphi*vr-Br*vphi)/radius
        tmp(ixo^s)=(wct(ixo^s,bphi_)*wct(ixo^s,mom(1)) &
                   -wct(ixo^s,br_)*wct(ixo^s,mom(3)))
        w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt*tmp(ixo^s)/x(ixo^s,1)
      end if
    case (spherical)
    {^nooned
      ! s[b2]=(vr*Btheta-vtheta*Br)/r
      !       + cot(theta)*psi/r
      tmp(ixo^s)= wct(ixo^s,mom(1))*wct(ixo^s,mag(2)) &
                 -wct(ixo^s,mom(2))*wct(ixo^s,mag(1))
      if (b0field) then
         tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(1))*block%b0(ixo^s,2,0) &
                    -wct(ixo^s,mom(2))*block%b0(ixo^s,1,0)
      end if
      ! Divide by radius and add to w
      w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
    }
      if(ndir==3) then
        ! s[b3]=(vr*Bphi-vphi*Br)/r
        ! -cot(theta)*(vphi*Btheta-vtheta*Bphi)/r
        tmp(ixo^s)=wct(ixo^s,mom(1))*wct(ixo^s,mag(3)) &
                -wct(ixo^s,mom(3))*wct(ixo^s,mag(1)){^nooned &
               -(wct(ixo^s,mom(3))*wct(ixo^s,mag(2)) &
                -wct(ixo^s,mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
                              /dsin(x(ixo^s,2)) }
        if (b0field) then
           tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(1))*block%b0(ixo^s,3,0) &
              -wct(ixo^s,mom(3))*block%b0(ixo^s,1,0){^nooned &
              -(wct(ixo^s,mom(3))*block%b0(ixo^s,2,0) &
               -wct(ixo^s,mom(2))*block%b0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
                              /dsin(x(ixo^s,2)) }
        end if
        ! Divide by radius and add to w
        w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
      end if
    end select
 
  end subroutine addgeometrymf
 
  !> Calculate idirmin and the idirmin:3 components of the common current array
  !> make sure that dxlevel(^D) is set correctly.
  subroutine get_current(w,ixI^L,ixO^L,idirmin,current)
    use mod_global_parameters
    use mod_geometry
 
    integer :: idirmin0
    integer :: ixO^L, idirmin, ixI^L
    double precision :: 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),bvec(ixI^S,1:ndir)
    integer :: idir
 
    idirmin0 = 7-2*ndir
 
    bvec(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))
 
    call curlvector(bvec,ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
 
    if(b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
        block%J0(ixo^s,idirmin0:3)
 
  end subroutine get_current
 
  !> Calculate div B within ixO
  subroutine get_divb(w,ixI^L,ixO^L,divb)
    use mod_geometry
    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)
 
    double precision                   :: bvec(ixI^S,1:ndir)
 
    bvec(ixi^s,:)=w(ixi^s,mag(:))
 
    select case(typediv)
    case("central")
      call divvector(bvec,ixi^l,ixo^l,divb)
    case("limited")
      call divvectors(bvec,ixi^l,ixo^l,divb)
    end select
  end subroutine get_divb
 
end module mod_magnetofriction