amun-code/sources/schemes.F90

!!******************************************************************************
!!
!!  This file is part of the AMUN source code, a program to perform
!!  Newtonian or relativistic magnetohydrodynamical simulations on uniform or
!!  adaptive mesh.
!!
!!  Copyright (C) 2008-2021 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License as published by
!!  the Free Software Foundation, either version 3 of the License, or
!!  (at your option) any later version.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!
!!******************************************************************************
!!
!! module: SCHEMES
!!
!!  The module provides and interface to numerical schemes, subroutines to
!!  calculate variable increment and one dimensional Riemann solvers.
!!
!!  If you implement a new set of equations, you have to add at least one
!!  corresponding update_flux subroutine, and one Riemann solver.
!!
!!******************************************************************************
!
module schemes

  implicit none

! Rieman solver and state vectors names
!
  character(len=255) , save :: name_sol = ""
  character(len=255) , save :: name_sts = "primitive"

! 4-vector reconstruction flag
!
  logical            , save :: states_4vec = .false.

! high order fluxes
!
  logical            , save :: high_order_flux = .false.

! interfaces for procedure pointers
!
  abstract interface
    subroutine update_flux_iface(dx, q, f)
      real(kind=8), dimension(:)        , intent(in)  :: dx
      real(kind=8), dimension(:,:,:,:)  , intent(in)  :: q
      real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
    end subroutine
    subroutine riemann_solver_iface(ql, qr, ul, ur, fl, fr, cl, cr, f)
      real(kind=8), dimension(:,:), intent(in)  :: ql, qr, ul, ur, fl, fr
      real(kind=8), dimension(:,:), intent(in)  :: cl, cr
      real(kind=8), dimension(:,:), intent(out) :: f
    end subroutine
  end interface

! pointer to the Riemann solver
!
  procedure(riemann_solver_iface), pointer, save :: riemann_solver => null()

! by default everything is private
!
  private

! declare public subroutines
!
  public :: initialize_schemes, finalize_schemes, print_schemes
  public :: update_flux

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
  contains
!
!===============================================================================
!!
!!***  PUBLIC SUBROUTINES  *****************************************************
!!
!===============================================================================
!
! subroutine INITIALIZE_SCHEMES:
! -----------------------------
!
!   Subroutine initiate the module by setting module parameters and subroutine
!   pointers.
!
!   Arguments:
!
!     verbose - a logical flag turning the information printing;
!     status  - an integer flag for error return value;
!
!===============================================================================
!
  subroutine initialize_schemes(verbose, status)

    use equations , only : eqsys, eos
    use parameters, only : get_parameter

    implicit none

    logical, intent(in)  :: verbose
    integer, intent(out) :: status

    character(len=64) :: solver = "HLL"
    character(len=64) :: statev = "primitive"
    character(len=64) :: flux   = "off"

!-------------------------------------------------------------------------------
!
    status = 0

! get the Riemann solver
!
    call get_parameter("riemann_solver" , solver)
    call get_parameter("state_variables", statev)

! depending on the system of equations initialize the module variables
!
    select case(trim(eqsys))

!--- HYDRODYNAMICS ---
!
    case("hd", "HD", "hydro", "HYDRO", "hydrodynamic", "HYDRODYNAMIC")

! depending on the equation of state complete the initialization
!
      select case(trim(eos))

      case("iso", "ISO", "isothermal", "ISOTHERMAL")

! select the Riemann solver
!
        select case(trim(solver))

        case("hll", "HLL")

! set the solver name
!
          name_sol =  "HLL"

! set pointers to subroutines
!
          riemann_solver => riemann_hll

        case("roe", "ROE")

! set the solver name
!
          name_sol =  "ROE"

! set pointers to subroutines
!
          riemann_solver => riemann_hd_iso_roe

        case default

          if (verbose) then
            write(*,*)
            write(*,"(1x,a)") "ERROR!"
            write(*,"(1x,a)") "The selected Riemann solver implemented "    // &
                              "for isothermal HD: '" // trim(solver) // "'."
            write(*,"(1x,a)") "Available Riemann solvers: 'hll', 'roe'."
          end if
          status = 1

        end select

      case("adi", "ADI", "adiabatic", "ADIABATIC")

! select the Riemann solver
!
        select case(trim(solver))

        case("hll", "HLL")

! set the solver name
!
          name_sol =  "HLL"

! set pointers to subroutines
!
          riemann_solver => riemann_hll

        case("hllc", "HLLC")

! set the solver name
!
          name_sol =  "HLLC"

! set pointers to subroutines
!
          riemann_solver => riemann_hd_hllc

        case("roe", "ROE")

! set the solver name
!
          name_sol =  "ROE"

! set pointers to subroutines
!
          riemann_solver => riemann_hd_adi_roe

        case default

          if (verbose) then
            write(*,*)
            write(*,"(1x,a)") "ERROR!"
            write(*,"(1x,a)") "The selected Riemann solver implemented "    // &
                              "for adiabatic HD: '" // trim(solver) // "'."
            write(*,"(1x,a)") "Available Riemann solvers: 'hll', 'hllc', 'roe'."
          end if
          status = 1

        end select

      end select

!--- MAGNETOHYDRODYNAMICS ---
!
    case("mhd", "MHD", "magnetohydrodynamic", "MAGNETOHYDRODYNAMIC")

! depending on the equation of state complete the initialization
!
      select case(trim(eos))

      case("iso", "ISO", "isothermal", "ISOTHERMAL")

! select the Riemann solver
!
        select case(trim(solver))

        case("hll", "HLL")

! set the solver name
!
          name_sol =  "HLL"

! set pointers to subroutines
!
          riemann_solver => riemann_hll

        case ("hlld", "HLLD")

! set the solver name
!
          name_sol =  "HLLD"

! set the solver pointer
!
          riemann_solver => riemann_mhd_iso_hlld

        case("roe", "ROE")

! set the solver name
!
          name_sol =  "ROE"

! set pointers to subroutines
!
          riemann_solver => riemann_mhd_iso_roe

        case default

          if (verbose) then
            write(*,*)
            write(*,"(1x,a)") "ERROR!"
            write(*,"(1x,a)") "The selected Riemann solver implemented "    // &
                              "for isothermal MHD: '" // trim(solver) // "'."
            write(*,"(1x,a)") "Available Riemann solvers: 'hll', 'hlld'"    // &
                              ", 'roe'."
          end if
          status = 1

        end select

      case("adi", "ADI", "adiabatic", "ADIABATIC")

! select the Riemann solver
!
        select case(trim(solver))

        case ("hll", "HLL")

! set the solver name
!
          name_sol =  "HLL"

! set pointers to subroutines
!
          riemann_solver => riemann_hll

        case ("hllc", "HLLC")

! set the solver name
!
          name_sol =  "HLLC"

! set the solver pointer
!
          riemann_solver => riemann_mhd_hllc

        case ("hlld", "HLLD")

! set the solver name
!
          name_sol =  "HLLD"

! set the solver pointer
!
          riemann_solver => riemann_mhd_adi_hlld

        case("roe", "ROE")

! set the solver name
!
          name_sol =  "ROE"

! set pointers to subroutines
!
          riemann_solver => riemann_mhd_adi_roe

        case default

          if (verbose) then
            write(*,*)
            write(*,"(1x,a)") "ERROR!"
            write(*,"(1x,a)") "The selected Riemann solver implemented "    // &
                              "for adiabatic MHD: '" // trim(solver) // "'."
            write(*,"(1x,a)") "Available Riemann solvers: 'hll', 'hllc'"    // &
                              ", 'hlld', 'roe'."
          end if
          status = 1

        end select

      end select

!--- SPECIAL RELATIVITY HYDRODYNAMICS ---
!
    case("srhd", "SRHD")

! depending on the equation of state complete the initialization
!
      select case(trim(eos))

      case("adi", "ADI", "adiabatic", "ADIABATIC")

! select the state reconstruction method
!
        select case(trim(statev))

        case("4vec", "4-vector", "4VEC", "4-VECTOR")

! set the state reconstruction name
!
          name_sts =  "4-vector"

! set 4-vector reconstruction flag
!
          states_4vec = .true.

! in the case of state variables, revert to primitive
!
        case default

! set the state reconstruction name
!
          name_sts =  "primitive"

        end select

! select the Riemann solver
!
        select case(trim(solver))

        case ("hll", "HLL")

! set the solver name
!
          name_sol =  "HLL"

! set pointers to subroutines
!
          riemann_solver => riemann_hll

        case("hllc", "HLLC", "hllcm", "HLLCM", "hllc-m", "HLLC-M")

! set the solver name
!
          name_sol =  "HLLC (Mignone & Bodo 2005)"

! set pointers to subroutines
!
          riemann_solver => riemann_srhd_hllc

        case default

          if (verbose) then
            write(*,*)
            write(*,"(1x,a)") "ERROR!"
            write(*,"(1x,a)") "The selected Riemann solver implemented "    // &
                              "for adiabatic SR-HD: '" // trim(solver) // "'."
            write(*,"(1x,a)") "Available Riemann solvers: 'hll', 'hllc'."
          end if
          status = 1

        end select

      end select

!--- SPECIAL RELATIVITY MAGNETOHYDRODYNAMICS ---
!
    case("srmhd", "SRMHD")

! depending on the equation of state complete the initialization
!
      select case(trim(eos))

      case("adi", "ADI", "adiabatic", "ADIABATIC")

! select the state reconstruction method
!
        select case(trim(statev))

        case("4vec", "4-vector", "4VEC", "4-VECTOR")

! set the state reconstruction name
!
          name_sts =  "4-vector"

! set 4-vector reconstruction flag
!
          states_4vec = .true.

! in the case of state variables, revert to primitive
!
        case default

! set the state reconstruction name
!
          name_sts =  "primitive"

        end select

! select the Riemann solver
!
        select case(trim(solver))

        case ("hll", "HLL")

! set the solver name
!
          name_sol =  "HLL"

! set pointers to subroutines
!
          riemann_solver => riemann_hll

        case("hllc", "HLLC", "hllcm", "HLLCM", "hllc-m", "HLLC-M")

! set the solver name
!
          name_sol =  "HLLC (Mignone & Bodo 2006)"

! set pointers to subroutines
!
          riemann_solver => riemann_srmhd_hllc

        case default

          if (verbose) then
            write(*,*)
            write(*,"(1x,a)") "ERROR!"
            write(*,"(1x,a)") "The selected Riemann solver implemented "    // &
                              "for adiabatic SR-MHD: '" // trim(solver) // "'."
            write(*,"(1x,a)") "Available Riemann solvers: 'hll', 'hllc'."
          end if
          status = 1

        end select

      end select

    end select

! flag for higher order flux correction
!
    call get_parameter("high_order_flux", flux)

    select case(trim(flux))
    case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
      high_order_flux = .true.
    case default
      high_order_flux = .false.
    end select

!-------------------------------------------------------------------------------
!
  end subroutine initialize_schemes
!
!===============================================================================
!
! subroutine FINALIZE_SCHEMES:
! ---------------------------
!
!   Subroutine releases memory used by the module.
!
!   Arguments:
!
!     status - an integer flag for error return value;
!
!===============================================================================
!
  subroutine finalize_schemes(status)

    implicit none

    integer, intent(out) :: status

!-------------------------------------------------------------------------------
!
    status = 0

    nullify(riemann_solver)

!-------------------------------------------------------------------------------
!
  end subroutine finalize_schemes
!
!===============================================================================
!
! subroutine PRINT_SCHEMES:
! ------------------------
!
!   Subroutine prints module parameters and settings.
!
!   Arguments:
!
!     verbose - a logical flag turning the information printing;
!
!===============================================================================
!
  subroutine print_schemes(verbose)

    use helpers, only : print_section, print_parameter

    implicit none

    logical, intent(in) :: verbose

!-------------------------------------------------------------------------------
!
    if (verbose) then
      call print_section(verbose, "Schemes")
      call print_parameter(verbose, "Riemann solver" , name_sol)
      call print_parameter(verbose, "state variables", name_sts)
      if (high_order_flux) then
        call print_parameter(verbose, "high order flux correction", "on" )
      else
        call print_parameter(verbose, "high order flux correction", "off")
      end if
    end if

!-------------------------------------------------------------------------------
!
  end subroutine print_schemes
!
!===============================================================================
!
! subroutine UPDATE_FLUX:
! ----------------------
!
!   Subroutine solves the Riemann problem along each direction and calculates
!   the numerical fluxes, which are used later to calculate the conserved
!   variable increments.
!
!   Arguments:
!
!     dx - the spatial step;
!     q  - the array of primitive variables;
!     s  - the state arrays for primitive variables;
!     f  - the array of numerical fluxes;
!
!===============================================================================
!
  subroutine update_flux(dx, q, s, f)

    use coordinates, only : nn => bcells, nbl, neu
    use equations  , only : relativistic
    use equations  , only : nf, ivars, ivx, ivz

    implicit none

    real(kind=8), dimension(:)          , intent(in)    :: dx
    real(kind=8), dimension(:,:,:,:)    , intent(in)    :: q
    real(kind=8), dimension(:,:,:,:,:,:), intent(inout) :: s
    real(kind=8), dimension(:,:,:,:,:)  , intent(out)   :: f

    integer      :: n, i, j, k = 1, l, p
    real(kind=8) :: vm

    real(kind=8), dimension(nf,nn,2) :: qi
    real(kind=8), dimension(nf,nn)   :: fi

!-------------------------------------------------------------------------------
!
! in the relativistic case, apply the reconstruction on variables
! using the four-velocities if requested
!
    if (relativistic .and. states_4vec) then

      f(:,:,:,:,1) = q(:,:,:,:)

#if NDIMS == 3
      do k = 1, nn
#endif /* NDIMS == 3 */
        do j = 1, nn
          do i = 1, nn

            vm = sqrt(1.0d+00 - sum(q(ivx:ivz,i,j,k)**2))

            f(ivx:ivz,i,j,k,1) = f(ivx:ivz,i,j,k,1) / vm

          end do ! i = 1, nn
        end do ! j = 1, nn
#if NDIMS == 3
      end do ! k = 1, nn
#endif /* NDIMS == 3 */

    call reconstruct_interfaces(dx(:), f(:,:,:,:,1), s(:,:,:,:,:,:))

! convert the states' four-velocities back to velocities
!
#if NDIMS == 3
      do k = 1, nn
#endif /* NDIMS == 3 */
        do j = 1, nn
          do i = 1, nn

            do l = 1, 2
              do p = 1, NDIMS

                vm = sqrt(1.0d+00 + sum(s(ivx:ivz,i,j,k,l,p)**2))

                s(ivx:ivz,i,j,k,l,p) = s(ivx:ivz,i,j,k,l,p) / vm

              end do ! p = 1, ndims
            end do ! l = 1, 2
          end do ! i = 1, nn
        end do ! j = 1, nn
#if NDIMS == 3
      end do ! k = 1, nn
#endif /* NDIMS == 3 */

    else

      call reconstruct_interfaces(dx(:), q(:,:,:,:), s(:,:,:,:,:,:))

    end if

    f(:,:,:,:,:) = 0.0d+00

!  calculate the flux along the X-direction
!
#if NDIMS == 3
    do k = nbl, neu
#endif /* NDIMS == 3 */
      do j = nbl, neu

! copy states to directional lines with proper vector component ordering
!
        do n = 1, nf
          qi(n,:,1:2) = s(ivars(1,n),:,j,k,1:2,1)
        end do

! call one dimensional Riemann solver in order to obtain numerical fluxes
!
        call numerical_flux(qi(:,:,:), fi(:,:))

! update the array of fluxes
!
        do n = 1, nf
          f(ivars(1,n),:,j,k,1) = fi(n,:)
        end do

      end do ! j = nbl, neu

!  calculate the flux along the Y direction
!
      do i = nbl, neu

! copy directional variable vectors to pass to the one dimensional solver
!
        do n = 1, nf
          qi(n,:,1:2) = s(ivars(2,n),i,:,k,1:2,2)
        end do

! call one dimensional Riemann solver in order to obtain numerical fluxes
!
        call numerical_flux(qi(:,:,:), fi(:,:))

! update the array of fluxes
!
        do n = 1, nf
          f(ivars(2,n),i,:,k,2) = fi(n,:)
        end do

      end do ! i = nbl, neu
#if NDIMS == 3
    end do ! k = nbl, neu
#endif /* NDIMS == 3 */

#if NDIMS == 3
!  calculate the flux along the Z direction
!
    do j = nbl, neu
      do i = nbl, neu

! copy directional variable vectors to pass to the one dimensional solver
!
        do n = 1, nf
          qi(n,:,1:2) = s(ivars(3,n),i,j,:,1:2,3)
        end do

! call one dimensional Riemann solver in order to obtain numerical fluxes
!
        call numerical_flux(qi(:,:,:), fi(:,:))

! update the array of fluxes
!
        do n = 1, nf
          f(ivars(3,n),i,j,:,3) = fi(n,:)
        end do

      end do ! i = nbl, neu
    end do ! j = nbl, neu
#endif /* NDIMS == 3 */

!-------------------------------------------------------------------------------
!
  end subroutine update_flux
!
!===============================================================================
!!
!!***  PRIVATE SUBROUTINES  ****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine RECONSTRUCT_INTERFACES:
! ---------------------------------
!
!   Subroutine reconstructs the Riemann states along all directions.
!
!   Arguments:
!
!     dx   - the spatial step;
!     q    - the array of primitive variables;
!     qi   - the array of reconstructed states (2 in each direction);
!
!===============================================================================
!
  subroutine reconstruct_interfaces(dx, q, qi)

! include external procedures
!
    use equations     , only : nf
    use equations     , only : positive
    use interpolations, only : interfaces

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:)          , intent(in)  :: dx
    real(kind=8), dimension(:,:,:,:)    , intent(in)  :: q
    real(kind=8), dimension(:,:,:,:,:,:), intent(out) :: qi

! local variables
!
    integer :: p
!
!-------------------------------------------------------------------------------
!
! interpolate interfaces for all flux variables
!
    do p = 1, nf
      call interfaces(positive(p), dx(1:NDIMS), q (p,:,:,:),                   &
                                                qi(p,:,:,:,1:2,1:NDIMS))
    end do

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_interfaces
!
!===============================================================================
!
! subroutine NUMERICAL_FLUX:
! -------------------------
!
!   Subroutine prepares Riemann states and calls the selected Riemann solver
!   in order to get the numerical flux. If requested, the resulting flux
!   is corrected by higher order correction terms.
!
!   Arguments:
!
!     q - the array of primitive variables at the Riemann states;
!     f - the output array of fluxes;
!
!===============================================================================
!
  subroutine numerical_flux(q, f)

    use equations, only : magnetized, cmax, ibx, ibp
    use equations, only : prim2cons, fluxspeed

    implicit none

    real(kind=8), dimension(:,:,:), intent(inout) :: q
    real(kind=8), dimension(:,:)  , intent(out)   :: f

    integer      :: i
    real(kind=8) :: bx, bp

    real(kind=8), dimension(size(q,1),size(q,2)) :: ql, qr, ul, ur, fl, fr
    real(kind=8), dimension(        2,size(q,2)) :: cl, cr

!-------------------------------------------------------------------------------
!
! copy the state vectors
!
    ql(:,:) = q(:,:,1)
    qr(:,:) = q(:,:,2)

! obtain the state values for Bx and Psi for the GLM-MHD equations
!
    if (magnetized) then
      do i = 1, size(q, 2)

        bx        = 0.5d+00 * ((qr(ibx,i) + ql(ibx,i))                         &
                                             - (qr(ibp,i) - ql(ibp,i)) / cmax)
        bp        = 0.5d+00 * ((qr(ibp,i) + ql(ibp,i))                         &
                                             - (qr(ibx,i) - ql(ibx,i)) * cmax)

        ql(ibx,i) = bx
        qr(ibx,i) = bx
        ql(ibp,i) = bp
        qr(ibp,i) = bp

      end do
    end if

! calculate corresponding conserved variables of the left and right states
!
    call prim2cons(ql(:,:), ul(:,:))
    call prim2cons(qr(:,:), ur(:,:))

! calculate the physical fluxes and speeds at the states
!
    call fluxspeed(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
    call fluxspeed(qr(:,:), ur(:,:), fr(:,:), cr(:,:))

! get the HLL fluxes
!
    call riemann_solver(ql(:,:), qr(:,:), ul(:,:), ur(:,:),                    &
                        fl(:,:), fr(:,:), cl(:,:), cr(:,:), f(:,:))

! higher order flux corrections
!
    call higher_order_flux_correction(f)

!-------------------------------------------------------------------------------
!
  end subroutine numerical_flux
!
!===============================================================================
!
! subroutine RIEMANN_HLL:
! ----------------------
!
!   Subroutine solves one dimensional general Riemann problem using
!   the Harten-Lax-van Leer (HLL) method.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the Riemann states;
!     ul, ur - the conservative variables of the Riemann states;
!     fl, fr - the physical fluxes of the Riemann states;
!     f      - the Riemann average flux;
!
!   References:
!
!     [1] Harten, A., Lax, P. D. & Van Leer, B.
!         "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic
!          Conservation Laws",
!         SIAM Review, 1983, Volume 25, Number 1, pp. 35-61
!
!===============================================================================
!
  subroutine riemann_hll(ql, qr, ul, ur, fl, fr, cl, cr, f)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: ql, qr, ul, ur, fl, fr, cl, cr
    real(kind=8), dimension(:,:), intent(out) :: f

! local variables
!
    integer      :: n, i
    real(kind=8) :: sl, sr, srml

! local arrays to store the states
!
    real(kind=8), dimension(size(ql,1)) :: wl, wr
!
!-------------------------------------------------------------------------------
!
! get the length of vector
!
    n = size(ql, 2)

! iterate over all points
!
    do i = 1, n

! estimate the minimum and maximum speeds
!
      sl = min(cl(1,i), cr(1,i))
      sr = max(cl(2,i), cr(2,i))

! calculate the HLL flux
!
      if (sl >= 0.0d+00) then

        f(:,i) = fl(:,i)

      else if (sr <= 0.0d+00) then

        f(:,i) = fr(:,i)

      else ! sl < 0 < sr

! calculate speed difference
!
        srml = sr - sl

! calculate vectors of the left and right-going waves
!
        wl(:)  = sl * ul(:,i) - fl(:,i)
        wr(:)  = sr * ur(:,i) - fr(:,i)

! calculate the fluxes for the intermediate state
!
        f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml

      end if ! sl < 0 < sr

    end do ! i = 1, n

!-------------------------------------------------------------------------------
!
  end subroutine riemann_hll
!
!===============================================================================
!
! subroutine RIEMANN_HD_HLLC:
! --------------------------
!
!   Subroutine solves one dimensional Riemann problem using the HLLC
!   method by Gurski or Li.  The HLLC and HLLC-C differ by definitions of
!   the tangential components of the velocity and magnetic field.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the Riemann states;
!     ul, ur - the conservative variables of the Riemann states;
!     fl, fr - the physical fluxes of the Riemann states;
!     f      - the Riemann average flux;
!
!   References:
!
!     [1] Gurski, K. F.,
!         "An HLLC-Type Approximate Riemann Solver for Ideal
!          Magnetohydrodynamics",
!         SIAM Journal on Scientific Computing, 2004, Volume 25, Issue 6,
!         pp. 2165–2187
!     [2] Li, S.,
!         "An HLLC Riemann solver for magneto-hydrodynamics",
!         Journal of Computational Physics, 2005, Volume 203, Issue 1,
!         pp. 344-357
!
!===============================================================================
!
  subroutine riemann_hd_hllc(ql, qr, ul, ur, fl, fr, cl, cr, f)

! include external procedures
!
    use equations, only : idn, ivy, ivz, imx, imy, imz, ien

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: ql, qr, ul, ur, fl, fr, cl, cr
    real(kind=8), dimension(:,:), intent(out) :: f

! local variables
!
    integer      :: i
    real(kind=8) :: dn, pr
    real(kind=8) :: sl, sr, sm, slmm, srmm

! local arrays to store the states
!
    real(kind=8), dimension(size(ql,1)) :: wl, wr, ui
!
!-------------------------------------------------------------------------------
!
! iterate over all points
!
    do i = 1, size(ql, 2)

! estimate the minimum and maximum speeds
!
      sl = min(cl(1,i), cr(1,i))
      sr = max(cl(2,i), cr(2,i))

! calculate the HLL flux
!
      if (sl >= 0.0d+00) then

        f(:,i) = fl(:,i)

      else if (sr <= 0.0d+00) then

        f(:,i) = fr(:,i)

      else ! sl < 0 < sr

! calculate vectors of the left and right-going waves
!
        wl(:) = sl * ul(:,i) - fl(:,i)
        wr(:) = sr * ur(:,i) - fr(:,i)

! the speed of contact discontinuity
!
        dn =  wr(idn) - wl(idn)
        sm = (wr(imx) - wl(imx)) / dn

! calculate the pressure of the intermediate state
!
        pr = wl(idn) * wr(imx) / dn - wr(idn) * wl(imx) / dn

! separate intermediate states depending on the sign of the advection speed
!
        if (sm > 0.0d+00) then ! sm > 0

! the left speed difference
!
          slmm    =  sl - sm

! conservative variables for the left intermediate state
!
          ui(idn) =  wl(idn) / slmm
          ui(imx) =  ui(idn) * sm
          ui(imy) =  ui(idn) * ql(ivy,i)
          ui(imz) =  ui(idn) * ql(ivz,i)
          ui(ien) = (wl(ien) + sm * pr) / slmm

! the left intermediate flux
!
          f(:,i)  = sl * ui(:) - wl(:)

        else if (sm < 0.0d+00) then ! sm < 0

! the right speed difference
!
          srmm    =  sr - sm

! conservative variables for the right intermediate state
!
          ui(idn) =  wr(idn) / srmm
          ui(imx) =  ui(idn) * sm
          ui(imy) =  ui(idn) * qr(ivy,i)
          ui(imz) =  ui(idn) * qr(ivz,i)
          ui(ien) = (wr(ien) + sm * pr) / srmm

! the right intermediate flux
!
          f(:,i)  = sr * ui(:) - wr(:)

        else ! sm = 0

! intermediate flux is constant across the contact discontinuity and all except
! the parallel momentum flux are zero
!
          f(idn,i) =   0.0d+00
          f(imx,i) = - wl(imx)
          f(imy,i) =   0.0d+00
          f(imz,i) =   0.0d+00
          f(ien,i) =   0.0d+00

        end if ! sm = 0

      end if ! sl < 0 < sr

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_hd_hllc
!
!===============================================================================
!
! subroutine RIEMANN_MHD_HLLC:
! ---------------------------
!
!   Subroutine solves one dimensional Riemann problem using the HLLC method,
!   by Toro.  In the HLLC method the tangential components of the velocity are
!   discontinuous actoss the contact dicontinuity.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the Riemann states;
!     ul, ur - the conservative variables of the Riemann states;
!     fl, fr - the physical fluxes of the Riemann states;
!     f      - the Riemann average flux;
!
!   References:
!
!     [1] Toro, E. F., Spruce, M., & Speares, W.
!         "Restoration of the contact surface in the HLL-Riemann solver",
!         Shock Waves, 1994, Volume 4, Issue 1, pp. 25-34
!
!===============================================================================
!
  subroutine riemann_mhd_hllc(ql, qr, ul, ur, fl, fr, cl, cr, f)

! include external procedures
!
    use equations, only : idn, ivy, ivz, imx, imy, imz, ien
    use equations, only : ibx, iby, ibz, ibp

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: ql, qr, ul, ur, fl, fr, cl, cr
    real(kind=8), dimension(:,:), intent(out) :: f

! local variables
!
    integer      :: i
    real(kind=8) :: dn, pt, vy, vz, bx, by, bz, vb, b2
    real(kind=8) :: sl, sr, sm, srml, slmm, srmm

! local arrays to store the states
!
    real(kind=8), dimension(size(ql,1)) :: wl, wr, ui
!
!-------------------------------------------------------------------------------
!
! iterate over all points
!
    do i = 1, size(ql, 2)

! estimate the minimum and maximum speeds
!
      sl = min(cl(1,i), cr(1,i))
      sr = max(cl(2,i), cr(2,i))

! calculate the HLLC flux
!
      if (sl >= 0.0d+00) then

        f(:,i) = fl(:,i)

      else if (sr <= 0.0d+00) then

        f(:,i) = fr(:,i)

      else ! sl < 0 < sr

! Bx and square of Bx, i.e. Bₓ²
!
        bx = ql(ibx,i)
        b2 = ql(ibx,i) * qr(ibx,i)

! speed difference
!
        srml = sr - sl

! calculate vectors of the left and right-going waves
!
        wl(:) = sl * ul(:,i) - fl(:,i)
        wr(:) = sr * ur(:,i) - fr(:,i)

! the speed of contact discontinuity
!
        dn =  wr(idn) - wl(idn)
        sm = (wr(imx) - wl(imx)) / dn

! the total pressure, constant across the contact discontinuity
!
        pt = wl(idn) * wr(imx) / dn - wr(idn) * wl(imx) / dn + b2

! separate the cases when Bₓ = 0 and Bₓ ≠ 0
!
        if (b2 > 0.0d+00) then

! constant intermediate state tangential components of velocity and
! magnetic field
!
          vy = (wr(imy) - wl(imy)) / dn
          vz = (wr(imz) - wl(imz)) / dn
          by = (wr(iby) - wl(iby)) / srml
          bz = (wr(ibz) - wl(ibz)) / srml

! scalar product of velocity and magnetic field for the intermediate states
!
          vb = sm * bx + vy * by + vz * bz

! separate intermediate states depending on the sign of the advection speed
!
          if (sm > 0.0d+00) then ! sm > 0

! the left speed difference
!
            slmm = sl - sm

! conservative variables for the left intermediate state
!
            ui(idn) =  wl(idn) / slmm
            ui(imx) =  ui(idn) * sm
            ui(imy) =  ui(idn) * vy
            ui(imz) =  ui(idn) * vz
            ui(ibx) =  bx
            ui(iby) =  by
            ui(ibz) =  bz
            ui(ibp) =  ql(ibp,i)
            ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm

! the left intermediate flux
!
            f(:,i)  = sl * ui(:) - wl(:)

          else if (sm < 0.0d+00) then ! sm < 0

! the right speed difference
!
            srmm = sr - sm

! conservative variables for the right intermediate state
!
            ui(idn) =  wr(idn) / srmm
            ui(imx) =  ui(idn) * sm
            ui(imy) =  ui(idn) * vy
            ui(imz) =  ui(idn) * vz
            ui(ibx) =  bx
            ui(iby) =  by
            ui(ibz) =  bz
            ui(ibp) =  qr(ibp,i)
            ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm

! the right intermediate flux
!
            f(:,i)  = sr * ui(:) - wr(:)

          else ! sm = 0

! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
            f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml

          end if ! sm = 0

        else ! Bₓ = 0

! separate intermediate states depending on the sign of the advection speed
!
          if (sm > 0.0d+00) then ! sm > 0

! the left speed difference
!
            slmm = sl - sm

! conservative variables for the left intermediate state
!
            ui(idn) =  wl(idn) / slmm
            ui(imx) =  ui(idn) * sm
            ui(imy) =  ui(idn) * ql(ivy,i)
            ui(imz) =  ui(idn) * ql(ivz,i)
            ui(ibx) =  0.0d+00
            ui(iby) =  wl(iby) / slmm
            ui(ibz) =  wl(ibz) / slmm
            ui(ibp) =  ql(ibp,i)
            ui(ien) = (wl(ien) + sm * pt) / slmm

! the left intermediate flux
!
            f(:,i)  = sl * ui(:) - wl(:)

          else if (sm < 0.0d+00) then ! sm < 0

! the right speed difference
!
            srmm = sr - sm

! conservative variables for the right intermediate state
!
            ui(idn) =  wr(idn) / srmm
            ui(imx) =  ui(idn) * sm
            ui(imy) =  ui(idn) * qr(ivy,i)
            ui(imz) =  ui(idn) * qr(ivz,i)
            ui(ibx) =  0.0d+00
            ui(iby) =  wr(iby) / srmm
            ui(ibz) =  wr(ibz) / srmm
            ui(ibp) =  qr(ibp,i)
            ui(ien) = (wr(ien) + sm * pt) / srmm

! the right intermediate flux
!
            f(:,i)  = sr * ui(:) - wr(:)

          else ! sm = 0

! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
            f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml

          end if ! sm = 0

        end if ! Bₓ = 0

      end if ! sl < 0 < sr

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_mhd_hllc
!
!===============================================================================
!
! subroutine RIEMANN_SRHD_HLLC:
! ----------------------------
!
!   Subroutine solves one dimensional Riemann problem using
!   the Harten-Lax-van Leer method with contact discontinuity resolution (HLLC)
!   by Mignone & Bodo.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the Riemann states;
!     ul, ur - the conservative variables of the Riemann states;
!     fl, fr - the physical fluxes of the Riemann states;
!     f      - the Riemann average flux;
!
!   References:
!
!     [1] Mignone, A. & Bodo, G.
!         "An HLLC Riemann solver for relativistic flows - I. Hydrodynamics",
!         Monthly Notices of the Royal Astronomical Society,
!         2005, Volume 364, Pages 126-136
!
!===============================================================================
!
  subroutine riemann_srhd_hllc(ql, qr, ul, ur, fl, fr, cl, cr, f)

! include external procedures
!
    use algebra  , only : quadratic
    use equations, only : idn, imx, imy, imz, ien

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: ql, qr, ul, ur, fl, fr, cl, cr
    real(kind=8), dimension(:,:), intent(out) :: f

! local variables
!
    integer      :: i, nr
    real(kind=8) :: sl, sr, srml, sm
    real(kind=8) :: pr, dv

! local arrays to store the states
!
    real(kind=8), dimension(size(ql,1)) :: wl, wr, uh, us, fh
    real(kind=8), dimension(3)          :: a
    real(kind=8), dimension(2)          :: x
!
!-------------------------------------------------------------------------------
!
! iterate over all points
!
    do i = 1, size(ql, 2)

! estimate the minimum and maximum speeds
!
      sl = min(cl(1,i), cr(1,i))
      sr = max(cl(2,i), cr(2,i))

! calculate the HLL flux
!
      if (sl >= 0.0d+00) then

        f(:,i) = fl(:,i)

      else if (sr <= 0.0d+00) then

        f(:,i) = fr(:,i)

      else ! sl < 0 < sr

! calculate the inverse of speed difference
!
        srml = sr - sl

! calculate vectors of the left and right-going waves
!
        wl(:)  = sl * ul(:,i) - fl(:,i)
        wr(:)  = sr * ur(:,i) - fr(:,i)

! calculate fluxes for the intermediate state
!
        uh(:)  = (     wr(:) -      wl(:)) / srml
        fh(:)  = (sl * wr(:) - sr * wl(:)) / srml

! correct the energy waves
!
        wl(ien)   = wl(ien) + wl(idn)
        wr(ien)   = wr(ien) + wr(idn)

! prepare the quadratic coefficients (eq. 18 in [1])
!
        a(1) = uh(imx)
        a(2) = - (fh(imx) + uh(ien) + uh(idn))
        a(3) = fh(ien) + fh(idn)

! solve the quadratic equation
!
        nr   = quadratic(a(1:3), x(1:2))

! if Δ < 0, just use the HLL flux
!
        if (nr < 1) then
          f(:,i) = fh(:)
        else

! get the contact dicontinuity speed
!
          sm = x(1)

! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
!
          if ((sm <= sl) .or. (sm >= sr)) then
            f(:,i) = fh(:)
          else

! calculate total pressure (eq. 17 in [1])
!
            pr = fh(imx) - (fh(ien) + fh(idn)) * sm

! if the pressure is negative, use the HLL flux
!
            if (pr <= 0.0d+00) then
              f(:,i) = fh(:)
            else

! depending in the sign of the contact dicontinuity speed, calculate the proper
! state and corresponding flux
!
              if (sm > 0.0d+00) then

! calculate the conserved variable vector (eqs. 16 in [1])
!
                dv      = sl - sm
                us(idn) = wl(idn) / dv
                us(imy) = wl(imy) / dv
                us(imz) = wl(imz) / dv
                us(ien) = (wl(ien) + pr * sm) / dv
                us(imx) = (us(ien) + pr) * sm
                us(ien) = us(ien) - us(idn)

! calculate the flux (eq. 14 in [1])
!
                f(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))

              else if (sm < 0.0d+00) then

! calculate the conserved variable vector (eqs. 16 in [1])
!
                dv      = sr - sm
                us(idn) = wr(idn) / dv
                us(imy) = wr(imy) / dv
                us(imz) = wr(imz) / dv
                us(ien) = (wr(ien) + pr * sm) / dv
                us(imx) = (us(ien) + pr) * sm
                us(ien) = us(ien) - us(idn)

! calculate the flux (eq. 14 in [1])
!
                f(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))

              else

! intermediate flux is constant across the contact discontinuity and all fluxes
! except the parallel momentum one are zero
!
                f(idn,i) = 0.0d+00
                f(imx,i) = pr
                f(imy,i) = 0.0d+00
                f(imz,i) = 0.0d+00
                f(ien,i) = 0.0d+00

              end if ! sm == 0

            end if ! p* < 0

          end if ! sl < sm < sr

        end if ! nr < 1

      end if ! sl < 0 < sr

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_srhd_hllc
!
!===============================================================================
!
! subroutine RIEMANN_SRMHD_HLLC:
! -----------------------------
!
!   Subroutine solves one dimensional Riemann problem using
!   the Harten-Lax-van Leer method with contact discontinuity resolution (HLLC)
!   by Mignone & Bodo.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the Riemann states;
!     ul, ur - the conservative variables of the Riemann states;
!     fl, fr - the physical fluxes of the Riemann states;
!     f      - the Riemann average flux;
!
!   References:
!
!     [1] Mignone, A. & Bodo, G.
!         "An HLLC Riemann solver for relativistic flows - II.
!          Magnetohydrodynamics",
!         Monthly Notices of the Royal Astronomical Society,
!         2006, Volume 368, Pages 1040-1054
!
!===============================================================================
!
  subroutine riemann_srmhd_hllc(ql, qr, ul, ur, fl, fr, cl, cr, f)

! include external procedures
!
    use algebra  , only : quadratic
    use equations, only : idn, ivx, imx, imy, imz, ien, ibx, iby, ibz, ibp

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: ql, qr, ul, ur, fl, fr, cl, cr
    real(kind=8), dimension(:,:), intent(out) :: f

! local variables
!
    integer      :: i, nr
    real(kind=8) :: sl, sr, srml, sm
    real(kind=8) :: bx2
    real(kind=8) :: pt, dv, fc, uu, ff, uf
    real(kind=8) :: vv, vb, gi

! local arrays to store the states
!
    real(kind=8), dimension(size(ql,1)) :: wl, wr, uh, us, fh
    real(kind=8), dimension(3)          :: a, vs
    real(kind=8), dimension(2)          :: x
!
!-------------------------------------------------------------------------------
!
! iterate over all points
!
    do i = 1, size(ql, 2)

! estimate the minimum and maximum speeds
!
      sl = min(cl(1,i), cr(1,i))
      sr = max(cl(2,i), cr(2,i))

! calculate the HLL flux
!
      if (sl >= 0.0d+00) then

        f(:,i) = fl(:,i)

      else if (sr <= 0.0d+00) then

        f(:,i) = fr(:,i)

      else ! sl < 0 < sr

! calculate the inverse of speed difference
!
        srml = sr - sl

! calculate vectors of the left and right-going waves
!
        wl(:)  = sl * ul(:,i) - fl(:,i)
        wr(:)  = sr * ur(:,i) - fr(:,i)

! calculate fluxes for the intermediate state
!
        uh(:)  = (     wr(:) -      wl(:)) / srml
        fh(:)  = (sl * wr(:) - sr * wl(:)) / srml

! correct the energy waves
!
        wl(ien)   = wl(ien) + wl(idn)
        wr(ien)   = wr(ien) + wr(idn)

! calculate Bₓ²
!
        bx2       = ql(ibx,i) * qr(ibx,i)

! calculate the contact dicontinuity speed solving eq. (41)
!
        if (bx2 >= 1.0d-16) then

! prepare the quadratic coefficients, (eq. 42 in [1])
!
          uu = sum(uh(iby:ibz) * uh(iby:ibz))
          ff = sum(fh(iby:ibz) * fh(iby:ibz))
          uf = sum(uh(iby:ibz) * fh(iby:ibz))

          a(1) = uh(imx) - uf
          a(2) = uu + ff - (fh(imx) + uh(ien) + uh(idn))
          a(3) = fh(ien) + fh(idn) - uf

! solve the quadratic equation
!
          nr   = quadratic(a(1:3), x(1:2))

! if Δ < 0, just use the HLL flux
!
          if (nr < 1) then
            f(:,i) = fh(:)
          else

! get the contact dicontinuity speed
!
            sm = x(1)

! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
!
            if ((sm <= sl) .or. (sm >= sr)) then
              f(:,i) = fh(:)
            else

! substitute magnetic field components from the HLL state (eq. 37 in [1])
!
              us(ibx) = ql(ibx,i)
              us(iby) = uh(iby)
              us(ibz) = uh(ibz)
              us(ibp) = ql(ibp,i)

! calculate velocity components without Bₓ normalization (eq. 38 in [1])
!
              vs(1)   = sm
              vs(2)   = (us(iby) * sm - fh(iby)) / us(ibx)
              vs(3)   = (us(ibz) * sm - fh(ibz)) / us(ibx)

! calculate v² and v.B
!
              vv      = sum(vs(1:3) * vs(  1:  3))
              vb      = sum(vs(1:3) * us(ibx:ibz))

! calculate inverse of Lorentz factor
!
              gi      = 1.0d+00 - vv

! calculate total pressure (eq. 40 in [1])
!
              pt = fh(imx) + gi * bx2 - (fh(ien) + fh(idn) - us(ibx) * vb) * sm

! if the pressure is negative, use the HLL flux
!
              if (pt <= 0.0d+00) then
                f(:,i) = fh(:)
              else

! depending in the sign of the contact dicontinuity speed, calculate the proper
! state and corresponding flux
!
                if (sm > 0.0d+00) then

! calculate the conserved variable vector (eqs. 43-46 in [1])
!
                  dv      = sl - sm
                  fc      = (sl - ql(ivx,i)) / dv
                  us(idn) = fc * ul(idn,i)
                  us(imy) = (sl * ul(imy,i) - fl(imy,i)                        &
                          - us(ibx) * (gi * us(iby) + vs(2) * vb)) / dv
                  us(imz) = (sl * ul(imz,i) - fl(imz,i)                        &
                          - us(ibx) * (gi * us(ibz) + vs(3) * vb)) / dv
                  us(ien) = (sl * (ul(ien,i) + ul(idn,i))                      &
                          - (fl(ien,i) + fl(idn,i))                            &
                          + pt * sm - us(ibx) * vb) / dv - us(idn)

! calculate normal component of momentum
!
                  us(imx) = (us(ien) + us(idn) + pt) * sm - us(ibx) * vb

! calculate the flux (eq. 14 in [1])
!
                  f(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))

                else if (sm < 0.0d+00) then

! calculate the conserved variable vector (eqs. 43-46 in [1])
!
                  dv      = sr - sm
                  fc      = (sr - qr(ivx,i)) / dv
                  us(idn) = fc * ur(idn,i)
                  us(imy) = (sr * ur(imy,i) - fr(imy,i)                        &
                          - us(ibx) * (gi * us(iby) + vs(2) * vb)) / dv
                  us(imz) = (sr * ur(imz,i) - fr(imz,i)                        &
                          - us(ibx) * (gi * us(ibz) + vs(3) * vb)) / dv
                  us(ien) = (sr * (ur(ien,i) + ur(idn,i))                      &
                          - (fr(ien,i) + fr(idn,i))                            &
                          + pt * sm - us(ibx) * vb) / dv - us(idn)

! calculate normal component of momentum
!
                  us(imx) = (us(ien) + us(idn) + pt) * sm - us(ibx) * vb

! calculate the flux (eq. 14 in [1])
!
                  f(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))

                else

! intermediate flux is constant across the contact discontinuity so use HLL flux
!
                  f(:,i) = fh(:)

                end if ! sm == 0

              end if ! p* < 0

            end if ! sl < sm < sr

          end if ! nr < 1

        else ! Bx² > 0

! prepare the quadratic coefficients (eq. 47 in [1])
!
          a(1) = uh(imx)
          a(2) = - (fh(imx) + uh(ien) + uh(idn))
          a(3) = fh(ien) + fh(idn)

! solve the quadratic equation
!
          nr   = quadratic(a(1:3), x(1:2))

! if Δ < 0, just use the HLL flux
!
          if (nr < 1) then
            f(:,i) = fh(:)
          else

! get the contact dicontinuity speed
!
            sm = x(1)

! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
!
            if ((sm <= sl) .or. (sm >= sr)) then
              f(:,i) = fh(:)
            else

! calculate total pressure (eq. 48 in [1])
!
              pt = fh(imx) - (fh(ien) + fh(idn)) * sm

! if the pressure is negative, use the HLL flux
!
              if (pt <= 0.0d+00) then
                f(:,i) = fh(:)
              else

! depending in the sign of the contact dicontinuity speed, calculate the proper
! state and corresponding flux
!
                if (sm > 0.0d+00) then

! calculate the conserved variable vector (eqs. 43-46 in [1])
!
                  dv      = sl - sm
                  us(idn) = wl(idn) / dv
                  us(imy) = wl(imy) / dv
                  us(imz) = wl(imz) / dv
                  us(ien) = (wl(ien) + pt * sm) / dv
                  us(imx) = (us(ien) + pt) * sm
                  us(ien) = us(ien) - us(idn)
                  us(ibx) = 0.0d+00
                  us(iby) = wl(iby) / dv
                  us(ibz) = wl(ibz) / dv
                  us(ibp) = ql(ibp,i)

! calculate the flux (eq. 27 in [1])
!
                  f(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))

                else if (sm < 0.0d+00) then

! calculate the conserved variable vector (eqs. 43-46 in [1])
!
                  dv      = sr - sm
                  us(idn) = wr(idn) / dv
                  us(imy) = wr(imy) / dv
                  us(imz) = wr(imz) / dv
                  us(ien) = (wr(ien) + pt * sm) / dv
                  us(imx) = (us(ien) + pt) * sm
                  us(ien) = us(ien) - us(idn)
                  us(ibx) = 0.0d+00
                  us(iby) = wr(iby) / dv
                  us(ibz) = wr(ibz) / dv
                  us(ibp) = qr(ibp,i)

! calculate the flux (eq. 27 in [1])
!
                  f(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))

                else

! intermediate flux is constant across the contact discontinuity so use HLL flux
!
                  f(:,i) = fh(:)

                end if ! sm == 0

              end if ! p* < 0

            end if ! sl < sm < sr

          end if ! nr < 1

        end if ! Bx² == 0 or > 0

      end if ! sl < 0 < sr

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_srmhd_hllc
!
!===============================================================================
!
! subroutine RIEMANN_ISO_MHD_HLLD:
! -------------------------------
!
!   Subroutine solves one dimensional Riemann problem using the isothermal HLLD
!   method described by Mignone.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the Riemann states;
!     ul, ur - the conservative variables of the Riemann states;
!     fl, fr - the physical fluxes of the Riemann states;
!     f      - the Riemann average flux;
!
!   References:
!
!     [1] Mignone, A.,
!         "A simple and accurate Riemann solver for isothermal MHD",
!         Journal of Computational Physics, 2007, 225, pp. 1427-1441
!
!===============================================================================
!
  subroutine riemann_mhd_iso_hlld(ql, qr, ul, ur, fl, fr, cl, cr, f)

! include external procedures
!
    use equations, only : idn, imx, imy, imz, ibx, iby, ibz, ibp

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: ql, qr, ul, ur, fl, fr, cl, cr
    real(kind=8), dimension(:,:), intent(out) :: f

! local variables
!
    integer      :: i
    real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm, sdif
    real(kind=8) :: bx, b2, dn, dvl, dvr, ca

! local arrays to store the states
!
    real(kind=8), dimension(size(ql,1)) :: wl, wr, wcl, wcr, ui
!
!-------------------------------------------------------------------------------
!
! iterate over all points
!
    do i = 1, size(ql, 2)

! estimate the minimum and maximum speeds
!
      sl = min(cl(1,i), cr(1,i))
      sr = max(cl(2,i), cr(2,i))

! calculate the HLL flux
!
      if (sl >= 0.0d+00) then

        f(:,i) = fl(:,i)

      else if (sr <= 0.0d+00) then

        f(:,i) = fr(:,i)

      else ! sl < 0 < sr

! square of Bₓ, i.e. Bₓ²
!
        bx = ql(ibx,i)
        b2 = ql(ibx,i) * qr(ibx,i)

! calculate speed difference
!
        srml = sr - sl

! calculate vectors of the left and right-going waves
!
        wl(:) = sl * ul(:,i) - fl(:,i)
        wr(:) = sr * ur(:,i) - fr(:,i)

! the advection speed in the intermediate states
!
        dn =  wr(idn) - wl(idn)
        sm = (wr(imx) - wl(imx)) / dn

! speed differences
!
        slmm = sl - sm
        srmm = sr - sm

! calculate density of the intermediate states
!
        dn = dn / srml

! left and right Alfvén speeds
!
        ca  = abs(bx) / sqrt(dn)
        sml = sm - ca
        smr = sm + ca

! calculate denominators in order to check degeneracy
!
        dvl = slmm * wl(idn) - b2
        dvr = srmm * wr(idn) - b2

! check degeneracy Sl* -> Sl or Sr* -> Sr
!
        if (sml > sl .and. smr < sr .and. min(dvl, dvr) > 0.0d+00) then ! no degeneracies

! choose the correct state depending on the speed signs
!
          if (sml > 0.0d+00) then ! sl* > 0

! conservative variables for the outer left intermediate state
!
            ui(idn) = dn
            ui(imx) = dn * sm
            ui(imy) = dn * (slmm * wl(imy) - bx * wl(iby)) / dvl
            ui(imz) = dn * (slmm * wl(imz) - bx * wl(ibz)) / dvl
            ui(ibx) = bx
            ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
            ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
            ui(ibp) = ul(ibp,i)

! the outer left intermediate flux
!
            f(:,i)  = sl * ui(:) - wl(:)

          else if (smr < 0.0d+00) then ! sr* < 0

! conservative variables for the outer right intermediate state
!
            ui(idn) = dn
            ui(imx) = dn * sm
            ui(imy) = dn * (srmm * wr(imy) - bx * wr(iby)) / dvr
            ui(imz) = dn * (srmm * wr(imz) - bx * wr(ibz)) / dvr
            ui(ibx) = bx
            ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
            ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
            ui(ibp) = ur(ibp,i)

! the outer right intermediate flux
!
            f(:,i)  = sr * ui(:) - wr(:)

          else ! sl* <= 0 <= sr*

! separate cases for non-zero and zero Bx
!
            if (b2 > 0.0d+00) then

! conservative variables for the outer left intermediate state
!
              ui(idn) = dn
              ui(imx) = dn * sm
              ui(imy) = dn * (slmm * wl(imy) - bx * wl(iby)) / dvl
              ui(imz) = dn * (slmm * wl(imz) - bx * wl(ibz)) / dvl
              ui(ibx) = bx
              ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
              ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
              ui(ibp) = ul(ibp,i)

! vector of the left-going Alfvén wave
!
              wcl(:)  = (sml - sl) * ui(:) + wl(:)

! conservative variables for the outer right intermediate state
!
              ui(idn) = dn
              ui(imx) = dn * sm
              ui(imy) = dn * (srmm * wr(imy) - bx * wr(iby)) / dvr
              ui(imz) = dn * (srmm * wr(imz) - bx * wr(ibz)) / dvr
              ui(ibx) = bx
              ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
              ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
              ui(ibp) = ur(ibp,i)

! vector of the right-going Alfvén wave
!
              wcr(:)  = (smr - sr) * ui(:) + wr(:)

! the flux corresponding to the middle intermediate state
!
              f(:,i)  = 0.5d+00 * (sml * wcr(:) / ca - smr * wcl(:) / ca)

            else ! Bx = 0 -> Sₘ = 0

! no Alfvén wave and Sₘ = 0, so revert to the HLL flux
!
              f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml

            end if

          end if ! sl* < 0 < sr*

        else ! some degeneracies

! separate degeneracies
!
          if (sml > sl .and. dvl > 0.0d+00) then ! sr* > sr

! conservative variables for the outer left intermediate state
!
            ui(idn) = dn
            ui(imx) = dn * sm
            ui(imy) = dn * (slmm * wl(imy) - bx * wl(iby)) / dvl
            ui(imz) = dn * (slmm * wl(imz) - bx * wl(ibz)) / dvl
            ui(ibx) = bx
            ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
            ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
            ui(ibp) = ul(ibp,i)

! choose the correct state depending on the speed signs
!
            if (sml > 0.0d+00) then ! sl* > 0

! the outer left intermediate flux
!
              f(:,i)  = sl * ui(:) - wl(:)

            else ! sl* <= 0

! separate cases for non-zero and zero Bx
!
              if (b2 > 0.0d+00) then

! vector of the left-going Alfvén wave
!
                wcl(:)  = (sml - sl) * ui(:) + wl(:)

! the flux corresponding to the middle intermediate state
!
                sdif    = sr - sml
                f(:,i)  = sml * wr(:) / sdif - sr * wcl(:) / sdif

              else ! Bx = 0 -> Sₘ = 0

! no Alfvén wave and Sₘ = 0, so revert to the HLL flux
!
                f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml

              end if

            end if ! sl* < 0

          else if (smr < sr .and. dvr > 0.0d+00) then ! sl* < sl

! conservative variables for the outer right intermediate state
!
            ui(idn) = dn
            ui(imx) = dn * sm
            ui(imy) = dn * (srmm * wr(imy) - bx * wr(iby)) / dvr
            ui(imz) = dn * (srmm * wr(imz) - bx * wr(ibz)) / dvr
            ui(ibx) = bx
            ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
            ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
            ui(ibp) = ur(ibp,i)

! choose the correct state depending on the speed signs
!
            if (smr < 0.0d+00) then ! sr* < 0

! the outer right intermediate flux
!
              f(:,i)  = sr * ui(:) - wr(:)

            else ! sr* >= 0

! separate cases for non-zero and zero Bx
!
              if (b2 > 0.0d+00) then

! vector of the right-going Alfvén wave
!
                wcr(:)  = (smr - sr) * ui(:) + wr(:)

! the flux corresponding to the middle intermediate state
!
                sdif    = smr - sl
                f(:,i)  = sl * wcr(:) / sdif - smr * wl(:) / sdif

              else ! Bx = 0 -> Sₘ = 0

! no Alfvén wave and Sₘ = 0, so revert to the HLL flux
!
                f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml

              end if

            end if ! sr* > 0

          else ! sl* < sl & sr* > sr

! both outer states are degenerate so revert to the HLL flux
!
            f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml

          end if ! sl* < sl & sr* > sr

        end if ! degeneracies

      end if ! sl < 0 < sr

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_mhd_iso_hlld
!
!===============================================================================
!
! subroutine RIEMANN_ADI_MHD_HLLD:
! -------------------------------
!
!   Subroutine solves one dimensional Riemann problem using the adiabatic HLLD
!   method described by Miyoshi & Kusano.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the Riemann states;
!     ul, ur - the conservative variables of the Riemann states;
!     fl, fr - the physical fluxes of the Riemann states;
!     f      - the Riemann average flux;
!
!   References:
!
!     [1] Miyoshi, T. & Kusano, K.,
!         "A multi-state HLL approximate Riemann solver for ideal
!          magnetohydrodynamics",
!         Journal of Computational Physics, 2005, 208, pp. 315-344
!
!===============================================================================
!
  subroutine riemann_mhd_adi_hlld(ql, qr, ul, ur, fl, fr, cl, cr, f)

! include external procedures
!
    use equations, only : idn, imx, imy, imz, ien, ibx, iby, ibz, ibp

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: ql, qr, ul, ur, fl, fr, cl, cr
    real(kind=8), dimension(:,:), intent(out) :: f

! local variables
!
    integer      :: i
    real(kind=8) :: sl, sr, sm, srml, slmm, srmm
    real(kind=8) :: dn, bx, b2, pt, vy, vz, by, bz, vb
    real(kind=8) :: dnl, dnr, cal, car, sml, smr
    real(kind=8) :: dv, dvl, dvr

! local arrays to store the states
!
    real(kind=8), dimension(size(ql,1)) :: wl, wr, wcl, wcr, ui
!
!-------------------------------------------------------------------------------
!
! iterate over all points
!
    do i = 1, size(ql, 2)

! estimate the minimum and maximum speeds
!
      sl = min(cl(1,i), cr(1,i))
      sr = max(cl(2,i), cr(2,i))

! calculate the HLLD flux
!
      if (sl >= 0.0d+00) then ! sl ≥ 0

        f(:,i) = fl(:,i)

      else if (sr <= 0.0d+00) then ! sr ≤ 0

        f(:,i) = fr(:,i)

      else ! sl < 0 < sr

! square of Bₓ, i.e. Bₓ²
!
        bx = ql(ibx,i)
        b2 = ql(ibx,i) * qr(ibx,i)

! speed difference
!
        srml = sr - sl

! calculate vectors of the left and right-going waves
!
        wl(:) = sl * ul(:,i) - fl(:,i)
        wr(:) = sr * ur(:,i) - fr(:,i)

! the speed of contact discontinuity
!
        dn =  wr(idn) - wl(idn)
        sm = (wr(imx) - wl(imx)) / dn

! speed differences
!
        slmm = sl - sm
        srmm = sr - sm

! left and right state densities
!
        dnl = wl(idn) / slmm
        dnr = wr(idn) / srmm

! the total pressure, constant across the contact discontinuity and Alfvén waves
!
        pt = wl(idn) * wr(imx) / dn - wr(idn) * wl(imx) / dn + b2

! left and right Alfvén speeds
!
        cal = abs(bx) / sqrt(dnl)
        car = abs(bx) / sqrt(dnr)
        sml = sm - cal
        smr = sm + car

! calculate division factors (also used to determine degeneracies)
!
        dvl = slmm * wl(idn) - b2
        dvr = srmm * wr(idn) - b2

! check degeneracy Sl* -> Sl or Sr* -> Sr
!
        if (sml > sl .and. smr < sr .and. min(dvl, dvr) > 0.0d+00) then ! no degeneracy

! choose the correct state depending on the speed signs
!
          if (sml > 0.0d+00) then ! sl* > 0

! primitive variables for the outer left intermediate state
!
            vy = (   slmm * wl(imy) - bx * wl(iby)) / dvl
            vz = (   slmm * wl(imz) - bx * wl(ibz)) / dvl
            by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
            bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
            vb = sm * bx + vy * by + vz * bz

! conservative variables for the outer left intermediate state
!
            ui(idn) = dnl
            ui(imx) = dnl * sm
            ui(imy) = dnl * vy
            ui(imz) = dnl * vz
            ui(ibx) = bx
            ui(iby) = by
            ui(ibz) = bz
            ui(ibp) = ul(ibp,i)
            ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm

! the outer left intermediate flux
!
            f(:,i) = sl * ui(:) - wl(:)

          else if (smr < 0.0d+00) then ! sr* < 0

! primitive variables for the outer right intermediate state
!
            vy = (   srmm * wr(imy) - bx * wr(iby)) / dvr
            vz = (   srmm * wr(imz) - bx * wr(ibz)) / dvr
            by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
            bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
            vb = sm * bx + vy * by + vz * bz

! conservative variables for the outer right intermediate state
!
            ui(idn) = dnr
            ui(imx) = dnr * sm
            ui(imy) = dnr * vy
            ui(imz) = dnr * vz
            ui(ibx) = bx
            ui(iby) = by
            ui(ibz) = bz
            ui(ibp) = ur(ibp,i)
            ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm

! the outer right intermediate flux
!
            f(:,i) = sr * ui(:) - wr(:)

          else ! sl* < 0 < sr*

! separate cases with non-zero and zero Bx
!
            if (b2 > 0.0d+00) then

! primitive variables for the outer left intermediate state
!
              vy = (   slmm * wl(imy) - bx * wl(iby)) / dvl
              vz = (   slmm * wl(imz) - bx * wl(ibz)) / dvl
              by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
              bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
              vb = sm * bx + vy * by + vz * bz

! conservative variables for the outer left intermediate state
!
              ui(idn) = dnl
              ui(imx) = dnl * sm
              ui(imy) = dnl * vy
              ui(imz) = dnl * vz
              ui(ibx) = bx
              ui(iby) = by
              ui(ibz) = bz
              ui(ibp) = ul(ibp,i)
              ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm

! vector of the left-going Alfvén wave
!
              wcl(:) = (sml - sl) * ui(:) + wl(:)

! primitive variables for the outer right intermediate state
!
              vy = (   srmm * wr(imy) - bx * wr(iby)) / dvr
              vz = (   srmm * wr(imz) - bx * wr(ibz)) / dvr
              by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
              bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
              vb = sm * bx + vy * by + vz * bz

! conservative variables for the outer right intermediate state
!
              ui(idn) = dnr
              ui(imx) = dnr * sm
              ui(imy) = dnr * vy
              ui(imz) = dnr * vz
              ui(ibx) = bx
              ui(iby) = by
              ui(ibz) = bz
              ui(ibp) = ur(ibp,i)
              ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm

! vector of the right-going Alfvén wave
!
              wcr(:) = (smr - sr) * ui(:) + wr(:)

! prepare constant primitive variables of the intermediate states
!
              dv = abs(bx) * (sqrt(dnr) + sqrt(dnl))
              vy = (wcr(imy) - wcl(imy)) / dv
              vz = (wcr(imz) - wcl(imz)) / dv
              dv = car       + cal
              by = (wcr(iby) - wcl(iby)) / dv
              bz = (wcr(ibz) - wcl(ibz)) / dv
              vb = sm * bx + vy * by + vz * bz

! choose the correct state depending on the sign of contact discontinuity
! advection speed
!
              if (sm > 0.0d+00) then ! sm > 0

! conservative variables for the inmost left intermediate state
!
                ui(idn) = dnl
                ui(imx) = dnl * sm
                ui(imy) = dnl * vy
                ui(imz) = dnl * vz
                ui(ibx) = bx
                ui(iby) = by
                ui(ibz) = bz
                ui(ibp) = ul(ibp,i)
                ui(ien) = - (wcl(ien) + sm * pt - bx * vb) / cal

! the inmost left intermediate flux
!
                f(:,i) = sml * ui(:) - wcl(:)

              else if (sm < 0.0d+00) then ! sm < 0

! conservative variables for the inmost right intermediate state
!
                ui(idn) = dnr
                ui(imx) = dnr * sm
                ui(imy) = dnr * vy
                ui(imz) = dnr * vz
                ui(ibx) = bx
                ui(iby) = by
                ui(ibz) = bz
                ui(ibp) = ur(ibp,i)
                ui(ien) =   (wcr(ien) + sm * pt - bx * vb) / car

! the inmost right intermediate flux
!
                f(:,i) = smr * ui(:) - wcr(:)

              else ! sm = 0

! in the case when Sₘ = 0 and Bₓ² > 0, all variables are continuous, therefore
! the flux can be averaged from the Alfvén waves using a simple HLL formula;
!
                f(:,i) = sml * wcr(:) / dv - smr * wcl(:) / dv

              end if ! sm = 0

            else ! Bx = 0

              if (sm > 0.0d+00) then

! primitive variables for the outer left intermediate state
!
                vy =    slmm * wl(imy) / dvl
                vz =    slmm * wl(imz) / dvl
                by = wl(idn) * wl(iby) / dvl
                bz = wl(idn) * wl(ibz) / dvl
                vb = vy * by + vz * bz

! conservative variables for the outer left intermediate state
!
                ui(idn) = dnl
                ui(imx) = dnl * sm
                ui(imy) = dnl * vy
                ui(imz) = dnl * vz
                ui(ibx) = 0.0d+00
                ui(iby) = by
                ui(ibz) = bz
                ui(ibp) = ul(ibp,i)
                ui(ien) = (wl(ien) + sm * pt) / slmm

! the outer left intermediate flux
!
                f(:,i) = sl * ui(:) - wl(:)

              else if (sm < 0.0d+00) then

! primitive variables for the outer right intermediate state
!
                vy = (   srmm * wr(imy)) / dvr
                vz = (   srmm * wr(imz)) / dvr
                by = (wr(idn) * wr(iby)) / dvr
                bz = (wr(idn) * wr(ibz)) / dvr
                vb = vy * by + vz * bz

! conservative variables for the outer right intermediate state
!
                ui(idn) = dnr
                ui(imx) = dnr * sm
                ui(imy) = dnr * vy
                ui(imz) = dnr * vz
                ui(ibx) = 0.0d+00
                ui(iby) = by
                ui(ibz) = bz
                ui(ibp) = ur(ibp,i)
                ui(ien) = (wr(ien) + sm * pt) / srmm

! the outer right intermediate flux
!
                f(:,i) = sr * ui(:) - wr(:)

              else ! Sm = 0

! since Bx = 0 and Sm = 0, then revert to the HLL flux
!
                f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml

              end if

            end if ! Bx = 0

          end if ! sl* < 0 < sr*

        else ! some degeneracies

! separate degeneracies
!
          if (sml > sl .and. dvl > 0.0d+00) then ! sr* > sr

! primitive variables for the outer left intermediate state
!
            vy = (   slmm * wl(imy) - bx * wl(iby)) / dvl
            vz = (   slmm * wl(imz) - bx * wl(ibz)) / dvl
            by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
            bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
            vb = sm * bx + vy * by + vz * bz

! conservative variables for the outer left intermediate state
!
            ui(idn) = dnl
            ui(imx) = dnl * sm
            ui(imy) = dnl * vy
            ui(imz) = dnl * vz
            ui(ibx) = bx
            ui(iby) = by
            ui(ibz) = bz
            ui(ibp) = ul(ibp,i)
            ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm

! choose the correct state depending on the speed signs
!
            if (sml > 0.0d+00) then ! sl* > 0

! the outer left intermediate flux
!
              f(:,i)  = sl * ui(:) - wl(:)

            else ! sl* <= 0

! separate cases with non-zero and zero Bx
!
              if (b2 > 0.0d+00) then

! vector of the left-going Alfvén wave
!
                wcl(:)  = (sml - sl) * ui(:) + wl(:)

! primitive variables for the inner left intermediate state
!
                dv = srmm * dnr + cal * dnl
                vy = (wr(imy) - wcl(imy)) / dv
                vz = (wr(imz) - wcl(imz)) / dv
                dv = sr - sml
                by = (wr(iby) - wcl(iby)) / dv
                bz = (wr(ibz) - wcl(ibz)) / dv
                vb = sm * bx + vy * by + vz * bz

! conservative variables for the inner left intermediate state
!
                ui(idn) = dnl
                ui(imx) = dnl * sm
                ui(imy) = dnl * vy
                ui(imz) = dnl * vz
                ui(ibx) = bx
                ui(iby) = by
                ui(ibz) = bz
                ui(ibp) = ul(ibp,i)
                ui(ien) = - (wcl(ien) + sm * pt - bx * vb) / cal

! choose the correct state depending on the sign of contact discontinuity
! advection speed
!
                if (sm >= 0.0d+00) then ! sm >= 0

! the inner left intermediate flux
!
                  f(:,i)  = sml * ui(:) - wcl(:)

                else ! sm < 0

! vector of the left-going Alfvén wave
!
                  wcr(:)  = (sm - sml) * ui(:) + wcl(:)

! calculate the average flux over the right inner intermediate state
!
                  f(:,i)  = sm * wr(:) / srmm - sr * wcr(:) / srmm

                end if ! sm < 0

              else ! bx = 0

! no Alfvén wave, so revert to the HLL flux
!
                f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml

              end if

            end if ! sl* < 0

          else if (smr < sr .and. dvr > 0.0d+00) then ! sl* < sl

! primitive variables for the outer right intermediate state
!
            vy = (   srmm * wr(imy) - bx * wr(iby)) / dvr
            vz = (   srmm * wr(imz) - bx * wr(ibz)) / dvr
            by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
            bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
            vb = sm * bx + vy * by + vz * bz

! conservative variables for the outer right intermediate state
!
            ui(idn) = dnr
            ui(imx) = dnr * sm
            ui(imy) = dnr * vy
            ui(imz) = dnr * vz
            ui(ibx) = bx
            ui(iby) = by
            ui(ibz) = bz
            ui(ibp) = ur(ibp,i)
            ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm

! choose the correct state depending on the speed signs
!
            if (smr < 0.0d+00) then ! sr* < 0

! the outer right intermediate flux
!
              f(:,i)  = sr * ui(:) - wr(:)

            else ! sr* > 0

! separate cases with non-zero and zero Bx
!
              if (b2 > 0.0d+00) then

! vector of the right-going Alfvén wave
!
                wcr(:)  = (smr - sr) * ui(:) + wr(:)

! primitive variables for the inner right intermediate state
!
                dv = slmm * dnl - car * dnr
                vy = (wl(imy) - wcr(imy)) / dv
                vz = (wl(imz) - wcr(imz)) / dv
                dv = sl - smr
                by = (wl(iby) - wcr(iby)) / dv
                bz = (wl(ibz) - wcr(ibz)) / dv
                vb = sm * bx + vy * by + vz * bz

! conservative variables for the inner left intermediate state
!
                ui(idn) = dnr
                ui(imx) = dnr * sm
                ui(imy) = dnr * vy
                ui(imz) = dnr * vz
                ui(ibx) = bx
                ui(iby) = by
                ui(ibz) = bz
                ui(ibp) = ur(ibp,i)
                ui(ien) =   (wcr(ien) + sm * pt - bx * vb) / car

! choose the correct state depending on the sign of contact discontinuity
! advection speed
!
                if (sm <= 0.0d+00) then ! sm <= 0

! the inner right intermediate flux
!
                  f(:,i)  = smr * ui(:) - wcr(:)

                else ! sm > 0

! vector of the right-going Alfvén wave
!
                    wcl(:)  = (sm - smr) * ui(:) + wcr(:)

! calculate the average flux over the left inner intermediate state
!
                    f(:,i)  = sm * wl(:) / slmm - sl * wcl(:) / slmm

                end if ! sm > 0

              else ! bx = 0

! no Alfvén wave, so revert to the HLL flux
!
                f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml

              end if

            end if ! sr* > 0

          else ! sl* < sl & sr* > sr

! so far we revert to HLL flux in the case of degeneracies
!
            f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml

          end if ! sl* < sl & sr* > sr

        end if ! deneneracies

      end if ! sl < 0 < sr

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_mhd_adi_hlld
!
!===============================================================================
!
! subroutine RIEMANN_HD_ISO_ROE:
! -----------------------------
!
!   Subroutine solves one dimensional Riemann problem using
!   the Roe's method.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the Riemann states;
!     ul, ur - the conservative variables of the Riemann states;
!     fl, fr - the physical fluxes of the Riemann states;
!     f      - the Riemann average flux;
!
!   References:
!
!     [1] Roe, P. L.
!         "Approximate Riemann Solvers, Parameter Vectors, and Difference
!          Schemes",
!         Journal of Computational Physics, 1981, 43, pp. 357-372
!     [2] Toro, E. F.,
!         "Riemann Solvers and Numerical Methods for Fluid dynamics",
!         Springer-Verlag Berlin Heidelberg, 2009
!
!===============================================================================
!
  subroutine riemann_hd_iso_roe(ql, qr, ul, ur, fl, fr, cl, cr, f)

! include external procedures
!
    use equations, only : idn, ivx, ivy, ivz
    use equations, only : eigensystem_roe

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: ql, qr, ul, ur, fl, fr, cl, cr
    real(kind=8), dimension(:,:), intent(out) :: f

! local variables
!
    integer      :: nf, i, p
    real(kind=8) :: sdl, sdr, sds

! local arrays to store the states
!
    real(kind=8), dimension(size(ql,1))            :: qi, ci, al
    real(kind=8), dimension(size(ql,1),size(ql,1)) :: li, ri
!
!-------------------------------------------------------------------------------
!
! get the number of fluxes
!
    nf = size(ql, 1)

! iterate over all points
!
    do i = 1, size(ql, 2)

! calculate variables for the Roe intermediate state
!
      sdl = sqrt(ql(idn,i))
      sdr = sqrt(qr(idn,i))
      sds = sdl + sdr

! prepare the Roe intermediate state vector (eq. 11.60 in [2])
!
      qi(idn)     = sdl * sdr
      qi(ivx:ivz) = sdl * ql(ivx:ivz,i) / sds + sdr * qr(ivx:ivz,i) / sds

! obtain eigenvalues and eigenvectors
!
      call eigensystem_roe(0.0d+00, 0.0d+00, qi(:), ci(:), ri(:,:), li(:,:))

! return upwind fluxes
!
      if (ci(1) >= 0.0d+00) then

        f(:,i) = fl(:,i)

      else if (ci(nf) <= 0.0d+00) then

        f(:,i) = fr(:,i)

      else

! calculate wave amplitudes α = L.ΔU
!
        al(:) = 0.0d+00
        do p = 1, nf
          al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i))
        end do

! calculate the flux average
!
        f(:,i) = 0.5d+00 * (fl(:,i) + fr(:,i))

! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K)
!
        if (qi(ivx) >= 0.0d+00) then
          do p = 1, nf
            f(:,i) = f(:,i) - 0.5d+00 * al(p) * abs(ci(p)) * ri(p,:)
          end do
        else
          do p = nf, 1, - 1
            f(:,i) = f(:,i) - 0.5d+00 * al(p) * abs(ci(p)) * ri(p,:)
          end do
        end if

      end if ! sl < 0 < sr

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_hd_iso_roe
!
!===============================================================================
!
! subroutine RIEMANN_HD_ISO_ROE:
! -----------------------------
!
!   Subroutine solves one dimensional Riemann problem using
!   the Roe's method.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the Riemann states;
!     ul, ur - the conservative variables of the Riemann states;
!     fl, fr - the physical fluxes of the Riemann states;
!     f      - the Riemann average flux;
!
!   References:
!
!     [1] Roe, P. L.
!         "Approximate Riemann Solvers, Parameter Vectors, and Difference
!          Schemes",
!         Journal of Computational Physics, 1981, 43, pp. 357-372
!     [2] Toro, E. F.,
!         "Riemann Solvers and Numerical Methods for Fluid dynamics",
!         Springer-Verlag Berlin Heidelberg, 2009
!
!===============================================================================
!
  subroutine riemann_hd_adi_roe(ql, qr, ul, ur, fl, fr, cl, cr, f)

! include external procedures
!
    use equations, only : idn, ivx, ivy, ivz, ipr, ien
    use equations, only : eigensystem_roe

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: ql, qr, ul, ur, fl, fr, cl, cr
    real(kind=8), dimension(:,:), intent(out) :: f

! local variables
!
    integer      :: nf, i, p
    real(kind=8) :: sdl, sdr, sds

! local arrays to store the states
!
    real(kind=8), dimension(size(ql,1))            :: qi, ci, al
    real(kind=8), dimension(size(ql,1),size(ql,1)) :: li, ri
!
!-------------------------------------------------------------------------------
!
! get the number of fluxes
!
    nf = size(ql, 1)

! iterate over all points
!
    do i = 1, size(ql, 2)

! calculate variables for the Roe intermediate state
!
      sdl = sqrt(ql(idn,i))
      sdr = sqrt(qr(idn,i))
      sds = sdl + sdr

! prepare the Roe intermediate state vector (eq. 11.60 in [2])
!
      qi(idn)     = sdl * sdr
      qi(ivx:ivz) = sdl * ql(ivx:ivz,i) / sds + sdr * qr(ivx:ivz,i) / sds
      qi(ipr)     = ((ul(ien,i) + ql(ipr,i)) / sdl                             &
                  +  (ur(ien,i) + qr(ipr,i)) / sdr) / sds

! obtain eigenvalues and eigenvectors
!
      call eigensystem_roe(0.0d+00, 0.0d+00, qi(:), ci(:), ri(:,:), li(:,:))

! return upwind fluxes
!
      if (ci(1) >= 0.0d+00) then

        f(:,i) = fl(:,i)

      else if (ci(nf) <= 0.0d+00) then

        f(:,i) = fr(:,i)

      else

! calculate wave amplitudes α = L.ΔU
!
        al(:) = 0.0d+00
        do p = 1, nf
          al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i))
        end do

! calculate the flux average
!
        f(:,i) = 0.5d+00 * (fl(:,i) + fr(:,i))

! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K)
!
        if (qi(ivx) >= 0.0d+00) then
          do p = 1, nf
            f(:,i) = f(:,i) - 0.5d+00 * al(p) * abs(ci(p)) * ri(p,:)
          end do
        else
          do p = nf, 1, - 1
            f(:,i) = f(:,i) - 0.5d+00 * al(p) * abs(ci(p)) * ri(p,:)
          end do
        end if

      end if ! sl < 0 < sr

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_hd_adi_roe
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ISO_ROE:
! ------------------------------
!
!   Subroutine solves one dimensional Riemann problem using
!   the Roe's method.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the Riemann states;
!     ul, ur - the conservative variables of the Riemann states;
!     fl, fr - the physical fluxes of the Riemann states;
!     f      - the Riemann average flux;
!
!   References:
!
!     [1] Stone, J. M. & Gardiner, T. A.,
!         "ATHENA: A New Code for Astrophysical MHD",
!         The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177
!     [2] Toro, E. F.,
!         "Riemann Solvers and Numerical Methods for Fluid dynamics",
!         Springer-Verlag Berlin Heidelberg, 2009
!
!===============================================================================
!
  subroutine riemann_mhd_iso_roe(ql, qr, ul, ur, fl, fr, cl, cr, f)

! include external procedures
!
    use equations, only : idn, ivx, ivy, ivz, imx, imy, imz, ibx, iby, ibz, ibp
    use equations, only : eigensystem_roe

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: ql, qr, ul, ur, fl, fr, cl, cr
    real(kind=8), dimension(:,:), intent(out) :: f

! local variables
!
    integer      :: nf, i, p
    real(kind=8) :: sdl, sdr, sds
    real(kind=8) :: xx, yy

! local arrays to store the states
!
    real(kind=8), dimension(size(ql,1))            :: qi, ci, al
    real(kind=8), dimension(size(ql,1),size(ql,1)) :: li, ri
!
!-------------------------------------------------------------------------------
!
! get the number of fluxes
!
    nf = size(ql, 1)

! iterate over all points
!
    do i = 1, size(ql, 2)

! calculate variables for the Roe intermediate state
!
      sdl = sqrt(ql(idn,i))
      sdr = sqrt(qr(idn,i))
      sds = sdl + sdr

! prepare the Roe intermediate state vector (eq. 11.60 in [2])
!
      qi(idn)     = sdl * sdr
      qi(ivx:ivz) = sdl * ql(ivx:ivz,i) / sds + sdr * qr(ivx:ivz,i) / sds
      qi(ibx)     = ql(ibx,i)
      qi(iby:ibz) = sdr * ql(iby:ibz,i) / sds + sdl * qr(iby:ibz,i) / sds
      qi(ibp)     = ql(ibp,i)

! prepare coefficients
!
      xx = 0.5d+00 * ((ql(iby,i) - qr(iby,i))**2                               &
                                        + (ql(ibz,i) - qr(ibz,i))**2) / sds**2
      yy = 0.5d+00 * (ql(idn,i) + qr(idn,i)) / qi(idn)

! obtain eigenvalues and eigenvectors
!
      call eigensystem_roe(xx, yy, qi(:), ci(:), ri(:,:), li(:,:))

! return upwind fluxes
!
      if (ci(1) >= 0.0d+00) then

        f(:,i) = fl(:,i)

      else if (ci(nf) <= 0.0d+00) then

        f(:,i) = fr(:,i)

      else

! prepare fluxes which do not change across the states
!
        f(ibx,i) = fl(ibx,i)
        f(ibp,i) = fl(ibp,i)

! calculate wave amplitudes α = L.ΔU
!
        al(:) = 0.0d+00
        do p = 1, nf
          al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i))
        end do

! calculate the flux average
!
        f(idn,i) = 0.5d+00 * (fl(idn,i) + fr(idn,i))
        f(imx,i) = 0.5d+00 * (fl(imx,i) + fr(imx,i))
        f(imy,i) = 0.5d+00 * (fl(imy,i) + fr(imy,i))
        f(imz,i) = 0.5d+00 * (fl(imz,i) + fr(imz,i))
        f(iby,i) = 0.5d+00 * (fl(iby,i) + fr(iby,i))
        f(ibz,i) = 0.5d+00 * (fl(ibz,i) + fr(ibz,i))

! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K)
!
        if (qi(ivx) >= 0.0d+00) then
          do p = 1, nf
            xx       = - 0.5d+00 * al(p) * abs(ci(p))

            f(idn,i) = f(idn,i) + xx * ri(p,idn)
            f(imx,i) = f(imx,i) + xx * ri(p,imx)
            f(imy,i) = f(imy,i) + xx * ri(p,imy)
            f(imz,i) = f(imz,i) + xx * ri(p,imz)
            f(iby,i) = f(iby,i) + xx * ri(p,iby)
            f(ibz,i) = f(ibz,i) + xx * ri(p,ibz)
          end do
        else
          do p = nf, 1, -1
            xx       = - 0.5d+00 * al(p) * abs(ci(p))

            f(idn,i) = f(idn,i) + xx * ri(p,idn)
            f(imx,i) = f(imx,i) + xx * ri(p,imx)
            f(imy,i) = f(imy,i) + xx * ri(p,imy)
            f(imz,i) = f(imz,i) + xx * ri(p,imz)
            f(iby,i) = f(iby,i) + xx * ri(p,iby)
            f(ibz,i) = f(ibz,i) + xx * ri(p,ibz)
          end do
        end if

      end if ! sl < 0 < sr

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_mhd_iso_roe
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ADI_ROE:
! ------------------------------
!
!   Subroutine solves one dimensional Riemann problem using
!   the Roe's method.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the Riemann states;
!     ul, ur - the conservative variables of the Riemann states;
!     fl, fr - the physical fluxes of the Riemann states;
!     f      - the Riemann average flux;
!
!   References:
!
!     [1] Stone, J. M. & Gardiner, T. A.,
!         "ATHENA: A New Code for Astrophysical MHD",
!         The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177
!     [2] Toro, E. F.,
!         "Riemann Solvers and Numerical Methods for Fluid dynamics",
!         Springer-Verlag Berlin Heidelberg, 2009
!
!===============================================================================
!
  subroutine riemann_mhd_adi_roe(ql, qr, ul, ur, fl, fr, cl, cr, f)

! include external procedures
!
    use equations, only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
    use equations, only : ibx, iby, ibz, ibp
    use equations, only : eigensystem_roe

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: ql, qr, ul, ur, fl, fr, cl, cr
    real(kind=8), dimension(:,:), intent(out) :: f

! local variables
!
    integer      :: nf, i, p
    real(kind=8) :: sdl, sdr, sds
    real(kind=8) :: xx, yy
    real(kind=8) :: pml, pmr

! local arrays to store the states
!
    real(kind=8), dimension(size(ql,1))            :: qi, ci, al
    real(kind=8), dimension(size(ql,1),size(ql,1)) :: li, ri
!
!-------------------------------------------------------------------------------
!
! get the number of fluxes
!
    nf = size(ql, 1)

! iterate over all points
!
    do i = 1, size(ql, 2)

! calculate variables for the Roe intermediate state
!
      sdl = sqrt(ql(idn,i))
      sdr = sqrt(qr(idn,i))
      sds = sdl + sdr

! prepare magnetic pressures
!
      pml = 0.5d+00 * sum(ql(ibx:ibz,i)**2)
      pmr = 0.5d+00 * sum(qr(ibx:ibz,i)**2)

! prepare the Roe intermediate state vector (eq. 11.60 in [2])
!
      qi(idn)     = sdl * sdr
      qi(ivx:ivz) = sdl * ql(ivx:ivz,i) / sds + sdr * qr(ivx:ivz,i) / sds
      qi(ipr)     = ((ul(ien,i) + ql(ipr,i) + pml) / sdl                       &
                  +  (ur(ien,i) + qr(ipr,i) + pmr) / sdr) / sds
      qi(ibx)     = ql(ibx,i)
      qi(iby:ibz) = sdr * ql(iby:ibz,i) / sds + sdl * qr(iby:ibz,i) / sds
      qi(ibp)     = ql(ibp,i)

! prepare coefficients
!
      xx = 0.5d+00 * ((ql(iby,i) - qr(iby,i))**2                               &
                                        + (ql(ibz,i) - qr(ibz,i))**2) / sds**2
      yy = 0.5d+00 * (ql(idn,i) + qr(idn,i)) / qi(idn)

! obtain eigenvalues and eigenvectors
!
      call eigensystem_roe(xx, yy, qi(:), ci(:), ri(:,:), li(:,:))

! return upwind fluxes
!
      if (ci(1) >= 0.0d+00) then

        f(:,i) = fl(:,i)

      else if (ci(nf) <= 0.0d+00) then

        f(:,i) = fr(:,i)

      else

! prepare fluxes which do not change across the states
!
        f(ibx,i) = fl(ibx,i)
        f(ibp,i) = fl(ibp,i)

! calculate wave amplitudes α = L.ΔU
!
        al(:) = 0.0d+00
        do p = 1, nf
          al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i))
        end do

! calculate the flux average
!
        f(idn,i) = 0.5d+00 * (fl(idn,i) + fr(idn,i))
        f(imx,i) = 0.5d+00 * (fl(imx,i) + fr(imx,i))
        f(imy,i) = 0.5d+00 * (fl(imy,i) + fr(imy,i))
        f(imz,i) = 0.5d+00 * (fl(imz,i) + fr(imz,i))
        f(ien,i) = 0.5d+00 * (fl(ien,i) + fr(ien,i))
        f(iby,i) = 0.5d+00 * (fl(iby,i) + fr(iby,i))
        f(ibz,i) = 0.5d+00 * (fl(ibz,i) + fr(ibz,i))

! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K)
!
        if (qi(ivx) >= 0.0d+00) then
          do p = 1, nf
            xx       = - 0.5d+00 * al(p) * abs(ci(p))

            f(idn,i) = f(idn,i) + xx * ri(p,idn)
            f(imx,i) = f(imx,i) + xx * ri(p,imx)
            f(imy,i) = f(imy,i) + xx * ri(p,imy)
            f(imz,i) = f(imz,i) + xx * ri(p,imz)
            f(ien,i) = f(ien,i) + xx * ri(p,ien)
            f(iby,i) = f(iby,i) + xx * ri(p,iby)
            f(ibz,i) = f(ibz,i) + xx * ri(p,ibz)
          end do
        else
          do p = nf, 1, -1
            xx       = - 0.5d+00 * al(p) * abs(ci(p))

            f(idn,i) = f(idn,i) + xx * ri(p,idn)
            f(imx,i) = f(imx,i) + xx * ri(p,imx)
            f(imy,i) = f(imy,i) + xx * ri(p,imy)
            f(imz,i) = f(imz,i) + xx * ri(p,imz)
            f(ien,i) = f(ien,i) + xx * ri(p,ien)
            f(iby,i) = f(iby,i) + xx * ri(p,iby)
            f(ibz,i) = f(ibz,i) + xx * ri(p,ibz)
          end do
        end if

      end if ! sl < 0 < sr

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_mhd_adi_roe
!
!===============================================================================
!
! subroutine HIGHER_ORDER_FLUX_CORRECTION:
! ---------------------------------------
!
!   Subroutine adds higher order corrections to the numerical fluxes.
!
!   Arguments:
!
!     f - the vector of numerical fluxes;
!
!===============================================================================
!
  subroutine higher_order_flux_correction(f)

! include external procedures
!
    use interpolations, only : order

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:,:), intent(inout)  :: f

! local variables
!
    integer :: n

! local arrays to store the states
!
    real(kind=8), dimension(size(f,1),size(f,2)) :: f2, f4
!
!-------------------------------------------------------------------------------
!
! depending on the scheme order calculate and add higher order corrections,
! if required
!
    if (.not. high_order_flux .or. order <= 2) return

! get the length of vector
!
    n = size(f, 2)

! 2nd order correction
!
    f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
    f2(:,1    ) = 0.0d+00
    f2(:,  n  ) = 0.0d+00

! 4th order correction
!
    if (order > 4) then

      f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n  )                         &
                  - 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1))                        &
                  + 6.0d+00 *  f(:,3:n-2)) / 6.4d+02
      f4(:,1    ) = 0.0d+00
      f4(:,2    ) = 0.0d+00
      f4(:,  n-1) = 0.0d+00
      f4(:,  n  ) = 0.0d+00

! correct the flux for 4th order
!
      f(:,:) = f(:,:) + f4(:,:)

    end if

! correct the flux for 2nd order
!
    f(:,:) =  f(:,:) + f2(:,:)

!-------------------------------------------------------------------------------
!
  end subroutine higher_order_flux_correction

!===============================================================================
!
end module schemes