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_integrator='onestep' ! time marching scheme, or 'twostep','threestep'
!>    flux_scheme=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.1 ! 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
  !> 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

  !> 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
  logical :: mf_advance
  logical :: fix_conserve_at_step = .true.

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_cc, mf_cy, mf_cy_max, mf_cdivb, mf_tvdlfeps_min

    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_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.1d0   ! divb cleaning coefficient controls diffusion speed of divb
    mf_tvdlfeps=1.d0 ! coefficient to control the TVDLF dissipation
    mf_tvdlfeps_min = mf_tvdlfeps ! minimum of the TVDLF dissipation coefficient

    call mf_params_read(par_files)

    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

    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
    integer :: i,iigrid, igrid, idims,ix^D,hxM^LL,fhmf,tmpit,i^D
    logical :: patchwi(ixg^t), stagger_flag

    ! 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=it
    if(snapshotini==-1 .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)-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)-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
    ! 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)
        pso(igrid)%w(ixg^t,mag(:))=ps(igrid)%w(ixg^t,mag(:))
        ^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

      i=i+1
      tmf=tmf+dtfff
      if(mod(i,10)==0) then
        ! calculate metrics
        call metrics
        call printlog_mf
      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
        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) call mpi_file_close(fhmf,ierrmpi)
    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
        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)
        character(len=800) :: filename,filehead
        character(len=2048) :: line,datastr
        logical, save :: logmfopened=.false.

        if(mype==0) then
          if(.not.logmfopened) then
            ! generate filename
            write(filename,"(a,a)") trim(base_filename), "_mflog.csv"

            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="  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
          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                   :: 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
    integer :: i,iigrid, igrid
    double precision :: vhatmax,vhatmax_pe,vhatmaxgrid

    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

    if(b0field) then
      tmp(ixo^s)=sum((w(ixo^s,mag(:))+block%b0(ixo^s,:,0))**2,dim=ndim+1)  ! |B|**2
    else
      tmp(ixo^s)=sum(w(ixo^s,mag(:))**2,dim=ndim+1)         ! |B|**2
    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\}
}
  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

    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

    ! copy w instead of wold because of potential use of dimsplit or sourcesplit
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      ps1(igrid)%w=ps(igrid)%w
    end do

    istep=0

    select case (time_integrator)
     case ("onestep")
       call advect1mf(flux_scheme,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_scheme,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_scheme,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_scheme,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_scheme,qdt,2.0d0/3.0d0,idim^lim,qt+qdt/2.0d0,ps2,qt+qdt/3.0d0,ps)
     case default
       write(unitterm,*) "time_integrator=",time_integrator
       write(unitterm,*) "Error in advectmf: Unknown time integration method"
       call mpistop("Correct time_integrator")
    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
    character(len=*), 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,pso(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,wold)
    ! This subroutine is equivalent to VAC's `advect1' for one grid
    use mod_global_parameters
    use mod_fix_conserve

    character(len=*), 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), wold(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);
    saveigrid=igrid
    fc=0.d0

    typelimiter=type_limiter(node(plevel_,igrid))
    typegradlimiter=type_gradient_limiter(node(plevel_,igrid))

    ixo^l=ixg^l^lsubnghostcells;
    select case (method)
     case ("cd4")
       !================================
       ! 4th order central difference
       !================================
       call centdiff4mf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,wold,fc,dx^d,ps(igrid)%x)
     case ("tvdlf")
       !================================
       ! TVDLF
       !================================
       call tvdlfmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,wold,fc,dx^d,ps(igrid)%x)
     case ('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 ("fd")
       !================================
       ! finite difference
       !================================
       call fdmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,wold,fc,dx^d,ps(igrid)%x)
    case default
       if(mype==0) write(unitterm,*)'Error in advect1_gridmf:',method,' is not there!'
       call mpistop("The scheme is not implemented.")
    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

    integer :: jxR^L, ixC^L, jxC^L, iw
    double precision   :: ldw(ixi^s), rdw(ixi^s), dwC(ixi^s)

    if (typelimiter == limiter_mp5) then
       call mp5limiter(ixi^l,ixl^l,idim,w,wlc,wrc)
    else if (typelimiter == 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,typelimiter,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,wold,fC,dx^D,x)
    ! method=='tvdlf'  --> 2nd order TVD-Lax-Friedrich scheme.
    ! method=='tvdlf1' --> 1st order TVD-Lax-Friedrich scheme.
    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)               :: wCT, wnew, wold
    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
       block%iw0=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

    !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

    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
       block%iw0=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

    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,wold,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, wold
    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 (typelimiter)
    case (limiter_mp5)
       call mp5limiterl(ixi^l,il^l,idims,w,wlc)
    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,typelimiter,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 (typelimiter)
    case (limiter_mp5)
       call mp5limiterr(ixi^l,il^l,idims,w,wrc)
    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,typelimiter,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,wold,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

    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), wold(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
       block%iw0=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

    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)
    integer :: idims, ix^L, ixp^L, i^D, iside
    double precision :: divb(ixi^s),graddivb(ixi^s),bdivb(ixi^s,1:ndir)

    ! 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)
    integer :: idir

    ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
    double precision :: current(ixi^s,7-2*ndir:3),bvec(ixi^s,1:ndir)

    idirmin0 = 7-2*ndir

    bvec(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))

    call curlvector(bvec,ixi^l,ixo^l,current,idirmin,idirmin0,ndir)

    if(b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
        block%J0(ixo^s,idirmin0:3)

  end subroutine get_current

  !> 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