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")
 
    wrp=0.d0
    wlp=0.d0
    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:nw_recon)=wprim(ixo^s,1:nw_recon)
      wlp(ixo^s,1:nw_recon)=wprim(ixo^s,1:nw_recon)
 
      ! 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,wprim,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
    double precision, dimension(ixI^S,1:nwflux,1:ndim)    :: fc
    double precision, dimension(ixI^S,sdim:3)             :: fe
    type(state)                                           :: sct, snew
 
    ! 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 :: idims, iw, ix^d, hx^d, 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)
 
    ! 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
 
       kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
       kxr^l=kxc^l+kr(idims,^d);
       ! 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.
       do iw=1,nw_recon
        {do ix^db=iximin^db,iximax^db\}
           ! fill all cells for averaging to fix small values
           wrp(ix^d,iw)=wprim(ix^d,iw)
           wlp(ix^d,iw)=wprim(ix^d,iw)
        {end do\}
         wrp(kxc^s,iw)=wprim(kxr^s,iw)
       end do
 
       hxo^l=ixo^l-kr(idims,^d);
       if(stagger_grid) then
         ! ct needs all transverse cells
         ixcmax^d=ixomax^d+nghostcells-nghostcells*kr(idims,^d);
         ixcmin^d=hxomin^d-nghostcells+nghostcells*kr(idims,^d);
       else
         ! ixC is centered index in the idims direction from ixOmin-1/2 to ixOmax+1/2
         ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
       end if
 
       ! 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
      if(local_timestep) then
        dxinv(1:ndim)=-dtfactor/dxs(1:ndim)
        do idims= idims^lim
          hx^d=kr(idims,^d)\
          hxomin^d=ixomin^d-hx^d\
          do iw=iwstart,nwflux
           {do ix^db=hxomin^db,ixomax^db\}
              fc(ix^d,iw,idims)=block%dt(ix^d)*dxinv(idims)*fc(ix^d,iw,idims)
           {end do\}
           {do ix^db=ixomin^db,ixomax^db\}
              wnew(ix^d,iw)=wnew(ix^d,iw)+fc(ix^d,iw,idims)-fc(ix^d-hx^d,iw,idims)
           {end do\}
          end do
          ! 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
      else
        dxinv(1:ndim)=-qdt/dxs(1:ndim)
        do idims= idims^lim
          hx^d=kr(idims,^d)\
          hxomin^d=ixomin^d-hx^d\
          do iw=iwstart,nwflux
           {do ix^db=hxomin^db,ixomax^db\}
              fc(ix^d,iw,idims)=dxinv(idims)*fc(ix^d,iw,idims)
           {end do\}
           {do ix^db=ixomin^db,ixomax^db\}
              wnew(ix^d,iw)=wnew(ix^d,iw)+fc(ix^d,iw,idims)-fc(ix^d-hx^d,iw,idims)
           {end do\}
          end do
          ! 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
      end if
    else
      inv_volume(ixo^s) = 1.d0/block%dvolume(ixo^s)
      if(local_timestep) then
        do idims= idims^lim
          hx^d=kr(idims,^d)\
          hxomin^d=ixomin^d-hx^d\
          do iw=iwstart,nwflux
           {do ix^db=hxomin^db,ixomax^db\}
              fc(ix^d,iw,idims)=-block%dt(ix^d)*dtfactor*fc(ix^d,iw,idims)*block%surfaceC(ix^d,idims)
           {end do\}
           {do ix^db=ixomin^db,ixomax^db\}
              wnew(ix^d,iw)=wnew(ix^d,iw)+(fc(ix^d,iw,idims)-fc(ix^d-hx^d,iw,idims))*inv_volume(ix^d)
           {end do\}
          end do
          ! 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
      else
        do idims= idims^lim
          hx^d=kr(idims,^d)\
          hxomin^d=ixomin^d-hx^d\
          do iw=iwstart,nwflux
           {do ix^db=hxomin^db,ixomax^db\}
             fc(ix^d,iw,idims)=-qdt*fc(ix^d,iw,idims)*block%surfaceC(ix^d,idims)
           {end do\}
           {do ix^db=ixomin^db,ixomax^db\}
              wnew(ix^d,iw)=wnew(ix^d,iw)+(fc(ix^d,iw,idims)-fc(ix^d-hx^d,iw,idims))*inv_volume(ix^d)
           {end do\}
          end do
        end do
      end if
    end if
 
    if (.not.slab.and.idimsmin==1) &
         call phys_add_source_geom(qdt,dtfactor,ixi^l,ixo^l,wct,wprim,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
        fc(ixc^s,iw,idims)=half*(flc(ixc^s,iw)+frc(ixc^s,iw))
      end do
    end subroutine get_riemann_flux_tvdmu
 
    subroutine get_riemann_flux_tvdlf(iws,iwe)
      integer, intent(in) :: iws,iwe
 
      integer :: ix^D,jx^D
      double precision :: fac(ixC^S),phi
 
      fac(ixc^s) = -0.5d0*tvdlfeps*cmaxc(ixc^s,ii)
      do iw=iws,iwe
         fc(ixc^s,iw,idims)=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
           if(flux_adaptive_diffusion) then
            {do ix^db=ixcmin^db,ixcmax^db\}
               jx^d=ix^d+kr(idims,^d)\
               !> adaptive diffusion from Rempel et al. 2009, see also Rempel et al. 2014
               !> the previous version is adopt
               if(((wrc(ix^d,iw)-wlc(ix^d,iw))*(sct%w(jx^d,iw)-sct%w(ix^d,iw))) .gt. 1.e-18) then
                 phi=min((wrc(ix^d,iw)-wlc(ix^d,iw))**2/((sct%w(jx^d,iw)-sct%w(ix^d,iw))**2+1.e-18),one)
               else
                 phi=1.d0
               end if
               fc(ix^d,iw,idims)=fc(ix^d,iw,idims)+fac(ix^d)*(wrc(ix^d,iw)-wlc(ix^d,iw))*phi
            {end do\}
           else
            {do ix^db=ixcmin^db,ixcmax^db\}
               fc(ix^d,iw,idims)=fc(ix^d,iw,idims)+fac(ix^d)*(wrc(ix^d,iw)-wlc(ix^d,iw))
            {end do\}
           end if
         end if
      end do
 
    end subroutine get_riemann_flux_tvdlf
 
    subroutine get_riemann_flux_hll(iws,iwe)
      integer, intent(in) :: iws,iwe
      integer :: ix^D
      double precision :: phi
 
      if(flux_adaptive_diffusion) then
        do iw=iws,iwe
          if(flux_type(idims, iw) == flux_tvdlf) then
            if(stagger_grid) then
              ! CT MHD set zero normal B flux
              fc(ixc^s,iw,idims)=0.d0
            else
              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))
            end if
          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
                !> reduced diffusion is from Wang et al. 2024
                phi=max(abs(cmaxc(ix^d,ii)),abs(cminc(ix^d,ii)))/(cmaxc(ix^d,ii)-cminc(ix^d,ii))
                fc(ix^d,iw,idims)=(cmaxc(ix^d,ii)*flc(ix^d, iw)-cminc(ix^d,ii)*frc(ix^d,iw)&
                      +phi*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\}
          end if
        end do
      else
        do iw=iws,iwe
          if(flux_type(idims, iw) == flux_tvdlf) then
            if(stagger_grid) then
              ! CT MHD set zero normal B flux
              fc(ixc^s,iw,idims)=0.d0
            else
              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))
            end if
          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
                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\}
          end if
        end do
      end if
    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, dimension(ixI^S) :: patchf
      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,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
      ! magnetic field from the right and the left reconstruction
      double precision, dimension(ixI^S,ndir) :: BR, BL
      integer :: ip1,ip2,ip3,idir,ix^D,^C&b^C_,^C&m^C_
      integer  :: rho_, p_, e_, mom(1:ndir), mag(1:ndir)
 
      associate(sr=>cmaxc,sl=>cminc)
 
      rho_=iw_rho
      ^c&mom(^c)=iw_mom(^c)\
      m^c_=mom(^c);
      ^c&mag(^c)=iw_mag(^c)\
      b^c_=mag(^c);
      e_ = iw_e 
      p_ = e_
 
      ip1=idims
      ip3=3
      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
     {!DEC$ VECTOR ALWAYS
      do ix^db=ixcmin^db,ixcmax^db\}
        ptr(ix^d)=wrp(ix^d,p_)+0.5d0*(^c&br(ix^d,^c)**2+)
        ptl(ix^d)=wlp(ix^d,p_)+0.5d0*(^c&bl(ix^d,^c)**2+)
        if(iw_equi_rho>0) then
          sur(ix^d)=(sr(ix^d,ii)-wrp(ix^d,mom(ip1)))*(wrc(ix^d,rho_)+block%equi_vars(ix^d,iw_equi_rho,ip1))
          sul(ix^d)=(sl(ix^d,ii)-wlp(ix^d,mom(ip1)))*(wlc(ix^d,rho_)+block%equi_vars(ix^d,iw_equi_rho,ip1))
        else
          sur(ix^d)=(sr(ix^d,ii)-wrp(ix^d,mom(ip1)))*wrc(ix^d,rho_)
          sul(ix^d)=(sl(ix^d,ii)-wlp(ix^d,mom(ip1)))*wlc(ix^d,rho_)
        end if
        ! Miyoshi equation (38) and Guo euqation (20)
        sm(ix^d)=(sur(ix^d)*wrp(ix^d,mom(ip1))-sul(ix^d)*wlp(ix^d,mom(ip1))-&
                   ptr(ix^d)+ptl(ix^d))/(sur(ix^d)-sul(ix^d))
        ! Miyoshi equation (39) and Guo euqation (28)
        w1r(ix^d,mom(ip1))=sm(ix^d)
        w1l(ix^d,mom(ip1))=sm(ix^d)
        w2r(ix^d,mom(ip1))=sm(ix^d)
        w2l(ix^d,mom(ip1))=sm(ix^d)
        ! Guo equation (22)
        w1r(ix^d,mag(ip1))=bx(ix^d)
        w1l(ix^d,mag(ip1))=bx(ix^d)
        if(b0field) then
          ptr(ix^d)=wrp(ix^d,p_)+0.5d0*(^c&wrc(ix^d,b^c_)**2+)
          ptl(ix^d)=wlp(ix^d,p_)+0.5d0*(^c&wlc(ix^d,b^c_)**2+)
        end if
        ! Miyoshi equation (43) and Guo equation (27)
        w1r(ix^d,rho_)=sur(ix^d)/(sr(ix^d,ii)-sm(ix^d))
        w1l(ix^d,rho_)=sul(ix^d)/(sl(ix^d,ii)-sm(ix^d))
        ip2=mod(ip1+1,ndir)
        if(ip2==0) ip2=ndir
        r1r(ix^d)=sur(ix^d)*(sr(ix^d,ii)-sm(ix^d))-bx(ix^d)**2
        if(r1r(ix^d)/=0.d0) r1r(ix^d)=1.d0/r1r(ix^d)
        r1l(ix^d)=sul(ix^d)*(sl(ix^d,ii)-sm(ix^d))-bx(ix^d)**2
        if(r1l(ix^d)/=0.d0) r1l(ix^d)=1.d0/r1l(ix^d)
        ! Miyoshi equation (44)
        w1r(ix^d,mom(ip2))=wrp(ix^d,mom(ip2))-bx(ix^d)*br(ix^d,ip2)*&
          (sm(ix^d)-wrp(ix^d,mom(ip1)))*r1r(ix^d)
        w1l(ix^d,mom(ip2))=wlp(ix^d,mom(ip2))-bx(ix^d)*bl(ix^d,ip2)*&
          (sm(ix^d)-wlp(ix^d,mom(ip1)))*r1l(ix^d)
        ! partial solution for later usage
        w1r(ix^d,mag(ip2))=(sur(ix^d)*(sr(ix^d,ii)-wrp(ix^d,mom(ip1)))-bx(ix^d)**2)*r1r(ix^d)
        w1l(ix^d,mag(ip2))=(sul(ix^d)*(sl(ix^d,ii)-wlp(ix^d,mom(ip1)))-bx(ix^d)**2)*r1l(ix^d)
        {^ifthreec
        ip3=mod(ip1+2,ndir)
        if(ip3==0) ip3=ndir
        ! Miyoshi equation (46)
        w1r(ix^d,mom(ip3))=wrp(ix^d,mom(ip3))-bx(ix^d)*br(ix^d,ip3)*&
          (sm(ix^d)-wrp(ix^d,mom(ip1)))*r1r(ix^d)
        w1l(ix^d,mom(ip3))=wlp(ix^d,mom(ip3))-bx(ix^d)*bl(ix^d,ip3)*&
          (sm(ix^d)-wlp(ix^d,mom(ip1)))*r1l(ix^d)
        ! Miyoshi equation (47)
        w1r(ix^d,mag(ip3))=br(ix^d,ip3)*w1r(ix^d,mag(ip2))
        w1l(ix^d,mag(ip3))=bl(ix^d,ip3)*w1l(ix^d,mag(ip2))
        }
        ! Miyoshi equation (45)
        w1r(ix^d,mag(ip2))=br(ix^d,ip2)*w1r(ix^d,mag(ip2))
        w1l(ix^d,mag(ip2))=bl(ix^d,ip2)*w1l(ix^d,mag(ip2))
        if(b0field) then
          ! Guo equation (26)
          ^c&w1r(ix^d,b^c_)=w1r(ix^d,b^c_)-block%B0(ix^d,^c,ip1)\
          ^c&w1l(ix^d,b^c_)=w1l(ix^d,b^c_)-block%B0(ix^d,^c,ip1)\
        end if
        ! equation (48)
        if(phys_energy) then
          ! Guo equation (25) equivalent to Miyoshi equation (41)
          w1r(ix^d,p_)=sur(ix^d)*(sm(ix^d)-wrp(ix^d,mom(ip1)))+ptr(ix^d)
          w1l(ix^d,p_)=w1r(ix^d,p_)
          if(b0field) then
            ! Guo equation (32)
            w1r(ix^d,p_)=w1r(ix^d,p_)+(^c&block%B0(ix^d,^c,ip1)*(wrc(ix^d,b^c_)-w1r(ix^d,b^c_))+)
            w1l(ix^d,p_)=w1l(ix^d,p_)+(^c&block%B0(ix^d,^c,ip1)*(wlc(ix^d,b^c_)-w1l(ix^d,b^c_))+)
          end if
          ! Miyoshi equation (48) and main part of Guo euqation (31)
          w1r(ix^d,e_)=((sr(ix^d,ii)-wrp(ix^d,mom(ip1)))*wrc(ix^d,e_)-ptr(ix^d)*wrp(ix^d,mom(ip1))+&
            w1r(ix^d,p_)*sm(ix^d)+bx(ix^d)*((^c&wrp(ix^d,m^c_)*wrc(ix^d,b^c_)+)-&
            (^c&w1r(ix^d,m^c_)*w1r(ix^d,b^c_)+)))/(sr(ix^d,ii)-sm(ix^d))
          w1l(ix^d,e_)=((sl(ix^d,ii)-wlp(ix^d,mom(ip1)))*wlc(ix^d,e_)-ptl(ix^d)*wlp(ix^d,mom(ip1))+&
            w1l(ix^d,p_)*sm(ix^d)+bx(ix^d)*((^c&wlp(ix^d,m^c_)*wlc(ix^d,b^c_)+)-&
            (^c&w1l(ix^d,m^c_)*w1l(ix^d,b^c_)+)))/(sl(ix^d,ii)-sm(ix^d))
          if(b0field) then
            ! Guo equation (31)
            w1r(ix^d,e_)=w1r(ix^d,e_)+((^c&w1r(ix^d,b^c_)*block%B0(ix^d,^c,ip1)+)*sm(ix^d)-&
                 (^c&wrc(ix^d,b^c_)*block%B0(ix^d,^c,ip1)+)*wrp(ix^d,mom(ip1)))/(sr(ix^d,ii)-sm(ix^d))
            w1l(ix^d,e_)=w1l(ix^d,e_)+((^c&w1l(ix^d,b^c_)*block%B0(ix^d,^c,ip1)+)*sm(ix^d)-&
                 (^c&wlc(ix^d,b^c_)*block%B0(ix^d,^c,ip1)+)*wlp(ix^d,mom(ip1)))/(sl(ix^d,ii)-sm(ix^d))
          end if
          if(iw_equi_p>0) then
            w1r(ix^d,e_)=w1r(ix^d,e_)+1d0/(phys_gamma-1)*block%equi_vars(ix^d,iw_equi_p,ip1)*&
               (sm(ix^d)-wrp(ix^d,mom(ip1)))/(sr(ix^d,ii)-sm(ix^d))
            w1l(ix^d,e_)=w1l(ix^d,e_)+1d0/(phys_gamma-1)*block%equi_vars(ix^d,iw_equi_p,ip1)*&
               (sm(ix^d)-wlp(ix^d,mom(ip1)))/(sl(ix^d,ii)-sm(ix^d))
          end if
        end if
 
        ! Miyoshi equation (49) and Guo equation (35)
        w2r(ix^d,rho_)=w1r(ix^d,rho_)
        w2l(ix^d,rho_)=w1l(ix^d,rho_)
        w2r(ix^d,mag(ip1))=w1r(ix^d,mag(ip1))
        w2l(ix^d,mag(ip1))=w1l(ix^d,mag(ip1))
        r1r(ix^d)=sqrt(w1r(ix^d,rho_))
        r1l(ix^d)=sqrt(w1l(ix^d,rho_))
        tmp(ix^d)=1.d0/(r1r(ix^d)+r1l(ix^d))
        signbx(ix^d)=sign(1.d0,bx(ix^d))
        ! Miyoshi equation (51) and Guo equation (33)
        s1r(ix^d)=sm(ix^d)+abs(bx(ix^d))/r1r(ix^d)
        s1l(ix^d)=sm(ix^d)-abs(bx(ix^d))/r1l(ix^d)
        ! Miyoshi equation (59) and Guo equation (41)
        w2r(ix^d,mom(ip2))=(r1l(ix^d)*w1l(ix^d,mom(ip2))+r1r(ix^d)*w1r(ix^d,mom(ip2))+&
            (w1r(ix^d,mag(ip2))-w1l(ix^d,mag(ip2)))*signbx(ix^d))*tmp(ix^d)
        w2l(ix^d,mom(ip2))=w2r(ix^d,mom(ip2))
        ! Miyoshi equation (61) and Guo equation (43)
        w2r(ix^d,mag(ip2))=(r1l(ix^d)*w1r(ix^d,mag(ip2))+r1r(ix^d)*w1l(ix^d,mag(ip2))+&
            r1l(ix^d)*r1r(ix^d)*(w1r(ix^d,mom(ip2))-w1l(ix^d,mom(ip2)))*signbx(ix^d))*tmp(ix^d)
        w2l(ix^d,mag(ip2))=w2r(ix^d,mag(ip2))
        {^ifthreec
        ! Miyoshi equation (60) and Guo equation (42)
        w2r(ix^d,mom(ip3))=(r1l(ix^d)*w1l(ix^d,mom(ip3))+r1r(ix^d)*w1r(ix^d,mom(ip3))+&
            (w1r(ix^d,mag(ip3))-w1l(ix^d,mag(ip3)))*signbx(ix^d))*tmp(ix^d)
        w2l(ix^d,mom(ip3))=w2r(ix^d,mom(ip3))
        ! Miyoshi equation (62) and Guo equation (44)
        w2r(ix^d,mag(ip3))=(r1l(ix^d)*w1r(ix^d,mag(ip3))+r1r(ix^d)*w1l(ix^d,mag(ip3))+&
            r1l(ix^d)*r1r(ix^d)*(w1r(ix^d,mom(ip3))-w1l(ix^d,mom(ip3)))*signbx(ix^d))*tmp(ix^d)
        w2l(ix^d,mag(ip3))=w2r(ix^d,mag(ip3))
        }
        ! Miyoshi equation (63) and Guo equation (45)
        if(phys_energy) then
          w2r(ix^d,e_)=w1r(ix^d,e_)+r1r(ix^d)*((^c&w1r(ix^d,m^c_)*w1r(ix^d,b^c_)+)-&
            (^c&w2r(ix^d,m^c_)*w2r(ix^d,b^c_)+))*signbx(ix^d)
          w2l(ix^d,e_)=w1l(ix^d,e_)-r1l(ix^d)*((^c&w1l(ix^d,m^c_)*w1l(ix^d,b^c_)+)-&
            (^c&w2l(ix^d,m^c_)*w2l(ix^d,b^c_)+))*signbx(ix^d)
        end if
 
        ! convert velocity to momentum
        ^c&w1r(ix^d,m^c_)=w1r(ix^d,m^c_)*w1r(ix^d,rho_)\
        ^c&w1l(ix^d,m^c_)=w1l(ix^d,m^c_)*w1l(ix^d,rho_)\
        ^c&w2r(ix^d,m^c_)=w2r(ix^d,m^c_)*w2r(ix^d,rho_)\
        ^c&w2l(ix^d,m^c_)=w2l(ix^d,m^c_)*w2l(ix^d,rho_)\
        if(iw_equi_rho>0) then
          w1r(ix^d,rho_)=w1r(ix^d,rho_)-block%equi_vars(ix^d,iw_equi_rho,ip1)
          w1l(ix^d,rho_)=w1l(ix^d,rho_)-block%equi_vars(ix^d,iw_equi_rho,ip1)
          w2r(ix^d,rho_)=w2r(ix^d,rho_)-block%equi_vars(ix^d,iw_equi_rho,ip1)
          w2l(ix^d,rho_)=w2l(ix^d,rho_)-block%equi_vars(ix^d,iw_equi_rho,ip1)
        end if
     {end do\}
 
      do iw=iws,iwe
        if(flux_type(idims, iw)==flux_special) then
          ! known flux (fLC=fRC) for normal B and psi_ in GLM method
         {!dir$ ivdep
          do ix^db=ixcmin^db,ixcmax^db\}
            fc(ix^d,iw,ip1)=flc(ix^d,iw)
         {end do\}
        else if(flux_type(idims, iw)==flux_hll) then
          ! using hll flux for tracers
         {!dir$ ivdep
          do ix^db=ixcmin^db,ixcmax^db\}
            fc(ix^d,iw,ip1)=(sr(ix^d,ii)*flc(ix^d,iw)-sl(ix^d,ii)*frc(ix^d,iw) &
              +sr(ix^d,ii)*sl(ix^d,ii)*(wrc(ix^d,iw)-wlc(ix^d,iw)))/(sr(ix^d,ii)-sl(ix^d,ii))
         {end do\}
        else
          ! construct hlld flux
          ! Miyoshi equation (66) and Guo equation (46)
         {!dir$ ivdep
          !DEC$ VECTOR ALWAYS
          do ix^db=ixcmin^db,ixcmax^db\}
            if(sl(ix^d,ii)>0.d0) then
              fc(ix^d,iw,ip1)=flc(ix^d,iw)
            else if(s1l(ix^d)>=0.d0) then
              fc(ix^d,iw,ip1)=flc(ix^d,iw)+sl(ix^d,ii)*(w1l(ix^d,iw)-wlc(ix^d,iw))
            else if(sm(ix^d)>=0.d0) then
              fc(ix^d,iw,ip1)=flc(ix^d,iw)+sl(ix^d,ii)*(w1l(ix^d,iw)-wlc(ix^d,iw))+&
                s1l(ix^d)*(w2l(ix^d,iw)-w1l(ix^d,iw))
            else if(s1r(ix^d)>=0.d0) then
              fc(ix^d,iw,ip1)=frc(ix^d,iw)+sr(ix^d,ii)*(w1r(ix^d,iw)-wrc(ix^d,iw))+&
                s1r(ix^d)*(w2r(ix^d,iw)-w1r(ix^d,iw))
            else if(sr(ix^d,ii)>=0.d0) then
              fc(ix^d,iw,ip1)=frc(ix^d,iw)+sr(ix^d,ii)*(w1r(ix^d,iw)-wrc(ix^d,iw))
            else if(sr(ix^d,ii)<0.d0) then
              fc(ix^d,iw,ip1)=frc(ix^d,iw)
            end if
         {end do\}
        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(iws,iwe)
      implicit none
      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
      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
 
      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=iws,iwe
        ! 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,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,index_v_mag)>=0.d0) then
          fc(ix^d,iws:iwe,ip1)=f1r(ix^d,iws:iwe)
        else if(sr(ix^d,index_v_mag)<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_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,nw_recon
          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:nw_recon)=wlp(ixl^s,1:nw_recon)
   wrc(ixr^s,1:nw_recon)=wrp(ixr^s,1:nw_recon)
   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,nw_recon+1:nw_recon+nwaux)=wlc(ixl^s,nw_recon+1:nw_recon+nwaux)
      wrp(ixr^s,nw_recon+1:nw_recon+nwaux)=wrc(ixr^s,nw_recon+1:nw_recon+nwaux)
   endif
 
  end subroutine reconstruct_lr
 
end module mod_finite_volume