mod_finite_volume.t

Go to the documentation of this file.

!> Module with finite volume methods for fluxes
module mod_finite_volume
#include "amrvac.h"
  implicit none
  private
 
  public :: finite_volume
  public :: hancock
  public :: reconstruct_lr
 
contains
 
  !> The non-conservative Hancock predictor for TVDLF
  !>
  !> 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
  subroutine hancock(qdt,dtfactor,ixI^L,ixO^L,idims^LIM,qtC,sCT,qt,snew,dxs,x)
    use mod_physics
    use mod_global_parameters
    use mod_source, only: addsource2
    use mod_comm_lib, only: mpistop
 
    integer, intent(in) :: ixi^l, ixo^l, idims^lim
    double precision, intent(in) :: qdt, dtfactor,qtc, qt, dxs(ndim), x(ixi^s,1:ndim)
    type(state) :: sct, snew
 
    double precision, dimension(ixI^S,1:nw) :: wprim, wlc, wrc
    ! left and right constructed status in primitive form, needed for better performance
    double precision, dimension(ixI^S,1:nw) :: wlp, wrp
    double precision, dimension(ixO^S)      :: inv_volume
    double precision :: flc(ixi^s, nwflux), frc(ixi^s, nwflux)
    double precision :: dxinv(1:ndim)
    integer :: idims, iw, ix^l, hxo^l
    logical :: active
 
    associate(wct=>sct%w,wnew=>snew%w)
    ! Expand limits in each idims direction in which fluxes are added
    ix^l=ixo^l;
    do idims= idims^lim
       ix^l=ix^l^laddkr(idims,^d);
    end do
    if (ixi^l^ltix^l|.or.|.or.) &
         call mpistop("Error in Hancock: Nonconforming input limits")
 
    wprim=wct
    call phys_to_primitive(ixi^l,ixi^l,wprim,x)
 
    dxinv=-qdt/dxs
    if(.not.slab_uniform) inv_volume = 1.d0/block%dvolume(ixo^s)
    do idims= idims^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);
 
      wrp(hxo^s,1:nwflux)=wprim(ixo^s,1:nwflux)
      wlp(ixo^s,1:nwflux)=wprim(ixo^s,1:nwflux)
 
      ! apply limited reconstruction for left and right status at cell interfaces
      call reconstruct_lr(ixi^l,ixo^l,hxo^l,idims,wprim,wlc,wrc,wlp,wrp,x,dxs(idims))
 
      ! Calculate the fLC and fRC fluxes
      call phys_get_flux(wrc,wrp,x,ixi^l,hxo^l,idims,frc)
      call phys_get_flux(wlc,wlp,x,ixi^l,ixo^l,idims,flc)
 
      ! Advect w(iw)
      if (slab_uniform) then
        if(local_timestep) then
          do iw=1,nwflux
            wnew(ixo^s,iw)=wnew(ixo^s,iw)-block%dt(ixo^s)*dtfactor/dxs(idims)* &
                 (flc(ixo^s, iw)-frc(hxo^s, iw))
          end do
        else  
          do iw=1,nwflux
            wnew(ixo^s,iw)=wnew(ixo^s,iw)+dxinv(idims)* &
                 (flc(ixo^s, iw)-frc(hxo^s, iw))
          end do
        endif
      else
        if(local_timestep) then
          do iw=1,nwflux
            wnew(ixo^s,iw)=wnew(ixo^s,iw) - block%dt(ixo^s)*dtfactor * inv_volume &
                 *(block%surfaceC(ixo^s,idims)*flc(ixo^s, iw) &
                 -block%surfaceC(hxo^s,idims)*frc(hxo^s, iw))
          end do
        else
          do iw=1,nwflux
            wnew(ixo^s,iw)=wnew(ixo^s,iw) - qdt * inv_volume &
                 *(block%surfaceC(ixo^s,idims)*flc(ixo^s, iw) &
                 -block%surfaceC(hxo^s,idims)*frc(hxo^s, iw))
          end do
        end if
      end if
    end do ! next idims
    b0i=0
 
    if (.not.slab.and.idimsmin==1) call phys_add_source_geom(qdt,dtfactor,ixi^l,ixo^l,wct,wnew,x)
 
    call addsource2(qdt*dble(idimsmax-idimsmin+1)/dble(ndim), &
          dtfactor*dble(idimsmax-idimsmin+1)/dble(ndim),& 
         ixi^l,ixo^l,1,nw,qtc,wct,wprim,qt,wnew,x,.false.,active)
 
    ! check and optionally correct unphysical values
    if(fix_small_values) then
       call phys_handle_small_values(.false.,wnew,x,ixi^l,ixo^l,'exit hancock finite_volume')
    endif
    end associate
  end subroutine hancock
 
  !> finite volume method
  subroutine finite_volume(method,qdt,dtfactor,ixI^L,ixO^L,idims^LIM, &
       qtC,sCT,qt,snew,fC,fE,dxs,x)
    use mod_physics
    use mod_variables
    use mod_global_parameters
    use mod_tvd, only:tvdlimit2
    use mod_source, only: addsource2
    use mod_usr_methods
    use mod_comm_lib, only: mpistop
 
    integer, intent(in)                                   :: method
    double precision, intent(in)                          :: qdt, dtfactor, qtc, qt, dxs(ndim)
    integer, intent(in)                                   :: ixi^l, ixo^l, idims^lim
    double precision, dimension(ixI^S,1:ndim), intent(in) :: x
    type(state)                                           :: sct, snew
    double precision, dimension(ixI^S,1:nwflux,1:ndim)    :: fc
    double precision, dimension(ixI^S,sdim:3)             :: fe
 
    ! primitive w at cell center
    double precision, dimension(ixI^S,1:nw) :: wprim
    ! left and right constructed status in conservative form
    double precision, dimension(ixI^S,1:nw) :: wlc, wrc
    ! left and right constructed status in primitive form, needed for better performance
    double precision, dimension(ixI^S,1:nw) :: wlp, wrp
    double precision, dimension(ixI^S,1:nwflux) :: flc, frc
    double precision, dimension(ixI^S,1:number_species)      :: cmaxc
    double precision, dimension(ixI^S,1:number_species)      :: cminc
    double precision, dimension(ixI^S)      :: hspeed
    double precision, dimension(ixO^S)      :: inv_volume
    double precision, dimension(1:ndim)     :: dxinv
    integer, dimension(ixI^S)               :: patchf
    integer :: idims, iw, ix^l, hxo^l, ixc^l, ixcr^l, kxc^l, kxr^l, ii
    logical :: active
    type(ct_velocity) :: vcts
 
    associate(wct=>sct%w, wnew=>snew%w)
 
    fc=0.d0
    flc=0.d0
    frc=0.d0
 
    ! The flux calculation contracts by one in the idims direction it is applied.
    ! The limiter contracts the same directions by one more, so expand ixO by 2.
    ix^l=ixo^l;
    do idims= idims^lim
       ix^l=ix^l^ladd2*kr(idims,^d);
    end do
    if (ixi^l^ltix^l|.or.|.or.) &
         call mpistop("Error in fv : Nonconforming input limits")
 
    wprim=wct
    call phys_to_primitive(ixi^l,ixi^l,wprim,x)
 
    do idims= idims^lim
       ! use interface value of w0 at idims
       b0i=idims
 
       hxo^l=ixo^l-kr(idims,^d);
 
       kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
       kxr^l=kxc^l+kr(idims,^d);
 
       if(stagger_grid) then
         ! ct needs all transverse cells
         ixcmax^d=ixomax^d+nghostcells-nghostcells*kr(idims,^d);
         ixcmin^d=hxomin^d-nghostcells+nghostcells*kr(idims,^d);
       else
         ! ixC is centered index in the idims direction from ixOmin-1/2 to ixOmax+1/2
         ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
       end if
 
       ! wRp and wLp are defined at the same locations, and will correspond to
       ! the left and right reconstructed values at a cell face. Their indexing
       ! is similar to cell-centered values, but in direction idims they are
       ! shifted half a cell towards the 'lower' direction.
       wrp(kxc^s,1:nw)=wprim(kxr^s,1:nw)
       wlp(kxc^s,1:nw)=wprim(kxc^s,1:nw)
 
       ! Determine stencil size
       {ixcrmin^d = max(ixcmin^d - phys_wider_stencil,ixglo^d)\}
       {ixcrmax^d = min(ixcmax^d + phys_wider_stencil,ixghi^d)\}
 
       ! apply limited reconstruction for left and right status at cell interfaces
       call reconstruct_lr(ixi^l,ixcr^l,ixcr^l,idims,wprim,wlc,wrc,wlp,wrp,x,dxs(idims))
 
       ! special modification of left and right status before flux evaluation
       call phys_modify_wlr(ixi^l,ixcr^l,qt,wlc,wrc,wlp,wrp,sct,idims)
 
       ! evaluate physical fluxes according to reconstructed status
       call phys_get_flux(wlc,wlp,x,ixi^l,ixc^l,idims,flc)
       call phys_get_flux(wrc,wrp,x,ixi^l,ixc^l,idims,frc)
 
       if(h_correction) then
         call phys_get_h_speed(wprim,x,ixi^l,ixo^l,idims,hspeed)
       end if
       ! estimating bounds for the minimum and maximum signal velocities
       if(method==fs_tvdlf.or.method==fs_tvdmu) then
         call phys_get_cbounds(wlc,wrc,wlp,wrp,x,ixi^l,ixc^l,idims,hspeed,cmaxc)
         ! index of var  velocity appears in the induction eq. 
         if(stagger_grid) call phys_get_ct_velocity(vcts,wlp,wrp,ixi^l,ixc^l,idims,cmaxc(ixi^s,index_v_mag))
       else
         call phys_get_cbounds(wlc,wrc,wlp,wrp,x,ixi^l,ixc^l,idims,hspeed,cmaxc,cminc)
         if(stagger_grid) call phys_get_ct_velocity(vcts,wlp,wrp,ixi^l,ixc^l,idims,cmaxc(ixi^s,index_v_mag),cminc(ixi^s,index_v_mag))
       end if
 
       ! use approximate Riemann solver to get flux at interfaces
       select case(method)
       case(fs_hll)
         do ii=1,number_species
           call get_riemann_flux_hll(start_indices(ii),stop_indices(ii))
         end do
       case(fs_hllc,fs_hllcd)
         do ii=1,number_species
           call get_riemann_flux_hllc(start_indices(ii),stop_indices(ii))
         end do
       case(fs_hlld)
         do ii=1,number_species
           if(ii==index_v_mag) then
             call get_riemann_flux_hlld(start_indices(ii),stop_indices(ii))
           else
             call get_riemann_flux_hll(start_indices(ii),stop_indices(ii))
           endif   
         end do
       case(fs_tvdlf)
         do ii=1,number_species
           call get_riemann_flux_tvdlf(start_indices(ii),stop_indices(ii))
         end do
       case(fs_tvdmu)
         call get_riemann_flux_tvdmu()
       case default
         call mpistop('unkown Riemann flux in finite volume')
       end select
 
    end do ! Next idims
    b0i=0
    if(stagger_grid) call phys_update_faces(ixi^l,ixo^l,qt,qdt,wprim,fc,fe,sct,snew,vcts)
    if(slab_uniform) then
      dxinv=-qdt/dxs
      do idims= idims^lim
        hxo^l=ixo^l-kr(idims,^d);
        ! TODO maybe put if outside loop idims: but too much code is copy pasted
        ! this is also done in hancock and fd, centdiff in mod_finite_difference
        if(local_timestep) then
          do iw=iwstart,nwflux
            fc(ixi^s,iw,idims)=-block%dt(ixi^s)*dtfactor/dxs(idims)*fc(ixi^s,iw,idims)
          end do
        else
          ! Multiply the fluxes by -dt/dx since Flux fixing expects this
          fc(ixi^s,1:nwflux,idims)=dxinv(idims)*fc(ixi^s,1:nwflux,idims)
        end if
 
        wnew(ixo^s,iwstart:nwflux)=wnew(ixo^s,iwstart:nwflux)+&
            (fc(ixo^s,iwstart:nwflux,idims)-fc(hxo^s,iwstart:nwflux,idims))
 
        ! For the MUSCL scheme apply the characteristic based limiter
        if(method==fs_tvdmu) &
           call tvdlimit2(method,qdt,ixi^l,ixc^l,ixo^l,idims,wlc,wrc,wnew,x,fc,dxs)
 
      end do ! Next idims
    else
      inv_volume = 1.d0/block%dvolume(ixo^s)
      do idims= idims^lim
        hxo^l=ixo^l-kr(idims,^d);
 
        if(local_timestep) then
          do iw=iwstart,nwflux
            fc(ixi^s,iw,idims)=-block%dt(ixi^s)*dtfactor*fc(ixi^s,iw,idims)*block%surfaceC(ixi^s,idims)
            wnew(ixo^s,iw)=wnew(ixo^s,iw) + (fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims)) * &
                inv_volume
          end do
        else
          do iw=iwstart,nwflux
            fc(ixi^s,iw,idims)=-qdt*fc(ixi^s,iw,idims)*block%surfaceC(ixi^s,idims)
            wnew(ixo^s,iw)=wnew(ixo^s,iw) + (fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims)) * &
                inv_volume
          end do
        end if 
        ! For the MUSCL scheme apply the characteristic based limiter
        if (method==fs_tvdmu) &
             call tvdlimit2(method,qdt,ixi^l,ixc^l,ixo^l,idims,wlc,wrc,wnew,x,fc,dxs)
 
      end do ! Next idims
    end if
 
    if (.not.slab.and.idimsmin==1) &
         call phys_add_source_geom(qdt,dtfactor,ixi^l,ixo^l,wct,wnew,x)
 
    if(stagger_grid) call phys_face_to_center(ixo^l,snew)
 
    ! check and optionally correct unphysical values
    if(fix_small_values) then
       call phys_handle_small_values(.false.,wnew,x,ixi^l,ixo^l,'multi-D finite_volume')
       if(crash) then
         ! replace erroneous values with values at previous step
         wnew=wct
         if(stagger_grid) snew%ws=sct%ws
       end if
    end if
 
    call addsource2(qdt*dble(idimsmax-idimsmin+1)/dble(ndim),& 
         dtfactor*dble(idimsmax-idimsmin+1)/dble(ndim),&
         ixi^l,ixo^l,1,nw,qtc,wct,wprim,qt,wnew,x,.false.,active)
 
  end associate
  contains
 
    subroutine get_riemann_flux_tvdmu()
      do iw=iwstart,nwflux
         ! To save memory we use fLC to store (F_L+F_R)/2=half*(fLC+fRC)
         flc(ixc^s, iw)=half*(flc(ixc^s, iw)+frc(ixc^s, iw))
         fc(ixc^s,iw,idims)=flc(ixc^s, iw)
      end do
    end subroutine get_riemann_flux_tvdmu
 
    subroutine get_riemann_flux_tvdlf(iws,iwe)
      integer, intent(in) :: iws,iwe
      double precision :: fac(ixC^S)
 
      fac(ixc^s) = -0.5d0*tvdlfeps*cmaxc(ixc^s,ii)
      ! Calculate fLC=f(uL_j+1/2) and fRC=f(uR_j+1/2) for each iw
      do iw=iws,iwe
         ! To save memory we use fLC to store (F_L+F_R)/2=half*(fLC+fRC)
         flc(ixc^s, iw)=0.5d0*(flc(ixc^s, iw)+frc(ixc^s, iw))
         ! Add TVDLF dissipation to the flux
         if (flux_type(idims, iw) /= flux_no_dissipation) then
            flc(ixc^s, iw)=flc(ixc^s, iw) + fac(ixc^s)*(wrc(ixc^s,iw)-wlc(ixc^s,iw))
         end if
         fc(ixc^s,iw,idims)=flc(ixc^s, iw)
      end do ! Next iw
 
    end subroutine get_riemann_flux_tvdlf
 
    subroutine get_riemann_flux_hll(iws,iwe)
      integer, intent(in) :: iws,iwe
      integer :: ix^D
 
      do iw=iws,iwe
        if(flux_type(idims, iw) == flux_tvdlf) then
          ! CT MHD does not need normal B flux
          if(stagger_grid) cycle
          fc(ixc^s,iw,idims) = -tvdlfeps*half*max(cmaxc(ixc^s,ii),dabs(cminc(ixc^s,ii))) * &
               (wrc(ixc^s,iw)-wlc(ixc^s,iw))
        else
         {do ix^db=ixcmin^db,ixcmax^db\}
           if(cminc(ix^d,ii) >= zero) then
             fc(ix^d,iw,idims)=flc(ix^d,iw)
           else if(cmaxc(ix^d,ii) <= zero) then
             fc(ix^d,iw,idims)=frc(ix^d,iw)
           else
             ! Add hll dissipation to the flux
             fc(ix^d,iw,idims)=(cmaxc(ix^d,ii)*flc(ix^d, iw)-cminc(ix^d,ii)*frc(ix^d,iw)&
                   +cminc(ix^d,ii)*cmaxc(ix^d,ii)*(wrc(ix^d,iw)-wlc(ix^d,iw)))&
                   /(cmaxc(ix^d,ii)-cminc(ix^d,ii))
           end if
         {end do\}
       endif 
      end do
 
    end subroutine get_riemann_flux_hll
 
    subroutine get_riemann_flux_hllc(iws,iwe)
      integer, intent(in) :: iws, iwe  
      double precision, dimension(ixI^S,1:nwflux)     :: whll, Fhll, fCD
      double precision, dimension(ixI^S)              :: lambdaCD
 
      integer  :: rho_, p_, e_, mom(1:ndir)
 
      rho_ = iw_rho
      if (allocated(iw_mom)) mom(:) = iw_mom(:)
      e_ = iw_e 
 
      if(associated(phys_hllc_init_species)) then
       call phys_hllc_init_species(ii, rho_, mom(:), e_)
      endif  
 
      p_ = e_
 
      patchf(ixc^s) =  1
      where(cminc(ixc^s,1) >= zero)
         patchf(ixc^s) = -2
      elsewhere(cmaxc(ixc^s,1) <= zero)
         patchf(ixc^s) =  2
      endwhere
      ! Use more diffusive scheme, is actually TVDLF and selected by patchf=4
      if(method==fs_hllcd) &
           call phys_diffuse_hllcd(ixi^l,ixc^l,idims,wlc,wrc,flc,frc,patchf)
 
      !---- calculate speed lambda at CD ----!
      if(any(patchf(ixc^s)==1)) &
           call phys_get_lcd(wlc,wrc,flc,frc,cminc(ixi^s,ii),cmaxc(ixi^s,ii),idims,ixi^l,ixc^l, &
           whll,fhll,lambdacd,patchf)
 
      ! now patchf may be -1 or 1 due to phys_get_lCD
      if(any(abs(patchf(ixc^s))== 1))then
         !======== flux at intermediate state ========!
         call phys_get_wcd(wlc,wrc,whll,frc,flc,fhll,patchf,lambdacd,&
              cminc(ixi^s,ii),cmaxc(ixi^s,ii),ixi^l,ixc^l,idims,fcd)
      endif ! Calculate the CD flux
 
      do iw=iws,iwe
         if (flux_type(idims, iw) == flux_tvdlf) then
            flc(ixc^s,iw)=-tvdlfeps*half*max(cmaxc(ixc^s,ii),abs(cminc(ixc^s,ii))) * &
                 (wrc(ixc^s,iw) - wlc(ixc^s,iw))
         else
            where(patchf(ixc^s)==-2)
               flc(ixc^s,iw)=flc(ixc^s,iw)
            elsewhere(abs(patchf(ixc^s))==1)
               flc(ixc^s,iw)=fcd(ixc^s,iw)
            elsewhere(patchf(ixc^s)==2)
               flc(ixc^s,iw)=frc(ixc^s,iw)
            elsewhere(patchf(ixc^s)==3)
               ! fallback option, reducing to HLL flux
               flc(ixc^s,iw)=fhll(ixc^s,iw)
            elsewhere(patchf(ixc^s)==4)
               ! fallback option, reducing to TVDLF flux
               flc(ixc^s,iw) = half*((flc(ixc^s,iw)+frc(ixc^s,iw)) &
                    -tvdlfeps * max(cmaxc(ixc^s,ii), dabs(cminc(ixc^s,ii))) * &
                    (wrc(ixc^s,iw)-wlc(ixc^s,iw)))
            endwhere
         end if
 
         fc(ixc^s,iw,idims)=flc(ixc^s,iw)
 
      end do ! Next iw
    end subroutine get_riemann_flux_hllc
 
    !> HLLD Riemann flux from Miyoshi 2005 JCP, 208, 315 and Guo 2016 JCP, 327, 543
    subroutine get_riemann_flux_hlld(iws,iwe)
      integer, intent(in) :: iws, iwe
      double precision, dimension(ixI^S,1:nwflux) :: w1R,w1L,f1R,f1L,f2R,f2L
      double precision, dimension(ixI^S,1:nwflux) :: w2R,w2L
      double precision, dimension(ixI^S) :: sm,s1R,s1L,suR,suL,Bx
      double precision, dimension(ixI^S) :: pts,ptR,ptL,signBx,r1L,r1R,tmp
      ! velocity from the right and the left reconstruction
      double precision, dimension(ixI^S,ndir) :: vRC, vLC
      ! magnetic field from the right and the left reconstruction
      double precision, dimension(ixI^S,ndir) :: BR, BL
      integer :: ip1,ip2,ip3,idir,ix^D
      integer  :: rho_, p_, e_, mom(1:ndir), mag(1:ndir)
 
      associate(sr=>cmaxc,sl=>cminc)
 
      rho_ = iw_rho
      mom(:) = iw_mom(:)
      mag(:) = iw_mag(:) 
      e_ = iw_e 
      p_ = e_
 
      f1r=0.d0
      f1l=0.d0
      f2r=0.d0
      f2l=0.d0
      w1l=0.d0
      w1r=0.d0
      w2l=0.d0
      w2r=0.d0
      ip1=idims
      ip3=3
      vrc(ixc^s,:)=wrp(ixc^s,mom(:))
      vlc(ixc^s,:)=wlp(ixc^s,mom(:))
      if(b0field) then
        br(ixc^s,:)=wrc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
        bl(ixc^s,:)=wlc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
      else
        br(ixc^s,:)=wrc(ixc^s,mag(:))
        bl(ixc^s,:)=wlc(ixc^s,mag(:))
      end if
      if(stagger_grid) then
        bx(ixc^s)=block%ws(ixc^s,ip1)
      else
        ! HLL estimation of normal magnetic field at cell interfaces
        ! Li, Shenghai, 2005 JCP, 203, 344, equation (33)
        bx(ixc^s)=(sr(ixc^s,ii)*br(ixc^s,ip1)-sl(ixc^s,ii)*bl(ixc^s,ip1))/(sr(ixc^s,ii)-sl(ixc^s,ii))
      end if
      ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(br(ixc^s,:)**2,dim=ndim+1)
      ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(bl(ixc^s,:)**2,dim=ndim+1)
      if(iw_equi_rho>0) then
        sur(ixc^s) = wrc(ixc^s,rho_)+ block%equi_vars(ixc^s,iw_equi_rho,ip1)
      else
        sur(ixc^s) = wrc(ixc^s,rho_)
      endif
      sur(ixc^s)=(sr(ixc^s,ii)-vrc(ixc^s,ip1))*sur(ixc^s)
      if(iw_equi_rho>0) then
        sul(ixc^s) = wlc(ixc^s,rho_)+ block%equi_vars(ixc^s,iw_equi_rho,ip1)
      else
        sul(ixc^s) = wlc(ixc^s,rho_)
      endif
      sul(ixc^s)=(sl(ixc^s,ii)-vlc(ixc^s,ip1))*sul(ixc^s)
      ! Miyoshi equation (38) and Guo euqation (20)
      sm(ixc^s)=(sur(ixc^s)*vrc(ixc^s,ip1)-sul(ixc^s)*vlc(ixc^s,ip1)-&
                 ptr(ixc^s)+ptl(ixc^s))/(sur(ixc^s)-sul(ixc^s))
      ! Miyoshi equation (39) and Guo euqation (28)
      w1r(ixc^s,mom(ip1))=sm(ixc^s)
      w1l(ixc^s,mom(ip1))=sm(ixc^s)
      w2r(ixc^s,mom(ip1))=sm(ixc^s)
      w2l(ixc^s,mom(ip1))=sm(ixc^s)
      ! Guo equation (22)
      w1r(ixc^s,mag(ip1))=bx(ixc^s)
      w1l(ixc^s,mag(ip1))=bx(ixc^s)
      if(b0field) then
        ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(wrc(ixc^s,mag(:))**2,dim=ndim+1)
        ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(wlc(ixc^s,mag(:))**2,dim=ndim+1)
      end if
 
      ! Miyoshi equation (43) and Guo equation (27)
      w1r(ixc^s,rho_)=sur(ixc^s)/(sr(ixc^s,ii)-sm(ixc^s))
      w1l(ixc^s,rho_)=sul(ixc^s)/(sl(ixc^s,ii)-sm(ixc^s))
 
      ip2=mod(ip1+1,ndir)
      if(ip2==0) ip2=ndir
      r1r(ixc^s)=sur(ixc^s)*(sr(ixc^s,ii)-sm(ixc^s))-bx(ixc^s)**2
      where(abs(r1r(ixc^s))>smalldouble)
        r1r(ixc^s)=1.d0/r1r(ixc^s)
      else where
        r1r(ixc^s)=0.d0
      end where
      r1l(ixc^s)=sul(ixc^s)*(sl(ixc^s,ii)-sm(ixc^s))-bx(ixc^s)**2
      where(abs(r1l(ixc^s))>smalldouble)
        r1l(ixc^s)=1.d0/r1l(ixc^s)
      else where
        r1l(ixc^s)=0.d0
      end where
      ! Miyoshi equation (44)
      w1r(ixc^s,mom(ip2))=vrc(ixc^s,ip2)-bx(ixc^s)*br(ixc^s,ip2)*&
        (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
      w1l(ixc^s,mom(ip2))=vlc(ixc^s,ip2)-bx(ixc^s)*bl(ixc^s,ip2)*&
        (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
      ! partial solution for later usage
      w1r(ixc^s,mag(ip2))=(sur(ixc^s)*(sr(ixc^s,ii)-vrc(ixc^s,ip1))-bx(ixc^s)**2)*r1r(ixc^s)
      w1l(ixc^s,mag(ip2))=(sul(ixc^s)*(sl(ixc^s,ii)-vlc(ixc^s,ip1))-bx(ixc^s)**2)*r1l(ixc^s)
      if(ndir==3) then
        ip3=mod(ip1+2,ndir)
        if(ip3==0) ip3=ndir
        ! Miyoshi equation (46)
        w1r(ixc^s,mom(ip3))=vrc(ixc^s,ip3)-bx(ixc^s)*br(ixc^s,ip3)*&
          (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
        w1l(ixc^s,mom(ip3))=vlc(ixc^s,ip3)-bx(ixc^s)*bl(ixc^s,ip3)*&
          (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
        ! Miyoshi equation (47)
        w1r(ixc^s,mag(ip3))=br(ixc^s,ip3)*w1r(ixc^s,mag(ip2))
        w1l(ixc^s,mag(ip3))=bl(ixc^s,ip3)*w1l(ixc^s,mag(ip2))
      end if
      ! Miyoshi equation (45)
      w1r(ixc^s,mag(ip2))=br(ixc^s,ip2)*w1r(ixc^s,mag(ip2))
      w1l(ixc^s,mag(ip2))=bl(ixc^s,ip2)*w1l(ixc^s,mag(ip2))
      if(b0field) then
        ! Guo equation (26)
        w1r(ixc^s,mag(:))=w1r(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
        w1l(ixc^s,mag(:))=w1l(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
      end if
      ! equation (48)
      if(phys_energy) then
        ! Guo equation (25) equivalent to Miyoshi equation (41)
        w1r(ixc^s,p_)=sur(ixc^s)*(sm(ixc^s)-vrc(ixc^s,ip1))+ptr(ixc^s)
        !w1L(ixC^S,p_)=suL(ixC^S)*(sm(ixC^S)-vLC(ixC^S,ip1))+ptL(ixC^S)
        w1l(ixc^s,p_)=w1r(ixc^s,p_)
        if(b0field) then
          ! Guo equation (32)
          w1r(ixc^s,p_)=w1r(ixc^s,p_)+sum(block%B0(ixc^s,:,ip1)*(wrc(ixc^s,mag(:))-w1r(ixc^s,mag(:))),dim=ndim+1)
          w1l(ixc^s,p_)=w1l(ixc^s,p_)+sum(block%B0(ixc^s,:,ip1)*(wlc(ixc^s,mag(:))-w1l(ixc^s,mag(:))),dim=ndim+1)
        end if
        ! Miyoshi equation (48) and main part of Guo euqation (31)
        w1r(ixc^s,e_)=((sr(ixc^s,ii)-vrc(ixc^s,ip1))*wrc(ixc^s,e_)-ptr(ixc^s)*vrc(ixc^s,ip1)+&
          w1r(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vrc(ixc^s,:)*wrc(ixc^s,mag(:)),dim=ndim+1)-&
          sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)))/(sr(ixc^s,ii)-sm(ixc^s))
        w1l(ixc^s,e_)=((sl(ixc^s,ii)-vlc(ixc^s,ip1))*wlc(ixc^s,e_)-ptl(ixc^s)*vlc(ixc^s,ip1)+&
          w1l(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vlc(ixc^s,:)*wlc(ixc^s,mag(:)),dim=ndim+1)-&
          sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)))/(sl(ixc^s,ii)-sm(ixc^s))
        if(b0field) then
          ! Guo equation (31)
          w1r(ixc^s,e_)=w1r(ixc^s,e_)+(sum(w1r(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
               sum(wrc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
          w1l(ixc^s,e_)=w1l(ixc^s,e_)+(sum(w1l(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
               sum(wlc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
        end if
        if(iw_equi_p>0) then
#if !defined(E_RM_W0) || E_RM_W0 == 1
          w1r(ixc^s,e_)= w1r(ixc^s,e_) + 1d0/(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
             (sm(ixc^s)-vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
          w1l(ixc^s,e_)= w1l(ixc^s,e_) + 1d0/(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
             (sm(ixc^s)-vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
#else
          w1r(ixc^s,e_)= w1r(ixc^s,e_) + phys_gamma /(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
             (sm(ixc^s)-vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
          w1l(ixc^s,e_)= w1l(ixc^s,e_) + phys_gamma /(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
             (sm(ixc^s)-vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
#endif
        endif
      end if
 
      ! Miyoshi equation (49) and Guo equation (35)
      w2r(ixc^s,rho_)=w1r(ixc^s,rho_)
      w2l(ixc^s,rho_)=w1l(ixc^s,rho_)
      w2r(ixc^s,mag(ip1))=w1r(ixc^s,mag(ip1))
      w2l(ixc^s,mag(ip1))=w1l(ixc^s,mag(ip1))
 
      r1r(ixc^s)=sqrt(w1r(ixc^s,rho_))
      r1l(ixc^s)=sqrt(w1l(ixc^s,rho_))
      tmp(ixc^s)=1.d0/(r1r(ixc^s)+r1l(ixc^s))
      signbx(ixc^s)=sign(1.d0,bx(ixc^s))
      ! Miyoshi equation (51) and Guo equation (33)
      s1r(ixc^s)=sm(ixc^s)+abs(bx(ixc^s))/r1r(ixc^s)
      s1l(ixc^s)=sm(ixc^s)-abs(bx(ixc^s))/r1l(ixc^s)
      ! Miyoshi equation (59) and Guo equation (41)
      w2r(ixc^s,mom(ip2))=(r1l(ixc^s)*w1l(ixc^s,mom(ip2))+r1r(ixc^s)*w1r(ixc^s,mom(ip2))+&
          (w1r(ixc^s,mag(ip2))-w1l(ixc^s,mag(ip2)))*signbx(ixc^s))*tmp(ixc^s)
      w2l(ixc^s,mom(ip2))=w2r(ixc^s,mom(ip2))
      ! Miyoshi equation (61) and Guo equation (43)
      w2r(ixc^s,mag(ip2))=(r1l(ixc^s)*w1r(ixc^s,mag(ip2))+r1r(ixc^s)*w1l(ixc^s,mag(ip2))+&
          r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip2))-w1l(ixc^s,mom(ip2)))*signbx(ixc^s))*tmp(ixc^s)
      w2l(ixc^s,mag(ip2))=w2r(ixc^s,mag(ip2))
      if(ndir==3) then
        ! Miyoshi equation (60) and Guo equation (42)
        w2r(ixc^s,mom(ip3))=(r1l(ixc^s)*w1l(ixc^s,mom(ip3))+r1r(ixc^s)*w1r(ixc^s,mom(ip3))+&
            (w1r(ixc^s,mag(ip3))-w1l(ixc^s,mag(ip3)))*signbx(ixc^s))*tmp(ixc^s)
        w2l(ixc^s,mom(ip3))=w2r(ixc^s,mom(ip3))
        ! Miyoshi equation (62) and Guo equation (44)
        w2r(ixc^s,mag(ip3))=(r1l(ixc^s)*w1r(ixc^s,mag(ip3))+r1r(ixc^s)*w1l(ixc^s,mag(ip3))+&
            r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip3))-w1l(ixc^s,mom(ip3)))*signbx(ixc^s))*tmp(ixc^s)
        w2l(ixc^s,mag(ip3))=w2r(ixc^s,mag(ip3))
      end if
      ! Miyoshi equation (63) and Guo equation (45)
      if(phys_energy) then
        w2r(ixc^s,e_)=w1r(ixc^s,e_)+r1r(ixc^s)*(sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)-&
          sum(w2r(ixc^s,mom(:))*w2r(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
        w2l(ixc^s,e_)=w1l(ixc^s,e_)-r1l(ixc^s)*(sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)-&
          sum(w2l(ixc^s,mom(:))*w2l(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
      end if
 
      ! convert velocity to momentum
      do idir=1,ndir
        w1r(ixc^s,mom(idir))=w1r(ixc^s,mom(idir))*w1r(ixc^s,rho_)
        w1l(ixc^s,mom(idir))=w1l(ixc^s,mom(idir))*w1l(ixc^s,rho_)
        w2r(ixc^s,mom(idir))=w2r(ixc^s,mom(idir))*w2r(ixc^s,rho_)
        w2l(ixc^s,mom(idir))=w2l(ixc^s,mom(idir))*w2l(ixc^s,rho_)
      end do
      if(iw_equi_rho>0) then
        w1r(ixc^s,rho_) = w1r(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
        w1l(ixc^s,rho_) = w1l(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
        w2r(ixc^s,rho_) = w2r(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
        w2l(ixc^s,rho_) = w2l(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
      endif
      ! get fluxes of intermedate states
      do iw=iws,iwe
        if(flux_type(idims, iw) == flux_special) then
          ! known flux (fLC=fRC) for normal B and psi_ in GLM method
          f1l(ixc^s,iw)=flc(ixc^s,iw)
          f1r(ixc^s,iw)=f1l(ixc^s,iw)
          f2l(ixc^s,iw)=f1l(ixc^s,iw)
          f2r(ixc^s,iw)=f1l(ixc^s,iw)
        else if(flux_type(idims, iw) == flux_hll) then
          ! using hll flux for tracers
          f1l(ixc^s,iw)=(sr(ixc^s,ii)*flc(ixc^s, iw)-sl(ixc^s,ii)*frc(ixc^s, iw) &
                    +sr(ixc^s,ii)*sl(ixc^s,ii)*(wrc(ixc^s,iw)-wlc(ixc^s,iw)))/(sr(ixc^s,ii)-sl(ixc^s,ii))
          f1r(ixc^s,iw)=f1l(ixc^s,iw)
          f2l(ixc^s,iw)=f1l(ixc^s,iw)
          f2r(ixc^s,iw)=f1l(ixc^s,iw)
        else
          ! construct hlld flux
          f1l(ixc^s,iw)=flc(ixc^s,iw)+sl(ixc^s,ii)*(w1l(ixc^s,iw)-wlc(ixc^s,iw))
          f1r(ixc^s,iw)=frc(ixc^s,iw)+sr(ixc^s,ii)*(w1r(ixc^s,iw)-wrc(ixc^s,iw))
          f2l(ixc^s,iw)=f1l(ixc^s,iw)+s1l(ixc^s)*(w2l(ixc^s,iw)-w1l(ixc^s,iw))
          f2r(ixc^s,iw)=f1r(ixc^s,iw)+s1r(ixc^s)*(w2r(ixc^s,iw)-w1r(ixc^s,iw))
        end if
      end do
 
      ! Miyoshi equation (66) and Guo equation (46)
     {do ix^db=ixcmin^db,ixcmax^db\}
        if(sl(ix^d,ii)>0.d0) then
          fc(ix^d,iws:iwe,ip1)=flc(ix^d,iws:iwe)
        else if(s1l(ix^d)>=0.d0) then
          fc(ix^d,iws:iwe,ip1)=f1l(ix^d,iws:iwe)
        else if(sm(ix^d)>=0.d0) then
          fc(ix^d,iws:iwe,ip1)=f2l(ix^d,iws:iwe)
        else if(s1r(ix^d)>=0.d0) then
          fc(ix^d,iws:iwe,ip1)=f2r(ix^d,iws:iwe)
        else if(sr(ix^d,ii)>=0.d0) then
          fc(ix^d,iws:iwe,ip1)=f1r(ix^d,iws:iwe)
        else if(sr(ix^d,ii)<0.d0) then
          fc(ix^d,iws:iwe,ip1)=frc(ix^d,iws:iwe)
        end if
     {end do\}
 
      end associate
    end subroutine get_riemann_flux_hlld
 
    !> HLLD Riemann flux from Miyoshi 2005 JCP, 208, 315 and Guo 2016 JCP, 327, 543
    !> https://arxiv.org/pdf/2108.04991.pdf
    subroutine get_riemann_flux_hlld_mag2()
      !use mod_mhd_phys
      use mod_variables
      use mod_physics
      implicit none
      double precision, dimension(ixI^S,1:nwflux) :: w1R,w1L,f1R,f1L,f2R,f2L
      double precision, dimension(ixI^S,1:nwflux) :: w2R,w2L
      double precision, dimension(ixI^S) :: sm,s1R,s1L,suR,suL,Bx
      double precision, dimension(ixI^S) :: pts,ptR,ptL,signBx,r1L,r1R,tmp
      ! velocity from the right and the left reconstruction
      double precision, dimension(ixI^S,ndir) :: vRC, vLC
      ! magnetic field from the right and the left reconstruction
      double precision, dimension(ixI^S,ndir) :: BR, BL
      integer :: ip1,ip2,ip3,idir,ix^D
      double precision :: phiPres, thetaSM, du, dv, dw
      integer :: ixV^L, ixVb^L, ixVc^L, ixVd^L, ixVe^L, ixVf^L
      integer  :: rho_, p_, e_, mom(1:ndir), mag(1:ndir)
      double precision, parameter :: aParam = 4d0
 
      rho_ = iw_rho
      mom(:) = iw_mom(:)
      mag(:) = iw_mag(:) 
      p_ = iw_e
      e_ = iw_e 
 
      associate(sr=>cmaxc,sl=>cminc)
 
      f1r=0.d0
      f1l=0.d0
      ip1=idims
      ip3=3
 
      vrc(ixc^s,:)=wrp(ixc^s,mom(:))
      vlc(ixc^s,:)=wlp(ixc^s,mom(:))
 
      ! reuse s1L s1R
      call get_hlld2_modif_c(wlp,x,ixi^l,ixo^l,s1l)
      call get_hlld2_modif_c(wrp,x,ixi^l,ixo^l,s1r)
      !phiPres = min(1, maxval(max(s1L(ixO^S),s1R(ixO^S))/cmaxC(ixO^S,1)))
      phipres = min(1d0, maxval(max(s1l(ixo^s),s1r(ixo^s)))/maxval(cmaxc(ixo^s,1)))
      phipres = phipres*(2d0 - phipres)  
 
     !we use here not reconstructed velocity: wprim?       
     ixv^l=ixo^l;
     !first dim
     ixvmin1=ixomin1+1  
     ixvmax1=ixomax1+1
     du = minval(wprim(ixv^s,mom(1))-wprim(ixo^s,mom(1)))
     if(du>0d0) du=0d0  
     dv = 0d0 
     dw = 0d0 
 
     {^nooned
     !second dim
     !i,j-1,k 
     ixv^l=ixo^l;
     ixvmin2=ixomin2-1  
     ixvmax2=ixomax2-1
 
     !i,j+1,k 
     ixvb^l=ixo^l;
     ixvbmin2=ixomin2+1  
     ixvbmax2=ixomax2+1
 
     !i+1,j,k 
     ixvc^l=ixo^l;
     ixvcmin1=ixomin1+1  
     ixvcmax1=ixomax1+1
 
     !i+1,j-1,k 
     ixvd^l=ixo^l;
     ixvdmin1=ixomin1+1  
     ixvdmax1=ixomax1+1
     ixvdmin2=ixomin2-1  
     ixvdmax2=ixomax2-1
 
     !i+1,j+1,k 
     ixve^l=ixo^l;
     ixvemin1=ixomin1+1  
     ixvemax1=ixomax1+1
     ixvemin2=ixomin2+1  
     ixvemax2=ixomax2+1
 
     dv = minval(min(wprim(ixo^s,mom(2))-wprim(ixv^s,mom(2)),&
                  wprim(ixvb^s,mom(2))-wprim(ixo^s,mom(2)),&
                  wprim(ixvc^s,mom(2))-wprim(ixvd^s,mom(2)),&
                  wprim(ixve^s,mom(2))-wprim(ixvc^s,mom(2))&
                ))
     if(dv>0d0) dv=0d0} 
 
      {^ifthreed
     !third dim
     !i,j,k-1 
     ixv^l=ixo^l;
     ixvmin3=ixomin3-1  
     ixvmax3=ixomax3-1
 
     !i,j,k+1 
     ixvb^l=ixo^l;
     ixvbmin3=ixomin3+1  
     ixvbmax3=ixomax3+1
 
     !i+1,j,k 
     ixvc^l=ixo^l;
     ixvcmin1=ixomin1+1  
     ixvcmax1=ixomax1+1
 
     !i+1,j,k-1 
     ixvd^l=ixo^l;
     ixvdmin1=ixomin1+1  
     ixvdmax1=ixomax1+1
     ixvdmin3=ixomin3-1  
     ixvdmax3=ixomax3-1
 
     !i+1,j,k+1 
     ixve^l=ixo^l;
     ixvemin1=ixomin1+1  
     ixvemax1=ixomax1+1
     ixvemin3=ixomin3+1  
     ixvemax3=ixomax3+1
     dw = minval(min(wprim(ixo^s,mom(3))-wprim(ixv^s,mom(3)),&
                  wprim(ixvb^s,mom(3))-wprim(ixo^s,mom(3)),&
                  wprim(ixvc^s,mom(3))-wprim(ixvd^s,mom(3)),&
                  wprim(ixve^s,mom(3))-wprim(ixvc^s,mom(3))&
                ))
     if(dw>0d0) dw=0d0}
     thetasm = maxval(cmaxc(ixo^s,1)) 
       
     thetasm = (min(1d0, (thetasm-du)/(thetasm-min(dv,dw))))**aparam
     !print*, "HLLD2 ", du,dv,dw, thetaSM, phiPres 
 
      if(b0field) then
        br(ixc^s,:)=wrc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
        bl(ixc^s,:)=wlc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
      else
        br(ixc^s,:)=wrc(ixc^s,mag(:))
        bl(ixc^s,:)=wlc(ixc^s,mag(:))
      end if
      ! HLL estimation of normal magnetic field at cell interfaces
      bx(ixc^s)=(sr(ixc^s,index_v_mag)*br(ixc^s,ip1)-sl(ixc^s,index_v_mag)*bl(ixc^s,ip1)-&
                 flc(ixc^s,mag(ip1))-frc(ixc^s,mag(ip1)))/(sr(ixc^s,index_v_mag)-sl(ixc^s,index_v_mag))
      ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(br(ixc^s,:)**2,dim=ndim+1)
      ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(bl(ixc^s,:)**2,dim=ndim+1)
      sur(ixc^s)=(sr(ixc^s,index_v_mag)-vrc(ixc^s,ip1))*wrc(ixc^s,rho_)
      sul(ixc^s)=(sl(ixc^s,index_v_mag)-vlc(ixc^s,ip1))*wlc(ixc^s,rho_)
      ! Miyoshi equation (38) and Guo euqation (20)
      sm(ixc^s)=(sur(ixc^s)*vrc(ixc^s,ip1)-sul(ixc^s)*vlc(ixc^s,ip1)-&
                 thetasm*(ptr(ixc^s)-ptl(ixc^s)) )/(sur(ixc^s)-sul(ixc^s))
      ! Miyoshi equation (39) and Guo euqation (28)
      w1r(ixc^s,mom(ip1))=sm(ixc^s)
      w1l(ixc^s,mom(ip1))=sm(ixc^s)
      w2r(ixc^s,mom(ip1))=sm(ixc^s)
      w2l(ixc^s,mom(ip1))=sm(ixc^s)
      ! Guo equation (22)
      w1r(ixc^s,mag(ip1))=bx(ixc^s)
      w1l(ixc^s,mag(ip1))=bx(ixc^s)
      if(b0field) then
        ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(wrc(ixc^s,mag(:))**2,dim=ndim+1)
        ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(wlc(ixc^s,mag(:))**2,dim=ndim+1)
      end if
 
      ! Miyoshi equation (43) and Guo equation (27)
      w1r(ixc^s,rho_)=sur(ixc^s)/(sr(ixc^s,index_v_mag)-sm(ixc^s))
      w1l(ixc^s,rho_)=sul(ixc^s)/(sl(ixc^s,index_v_mag)-sm(ixc^s))
 
      ip2=mod(ip1+1,ndir)
      if(ip2==0) ip2=ndir
      r1r(ixc^s)=sur(ixc^s)*(sr(ixc^s,index_v_mag)-sm(ixc^s))-bx(ixc^s)**2
      where(r1r(ixc^s)/=0.d0)
        r1r(ixc^s)=1.d0/r1r(ixc^s)
      endwhere
      r1l(ixc^s)=sul(ixc^s)*(sl(ixc^s,index_v_mag)-sm(ixc^s))-bx(ixc^s)**2
      where(r1l(ixc^s)/=0.d0)
        r1l(ixc^s)=1.d0/r1l(ixc^s)
      endwhere
      ! Miyoshi equation (44)
      w1r(ixc^s,mom(ip2))=vrc(ixc^s,ip2)-bx(ixc^s)*br(ixc^s,ip2)*&
        (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
      w1l(ixc^s,mom(ip2))=vlc(ixc^s,ip2)-bx(ixc^s)*bl(ixc^s,ip2)*&
        (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
      ! partial solution for later usage
      w1r(ixc^s,mag(ip2))=(sur(ixc^s)*(sr(ixc^s,index_v_mag)-vrc(ixc^s,ip1))-bx(ixc^s)**2)*r1r(ixc^s)
      w1l(ixc^s,mag(ip2))=(sul(ixc^s)*(sl(ixc^s,index_v_mag)-vlc(ixc^s,ip1))-bx(ixc^s)**2)*r1l(ixc^s)
      if(ndir==3) then
        ip3=mod(ip1+2,ndir)
        if(ip3==0) ip3=ndir
        ! Miyoshi equation (46)
        w1r(ixc^s,mom(ip3))=vrc(ixc^s,ip3)-bx(ixc^s)*br(ixc^s,ip3)*&
          (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
        w1l(ixc^s,mom(ip3))=vlc(ixc^s,ip3)-bx(ixc^s)*bl(ixc^s,ip3)*&
          (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
        ! Miyoshi equation (47)
        w1r(ixc^s,mag(ip3))=br(ixc^s,ip3)*w1r(ixc^s,mag(ip2))
        w1l(ixc^s,mag(ip3))=bl(ixc^s,ip3)*w1l(ixc^s,mag(ip2))
      end if
      ! Miyoshi equation (45)
      w1r(ixc^s,mag(ip2))=br(ixc^s,ip2)*w1r(ixc^s,mag(ip2))
      w1l(ixc^s,mag(ip2))=bl(ixc^s,ip2)*w1l(ixc^s,mag(ip2))
      if(b0field) then
        ! Guo equation (26)
        w1r(ixc^s,mag(:))=w1r(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
        w1l(ixc^s,mag(:))=w1l(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
      end if
      ! equation (48)
      if(phys_energy) then
        ! Guo equation (25) equivalent to Miyoshi equation (41)
        w1r(ixc^s,p_)=(sur(ixc^s)*ptl(ixc^s) - sul(ixc^s)*ptr(ixc^s) +&
          phipres * sur(ixc^s)*sul(ixc^s)*(vrc(ixc^s,ip1)-vlc(ixc^s,ip1)))/&
          (sur(ixc^s)-sul(ixc^s))
        w1l(ixc^s,p_)=w1r(ixc^s,p_)
        if(b0field) then
          ! Guo equation (32)
          w1r(ixc^s,p_)=w1r(ixc^s,p_)+sum(block%B0(ixc^s,:,ip1)*(wrc(ixc^s,mag(:))-w1r(ixc^s,mag(:))),dim=ndim+1)
          w1l(ixc^s,p_)=w1l(ixc^s,p_)+sum(block%B0(ixc^s,:,ip1)*(wlc(ixc^s,mag(:))-w1l(ixc^s,mag(:))),dim=ndim+1)
        end if
        ! Miyoshi equation (48) and main part of Guo euqation (31)
        w1r(ixc^s,e_)=((sr(ixc^s,index_v_mag)-vrc(ixc^s,ip1))*wrc(ixc^s,e_)-ptr(ixc^s)*vrc(ixc^s,ip1)+&
          w1r(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vrc(ixc^s,:)*wrc(ixc^s,mag(:)),dim=ndim+1)-&
          sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)))/(sr(ixc^s,index_v_mag)-sm(ixc^s))
        w1l(ixc^s,e_)=((sl(ixc^s,index_v_mag)-vlc(ixc^s,ip1))*wlc(ixc^s,e_)-ptl(ixc^s)*vlc(ixc^s,ip1)+&
          w1l(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vlc(ixc^s,:)*wlc(ixc^s,mag(:)),dim=ndim+1)-&
          sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)))/(sl(ixc^s,index_v_mag)-sm(ixc^s))
        if(b0field) then
          ! Guo equation (31)
          w1r(ixc^s,e_)=w1r(ixc^s,e_)+(sum(w1r(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
               sum(wrc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vrc(ixc^s,ip1))/(sr(ixc^s,index_v_mag)-sm(ixc^s))
          w1l(ixc^s,e_)=w1l(ixc^s,e_)+(sum(w1l(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
               sum(wlc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vlc(ixc^s,ip1))/(sl(ixc^s,index_v_mag)-sm(ixc^s))
        end if
      end if
 
      ! Miyoshi equation (49) and Guo equation (35)
      w2r(ixc^s,rho_)=w1r(ixc^s,rho_)
      w2l(ixc^s,rho_)=w1l(ixc^s,rho_)
      w2r(ixc^s,mag(ip1))=w1r(ixc^s,mag(ip1))
      w2l(ixc^s,mag(ip1))=w1l(ixc^s,mag(ip1))
 
      r1r(ixc^s)=sqrt(w1r(ixc^s,rho_))
      r1l(ixc^s)=sqrt(w1l(ixc^s,rho_))
      tmp(ixc^s)=1.d0/(r1r(ixc^s)+r1l(ixc^s))
      signbx(ixc^s)=sign(1.d0,bx(ixc^s))
      ! Miyoshi equation (51) and Guo equation (33)
      s1r(ixc^s)=sm(ixc^s)+abs(bx(ixc^s))/r1r(ixc^s)
      s1l(ixc^s)=sm(ixc^s)-abs(bx(ixc^s))/r1l(ixc^s)
      ! Miyoshi equation (59) and Guo equation (41)
      w2r(ixc^s,mom(ip2))=(r1l(ixc^s)*w1l(ixc^s,mom(ip2))+r1r(ixc^s)*w1r(ixc^s,mom(ip2))+&
          (w1r(ixc^s,mag(ip2))-w1l(ixc^s,mag(ip2)))*signbx(ixc^s))*tmp(ixc^s)
      w2l(ixc^s,mom(ip2))=w2r(ixc^s,mom(ip2))
      ! Miyoshi equation (61) and Guo equation (43)
      w2r(ixc^s,mag(ip2))=(r1l(ixc^s)*w1r(ixc^s,mag(ip2))+r1r(ixc^s)*w1l(ixc^s,mag(ip2))+&
          r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip2))-w1l(ixc^s,mom(ip2)))*signbx(ixc^s))*tmp(ixc^s)
      w2l(ixc^s,mag(ip2))=w2r(ixc^s,mag(ip2))
      if(ndir==3) then
        ! Miyoshi equation (60) and Guo equation (42)
        w2r(ixc^s,mom(ip3))=(r1l(ixc^s)*w1l(ixc^s,mom(ip3))+r1r(ixc^s)*w1r(ixc^s,mom(ip3))+&
            (w1r(ixc^s,mag(ip3))-w1l(ixc^s,mag(ip3)))*signbx(ixc^s))*tmp(ixc^s)
        w2l(ixc^s,mom(ip3))=w2r(ixc^s,mom(ip3))
        ! Miyoshi equation (62) and Guo equation (44)
        w2r(ixc^s,mag(ip3))=(r1l(ixc^s)*w1r(ixc^s,mag(ip3))+r1r(ixc^s)*w1l(ixc^s,mag(ip3))+&
            r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip3))-w1l(ixc^s,mom(ip3)))*signbx(ixc^s))*tmp(ixc^s)
        w2l(ixc^s,mag(ip3))=w2r(ixc^s,mag(ip3))
      end if
      ! Miyoshi equation (63) and Guo equation (45)
      if(phys_energy) then
        w2r(ixc^s,e_)=w1r(ixc^s,e_)+r1r(ixc^s)*(sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)-&
          sum(w2r(ixc^s,mom(:))*w2r(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
        w2l(ixc^s,e_)=w1l(ixc^s,e_)-r1l(ixc^s)*(sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)-&
          sum(w2l(ixc^s,mom(:))*w2l(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
      end if
 
      ! convert velocity to momentum
      do idir=1,ndir
        w1r(ixc^s,mom(idir))=w1r(ixc^s,mom(idir))*w1r(ixc^s,rho_)
        w1l(ixc^s,mom(idir))=w1l(ixc^s,mom(idir))*w1l(ixc^s,rho_)
        w2r(ixc^s,mom(idir))=w2r(ixc^s,mom(idir))*w2r(ixc^s,rho_)
        w2l(ixc^s,mom(idir))=w2l(ixc^s,mom(idir))*w2l(ixc^s,rho_)
      end do
 
      ! get fluxes of intermedate states
      do iw=1,nwflux
        ! CT MHD does not need normal B flux
        if(stagger_grid .and. flux_type(idims, iw) == flux_tvdlf) cycle
        if(flux_type(idims, iw) == flux_special) then
          ! known flux (fLC=fRC) for normal B and psi_ in GLM method
          f1l(ixc^s,iw)=flc(ixc^s,iw)
          f1r(ixc^s,iw)=f1l(ixc^s,iw)
          f2l(ixc^s,iw)=f1l(ixc^s,iw)
          f2r(ixc^s,iw)=f1l(ixc^s,iw)
        else if(flux_type(idims, iw) == flux_hll) then
          ! using hll flux for tracers
          f1l(ixc^s,iw)=(sr(ixc^s,index_v_mag)*flc(ixc^s, iw)-sl(ixc^s,index_v_mag)*frc(ixc^s, iw) &
                    +sr(ixc^s,index_v_mag)*sl(ixc^s,index_v_mag)*(wrc(ixc^s,iw)-wlc(ixc^s,iw)))/(sr(ixc^s,index_v_mag)-sl(ixc^s,index_v_mag))
          f1r(ixc^s,iw)=f1l(ixc^s,iw)
          f2l(ixc^s,iw)=f1l(ixc^s,iw)
          f2r(ixc^s,iw)=f1l(ixc^s,iw)
        else
          f1l(ixc^s,iw)=flc(ixc^s,iw)+sl(ixc^s,index_v_mag)*(w1l(ixc^s,iw)-wlc(ixc^s,iw))
          f1r(ixc^s,iw)=frc(ixc^s,iw)+sr(ixc^s,index_v_mag)*(w1r(ixc^s,iw)-wrc(ixc^s,iw))
          f2l(ixc^s,iw)=f1l(ixc^s,iw)+s1l(ixc^s)*(w2l(ixc^s,iw)-w1l(ixc^s,iw))
          f2r(ixc^s,iw)=f1r(ixc^s,iw)+s1r(ixc^s)*(w2r(ixc^s,iw)-w1r(ixc^s,iw))
        end if
      end do
 
      ! Miyoshi equation (66) and Guo equation (46)
     {do ix^db=ixcmin^db,ixcmax^db\}
        if(sl(ix^d,index_v_mag)>0.d0) then
          fc(ix^d,1:nwflux,ip1)=flc(ix^d,1:nwflux)
        else if(s1l(ix^d)>=0.d0) then
          fc(ix^d,1:nwflux,ip1)=f1l(ix^d,1:nwflux)
        else if(sm(ix^d)>=0.d0) then
          fc(ix^d,1:nwflux,ip1)=f2l(ix^d,1:nwflux)
        else if(s1r(ix^d)>=0.d0) then
          fc(ix^d,1:nwflux,ip1)=f2r(ix^d,1:nwflux)
        else if(sr(ix^d,index_v_mag)>=0.d0) then
          fc(ix^d,1:nwflux,ip1)=f1r(ix^d,1:nwflux)
        else if(sr(ix^d,index_v_mag)<0.d0) then
          fc(ix^d,1:nwflux,ip1)=frc(ix^d,1:nwflux)
        end if
     {end do\}
 
      end associate
    end subroutine get_riemann_flux_hlld_mag2
 
    !> Calculate fast magnetosonic wave speed
    subroutine get_hlld2_modif_c(w,x,ixI^L,ixO^L,csound)
      use mod_global_parameters
 
      integer, intent(in)          :: ixI^L, ixO^L
      double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
      double precision, intent(out):: csound(ixI^S)
      double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
      double precision :: inv_rho(ixO^S), gamma_A2(ixO^S)
      integer  :: rho_, p_, e_, mom(1:ndir), mag(1:ndir)
 
        rho_ = iw_rho
        mom(:) = iw_mom(:)
        mag(:) = iw_mag(:) 
        p_ = iw_e
        e_ = iw_e 
 
      inv_rho=1.d0/w(ixo^s,rho_)
 
      ! store |B|^2 in v
 
      if (b0field) then
        b2(ixo^s) = sum((w(ixo^s, mag(:))+block%B0(ixo^s,:,b0i))**2, dim=ndim+1)
      else
        b2(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=ndim+1)
      end if
 
 
      if (b0field) then
        avmincs2= w(ixo^s, mag(idims))+block%B0(ixo^s,idims,b0i)
      else
        avmincs2= w(ixo^s, mag(idims))
      end if
 
 
      csound(ixo^s) = sum(w(ixo^s, mom(:))**2, dim=ndim+1)
 
      cfast2(ixo^s)   = b2(ixo^s) * inv_rho+csound(ixo^s)
      avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
           * avmincs2(ixo^s)**2 &
           * inv_rho 
 
      where(avmincs2(ixo^s)<zero)
         avmincs2(ixo^s)=zero
      end where
 
      avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
 
     csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
 
    end subroutine get_hlld2_modif_c
 
  end subroutine finite_volume
 
  !> Determine the upwinded wLC(ixL) and wRC(ixR) from w.
  !> the wCT is only used when PPM is exploited.
  subroutine reconstruct_lr(ixI^L,ixL^L,ixR^L,idims,w,wLC,wRC,wLp,wRp,x,dxdim)
    use mod_physics
    use mod_global_parameters
    use mod_limiter
    use mod_comm_lib, only: mpistop
 
    integer, intent(in) :: ixi^l, ixl^l, ixr^l, idims
    double precision, intent(in) :: dxdim
    ! cell center w in primitive form
    double precision, dimension(ixI^S,1:nw) :: w
    ! left and right constructed status in conservative form
    double precision, dimension(ixI^S,1:nw) :: wlc, wrc
    ! left and right constructed status in primitive form
    double precision, dimension(ixI^S,1:nw) :: wlp, wrp
    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)
    double precision   :: a2max
 
    select case (type_limiter(block%level))
    case (limiter_mp5)
       call mp5limiter(ixi^l,ixl^l,idims,w,wlp,wrp)
    case (limiter_weno3)
       call weno3limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,1)
    case (limiter_wenoyc3)
       call weno3limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,2)
    case (limiter_weno5)
       call weno5limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,1)
    case (limiter_weno5nm)
       call weno5nmlimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,1)
    case (limiter_wenoz5)
       call weno5limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,2)
    case (limiter_wenoz5nm)
       call weno5nmlimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,2)
    case (limiter_wenozp5)
       call weno5limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,3)
    case (limiter_wenozp5nm)
       call weno5nmlimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,3)
    case (limiter_weno5cu6)
       call weno5cu6limiter(ixi^l,ixl^l,idims,w,wlp,wrp)
    case (limiter_teno5ad)
       call teno5adlimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp)
    case (limiter_weno7)
       call weno7limiter(ixi^l,ixl^l,idims,w,wlp,wrp,1)
    case (limiter_mpweno7)
       call weno7limiter(ixi^l,ixl^l,idims,w,wlp,wrp,2)
    case (limiter_venk)
       call venklimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp) 
       if(fix_small_values) then
          call phys_handle_small_values(.true.,wlp,x,ixi^l,ixl^l,'reconstruct left')
          call phys_handle_small_values(.true.,wrp,x,ixi^l,ixr^l,'reconstruct right')
       end if
    case (limiter_ppm)
       call ppmlimiter(ixi^l,ixm^ll,idims,w,w,wlp,wrp)
       if(fix_small_values) then
          call phys_handle_small_values(.true.,wlp,x,ixi^l,ixl^l,'reconstruct left')
          call phys_handle_small_values(.true.,wrp,x,ixi^l,ixr^l,'reconstruct right')
       end if
    case default
       jxr^l=ixr^l+kr(idims,^d);
       ixcmax^d=jxrmax^d; ixcmin^d=ixlmin^d-kr(idims,^d);
       jxc^l=ixc^l+kr(idims,^d);
       do iw=1,nwflux
          if (loglimit(iw)) then
             w(ixcmin^d:jxcmax^d,iw)=dlog10(w(ixcmin^d:jxcmax^d,iw))
             wlp(ixl^s,iw)=dlog10(wlp(ixl^s,iw))
             wrp(ixr^s,iw)=dlog10(wrp(ixr^s,iw))
          end if
 
          dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
          if(need_global_a2max) then 
            a2max=a2max_global(idims)
          else
            select case(idims)
            case(1)
              a2max=schmid_rad1
            {^iftwod
            case(2)
              a2max=schmid_rad2}
            {^ifthreed
            case(2)
              a2max=schmid_rad2
            case(3)
              a2max=schmid_rad3}
            case default
              call mpistop("idims is wrong in mod_limiter")
            end select
          end if
 
          ! limit flux from left and/or right
          call dwlimiter2(dwc,ixi^l,ixc^l,idims,type_limiter(block%level),ldw,rdw,a2max=a2max)
          wlp(ixl^s,iw)=wlp(ixl^s,iw)+half*ldw(ixl^s)
          wrp(ixr^s,iw)=wrp(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)
             wlp(ixl^s,iw)=10.0d0**wlp(ixl^s,iw)
             wrp(ixr^s,iw)=10.0d0**wrp(ixr^s,iw)
          end if
       end do
       if(fix_small_values) then
          call phys_handle_small_values(.true.,wlp,x,ixi^l,ixl^l,'reconstruct left')
          call phys_handle_small_values(.true.,wrp,x,ixi^l,ixr^l,'reconstruct right')
       end if
    end select
 
   wlc(ixl^s,1:nwflux) = wlp(ixl^s,1:nwflux)
   wrc(ixr^s,1:nwflux) = wrp(ixr^s,1:nwflux)
   call phys_to_conserved(ixi^l,ixl^l,wlc,x)
   call phys_to_conserved(ixi^l,ixr^l,wrc,x)
   if(nwaux>0)then
      wlp(ixl^s,nwflux+1:nwflux+nwaux) = wlc(ixl^s,nwflux+1:nwflux+nwaux)
      wrp(ixr^s,nwflux+1:nwflux+nwaux) = wrc(ixr^s,nwflux+1:nwflux+nwaux)
   endif
 
  end subroutine reconstruct_lr
 
end module mod_finite_volume