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-2022 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"

! KEPES solver
!
  logical            , save :: kepes = .false.

! 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, fi)
      real(kind=8), dimension(:,:), intent(in)  :: ql, qr
      real(kind=8), dimension(:,:), intent(out) :: fi
    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 helpers   , only : print_message, uppercase
    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"


    character(len=256) :: msg

    character(len=*), parameter :: fmt = "('The selected Riemann solver " // &
                             "is not implemented for the ',a,': ',a,'.'" // &
                             ",a,4x,'Available Riemann solvers are: ', a, '.')"
    character(len=*), parameter :: loc = 'SCHEMES::initialize_schemes()'

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

    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")

      select case(trim(eos))

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

        select case(trim(solver))

        case("hll", "HLL")
          name_sol =  "HLL"
          riemann_solver => riemann_hll

        case("roe", "ROE")
          name_sol =  "ROE"
          riemann_solver => riemann_hd_iso_roe

        case("kepes", "KEPES")
          name_sol =  "KEPES"
          kepes    = .true.
          riemann_solver => riemann_hd_iso_kepes

        case default

          if (verbose) then
            write(msg,fmt) 'isothermal hydrodynamics', &
                           "'" // uppercase(trim(solver)) // "'", &
                           new_line('a'), &
                           "'HLL', 'ROE', 'KEPES'"
            call print_message(loc, msg)
          end if
          status = 1

        end select

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

        select case(trim(solver))

        case("hll", "HLL")
          name_sol =  "HLL"
          riemann_solver => riemann_hll

        case("hllc", "HLLC")
          name_sol =  "HLLC"
          riemann_solver => riemann_hd_hllc

        case("roe", "ROE")
          name_sol =  "ROE"
          riemann_solver => riemann_hd_adi_roe

        case("kepes", "KEPES")
          name_sol =  "KEPES"
          kepes    = .true.
          riemann_solver => riemann_hd_adi_kepes

        case default

          if (verbose) then
            write(msg,fmt) 'adiabatic hydrodynamics', &
                           "'" // uppercase(trim(solver)) // "'", &
                           new_line('a'), &
                           "'HLL', 'HLLC', 'ROE', 'KEPES'"
            call print_message(loc, msg)
          end if
          status = 1

        end select

      end select

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

      select case(trim(eos))

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

        select case(trim(solver))

        case("hll", "HLL")
          name_sol =  "HLL"
          riemann_solver => riemann_hll

        case ("hlld", "HLLD")
          name_sol =  "HLLD"
          riemann_solver => riemann_mhd_iso_hlld

        case("roe", "ROE")
          name_sol =  "ROE"
          riemann_solver => riemann_mhd_iso_roe

        case("kepes", "KEPES")
          name_sol =  "KEPES"
          kepes    = .true.
          riemann_solver => riemann_mhd_iso_kepes

        case default

          if (verbose) then
            write(msg,fmt) 'isothermal magnetohydrodynamics', &
                           "'" // uppercase(trim(solver)) // "'", &
                           new_line('a'), &
                           "'HLL', 'HLLD', 'ROE', 'KEPES'"
            call print_message(loc, msg)
          end if
          status = 1

        end select

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

        select case(trim(solver))

        case ("hll", "HLL")
          name_sol =  "HLL"
          riemann_solver => riemann_hll

        case ("hllc", "HLLC")
          name_sol =  "HLLC"
          riemann_solver => riemann_mhd_hllc

        case ("hlld", "HLLD")
          name_sol =  "HLLD"
          riemann_solver => riemann_mhd_adi_hlld

        case("roe", "ROE")
          name_sol =  "ROE"
          riemann_solver => riemann_mhd_adi_roe

        case("kepes", "KEPES")
          name_sol =  "KEPES"
          kepes    = .true.
          riemann_solver => riemann_mhd_adi_kepes

        case default

          if (verbose) then
            write(msg,fmt) 'adiabatic magnetohydrodynamics', &
                           "'" // uppercase(trim(solver)) // "'", &
                           new_line('a'), &
                           "'HLL', 'HLLC', 'HLLD', 'ROE', 'KEPES'"
            call print_message(loc, msg)
          end if
          status = 1

        end select

      end select

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

      select case(trim(eos))

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

        select case(trim(statev))

        case("4vec", "4-vector", "4VEC", "4-VECTOR")
          name_sts =  "4-vector"
          states_4vec = .true.

        case default
          name_sts =  "primitive"

        end select

        select case(trim(solver))

        case ("hll", "HLL")
          name_sol =  "HLL"
          riemann_solver => riemann_hll

        case("hllc", "HLLC", "hllcm", "HLLCM", "hllc-m", "HLLC-M")
          name_sol =  "HLLC (Mignone & Bodo 2005)"
          riemann_solver => riemann_srhd_hllc

        case default

          if (verbose) then
            write(msg,fmt) 'adiabatic special relativity hydrodynamics', &
                           "'" // uppercase(trim(solver)) // "'", &
                           new_line('a'), &
                           "'HLL', 'HLLC'"
            call print_message(loc, msg)
          end if
          status = 1

        end select

      end select

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

      select case(trim(eos))

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

        select case(trim(statev))

        case("4vec", "4-vector", "4VEC", "4-VECTOR")
          name_sts =  "4-vector"
          states_4vec = .true.

        case default
          name_sts =  "primitive"

        end select

        select case(trim(solver))

        case ("hll", "HLL")
          name_sol =  "HLL"
          riemann_solver => riemann_hll

        case("hllc", "HLLC", "hllcm", "HLLCM", "hllc-m", "HLLC-M")
          name_sol =  "HLLC (Mignone & Bodo 2006)"
          riemann_solver => riemann_srmhd_hllc

        case default

          if (verbose) then
            write(msg,fmt) 'adiabatic special relativity magnetohydrodynamics',&
                           "'" // uppercase(trim(solver)) // "'", &
                           new_line('a'), &
                           "'HLL', 'HLLC'"
            call print_message(loc, msg)
          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, l, p
    real(kind=8) :: vm

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

!-------------------------------------------------------------------------------
!
#if NDIMS == 2
    k = 1
#endif /* NDIMS == 2 */

! 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 coordinates, only : nn => bcells
    use equations  , only : nf, ibx, ibp
    use equations  , only : magnetized, cmax

    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(nf,nn) :: ql, qr

!-------------------------------------------------------------------------------
!
    if (kepes) then
      call riemann_solver(q(:,:,1), q(:,:,2), f(:,:))
    else

! 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        = 5.0d-01 * ((qr(ibx,i) + ql(ibx,i)) &
                                             - (qr(ibp,i) - ql(ibp,i)) / cmax)
          bp        = 5.0d-01 * ((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

! get the numerical fluxes
!
      call riemann_solver(ql(:,:), qr(:,:), f(:,:))

! higher order flux corrections
!
      call higher_order_flux_correction(f)

    end if

!-------------------------------------------------------------------------------
!
  end subroutine numerical_flux
!
!*******************************************************************************
!
! GENERIC RIEMANN SOLVERS
!
!*******************************************************************************
!
!===============================================================================
!
! 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 left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   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,
!         https://doi.org/10.1137/1025002
!
!===============================================================================
!
  subroutine riemann_hll(ql, qr, fi)

    use coordinates, only : nn => bcells
    use equations  , only : nf
    use equations  , only : prim2cons, fluxspeed

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

    integer      :: i
    real(kind=8) :: sl, sr

    real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
    real(kind=8), dimension( 2,nn) :: cl, cr
    real(kind=8), dimension(nf   ) :: wl, wr

!-------------------------------------------------------------------------------
!
    call prim2cons(ql, ul)
    call prim2cons(qr, ur)

    call fluxspeed(ql, ul, fl, cl)
    call fluxspeed(qr, ur, fr, cr)

    do i = 1, nn

      sl = min(cl(1,i), cr(1,i))
      sr = max(cl(2,i), cr(2,i))

      if (sl >= 0.0d+00) then

        fi(:,i) = fl(:,i)

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

        fi(:,i) = fr(:,i)

      else ! sl < 0 < sr

        wl = sl * ul(:,i) - fl(:,i)
        wr = sr * ur(:,i) - fr(:,i)

        fi(:,i) = (sl * wr - sr * wl) / (sr - sl)

      end if ! sl < 0 < sr

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_hll
!
!*******************************************************************************
!
! ISOTHERMAL HYDRODYNAMIC RIEMANN SOLVERS
!
!*******************************************************************************
!
!===============================================================================
!
! subroutine RIEMANN_HD_ISO_ROE:
! -----------------------------
!
!   Subroutine solves one dimensional Riemann problem using
!   the Roe's method.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   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, fi)

    use coordinates, only : nn => bcells
    use equations  , only : nf, csnd
    use equations  , only : idn, ivx, ivy, ivz
    use equations  , only : prim2cons, fluxspeed

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

    integer      :: i
    real(kind=8) :: sdl, sdr, sds, cs, ch

    real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
    real(kind=8), dimension(nf   ) :: qs, lm, al, df

    logical                     , save :: first = .true.
    real(kind=8)                , save :: chi
    real(kind=8), dimension(4,4), save :: rvec, lvec
!$omp threadprivate(first, chi, rvec, lvec)

!-------------------------------------------------------------------------------
!
    if (first) then
      chi = 5.0d-01 / csnd

      rvec(:,1) = [ 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
      rvec(:,2) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
      rvec(:,3) = [ 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
      rvec(:,4) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]

      lvec(:,1) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
      lvec(:,2) = [   chi  , 0.0d+00, 0.0d+00, - chi   ]
      lvec(:,3) = [ 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
      lvec(:,4) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]

      first = .false.
    end if

    call prim2cons(ql, ul)
    call prim2cons(qr, ur)

    call fluxspeed(ql, ul, fl)
    call fluxspeed(qr, ur, fr)

    do i = 1, nn

      sdl = sqrt(ql(idn,i))
      sdr = sqrt(qr(idn,i))
      sds = sdl + sdr

      qs(idn)     = sdl * sdr
      qs(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds

      cs = sign(csnd, qs(ivx))
      ch = sign(chi , qs(ivx))

      lm(1) = qs(ivx) + cs
      lm(2) = qs(ivx)
      lm(3) = qs(ivx)
      lm(4) = qs(ivx) - cs

      rvec(1,2) = lm(1)
      rvec(4,2) = lm(4)
      rvec(1,3) = qs(ivy)
      rvec(4,3) = qs(ivy)
      rvec(1,4) = qs(ivz)
      rvec(4,4) = qs(ivz)

      lvec(1,1) = - ch * lm(4)
      lvec(2,1) = - qs(ivy)
      lvec(3,1) = - qs(ivz)
      lvec(4,1) =   ch * lm(1)
      lvec(1,2) =   ch
      lvec(4,2) = - ch

      al = abs(lm) * matmul(lvec, ur(:,i) - ul(:,i))
      df = matmul(al, rvec)
      fi(:,i) = 5.0d-01 * ((fl(:,i) + fr(:,i)) - df)

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_hd_iso_roe
!
!===============================================================================
!
! subroutine RIEMANN_HD_ISO_KEPES:
! -------------------------------
!
!   Subroutine solves one dimensional isothermal hydrodynamic Riemann problem
!   using the entropy stable KEPES method.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   References:
!
!     [1] Winters, A. R., Czernik, C., Schily, M. B., Gassner, G. J.,
!         "Entropy stable numerical approximations for the isothermal and
!          polytropic Euler equations",
!         BIT Numerical Mathematics (2020) 60:791–824,
!         https://doi.org/10.1007/s10543-019-00789-w
!
!===============================================================================
!
  subroutine riemann_hd_iso_kepes(ql, qr, fi)

    use coordinates, only : nn => bcells
    use equations  , only : nf, csnd, csnd2
    use equations  , only : idn, ivx, ivy, ivz, imx, imy, imz

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

    integer      :: i
    real(kind=8) :: v2l, v2r, cs
    real(kind=8) :: dnl, dna, vxa, vya, vza

    real(kind=8), dimension(nf) :: v, lm, tm

    logical                     , save :: first = .true.
    real(kind=8), dimension(4,4), save :: rm
!$omp threadprivate(first, rm)

!-------------------------------------------------------------------------------
!
    if (first) then
      rm(:,1) = [ 1.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
      rm(:,2) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
      rm(:,3) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
      rm(:,4) = [ 1.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]

      first = .false.
    end if

    do i = 1, nn

      v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i))
      v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i))

      dnl = lmean(ql(idn,i), qr(idn,i))
      dna = amean(ql(idn,i), qr(idn,i))
      vxa = amean(ql(ivx,i), qr(ivx,i))
      vya = amean(ql(ivy,i), qr(ivy,i))
      vza = amean(ql(ivz,i), qr(ivz,i))

      v(idn) = 5.0d-01 * (v2l - v2r)
      if (qr(idn,i) > ql(idn,i)) then
        v(idn) = v(idn) + csnd2 * log(qr(idn,i) / ql(idn,i))
      else if (ql(idn,i) > qr(idn,i)) then
        v(idn) = v(idn) - csnd2 * log(ql(idn,i) / qr(idn,i))
      end if
      v(ivx) = qr(ivx,i) - ql(ivx,i)
      v(ivy) = qr(ivy,i) - ql(ivy,i)
      v(ivz) = qr(ivz,i) - ql(ivz,i)

      cs = sign(csnd, vxa)

      lm(:) = [ vxa + cs, vxa, vxa, vxa - cs ]

      tm(1) = 5.0d-01 * dnl / csnd2
      tm(2) = dnl
      tm(3) = dnl
      tm(4) = tm(1)

      rm(2:4,1) = [ lm(1), vya, vza ]
      rm(2:4,4) = [ lm(4), vya, vza ]

      fi(idn,i) = dnl * vxa
      fi(imx,i) = fi(idn,i) * vxa + csnd2 * dna
      fi(imy,i) = fi(idn,i) * vya
      fi(imz,i) = fi(idn,i) * vza

      fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_hd_iso_kepes
!
!*******************************************************************************
!
! ADIABATIC HYDRODYNAMIC RIEMANN SOLVERS
!
!*******************************************************************************
!
!===============================================================================
!
! 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 left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   References:
!
!     [1] Toro, E. F.,
!         "The HLLC Riemann solver",
!         Shock Waves, 2019, 29, 8, 1065-1082,
!         https://doi.org/10.1007/s00193-019-00912-4
!
!===============================================================================
!
  subroutine riemann_hd_hllc(ql, qr, fi)

    use coordinates, only : nn => bcells
    use equations  , only : nf
    use equations  , only : idn, ivy, ivz, imx, imy, imz, ien
    use equations  , only : prim2cons, fluxspeed

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

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

    real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
    real(kind=8), dimension( 2,nn) :: cl, cr
    real(kind=8), dimension(nf   ) :: wl, wr, ui

!-------------------------------------------------------------------------------
!
    call prim2cons(ql, ul)
    call prim2cons(qr, ur)

    call fluxspeed(ql, ul, fl, cl)
    call fluxspeed(qr, ur, fr, cr)

    do i = 1, nn

      sl = min(cl(1,i), cr(1,i))
      sr = max(cl(2,i), cr(2,i))

      if (sl >= 0.0d+00) then

        fi(:,i) = fl(:,i)

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

        fi(:,i) = fr(:,i)

      else ! sl < 0 < sr

        wl = sl * ul(:,i) - fl(:,i)
        wr = sr * ur(:,i) - fr(:,i)

        dn =  wr(idn) - wl(idn)
        sm = (wr(imx) - wl(imx)) / dn

        pr = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn

        if (sm > 0.0d+00) then

          slmm = sl - sm

          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

          fi(:,i) = sl * ui - wl

        else if (sm < 0.0d+00) then

          srmm = sr - sm

          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

          fi(:,i) = sr * ui - wr

        else

          fi(:,i) = (sl * wr - sr * wl) / (sr - sl)

        end if

      end if

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_hd_hllc
!
!===============================================================================
!
! subroutine RIEMANN_HD_ADI_ROE:
! -----------------------------
!
!   Subroutine solves one dimensional Riemann problem using
!   the Roe's method.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   References:
!
!     [1] Roe, P. L.
!         "Approximate Riemann Solvers, Parameter Vectors, and Difference
!          Schemes",
!         Journal of Computational Physics, 1981, 43, pp. 357-372
!         https://doi.org/10.1016/0021-9991(81)90128-5
!     [2] Toro, E. F.,
!         "Riemann Solvers and Numerical Methods for Fluid dynamics",
!         Springer-Verlag Berlin Heidelberg, 2009
!
!===============================================================================
!
  subroutine riemann_hd_adi_roe(ql, qr, fi)

    use coordinates, only : nn => bcells
    use equations  , only : nf, adiabatic_index
    use equations  , only : idn, ivx, ivy, ivz, ipr, ien
    use equations  , only : prim2cons, fluxspeed

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

    integer      :: i
    real(kind=8) :: sdl, sdr, sds, vv, vh, c2, cs, vc, fa, fb, fc, f1, f2, f3

    real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
    real(kind=8), dimension(nf   ) :: qs, lm, al, df

    logical                     , save :: first = .true.
    real(kind=8)                , save :: adi_m1, adi_m1x
    real(kind=8), dimension(5,5), save :: rvec, lvec
!$omp threadprivate(first, adi_m1, adi_m1x, rvec, lvec)

!-------------------------------------------------------------------------------
!
    if (first) then

      adi_m1  = adiabatic_index - 1.0d+00
      adi_m1x = adiabatic_index / adi_m1

      rvec(:,1) = [ 1.0d+00, 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
      rvec(:,2) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
      rvec(:,3) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
      rvec(:,4) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
      rvec(:,5) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]

      lvec(:,1) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
      lvec(:,2) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
      lvec(:,3) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
      lvec(:,4) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
      lvec(:,5) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]

      first = .false.
    end if

    call prim2cons(ql, ul)
    call prim2cons(qr, ur)

    call fluxspeed(ql, ul, fl)
    call fluxspeed(qr, ur, fr)

    do i = 1, nn

      sdl = sqrt(ql(idn,i))
      sdr = sqrt(qr(idn,i))
      sds = sdl + sdr

      qs(idn)     = sdl * sdr
      qs(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds
      qs(ien)     = ((ul(ien,i) + ql(ipr,i)) / sdl &
                  +  (ur(ien,i) + qr(ipr,i)) / sdr) / sds ! enthalpy

      vv = sum(qs(ivx:ivz) * qs(ivx:ivz))
      vh = 5.0d-01 * vv
      c2 = adi_m1 * (qs(ien) - vh)
      cs = sign(sqrt(c2), qs(ivx))
      fa = adi_m1 / c2
      fb = 5.0d-01 * fa
      fc = 5.0d-01 / cs
      vc =   cs * qs(ivx)
      f1 =   fb * vh
      f2 =   fc * qs(ivx)
      f3 = - fb * qs(ivx)

      lm(1) = qs(ivx) + cs
      lm(2) = qs(ivx)
      lm(3) = qs(ivx)
      lm(4) = qs(ivx)
      lm(5) = qs(ivx) - cs

      rvec(1,2) = lm(1)
      rvec(1,3) = qs(ivy)
      rvec(1,4) = qs(ivz)
      rvec(1,5) = qs(ien) + vc
      rvec(2,2) = qs(ivx)
      rvec(2,3) = qs(ivy)
      rvec(2,4) = qs(ivz)
      rvec(2,5) = vh
      rvec(3,5) = qs(ivy)
      rvec(4,5) = qs(ivz)
      rvec(5,2) = lm(5)
      rvec(5,3) = qs(ivy)
      rvec(5,4) = qs(ivz)
      rvec(5,5) = qs(ien) - vc

      lvec(1,1) =   f1 - f2
      lvec(1,2) =   f3 + fc
      lvec(1,3) = - fb * qs(ivy)
      lvec(1,4) = - fb * qs(ivz)
      lvec(1,5) =   fb
      lvec(2,1) = - fa * vh + 1.0d+00
      lvec(2,2) =   fa * qs(ivx)
      lvec(2,3) =   fa * qs(ivy)
      lvec(2,4) =   fa * qs(ivz)
      lvec(2,5) = - fa
      lvec(3,1) = - qs(ivy)
      lvec(4,1) = - qs(ivz)
      lvec(5,1) =   f1 + f2
      lvec(5,2) =   f3 - fc
      lvec(5,3) =   lvec(1,3)
      lvec(5,4) =   lvec(1,4)
      lvec(5,5) =   fb

      al = abs(lm) * matmul(lvec, ur(:,i) - ul(:,i))
      df = matmul(al, rvec)
      fi(:,i) = 5.0d-01 * ((fl(:,i) + fr(:,i)) - df)

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_hd_adi_roe
!
!===============================================================================
!
! subroutine RIEMANN_HD_ADI_KEPES:
! -------------------------------
!
!   Subroutine solves one dimensional adiabatic hydrodynamic Riemann problem
!   using the entropy stable KEPES method.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   References:
!
!     [1] Winters, A. R., Derigs, D., Gassner, G. J., Walch, S.
!         "A uniquely defined entropy stable matrix dissipation operator
!          for high Mach number ideal MHD and compressible Euler
!          simulations",
!         Journal of Computational Physics, 2017, 332, pp. 274-289
!         https://doi.org/10.1016/j.jcp.2016.12.006
!
!===============================================================================
!
  subroutine riemann_hd_adi_kepes(ql, qr, fi)

    use coordinates, only : nn => bcells
    use equations  , only : nf, adiabatic_index
    use equations  , only : idn, ivx, ivy, ivz, ipr, imx, imy, imz, ien

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

    integer      :: i
    real(kind=8) :: bl, br, v2l, v2r, cs, vc
    real(kind=8) :: btl, dnl, prl, bta, dna, vxa, vya, vza, vva, pra, ent

    real(kind=8), dimension(nf) :: v, lm, tm

    logical                     , save :: first = .true.
    real(kind=8)                , save :: adi_m1, adi_m1x
    real(kind=8), dimension(5,5), save :: rm
!$omp threadprivate(first, adi_m1, adi_m1x, rm)

!-------------------------------------------------------------------------------
!
    if (first) then
      adi_m1  = adiabatic_index - 1.0d+00
      adi_m1x = adi_m1 / adiabatic_index

      rm(1,:) = [ 1.0d+00, 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
      rm(2,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
      rm(3,:) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
      rm(4,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
      rm(5,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]

      first = .false.
    end if

    do i = 1, nn

      bl  = ql(idn,i) / ql(ipr,i)
      br  = qr(idn,i) / qr(ipr,i)
      v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i))
      v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i))

      btl = lmean(bl, br)
      bta = amean(bl, br)
      dnl = lmean(ql(idn,i), qr(idn,i))
      dna = amean(ql(idn,i), qr(idn,i))
      vxa = amean(ql(ivx,i), qr(ivx,i))
      vya = amean(ql(ivy,i), qr(ivy,i))
      vza = amean(ql(ivz,i), qr(ivz,i))
      pra = dna / bta
      prl = dnl / btl
      vva = 5.0d-01 * sum(ql(ivx:ivz,i) * qr(ivx:ivz,i))
      ent = 1.0d+00 / adi_m1x / btl + vva

      v(idn) = 5.0d-01 * (bl * v2l - br * v2r)
      if (qr(idn,i) > ql(idn,i)) then
        v(idn) = v(idn) + log(qr(idn,i) / ql(idn,i))
      else if (ql(idn,i) > qr(idn,i)) then
        v(idn) = v(idn) - log(ql(idn,i) / qr(idn,i))
      end if
      if (br > bl) then
        v(idn) = v(idn) + log(br / bl) / adi_m1
      else if (bl > br) then
        v(idn) = v(idn) - log(bl / br) / adi_m1
      end if
      v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i)
      v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i)
      v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i)
      v(ipr) = bl - br

      cs = sign(sqrt(adiabatic_index * pra / dnl), vxa)
      vc = vxa * cs

      lm(1) = vxa + cs
      lm(2) = vxa
      lm(3) = vxa
      lm(4) = vxa
      lm(5) = vxa - cs

      tm(1) = 5.0d-01 * dnl / adiabatic_index
      tm(2) = dnl * adi_m1x
      tm(3) = pra
      tm(4) = pra
      tm(5) = tm(1)

      rm(2,1) = lm(1)
      rm(3,1) = vya
      rm(4,1) = vza
      rm(2,2) = vxa
      rm(3,2) = vya
      rm(4,2) = vza
      rm(2,5) = lm(5)
      rm(3,5) = vya
      rm(4,5) = vza
      rm(5,:) = [ ent + vc, vva, vya, vza, ent - vc ]

      fi(idn,i) = dnl * vxa
      fi(imx,i) = fi(idn,i) * vxa + pra
      fi(imy,i) = fi(idn,i) * vya
      fi(imz,i) = fi(idn,i) * vza
      fi(ien,i) = (prl / adi_m1 + pra + dnl * vva) * vxa

      fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_hd_adi_kepes
!
!*******************************************************************************
!
! RELATIVISTIC ADIABATIC HYDRODYNAMIC RIEMANN SOLVERS
!
!*******************************************************************************
!
!===============================================================================
!
! 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 left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   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, fi)

    use algebra    , only : quadratic
    use coordinates, only : nn => bcells
    use equations  , only : nf, idn, imx, imy, imz, ien
    use equations  , only : prim2cons, fluxspeed

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

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

    real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
    real(kind=8), dimension( 2,nn) :: cl, cr
    real(kind=8), dimension(nf   ) :: wl, wr, uh, us, fh
    real(kind=8), dimension(3)     :: a
    real(kind=8), dimension(2)     :: x

!-------------------------------------------------------------------------------
!
    call prim2cons(ql, ul)
    call prim2cons(qr, ur)

    call fluxspeed(ql, ul, fl, cl)
    call fluxspeed(qr, ur, fr, cr)

    do i = 1, nn

      sl = min(cl(1,i), cr(1,i))
      sr = max(cl(2,i), cr(2,i))

      if (sl >= 0.0d+00) then

        fi(:,i) = fl(:,i)

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

        fi(:,i) = fr(:,i)

      else ! sl < 0 < sr

        srml = sr - sl

        wl(:)  = sl * ul(:,i) - fl(:,i)
        wr(:)  = sr * ur(:,i) - fr(:,i)

        uh(:)  = (     wr(:) -      wl(:)) / srml
        fh(:)  = (sl * wr(:) - sr * wl(:)) / srml

        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
          fi(:,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
            fi(:,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
              fi(:,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])
!
                fi(:,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])
!
                fi(:,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
!
                fi(idn,i) = 0.0d+00
                fi(imx,i) = pr
                fi(imy,i) = 0.0d+00
                fi(imz,i) = 0.0d+00
                fi(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
!
!*******************************************************************************
!
! ISOTHERMAL MAGNETOHYDRODYNAMIC RIEMANN SOLVERS
!
!*******************************************************************************
!
!===============================================================================
!
! subroutine RIEMANN_ISO_MHD_HLLD:
! -------------------------------
!
!   Subroutine solves one dimensional Riemann problem using the modified
!   isothermal HLLD method by Mignone. The difference is in allowing a density
!   jumps across the Alfvén waves, due to the inner intermediate states being
!   an average of the states between two Alfvén waves, which also include
!   two slow magnetosonic waves. Moreover, due to the different left and right
!   intermediate state densities the left and right Alfvén speeds might be
!   different too.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   References:
!
!     [1] Mignone, A.,
!         "A simple and accurate Riemann solver for isothermal MHD",
!         Journal of Computational Physics, 2007, 225, pp. 1427-1441,
!         https://doi.org/10.1016/j.jcp.2007.01.033
!
!===============================================================================
!
  subroutine riemann_mhd_iso_hlld(ql, qr, fi)

    use coordinates, only : nn => bcells
    use equations  , only : nf, idn, imx, imy, imz, ibx, iby, ibz, ibp
    use equations  , only : prim2cons, fluxspeed

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

    integer      :: i
    real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm
    real(kind=8) :: bx, b2, dvl, dvr, dnl, dnr

    real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
    real(kind=8), dimension( 2,nn) :: cl, cr
    real(kind=8), dimension(nf   ) :: wl, wr, wcl, wcr, uil, uir

!-------------------------------------------------------------------------------
!
    call prim2cons(ql, ul)
    call prim2cons(qr, ur)

    call fluxspeed(ql, ul, fl, cl)
    call fluxspeed(qr, ur, fr, cr)

    do i = 1, nn

      sl = min(cl(1,i), cr(1,i))
      sr = max(cl(2,i), cr(2,i))

      if (sl >= 0.0d+00) then

        fi(:,i) = fl(:,i)

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

        fi(:,i) = fr(:,i)

      else

        bx = ql(ibx,i)
        b2 = ql(ibx,i) * qr(ibx,i)

        wl(:) = sl * ul(:,i) - fl(:,i)
        wr(:) = sr * ur(:,i) - fr(:,i)

        sm = (wr(imx) - wl(imx)) / (wr(idn) - wl(idn))

        srml = sr - sl
        slmm = sl - sm
        srmm = sr - sm

        uil(idn) = wl(idn) / slmm
        uir(idn) = wr(idn) / srmm
        uil(imx) = uil(idn) * sm
        uir(imx) = uir(idn) * sm
        uil(ibx) = bx
        uir(ibx) = bx
        uil(ibp) = ql(ibp,i)
        uir(ibp) = qr(ibp,i)

        sml = sm - abs(bx) / sqrt(uil(idn))
        smr = sm + abs(bx) / sqrt(uir(idn))

        dvl = slmm * wl(idn) - b2
        dvr = srmm * wr(idn) - b2

        dnr = uir(idn) / dvr
        dnl = uil(idn) / dvl

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

          if (sml > 0.0d+00) then ! sl* > 0

            uil(imy) = dnl * (slmm * wl(imy) - bx * wl(iby))
            uil(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz))
            uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
            uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl

            fi(:,i)  = sl * uil(:) - wl(:)

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

            uir(imy) = dnr * (srmm * wr(imy) - bx * wr(iby))
            uir(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz))
            uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
            uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr

            fi(:,i)  = sr * uir(:) - wr(:)

          else ! sl* <= 0 <= sr*

            if (smr > sml) then

              uil(imy) = dnl * (slmm * wl(imy) - bx * wl(iby))
              uil(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz))
              uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
              uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl

              uir(imy) = dnr * (srmm * wr(imy) - bx * wr(iby))
              uir(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz))
              uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
              uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr

              wcl(:) = (sml - sl) * uil(:) + wl(:)
              wcr(:) = (smr - sr) * uir(:) + wr(:)

              dvl = smr - sml
              fi(:,i) = (sml * wcr(:) - smr * wcl(:)) / dvl

            else ! Bx = 0 -> Sₘ = 0

              fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml

            end if

          end if ! sl* <= 0 <= sr*

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

          uil(imy) = dnl * (slmm * wl(imy) - bx * wl(iby))
          uil(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz))
          uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
          uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl

          if (sml > 0.0d+00) then ! sl* > 0

            fi(:,i) = sl * uil(:) - wl(:)

          else if (sr > sml) then

            wcl(:) = (sml - sl) * uil(:) + wl(:)

            dvl = sr - sml
            fi(:,i) = (sml * wr(:) - sr * wcl(:)) / dvl

          else ! Bx = 0 -> Sₘ = 0

            fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml

          end if ! sl* < 0

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

          uir(imy) = dnr * (srmm * wr(imy) - bx * wr(iby))
          uir(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz))
          uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
          uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr

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

            fi(:,i) = sr * uir(:) - wr(:)

          else if (smr > sl) then

            wcr(:) = (smr - sr) * uir(:) + wr(:)

            dvr = smr - sl
            fi(:,i) = (sl * wcr(:) - smr * wl(:)) / dvr

          else ! Bx = 0 -> Sₘ = 0

            fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml

          end if ! sr* > 0

        else ! sl* < sl & sr* > sr

          fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml

        end if

      end if

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_mhd_iso_hlld
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ISO_ROE:
! ------------------------------
!
!   Subroutine solves one dimensional Riemann problem using
!   the Roe's method.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   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, fi)

    use coordinates, only : nn => bcells
    use equations  , only : nf, csnd2
    use equations  , only : idn, imx, imy, imz, ivx, ivy, ivz, &
                            ibx, iby, ibz, ibp
    use equations  , only : prim2cons, fluxspeed

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

    integer      :: i
    real(kind=8) :: sdl, sdr, sds, xx, yy, bty, btz, br, br2
    real(kind=8) :: cc, ca, cs, cf, cc2, ca2, cs2, cf2, alf, als
    real(kind=8) :: vqstr, sqrtd, qs, qf, norm, css, cff
    real(kind=8) :: as_prime, af_prime, as, af, aspb, afpb
    real(kind=8) :: f1, f2, f3, f4

    real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
    real(kind=8), dimension(nf   ) :: qi
    real(kind=8), dimension(6    ) :: lm, al, df

    logical                     , save :: first = .true.
    integer     , dimension(6)  , save :: ivs
    real(kind=8), dimension(6,6), save :: rvec, lvec
!$omp threadprivate(first, ivs, rvec, lvec)

!-------------------------------------------------------------------------------
!
    if (first) then

      rvec(:,:) = 0.0d+00
      lvec(:,:) = 0.0d+00

      ivs = [ idn, imx, imy, imz, iby, ibz ]

      first = .false.
    end if

    call prim2cons(ql, ul)
    call prim2cons(qr, ur)

    call fluxspeed(ql, ul, fl)
    call fluxspeed(qr, ur, fr)

    do i = 1, nn

      sdl = sqrt(ql(idn,i))
      sdr = sqrt(qr(idn,i))
      sds = sdl + sdr

      qi(idn)     = sdl * sdr
      qi(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds
      qi(iby:ibz) = (sdr * ql(iby:ibz,i) + sdl * qr(iby:ibz,i)) / sds
      qi(ibx)     = ql(ibx,i)
      qi(ibp)     = ql(ibp,i)

      f1 = ql(iby,i) - qr(iby,i)
      f2 = ql(ibz,i) - qr(ibz,i)
      xx = 5.0d-01 * (f1 * f1 + f2 * f2) / (sds * sds)
      yy = 5.0d-01 * (ql(idn,i) + qr(idn,i)) / qi(idn)

      br2  = qi(iby) * qi(iby) + qi(ibz) * qi(ibz)
      br   = sqrt(br2)
      if (br2 > 0.0d+00) then
        bty  = qi(iby) / br
        btz  = qi(ibz) / br
      else
        bty  = 0.0d+00
        btz  = 0.0d+00
      end if

      cc2 = csnd2 + xx
      ca2 = qi(ibx) * qi(ibx) / qi(idn)
      f1  = br2 / qi(idn)
      f2  = ca2 + f1
      f3  = f2 - cc2
      f4  = sqrt(f3 * f3 + 4.0d+00 * cc2 * f1)
      cf2 = 5.0d-01 * (f2 + cc2 + f4)
      cs2 = cc2 * ca2 / cf2

      if (cs2 < cc2 .and. cc2 < cf2) then
        f1  = (cc2 - cs2) / (cf2 - cs2)
        alf = sqrt(f1)
        als = sqrt(1.0d+00 - f1)
      else if (cc2 >= cf2) then
        alf = 1.0d+00
        als = 0.0d+00
      else
        alf = 0.0d+00
        als = 1.0d+00
      end if

      cc =      sqrt(cc2)
      ca =      sqrt(ca2)
      cf = sign(sqrt(cf2), qi(ivx))
      cs = sign(sqrt(cs2), qi(ivx))

      lm(1) = qi(ivx) + cf
      lm(2) = qi(ivx) + ca
      lm(3) = qi(ivx) + cs
      lm(4) = qi(ivx) - cs
      lm(5) = qi(ivx) - ca
      lm(6) = qi(ivx) - cf

      sqrtd    = sqrt(qi(idn))
      f1       = sqrt(yy)
      qf       = cf * sign(alf, qi(ibx)) / f1
      qs       = cs * sign(als, qi(ibx)) / f1
      f1       = cc / (f1 * sqrtd)
      af_prime = alf * f1
      as_prime = als * f1

      rvec(1,1) =   alf
      rvec(1,2) =   alf * lm(1)
      rvec(1,3) =   alf * qi(ivy) - qs * bty
      rvec(1,4) =   alf * qi(ivz) - qs * btz
      rvec(1,5) =   as_prime * bty
      rvec(1,6) =   as_prime * btz

      rvec(6,1) =   alf
      rvec(6,2) =   alf * lm(6)
      rvec(6,3) =   alf * qi(ivy) + qs * bty
      rvec(6,4) =   alf * qi(ivz) + qs * btz
      rvec(6,5) =   rvec(1,5)
      rvec(6,6) =   rvec(1,6)

      rvec(3,1) =   als
      rvec(3,2) =   als * lm(3)
      rvec(3,3) =   als * qi(ivy) + qf * bty
      rvec(3,4) =   als * qi(ivz) + qf * btz
      rvec(3,5) = - af_prime * bty
      rvec(3,6) = - af_prime * btz

      rvec(4,1) =   als
      rvec(4,2) =   als * lm(4)
      rvec(4,3) =   als * qi(ivy) - qf * bty
      rvec(4,4) =   als * qi(ivz) - qf * btz
      rvec(4,5) =   rvec(3,5)
      rvec(4,6) =   rvec(3,6)

      f1 = sign(1.0d+00 / sqrtd, qi(ibx))
      rvec(2,3) =   btz
      rvec(2,4) = - bty
      rvec(2,5) = - btz * f1
      rvec(2,6) =   bty * f1

      rvec(5,3) = - rvec(2,3)
      rvec(5,4) = - rvec(2,4)
      rvec(5,5) =   rvec(2,5)
      rvec(5,6) =   rvec(2,6)

      norm  = 2.0d+00 * cc2
      cff   = alf * cf / norm
      css   = als * cs / norm
      norm  = norm / yy
      qf    = qf / norm
      qs    = qs / norm
      f1    = qi(idn) / norm
      af    = af_prime * f1
      as    = as_prime * f1
      f1    = br / norm
      afpb  = af_prime * f1
      aspb  = as_prime * f1
      vqstr = (qi(ivy) * bty + qi(ivz) * btz)

      f1 = qs * vqstr
      lvec(1,1) = - cff * lm(6) + f1 - aspb
      lvec(1,2) =   cff
      lvec(1,3) = - qs * bty
      lvec(1,4) = - qs * btz
      lvec(1,5) =   as * bty
      lvec(1,6) =   as * btz

      lvec(6,1) =   cff * lm(1) - f1 - aspb
      lvec(6,2) = - lvec(1,2)
      lvec(6,3) = - lvec(1,3)
      lvec(6,4) = - lvec(1,4)
      lvec(6,5) =   lvec(1,5)
      lvec(6,6) =   lvec(1,6)

      f1 = qf * vqstr
      lvec(3,1) = - css * lm(4) - f1 + afpb
      lvec(3,2) =   css
      lvec(3,3) =   qf * bty
      lvec(3,4) =   qf * btz
      lvec(3,5) = - af * bty
      lvec(3,6) = - af * btz

      lvec(4,1) =   css * lm(3) + f1 + afpb
      lvec(4,2) = - lvec(3,2)
      lvec(4,3) = - lvec(3,3)
      lvec(4,4) = - lvec(3,4)
      lvec(4,5) =   lvec(3,5)
      lvec(4,6) =   lvec(3,6)

      f1 = sign(5.0d-01 * sqrtd, qi(ibx))
      lvec(2,1) = - 5.0d-01 * (qi(ivy) * btz - qi(ivz) * bty)
      lvec(2,3) =   5.0d-01 * btz
      lvec(2,4) = - 5.0d-01 * bty
      lvec(2,5) = - f1 * btz
      lvec(2,6) =   f1 * bty

      lvec(5,1) = - lvec(2,1)
      lvec(5,3) = - lvec(2,3)
      lvec(5,4) = - lvec(2,4)
      lvec(5,5) =   lvec(2,5)
      lvec(5,6) =   lvec(2,6)

      al = abs(lm) * matmul(lvec, ur(ivs,i) - ul(ivs,i))
      df = matmul(al, rvec)
      fi(ivs,i) = 5.0d-01 * ((fl(ivs,i) + fr(ivs,i)) - df(:))
      fi(ibx,i) = fl(ibx,i)
      fi(ibp,i) = fl(ibp,i)

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_mhd_iso_roe
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ISO_KEPES:
! --------------------------------
!
!   Subroutine solves one dimensional isothermal magnetohydrodynamic
!   Riemann problem using the entropy stable KEPES method. The method is
!   a modification of the method described in [1] for the isothermal equation
!   of state.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   References:
!
!     [1] Derigs, D., Winters, A. R., Gassner, G. J., Walch, S., Bohm, M.,
!         "Ideal GLM-MHD: About the entropy consistent nine-wave magnetic
!          field divergence diminishing ideal magnetohydrodynamics equations",
!         Journal of Computational Physics, 2018, 364, pp. 420-467,
!         https://doi.org/10.1016/j.jcp.2018.03.002
!     [2] Winters, A. R., Czernik, C., Schily, M. B., Gassner, G. J.,
!         "Entropy stable numerical approximations for the isothermal and
!          polytropic Euler equations",
!         BIT Numerical Mathematics (2020) 60:791–824,
!         https://doi.org/10.1007/s10543-019-00789-w
!
!===============================================================================
!
  subroutine riemann_mhd_iso_kepes(ql, qr, fi)

    use coordinates, only : nn => bcells
    use equations  , only : nf, csnd, csnd2, ch => cglm
    use equations  , only : idn, imx, imy, imz, ivx, ivy, ivz, &
                            ibx, iby, ibz, ibp

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

    integer      :: i
    real(kind=8) :: dna, pta, vxa, vya, vza, bxa, bya, bza, bpa
    real(kind=8) :: dnl, v2l, v2r, b2l, b2r, bp2, das, daf, csq
    real(kind=8) :: ca, cs, cf, ca2, cs2, cf2, x2, x3
    real(kind=8) :: alf, als, f1, f2, f3, f4

    real(kind=8), dimension(nf) :: v, lm, tm

    logical                     , save :: first = .true.
    real(kind=8), dimension(8,8), save :: rm
!$omp threadprivate(first, rm)

!-------------------------------------------------------------------------------
!
    if (first) then
      rm(:,:) =   0.0d+00
      rm(5,4) =   1.0d+00
      rm(8,4) =   1.0d+00
      rm(5,5) =   1.0d+00
      rm(8,5) = - 1.0d+00

      first = .false.
    end if

    do i = 1, nn

      v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i))
      v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i))
      b2l = sum(ql(ibx:ibz,i) * ql(ibx:ibz,i))
      b2r = sum(qr(ibx:ibz,i) * qr(ibx:ibz,i))

      dnl = lmean(ql(idn,i), qr(idn,i))
      dna = amean(ql(idn,i), qr(idn,i))
      vxa = amean(ql(ivx,i), qr(ivx,i))
      vya = amean(ql(ivy,i), qr(ivy,i))
      vza = amean(ql(ivz,i), qr(ivz,i))
      bxa = amean(ql(ibx,i), qr(ibx,i))
      bya = amean(ql(iby,i), qr(iby,i))
      bza = amean(ql(ibz,i), qr(ibz,i))
      bpa = amean(ql(ibp,i), qr(ibp,i))
      pta = 5.0d-01 * amean(b2l, b2r) + csnd2 * dna

      v(idn) = 5.0d-01 * (v2l - v2r)
      if (qr(idn,i) > ql(idn,i)) then
        v(idn) = v(idn) + csnd2 * log(qr(idn,i) / ql(idn,i))
      else if (ql(idn,i) > qr(idn,i)) then
        v(idn) = v(idn) - csnd2 * log(ql(idn,i) / qr(idn,i))
      end if
      v(ivx) = qr(ivx,i) - ql(ivx,i)
      v(ivy) = qr(ivy,i) - ql(ivy,i)
      v(ivz) = qr(ivz,i) - ql(ivz,i)
      v(ibx) = qr(ibx,i) - ql(ibx,i)
      v(iby) = qr(iby,i) - ql(iby,i)
      v(ibz) = qr(ibz,i) - ql(ibz,i)
      v(ibp) = qr(ibp,i) - ql(ibp,i)

      bp2 = bya * bya + bza * bza
      if (bp2 > 0.0d+00) then
        f1 = sqrt(bp2)
        x2 = bya / f1
        x3 = bza / f1
      else
        x2 = 0.0d+00
        x3 = 0.0d+00
      end if

      ca2 = bxa * bxa / dnl
      f1  = bp2 / dnl
      f2  = ca2 + f1
      f3  = f2 - csnd2
      f4  = sqrt(f3 * f3 + 4.0d+00 * csnd2 * f1)
      cf2 = 5.0d-01 * (f2 + csnd2 + f4)
      cs2 = csnd2 * ca2 / cf2

      if (cs2 < csnd2 .and. csnd2 < cf2) then
        f1  = (csnd2 - cs2) / (cf2 - cs2)
        alf = sqrt(f1)
        als = sqrt(1.0d+00 - f1)
      else if (csnd2 >= cf2) then
        alf = 1.0d+00
        als = 0.0d+00
      else
        alf = 0.0d+00
        als = 1.0d+00
      end if

      cf = sign(sqrt(cf2), vxa)
      ca =      sqrt(ca2)
      cs = sign(sqrt(cs2), vxa)

      lm(1) = vxa + cf
      lm(2) = vxa + ca
      lm(3) = vxa + cs
      lm(4) = vxa + ch
      lm(5) = vxa - ch
      lm(6) = vxa - cs
      lm(7) = vxa - ca
      lm(8) = vxa - cf

      tm(1) = 5.0d-01 / (dnl * csnd2)
      tm(2) = 5.0d-01 / dnl**2
      tm(3) = tm(1)
      tm(4) = 5.0d-01
      tm(5) = tm(4)
      tm(6) = tm(1)
      tm(7) = tm(2)
      tm(8) = tm(1)

      das = dnl * als
      daf = dnl * alf
      csq = csnd * sqrt(dnl)

      f1 = daf
      f2 = sign(das, bxa) * cs
      f3 = als * csq
      rm(1,1) = f1
      rm(2,1) = f1 * (vxa + cf)
      rm(3,1) = f1 * vya - f2 * x2
      rm(4,1) = f1 * vza - f2 * x3
      rm(6,1) = f3 * x2
      rm(7,1) = f3 * x3

      rm(1,8) = f1
      rm(2,8) = f1 * (vxa - cf)
      rm(3,8) = f1 * vya + f2 * x2
      rm(4,8) = f1 * vza + f2 * x3
      rm(6,8) = rm(6,1)
      rm(7,8) = rm(7,1)

      f1 = dnl * sqrt(dna)
      rm(3,2) =    f1 * x3
      rm(4,2) = -  f1 * x2
      rm(6,2) = - dnl * x3
      rm(7,2) =   dnl * x2

      rm(3,7) = - rm(3,2)
      rm(4,7) = - rm(4,2)
      rm(6,7) =   rm(6,2)
      rm(7,7) =   rm(7,2)

      f1 = das
      f2 = sign(daf, bxa) * cf
      f3 = - alf * csq
      rm(1,3) = f1
      rm(2,3) = f1 * (vxa + cs)
      rm(3,3) = f1 * vya + f2 * x2
      rm(4,3) = f1 * vza + f2 * x3
      rm(6,3) = f3 * x2
      rm(7,3) = f3 * x3

      rm(1,6) = f1
      rm(2,6) = f1 * (vxa - cs)
      rm(3,6) = f1 * vya - f2 * x2
      rm(4,6) = f1 * vza - f2 * x3
      rm(6,6) = rm(6,3)
      rm(7,6) = rm(7,3)

      fi(idn,i) = dnl * vxa
      fi(imx,i) = fi(idn,i) * vxa - bxa * bxa + pta
      fi(imy,i) = fi(idn,i) * vya - bxa * bya
      fi(imz,i) = fi(idn,i) * vza - bxa * bza
      fi(ibx,i) = ch  * bpa
      fi(iby,i) = vxa * bya - bxa * vya
      fi(ibz,i) = vxa * bza - bxa * vza
      fi(ibp,i) = ch  * bxa

      fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_mhd_iso_kepes
!
!*******************************************************************************
!
! ADIABATIC MAGNETOHYDRODYNAMIC RIEMANN SOLVERS
!
!*******************************************************************************
!
!===============================================================================
!
! 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 left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   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, fi)

    use coordinates, only : nn => bcells
    use equations  , only : nf, idn, ivy, ivz, imx, imy, imz, ien, &
                                ibx, iby, ibz, ibp
    use equations  , only : prim2cons, fluxspeed

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

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

    real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
    real(kind=8), dimension( 2,nn) :: cl, cr
    real(kind=8), dimension(nf   ) :: wl, wr, ui

!-------------------------------------------------------------------------------
!
    call prim2cons(ql, ul)
    call prim2cons(qr, ur)

    call fluxspeed(ql, ul, fl, cl)
    call fluxspeed(qr, ur, fr, cr)

    do i = 1, nn

      sl = min(cl(1,i), cr(1,i))
      sr = max(cl(2,i), cr(2,i))

      if (sl >= 0.0d+00) then

        fi(:,i) = fl(:,i)

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

        fi(:,i) = fr(:,i)

      else ! sl < 0 < sr

        bx = ql(ibx,i)
        b2 = ql(ibx,i) * qr(ibx,i)

        srml = sr - sl

        wl(:) = sl * ul(:,i) - fl(:,i)
        wr(:) = sr * ur(:,i) - fr(:,i)

        dn =  wr(idn) - wl(idn)
        sm = (wr(imx) - wl(imx)) / dn

        pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2

        if (b2 > 0.0d+00) then

          vy = (wr(imy) - wl(imy)) / dn
          vz = (wr(imz) - wl(imz)) / dn
          by = (wr(iby) - wl(iby)) / srml
          bz = (wr(ibz) - wl(ibz)) / srml
          vb = sm * bx + vy * by + vz * bz

          if (sm > 0.0d+00) then ! sm > 0

            slmm = sl - sm

            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

            fi(:,i)  = sl * ui(:) - wl(:)

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

            srmm = sr - sm

            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

            fi(:,i)  = sr * ui(:) - wr(:)

          else ! sm = 0

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

          end if ! sm = 0

        else ! Bₓ = 0

          if (sm > 0.0d+00) then ! sm > 0

            slmm = sl - sm

            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

            fi(:,i)  = sl * ui(:) - wl(:)

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

            srmm = sr - sm

            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

            fi(:,i)  = sr * ui(:) - wr(:)

          else ! sm = 0

! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
            fi(:,i) = (sl * wr(:) - 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_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 left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   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, fi)

    use coordinates, only : nn => bcells
    use equations  , only : nf, idn, imx, imy, imz, ien, ibx, iby, ibz, ibp
    use equations  , only : prim2cons, fluxspeed

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

    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

    real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
    real(kind=8), dimension( 2,nn) :: cl, cr
    real(kind=8), dimension(nf   ) :: wl, wr, wcl, wcr, ui

!-------------------------------------------------------------------------------
!
    call prim2cons(ql, ul)
    call prim2cons(qr, ur)

    call fluxspeed(ql, ul, fl, cl)
    call fluxspeed(qr, ur, fr, cr)

    do i = 1, nn

      sl = min(cl(1,i), cr(1,i))
      sr = max(cl(2,i), cr(2,i))

      if (sl >= 0.0d+00) then ! sl ≥ 0

        fi(:,i) = fl(:,i)

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

        fi(:,i) = fr(:,i)

      else ! sl < 0 < sr

        bx = ql(ibx,i)
        b2 = ql(ibx,i) * qr(ibx,i)

        srml = sr - sl

        wl(:) = sl * ul(:,i) - fl(:,i)
        wr(:) = sr * ur(:,i) - fr(:,i)

        dn =  wr(idn) - wl(idn)
        sm = (wr(imx) - wl(imx)) / dn

        slmm = sl - sm
        srmm = sr - sm

        dnl = wl(idn) / slmm
        dnr = wr(idn) / srmm

        pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2

        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
!
            fi(:,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
!
            fi(:,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
!
                fi(:,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
!
                fi(:,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;
!
                fi(:,i) = (sml * wcr(:) - 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
!
                fi(:,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
!
                fi(:,i) = sr * ui(:) - wr(:)

              else ! Sm = 0

! since Bx = 0 and Sm = 0, then revert to the HLL flux
!
                fi(:,i) = (sl * wr(:) - 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
!
              fi(:,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
!
                  fi(:,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
!
                  fi(:,i)  = (sm * wr(:) - sr * wcr(:)) / srmm

                end if ! sm < 0

              else ! bx = 0

! no Alfvén wave, so revert to the HLL flux
!
                fi(:,i) = (sl * wr(:) - 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
!
              fi(:,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
!
                  fi(:,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
!
                  fi(:,i) = (sm * wl(:) - sl * wcl(:)) / slmm

                end if ! sm > 0

              else ! bx = 0

! no Alfvén wave, so revert to the HLL flux
!
                fi(:,i) = (sl * wr(:) - 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
!
            fi(:,i) = (sl * wr(:) - 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_MHD_ADI_ROE:
! ------------------------------
!
!   Subroutine solves one dimensional Riemann problem using
!   the Roe's method.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   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, fi)

    use coordinates, only : nn => bcells
    use equations  , only : nf, adiabatic_index
    use equations  , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien, &
                            ibx, iby, ibz, ibp
    use equations  , only : prim2cons, fluxspeed

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

    integer      :: i
    real(kind=8) :: sdl, sdr, sds
    real(kind=8) :: xx, yy, pml, pmr

    real(kind=8) :: cc2, ca2, cf2, cs2, cc, ca, cf, cs
    real(kind=8) :: v2, v2h, br2, br, hp, ayy, sqty
    real(kind=8) :: bty, btz, qf, qs, sqrtd, vqstr, vbet, norm
    real(kind=8) :: alf, als, af_prime, as_prime, afpbb, aspbb, afpb, aspb
    real(kind=8) :: f1, f2, f3, f4, f5

    real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
    real(kind=8), dimension(nf   ) :: qi
    real(kind=8), dimension(7    ) :: lm, al, df

    logical                     , save :: first = .true.
    integer     , dimension(7)  , save :: ivs
    real(kind=8)                , save :: adi_m1, adi_m2, adi_m2d1
    real(kind=8), dimension(7,7), save :: rvec, lvec
!$omp threadprivate(first, ivs, adi_m1, adi_m2, adi_m2d1, rvec, lvec)

!-------------------------------------------------------------------------------
!
    if (first) then
      adi_m1   = adiabatic_index - 1.0d+00
      adi_m2   = adiabatic_index - 2.0d+00
      adi_m2d1 = adi_m2 / adi_m1

      rvec(:,:) = 0.0d+00
      lvec(:,:) = 0.0d+00
      rvec(4,1) = 1.0d+00

      ivs = [ idn, imx, imy, imz, ipr, iby, ibz ]

      first = .false.
    end if

    call prim2cons(ql, ul)
    call prim2cons(qr, ur)

    call fluxspeed(ql, ul, fl)
    call fluxspeed(qr, ur, fr)

    do i = 1, nn

      pml = 5.0d-01 * sum(ql(ibx:ibz,i) * ql(ibx:ibz,i))
      pmr = 5.0d-01 * sum(qr(ibx:ibz,i) * qr(ibx:ibz,i))
      sdl = sqrt(ql(idn,i))
      sdr = sqrt(qr(idn,i))
      sds = sdl + sdr

      qi(idn)     = sdl * sdr
      qi(ivx:ivz) = (sdl * ql(ivx:ivz,i) + 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(iby:ibz) = (sdr * ql(iby:ibz,i) + sdl * qr(iby:ibz,i)) / sds
      qi(ibx)     = ql(ibx,i)
      qi(ibp)     = ql(ibp,i)

      f1 = qr(iby,i) - ql(iby,i)
      f2 = qr(ibz,i) - ql(ibz,i)
      xx = 5.0d-01 * (f1 * f1 + f2 * f2) / (sds * sds)
      yy = 5.0d-01 * (ql(idn,i) + qr(idn,i)) / qi(idn)

      br2 = qi(iby) * qi(iby) + qi(ibz) * qi(ibz)
      if (br2 > 0.0d+00) then
        br  = sqrt(br2)
        bty = qi(iby) / br
        btz = qi(ibz) / br
      else
        br  = 0.0d+00
        bty = 0.0d+00
        btz = 0.0d+00
      end if

      v2   = sum(qi(ivx:ivz) * qi(ivx:ivz))
      v2h  = 5.0d-01 * v2
      ayy  = adi_m1 - adi_m2 * yy
      sqty = sqrt(ayy)
      ca2  = qi(ibx) * qi(ibx)
      hp   = qi(ien) - (ca2 + br2) / qi(idn)
      ca2  = ca2 / qi(idn)
      f1   = adi_m1 * (hp - v2h)
      f2   = adi_m2 * xx
      if (f1 > f2) then
        cc2 = f1 - f2
        cc  = sqrt(cc2)
        f1  = ayy * br2 / qi(idn)
        f2  = ca2 + f1
        f3  = f2 + cc2
        f4  = f2 - cc2
        f5  = sqrt(f4 * f4 + 4.0d+00 * cc2 * f1)
        cf2 = 5.0d-01 * (f3 + f5)
        cs2 = cc2 * ca2 / cf2
      else
        cf2 = ca2 + ayy * br2 / qi(idn)
        cc2 = 0.0d+00
        cc  = 0.0d+00
        cs2 = 0.0d+00
      end if

      if (cs2 < cc2 .and. cc2 < cf2) then
        f1  = (cc2 - cs2) / (cf2 - cs2)
        alf = sqrt(f1)
        als = sqrt(1.0d+00 - f1)
      else if (cc2 >= cf2) then
        alf = 1.0d+00
        als = 0.0d+00
      else
        alf = 0.0d+00
        als = 1.0d+00
      end if

      cf = sign(sqrt(cf2), qi(ivx))
      ca = sign(sqrt(ca2), qi(ivx))
      cs = sign(sqrt(cs2), qi(ivx))

      lm(1) = qi(ivx) + cf
      lm(2) = qi(ivx) + ca
      lm(3) = qi(ivx) + cs
      lm(4) = qi(ivx)
      lm(5) = qi(ivx) - cs
      lm(6) = qi(ivx) - ca
      lm(7) = qi(ivx) - cf

      vbet = (qi(ivy) * bty + qi(ivz) * btz) / sqty

      sqrtd    = sqrt(qi(idn))
      qf       = cf * sign(alf, qi(ibx))
      qs       = cs * sign(als, qi(ibx))
      f1       = cc / sqrtd
      af_prime = f1 * alf
      as_prime = f1 * als
      f1       = br / sqty
      afpbb    = af_prime * f1
      aspbb    = as_prime * f1

      f1 = qs / sqty
      f2 = qs * vbet
      f3 = as_prime / sqty
      rvec(1,1) = alf
      rvec(1,2) = alf * lm(1)
      rvec(1,3) = alf * qi(ivy) - f1 * bty
      rvec(1,4) = alf * qi(ivz) - f1 * btz
      rvec(1,5) = alf * (hp + qi(ivx) * cf) - f2 + aspbb
      rvec(1,6) = f3 * bty
      rvec(1,7) = f3 * btz

      rvec(7,1) = alf
      rvec(7,2) = alf * lm(7)
      rvec(7,3) = alf * qi(ivy) + f1 * bty
      rvec(7,4) = alf * qi(ivz) + f1 * btz
      rvec(7,5) = alf * (hp - qi(ivx) * cf) + f2 + aspbb
      rvec(7,6) = rvec(1,6)
      rvec(7,7) = rvec(1,7)

      f1 = qf / sqty
      f2 = qf * vbet
      f3 = - af_prime / sqty
      rvec(3,1) = als
      rvec(3,2) = als * lm(3)
      rvec(3,3) = als * qi(ivy) + f1 * bty
      rvec(3,4) = als * qi(ivz) + f1 * btz
      rvec(3,5) = als * (hp + qi(ivx) * cs) + f2 - afpbb
      rvec(3,6) = f3 * bty
      rvec(3,7) = f3 * btz

      rvec(5,1) = als
      rvec(5,2) = als * lm(5)
      rvec(5,3) = als * qi(ivy) - f1 * bty
      rvec(5,4) = als * qi(ivz) - f1 * btz
      rvec(5,5) = als * (hp - qi(ivx) * cs) - f2 - afpbb
      rvec(5,6) = rvec(3,6)
      rvec(5,7) = rvec(3,7)

      f1 = sign(1.0d+00, qi(ivx))
      f2 = sign(1.0d+00 / sqrtd, qi(ibx))
      rvec(2,3) =   f1 * btz
      rvec(2,4) = - f1 * bty
      rvec(2,5) =   f1 * (qi(ivy) * btz - qi(ivz) * bty)
      rvec(2,6) = - f2 * btz
      rvec(2,7) =   f2 * bty

      rvec(6,3) = - rvec(2,3)
      rvec(6,4) = - rvec(2,4)
      rvec(6,5) = - rvec(2,5)
      rvec(6,6) =   rvec(2,6)
      rvec(6,7) =   rvec(2,7)

      rvec(4,2) = qi(ivx)
      rvec(4,3) = qi(ivy)
      rvec(4,4) = qi(ivz)
      rvec(4,5) = v2h + adi_m2d1 * xx

      norm  = 2.0d+00 * cc2
      f1    = sqty * br / norm
      afpb  = af_prime * f1
      aspb  = as_prime * f1
      sqty  = sqty / norm
      vqstr = (qi(ivy) * bty + qi(ivz) * btz) * sqty

      f1 = adi_m1 * alf / norm
      f2 = qs * sqty
      f3 = as_prime * qi(idn) * sqty
      f4 = alf * cf / norm
      f5 = qs * vqstr
      lvec(1,1) =   f1 * (v2 - hp) - f4 * lm(7) + f5 - aspb
      lvec(1,5) =   f1
      lvec(1,2) = - f1 * qi(ivx) + f4
      lvec(1,3) = - f1 * qi(ivy) - f2 * bty
      lvec(1,4) = - f1 * qi(ivz) - f2 * btz
      lvec(1,6) =   f3 * bty - f1 * qi(iby)
      lvec(1,7) =   f3 * btz - f1 * qi(ibz)

      lvec(7,1) =   f1 * (v2 - hp) + f4 * lm(1) - f5 - aspb
      lvec(7,5) =   f1
      lvec(7,2) = - f1 * qi(ivx) - f4
      lvec(7,3) = - f1 * qi(ivy) + f2 * bty
      lvec(7,4) = - f1 * qi(ivz) + f2 * btz
      lvec(7,6) =   lvec(1,6)
      lvec(7,7) =   lvec(1,7)

      f1 = adi_m1 * als / norm
      f2 = qf * sqty
      f3 = af_prime * qi(idn) * sqty
      f4 = als * cs / norm
      f5 = qf * vqstr
      lvec(3,1) =   f1 * (v2 - hp) - f4 * lm(5) - f5 + afpb
      lvec(3,5) =   f1
      lvec(3,2) = - f1 * qi(ivx) + f4
      lvec(3,3) = - f1 * qi(ivy) + f2 * bty
      lvec(3,4) = - f1 * qi(ivz) + f2 * btz
      lvec(3,6) = - f3 * bty - f1 * qi(iby)
      lvec(3,6) = - f3 * btz - f1 * qi(ibz)

      lvec(5,1) =   f1 * (v2 - hp) + f4 * lm(3) + f5 + afpb
      lvec(5,5) =   f1
      lvec(5,2) = - f1 * qi(ivx) - f4
      lvec(5,3) = - f1 * qi(ivy) - f2 * bty
      lvec(5,4) = - f1 * qi(ivz) - f2 * btz
      lvec(5,6) =   lvec(3,6)
      lvec(5,7) =   lvec(3,7)

      f1 = sign(5.0d-01, qi(ivx)) * bty
      f2 = sign(5.0d-01, qi(ivx)) * btz
      f3 = sign(1.0d+00, qi(ivx)) * sign(sqrtd, qi(ibx))
      lvec(2,1) =   qi(ivz) * f1 - qi(ivy) * f2
      lvec(2,3) =   f2
      lvec(2,4) = - f1
      lvec(2,6) = - f2 * f3
      lvec(2,7) =   f1 * f3

      lvec(6,1) = - lvec(2,1)
      lvec(6,3) = - lvec(2,3)
      lvec(6,4) = - lvec(2,4)
      lvec(6,6) =   lvec(2,6)
      lvec(6,7) =   lvec(2,7)

      f1 = 1.0d+00 / cc2
      lvec(4,1) =   1.0d+00 - f1 * (v2h - adi_m2d1 * xx)
      lvec(4,5) = - f1
      lvec(4,2) =   f1 * qi(ivx)
      lvec(4,3) =   f1 * qi(ivy)
      lvec(4,4) =   f1 * qi(ivz)
      lvec(4,6) =   f1 * qi(iby)
      lvec(4,7) =   f1 * qi(ibz)

      al = abs(lm) * matmul(lvec, ur(ivs,i) - ul(ivs,i))
      df = matmul(al, rvec)
      fi(ivs,i) = 5.0d-01 * ((fl(ivs,i) + fr(ivs,i)) - df(:))
      fi(ibx,i) = fl(ibx,i)
      fi(ibp,i) = fl(ibp,i)

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_mhd_adi_roe
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ADI_KEPES:
! --------------------------------
!
!   Subroutine solves one dimensional adiabatic magnetohydrodynamic
!   Riemann problem using the entropy stable KEPES method.
!
!   Arguments:
!
!     ql, qr - the primitive variables of the left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   References:
!
!     [1] Derigs, D., Winters, A. R., Gassner, G. J., Walch, S., Bohm, M.,
!         "Ideal GLM-MHD: About the entropy consistent nine-wave magnetic
!          field divergence diminishing ideal magnetohydrodynamics equations",
!         Journal of Computational Physics, 2018, 364, pp. 420-467,
!         https://doi.org/10.1016/j.jcp.2018.03.002
!
!===============================================================================
!
  subroutine riemann_mhd_adi_kepes(ql, qr, fi)

    use coordinates, only : nn => bcells
    use equations  , only : nf, adiabatic_index, cglm
    use equations  , only : idn, imx, imy, imz, ivx, ivy, ivz, &
                            ibx, iby, ibz, ibp, ipr, ien

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

    integer      :: i
    real(kind=8) :: dna, pra, vxa, vya, vza, bxa, bya, bza, bpa, dnl, prl, pta
    real(kind=8) :: v2l, v2r, b2l, b2r, bl, br, bta, btl, bp2, eka, ema, vva
    real(kind=8) :: ch, ca, cs, cf, ca2, cs2, cf2, aa2, x2, x3, ub2, uba
    real(kind=8) :: alf, als, das, daf, csq, wps, wms, wpf, wmf
    real(kind=8) :: f1, f2, f3, f4, f5

    real(kind=8), dimension(nf) :: v, lm, tm

    logical                     , save :: first = .true.
    real(kind=8)                , save :: adi_m1, adi_m1x, adi_m1xi
    real(kind=8), dimension(9,9), save :: rm
!$omp threadprivate(first, adi_m1, adi_m1x, adi_m1xi, rm)

!-------------------------------------------------------------------------------
!
    if (first) then
      adi_m1   = adiabatic_index - 1.0d+00
      adi_m1x  = adi_m1 / adiabatic_index
      adi_m1xi = adiabatic_index / adi_m1

      rm(:,:) =   0.0d+00
      rm(1,5) =   1.0d+00
      rm(6,4) =   1.0d+00
      rm(6,6) =   1.0d+00

      first = .false.
    end if

    do i = 1, nn

      v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i))
      v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i))
      b2l = sum(ql(ibx:ibz,i) * ql(ibx:ibz,i))
      b2r = sum(qr(ibx:ibz,i) * qr(ibx:ibz,i))
      bl  = ql(idn,i) / ql(ipr,i)
      br  = qr(idn,i) / qr(ipr,i)

      btl = lmean(bl, br)
      bta = amean(bl, br)
      dnl = lmean(ql(idn,i), qr(idn,i))
      dna = amean(ql(idn,i), qr(idn,i))
      vxa = amean(ql(ivx,i), qr(ivx,i))
      vya = amean(ql(ivy,i), qr(ivy,i))
      vza = amean(ql(ivz,i), qr(ivz,i))
      bxa = amean(ql(ibx,i), qr(ibx,i))
      bya = amean(ql(iby,i), qr(iby,i))
      bza = amean(ql(ibz,i), qr(ibz,i))
      bpa = amean(ql(ibp,i), qr(ibp,i))
      prl = lmean(ql(ipr,i), qr(ipr,i))
      pra = dna / bta
      eka = 5.0d-01 * amean(v2l, v2r)
      ema = 5.0d-01 * amean(b2l, b2r)
      pta = pra + ema
      vva = 5.0d-01 * sum(ql(ivx:ivz,i) * qr(ivx:ivz,i))
      ub2 = 5.0d-01 * amean(ql(ivx,i) * b2l, qr(ivx,i) * b2r)
      uba = amean(sum(ql(ivx:ivz,i) * ql(ibx:ibz,i)), &
                  sum(qr(ivx:ivz,i) * qr(ibx:ibz,i)))

      v(idn) = 5.0d-01 * (bl * v2l - br * v2r)
      if (qr(idn,i) > ql(idn,i)) then
        v(idn) = v(idn) + log(qr(idn,i) / ql(idn,i))
      else if (ql(idn,i) > qr(idn,i)) then
        v(idn) = v(idn) - log(ql(idn,i) / qr(idn,i))
      end if
      if (br > bl) then
        v(idn) = v(idn) + log(br / bl) / adi_m1
      else if (bl > br) then
        v(idn) = v(idn) - log(bl / br) / adi_m1
      end if
      v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i)
      v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i)
      v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i)
      v(ipr) = bl - br
      v(ibx) = br * qr(ibx,i) - bl * ql(ibx,i)
      v(iby) = br * qr(iby,i) - bl * ql(iby,i)
      v(ibz) = br * qr(ibz,i) - bl * ql(ibz,i)
      v(ibp) = br * qr(ibp,i) - bl * ql(ibp,i)

      bp2 = bya * bya + bza * bza
      if (bp2 > 0.0d+00) then
        f1 = sqrt(bp2)
        x2 = bya / f1
        x3 = bza / f1
      else
        x2 = 0.0d+00
        x3 = 0.0d+00
      end if

      aa2 = adiabatic_index * pra / dnl
      ca2 = bxa * bxa / dnl
      f1  = bp2 / dnl
      f2  = ca2 + f1
      f3  = f2 - aa2
      f4  = sqrt(f3 * f3 + 4.0d+00 * aa2 * f1)
      cf2 = 5.0d-01 * (f2 + aa2 + f4)
      cs2 = aa2 * ca2 / cf2

      if (cs2 < aa2 .and. aa2 < cf2) then
        f1  = (aa2 - cs2) / (cf2 - cs2)
        alf = sqrt(f1)
        als = sqrt(1.0d+00 - f1)
      else if (aa2 >= cf2) then
        alf = 1.0d+00
        als = 0.0d+00
      else
        alf = 0.0d+00
        als = 1.0d+00
      end if

      cf = sign(sqrt(cf2), vxa)
      cs = sign(sqrt(cs2), vxa)
      ca = sign(sqrt(ca2), vxa)
      ch = sign(     cglm, vxa)

      lm(1) = vxa + cf
      lm(2) = vxa + ca
      lm(3) = vxa + cs
      lm(4) = vxa + ch
      lm(5) = vxa
      lm(6) = vxa - ch
      lm(7) = vxa - cs
      lm(8) = vxa - ca
      lm(9) = vxa - cf

      f1 = 5.0d-01 / bta
      f2 = adiabatic_index * dnl
      tm(1) = 5.0d-01 / f2
      tm(2) = f1 / dnl**2
      tm(3) = tm(1)
      tm(4) = f1
      tm(5) = dnl * adi_m1x
      tm(6) = tm(4)
      tm(7) = tm(1)
      tm(8) = tm(2)
      tm(9) = tm(1)

      das = dnl * als
      daf = dnl * alf
      csq = sqrt(f2 / bta)

      f3  = vva + adi_m1xi * prl / dnl
      f4  = vya * x2 + vza * x3
      f5  = sqrt(adiabatic_index * bp2 / (dnl * bta))
      f1  = dnl * (als * f3 - alf * f5)
      f2  = dnl * (als * cs * vxa + sign(alf, bxa) * cf * f4)
      wps = f1 + f2
      wms = f1 - f2
      f1  = dnl * (alf * f3 + als * f5)
      f2  = dnl * (alf * cf * vxa - sign(als, bxa) * cs * f4)
      wpf = f1 + f2
      wmf = f1 - f2

      f1 =   daf
      f2 =   sign(das, bxa) * cs
      f3 =   als * csq
      rm(1,1) = f1
      rm(2,1) = f1 * lm(1)
      rm(3,1) = f1 * vya - f2 * x2
      rm(4,1) = f1 * vza - f2 * x3
      rm(5,1) = wpf
      rm(7,1) = f3 * x2
      rm(8,1) = f3 * x3

      rm(1,9) = f1
      rm(2,9) = f1 * lm(9)
      rm(3,9) = f1 * vya + f2 * x2
      rm(4,9) = f1 * vza + f2 * x3
      rm(5,9) = wmf
      rm(7,9) = rm(7,1)
      rm(8,9) = rm(8,1)

      f1 = sign(dnl * sqrt(dna), vxa)
      rm(3,2) =   f1 * x3
      rm(4,2) = - f1 * x2
      rm(5,2) = - f1 * (x2 * vza - x3 * vya)
      rm(7,2) = - dnl * x3
      rm(8,2) =   dnl * x2

      rm(3,8) = - rm(3,2)
      rm(4,8) = - rm(4,2)
      rm(5,8) = - rm(5,2)
      rm(7,8) =   rm(7,2)
      rm(8,8) =   rm(8,2)

      f1 =   das
      f2 =   sign(daf, bxa) * cf
      f3 = - alf * csq
      rm(1,3) = f1
      rm(2,3) = f1 * lm(3)
      rm(3,3) = f1 * vya + f2 * x2
      rm(4,3) = f1 * vza + f2 * x3
      rm(5,3) = wps
      rm(7,3) = f3 * x2
      rm(8,3) = f3 * x3

      rm(1,7) = f1
      rm(2,7) = f1 * lm(7)
      rm(3,7) = f1 * vya - f2 * x2
      rm(4,7) = f1 * vza - f2 * x3
      rm(5,7) = wms
      rm(7,7) = rm(7,3)
      rm(8,7) = rm(8,3)

      f1 = sign(1.0d+00, vxa)
      rm(5,4) = bxa + bpa * f1
      rm(5,6) = bxa - bpa * f1
      rm(9,4) =   f1
      rm(9,6) = - f1

      rm(2,5) = vxa
      rm(3,5) = vya
      rm(4,5) = vza
      rm(5,5) = vva

      fi(idn,i) = dnl * vxa
      fi(imx,i) = fi(idn,i) * vxa - bxa * bxa + pta
      fi(imy,i) = fi(idn,i) * vya - bxa * bya
      fi(imz,i) = fi(idn,i) * vza - bxa * bza
      fi(ibx,i) = cglm * bpa
      fi(iby,i) = vxa * bya - bxa * vya
      fi(ibz,i) = vxa * bza - bxa * vza
      fi(ibp,i) = cglm * bxa
      fi(ien,i) = fi(idn,i) * (1.0d+00 / (adi_m1 * btl) - eka) &
               + (fi(imx,i) * vxa + fi(imy,i) * vya + fi(imz,i) * vza) &
               + (fi(ibx,i) * bxa + fi(iby,i) * bya + fi(ibz,i) * bza) &
               +  fi(ibp,i) * bpa - ub2 + bxa * uba &
               - cglm * amean(ql(ibx,i) * ql(ibp,i), qr(ibx,i) * qr(ibp,i))

      fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))

    end do

!-------------------------------------------------------------------------------
!
  end subroutine riemann_mhd_adi_kepes
!
!*******************************************************************************
!
! RELATIVISTIC ADIABATIC MAGNETOHYDRODYNAMIC RIEMANN SOLVERS
!
!*******************************************************************************
!
!===============================================================================
!
! 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 left and right Riemann states;
!     fi     - the numerical flux at the cell interface;
!
!   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, fi)

    use algebra    , only : quadratic
    use coordinates, only : nn => bcells
    use equations  , only : nf, idn, ivx, imx, imy, imz, ien, ibx, iby, ibz, ibp
    use equations  , only : prim2cons, fluxspeed

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: ql, qr
    real(kind=8), dimension(:,:), intent(out) :: fi

    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

    real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
    real(kind=8), dimension( 2,nn) :: cl, cr
    real(kind=8), dimension(nf   ) :: wl, wr, uh, us, fh
    real(kind=8), dimension(3)     :: a, vs
    real(kind=8), dimension(2)     :: x

!-------------------------------------------------------------------------------
!
    call prim2cons(ql, ul)
    call prim2cons(qr, ur)

    call fluxspeed(ql, ul, fl, cl)
    call fluxspeed(qr, ur, fr, cr)

    do i = 1, nn

      sl = min(cl(1,i), cr(1,i))
      sr = max(cl(2,i), cr(2,i))

      if (sl >= 0.0d+00) then

        fi(:,i) = fl(:,i)

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

        fi(:,i) = fr(:,i)

      else ! sl < 0 < sr

        srml = sr - sl

        wl(:) = sl * ul(:,i) - fl(:,i)
        wr(:) = sr * ur(:,i) - fr(:,i)

        uh(:) = (     wr(:) -      wl(:)) / srml
        fh(:) = (sl * wr(:) - sr * wl(:)) / srml

        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
            fi(:,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
              fi(:,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
                fi(:,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])
!
                  fi(:,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])
!
                  fi(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))

                else

! intermediate flux is constant across the contact discontinuity so use HLL flux
!
                  fi(:,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
            fi(:,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
              fi(:,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
                fi(:,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])
!
                  fi(:,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])
!
                  fi(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))

                else

! intermediate flux is constant across the contact discontinuity so use HLL flux
!
                  fi(:,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
!
!*******************************************************************************
!
! SUPPORTING SUBROUTINES/FUNCTIONS
!
!*******************************************************************************
!
!===============================================================================
!
! 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
!
!===============================================================================
!
! function LMEAN:
! --------------
!
!   Function calculates the logarithmic mean using the optimized algorithm
!   by Ranocha et al.
!
!   Arguments:
!
!     l, r - the values to calculate the mean of;
!
!   References:
!
!     [1] Ranocha et al.,
!         "Efficient implementation of modern entropy stable and kinetic energy
!          preserving discontinuous Galerkin methods for conservation laws",
!         http://arxiv.org/abs/2112.10517v1
!
!===============================================================================
!
  real(kind=8) function lmean(l, r)

    implicit none

    real(kind=8), intent(in) :: l, r

    real(kind=8) :: u, d, s

    real(kind=8), parameter :: c1 = 2.0d+00, c2 = 2.0d+00 / 3.0d+00, &
                               c3 = 4.0d-01, c4 = 2.0d+00 / 7.0d+00

!-------------------------------------------------------------------------------
!
    d = r - l
    s = r + l
    u = (d / s)**2
    if (u < 1.0d-04) then
      lmean = s / (c1 + u * (c2 + u * (c3 + u * c4)))
    else
      if (d >= 0.0d+00) then
        s = log(r / l)
      else
        s = log(l / r)
      end if
      lmean = abs(d) / s
    end if

!-------------------------------------------------------------------------------
!
  end function lmean
!
!===============================================================================
!
! function AMEAN:
! --------------
!
!   Function calculate the arithmetic mean.
!
!   Arguments:
!
!     l, r - the values to calculate the mean of;
!
!===============================================================================
!
  real(kind=8) function amean(l, r)

    implicit none

    real(kind=8), intent(in) :: l, r

!-------------------------------------------------------------------------------
!
    amean = 5.0d-01 * (l + r)

!-------------------------------------------------------------------------------
!
  end function amean

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