amun-code/sources/equations.F90

!!******************************************************************************
!!
!!  This file is part of the AMUN source code, a program to perform
!!  Newtonian or relativistic magnetohydrodynamical simulations on uniform or
!!  adaptive mesh.
!!
!!  Copyright (C) 2008-2021 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License as published by
!!  the Free Software Foundation, either version 3 of the License, or
!!  (at your option) any later version.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!
!!******************************************************************************
!!
!! module: EQUATIONS
!!
!!  This module provides interface for the systems of equations.  Any set of
!!  equations gives us some basic informations, such as the number of variables,
!!  the primitive and conservative variable definitions, the conversion between
!!  those variables, the flux and characteristic speeds defined in terms of
!!  primitive variables.  All this information is provided by this module.
!!
!!  In order to implement a new set of equations, we need to:
!!
!!  1) define the number of independent variables (or equations) nv;
!!  2) define the variable indices and names (both primitive and conservative);
!!  3) provide subroutines for primitive-conservative variable conversion and
!!     point them to corresponding pointers;
!!  4) provide a subroutine to calculate physical fluxes and characteristic
!!     speeds;
!!  5) provide a subroutine to calculate the maximum speed;
!!  6) optionally, define and read all physical constants related to a given
!!     system;
!!
!!******************************************************************************
!
module equations

  implicit none

! the number of variables, fluxes and passive scalars
!
  integer(kind=4), save :: nv = 0
  integer(kind=4), save :: nf = 0
  integer(kind=4), save :: ns = 0

! interfaces for procedure pointers
!
  interface
    subroutine prim2cons_iface(q, u, s)
      real(kind=8), dimension(:,:), intent(in)  :: q
      real(kind=8), dimension(:,:), intent(out) :: u
      logical     , optional      , intent(in)  :: s
    end subroutine
    subroutine cons2prim_iface(u, q, s)
      real(kind=8), dimension(:,:), intent(in)  :: u
      real(kind=8), dimension(:,:), intent(out) :: q
      integer                     , intent(out) :: s
    end subroutine
    subroutine fluxspeed_iface(q, u, f, c)
      real(kind=8), dimension(:,:)          , intent(in)  :: q, u
      real(kind=8), dimension(:,:)          , intent(out) :: f
      real(kind=8), dimension(:,:), optional, intent(out) :: c
    end subroutine
    function maxspeed_iface(qq) result(maxspeed)
      real(kind=8), dimension(:,:,:,:), intent(in) :: qq
      real(kind=8) :: maxspeed
    end function
    subroutine get_maximum_speeds_iface(qq, um, cm)
      real(kind=8), dimension(:,:,:,:), intent(in) :: qq
      real(kind=8)                    , intent(out):: um, cm
    end subroutine
    subroutine esystem_roe_iface(x, y, q, c, r, l)
      real(kind=8)                , intent(in)    :: x, y
      real(kind=8), dimension(:)  , intent(in)    :: q
      real(kind=8), dimension(:)  , intent(inout) :: c
      real(kind=8), dimension(:,:), intent(inout) :: l, r
    end subroutine
    subroutine nr_iterate_iface(mm, bb, mb, en, dn, w, vv, info)
      real(kind=8), intent(in)    :: mm, bb, mb, en, dn
      real(kind=8), intent(inout) :: w, vv
      logical     , intent(out)   :: info
    end subroutine
  end interface

! pointers to the conversion procedures
!
  procedure(prim2cons_iface)  , pointer, save :: prim2cons       => null()
  procedure(cons2prim_iface)  , pointer, save :: cons2prim       => null()

! pointer to the flux procedure
!
  procedure(fluxspeed_iface)  , pointer, save :: fluxspeed       => null()

! pointer to the maxspeed procedure
!
  procedure(maxspeed_iface)   , pointer, save :: maxspeed        => null()

  procedure(get_maximum_speeds_iface), pointer, save :: &
                                                 get_maximum_speeds => null()

! pointer to the Roe eigensystem procedure
!
  procedure(esystem_roe_iface), pointer, save :: eigensystem_roe => null()

! pointer to the variable conversion method
!
  procedure(nr_iterate_iface) , pointer, save :: nr_iterate => null()

! the system of equations and the equation of state
!
  character(len=32), save :: eqsys = "hydrodynamic"
  character(len=32), save :: eos   = "adiabatic"

! the flag indicating if the set of equations is relativistic or magnetized
!
  logical          , save :: relativistic = .false.
  logical          , save :: magnetized   = .false.

! the variable conversion method
!
  character(len=32), save :: c2p = "1Dw"

! the names of equations and methods
!
  character(len=80), save :: name_eqsys = ""
  character(len=80), save :: name_eos   = ""
  character(len=80), save :: name_c2p   = ""

! direction indices
!
  integer(kind=4)  , save :: inx =  1, iny =  2, inz =  3

! variable indices
!
  integer(kind=4)  , save :: idn = -1
  integer(kind=4)  , save :: ivx = -1, ivy = -1, ivz = -1
  integer(kind=4)  , save :: imx = -1, imy = -1, imz = -1
  integer(kind=4)  , save :: ibx = -1, iby = -1, ibz = -1
  integer(kind=4)  , save :: ibp = -1
  integer(kind=4)  , save :: ipr = -1, ien = -1
  integer(kind=4)  , save :: isl = -1, isu = -1

! variable names
!
  character(len=4), dimension(:), allocatable, save :: pvars, cvars

! variable boundary values
!
  real(kind=8), dimension(:,:,:), allocatable, save :: qpbnd

! adiabatic heat ratio
!
  real(kind=8), save :: adiabatic_index = 5.0d+00 / 3.0d+00

! additional adiabatic parameters
!
  real(kind=8), save :: gammam1 = 2.0d+00 / 3.0d+00, gammam1i = 1.5d+00
  real(kind=8), save :: gammaxi = 2.0d+00 / 5.0d+00

! isothermal speed of sound and its second power
!
  real(kind=8), save :: csnd    = 1.0d+00, csnd2   = 1.0d+00

! maximum speeds in the system
!
  real(kind=8), save :: cmax    = 0.0d+00, cmax2   = 0.0d+00
  real(kind=8), save :: umax    = 0.0d+00, cglm    = 0.0d+00

! the lower limits for density and pressure to be treated as physical
!
  real(kind=8), save :: dmin    = 1.0d-16, pmin    = 1.0d-16

! the upper limits for the Lorentz factor and corresponding |v|²
!
  real(kind=8), save :: lmax    = 1.0d+06
  real(kind=8), save :: vmax    = 0.999999999999d+00

! the upper bound for the sonic Mach number
!
  real(kind=8), save :: msmax   = 1.0d+03
  real(kind=8), save :: msfac   = 3.0d-06 / 5.0d+00

! the tolerance for Newton-Raphson interative method, the maximum number of
! iterations and the number of extra iterations for polishing
!
  real(kind=8), save :: tol     = 1.0d-10
  integer     , save :: nrmax   = 100
  integer     , save :: nrext   = 2

! flag for unphysical cells correction, the maximum distance of neighbors for
! averaging region, and the minimum number of cells for averaging
!
  logical     , save :: fix_unphysical_cells = .false.
  integer     , save :: ngavg   = 2
  integer     , save :: npavg   = 4

! the variable indices for update_flux() from module SCHEMES
!
  integer, dimension(:,:), allocatable :: ivars

! the variable positivity indicator
!
  logical, dimension(:), allocatable :: positive

! the array of maximum errors for conservative variables
!
  real(kind=8), dimension(:), allocatable :: errors

! by default everything is private
!
  private

! declare public variables and subroutines
!
  public :: initialize_equations, finalize_equations, print_equations
  public :: cons2prim, prim2cons, fluxspeed
  public :: prim2cons_hd_iso, fluxspeed_hd_iso
  public :: prim2cons_hd_adi, fluxspeed_hd_adi
  public :: prim2cons_mhd_iso, fluxspeed_mhd_iso
  public :: prim2cons_mhd_adi, fluxspeed_mhd_adi
  public :: prim2cons_srhd_adi, fluxspeed_srhd_adi
  public :: prim2cons_srmhd_adi, fluxspeed_srmhd_adi
  public :: maxspeed, reset_maxspeed, get_maximum_speeds
  public :: eigensystem_roe
  public :: update_primitive_variables
  public :: fix_unphysical_cells, correct_unphysical_states
  public :: adiabatic_index, relativistic, magnetized
  public :: csnd, csnd2
  public :: cmax, cmax2, umax, cglm
  public :: nv, nf, ns
  public :: inx, iny, inz
  public :: idn, ivx, ivy, ivz, imx, imy, imz
  public :: ibx, iby, ibz, ibp, ipr, ien
  public :: isl, isu
  public :: eqsys, eos
  public :: pvars, cvars
  public :: qpbnd
  public :: ivars, positive, errors

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

! include external procedures and variables
!
    use parameters, only : get_parameter

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    character(len=32), intent(in)  :: system, state
    logical          , intent(in)  :: verbose
    integer          , intent(out) :: status

! local variables
!
    integer           :: p
    character(len=32) :: unphysical_fix = "off"
!
!-------------------------------------------------------------------------------
!
    status = 0

! set the system of equations
!
    eqsys = system

! set the equation of state
!
    eos   = state

! get the number of passive scalars
!
    call get_parameter("nscalars"        , ns )
    ns = min(ns, 100)

! get the primitive variable solver
!
    call get_parameter("primitive_solver", c2p)

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

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

! the name of equation system
!
      name_eqsys = "HD"
      eqsys      = "hd"

! initialize the indices of density, and velocity and momenta components
!
      idn = 1
      ivx = 2
      ivy = 3
      ivz = 4
      imx = 2
      imy = 3
      imz = 4

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

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

! the type of equation of state
!
        name_eos = "isothermal"
        eos      = "iso"

! initialize the number of variables:
!
!   - 1 component of density
!   - 3 components of velocity
!
        nv = 4

! allocate the variable indices
!
        allocate(ivars(NDIMS,nv), stat = status)

! fill up the variable indices
!
        ivars(1,:) = (/ idn, ivx, ivy, ivz /)
        ivars(2,:) = (/ idn, ivy, ivz, ivx /)
#if NDIMS == 3
        ivars(3,:) = (/ idn, ivz, ivx, ivy /)
#endif /* NDIMS == 3 */

! set pointers to subroutines
!
        prim2cons          => prim2cons_hd_iso
        cons2prim          => cons2prim_hd_iso
        fluxspeed          => fluxspeed_hd_iso
        maxspeed           => maxspeed_hd_iso
        eigensystem_roe    => esystem_roe_hd_iso
        get_maximum_speeds => get_maximum_speeds_hd_iso

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

! the type of equation of state
!
        name_eos  = "adiabatic"
        eos       = "adi"

! initialize the indices of pressure and total energy
!
        ipr = 5
        ien = 5

! initialize the number of variables:
!
!   - 1 component of density
!   - 3 components of velocity
!   - 1 component of pressure
!
        nv = 5

! allocate the variable indices
!
        allocate(ivars(NDIMS,nv), stat = status)

! fill up the variable indices
!
        ivars(1,:) = (/ idn, ivx, ivy, ivz, ipr /)
        ivars(2,:) = (/ idn, ivy, ivz, ivx, ipr /)
#if NDIMS == 3
        ivars(3,:) = (/ idn, ivz, ivx, ivy, ipr /)
#endif /* NDIMS == 3 */

! set pointers to subroutines
!
        prim2cons          => prim2cons_hd_adi
        cons2prim          => cons2prim_hd_adi
        fluxspeed          => fluxspeed_hd_adi
        maxspeed           => maxspeed_hd_adi
        eigensystem_roe    => esystem_roe_hd_adi
        get_maximum_speeds => get_maximum_speeds_hd_adi

! warn about the unimplemented equation of state
!
      case default

        if (verbose) then
          write(*,*)
          write(*,"(1x,a)") "ERROR!"
          write(*,"(1x,a)") "The selected equation of state is not " //        &
                            "implemented for HD: '" // trim(eos) // "'."
          write(*,"(1x,a)") "Available equations of state: 'iso', 'adi'."
        end if
        status = 1

      end select

! proceed if everything is fine
!
      if (status == 0) then

! initialize the number of fluxes
!
        nf = nv

! set passive scalars index limits and update the number of variables
!
        if (ns > 0) then
          isl = nv + 1
          isu = nv + ns
          nv  = isu
        end if

! allocate arrays for variable names
!
        allocate(pvars(nv), cvars(nv), stat = status)

        if (status == 0) then

! fill in the primitive variable names
!
          pvars(idn) = 'dens'
          pvars(ivx) = 'velx'
          pvars(ivy) = 'vely'
          pvars(ivz) = 'velz'
          if (ipr > 0) pvars(ipr) = 'pres'

! fill in the conservative variable names
!
          cvars(idn) = 'dens'
          cvars(imx) = 'momx'
          cvars(imy) = 'momy'
          cvars(imz) = 'momz'
          if (ien > 0) cvars(ien) = 'ener'

        end if ! status
      end if ! status

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

! the name of equation system
!
      name_eqsys = "MHD"
      eqsys      = "mhd"

! set magnetized flag
!
      magnetized = .true.

! initialize the indices of density, and velocity and momenta components
!
      idn = 1
      ivx = 2
      ivy = 3
      ivz = 4
      imx = 2
      imy = 3
      imz = 4

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

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

! the type of equation of state
!
        name_eos = "isothermal"
        eos      = "iso"

! initialize the indices of magnetic field components and divergence potential
!
        ibx = 5
        iby = 6
        ibz = 7
        ibp = 8

! initialize the number of variables:
!
!   - 1 component of density
!   - 3 components of velocity
!   - 3 components of magnetic field
!   - 1 component of divergence potential
!
        nv  = 8

! allocate the variable indices
!
        allocate(ivars(NDIMS,nv), stat = status)

! fill up the variable indices
!
        ivars(1,:) = (/ idn, ivx, ivy, ivz, ibx, iby, ibz, ibp /)
        ivars(2,:) = (/ idn, ivy, ivz, ivx, iby, ibz, ibx, ibp /)
#if NDIMS == 3
        ivars(3,:) = (/ idn, ivz, ivx, ivy, ibz, ibx, iby, ibp /)
#endif /* NDIMS == 3 */

! set pointers to the subroutines
!
        prim2cons          => prim2cons_mhd_iso
        cons2prim          => cons2prim_mhd_iso
        fluxspeed          => fluxspeed_mhd_iso
        maxspeed           => maxspeed_mhd_iso
        eigensystem_roe    => esystem_roe_mhd_iso
        get_maximum_speeds => get_maximum_speeds_mhd_iso

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

! the type of equation of state
!
        name_eos  = "adiabatic"
        eos       = "adi"

! initialize the indices of pressure, total energy, magnetic field components
! and divergence potential
!
        ipr = 5
        ien = 5
        ibx = 6
        iby = 7
        ibz = 8
        ibp = 9

! initialize the number of variables:
!
!   - 1 component of density
!   - 3 components of velocity
!   - 1 component of pressure
!   - 3 components of magnetic field
!   - 1 component of magnetic divergence potential
!
        nv  = 9

! allocate the variable indices
!
        allocate(ivars(NDIMS,nv), stat = status)

! fill up the variable indices
!
        ivars(1,:) = (/ idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp /)
        ivars(2,:) = (/ idn, ivy, ivz, ivx, ipr, iby, ibz, ibx, ibp /)
#if NDIMS == 3
        ivars(3,:) = (/ idn, ivz, ivx, ivy, ipr, ibz, ibx, iby, ibp /)
#endif /* NDIMS == 3 */

! set pointers to subroutines
!
        prim2cons          => prim2cons_mhd_adi
        cons2prim          => cons2prim_mhd_adi
        fluxspeed          => fluxspeed_mhd_adi
        maxspeed           => maxspeed_mhd_adi
        eigensystem_roe    => esystem_roe_mhd_adi
        get_maximum_speeds => get_maximum_speeds_mhd_adi

      case default

        if (verbose) then
          write(*,*)
          write(*,"(1x,a)") "ERROR!"
          write(*,"(1x,a)") "The selected equation of state is not " //        &
                            "implemented for MHD: '" // trim(eos) // "'."
          write(*,"(1x,a)") "Available equations of state: 'iso', 'adi'."
        end if
        status = 1

      end select

! proceed if everything is fine
!
      if (status == 0) then

! initialize the number of fluxes
!
        nf = nv

! set passive scalars index limits and update the number of variables
!
        if (ns > 0) then
          isl = nv + 1
          isu = nv + ns
          nv  = isu
        end if

! allocate arrays for variable names
!
        allocate(pvars(nv), cvars(nv), stat = status)

        if (status == 0) then

! fill in the primitive variable names
!
          pvars(idn) = 'dens'
          pvars(ivx) = 'velx'
          pvars(ivy) = 'vely'
          pvars(ivz) = 'velz'
          if (ipr > 0) pvars(ipr) = 'pres'
          pvars(ibx) = 'magx'
          pvars(iby) = 'magy'
          pvars(ibz) = 'magz'
          pvars(ibp) = 'bpsi'

! fill in the conservative variable names
!
          cvars(idn) = 'dens'
          cvars(imx) = 'momx'
          cvars(imy) = 'momy'
          cvars(imz) = 'momz'
          if (ien > 0) cvars(ien) = 'ener'
          cvars(ibx) = 'magx'
          cvars(iby) = 'magy'
          cvars(ibz) = 'magz'
          cvars(ibp) = 'bpsi'

        end if ! status
      end if ! status

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

! the name of equation system
!
      name_eqsys = "Special Relativity HD"
      eqsys      = "srhd"

! set relativistic flag
!
      relativistic = .true.

! initialize the indices of density, and velocity and momenta components
!
      idn = 1
      ivx = 2
      ivy = 3
      ivz = 4
      imx = 2
      imy = 3
      imz = 4

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

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

! the type of equation of state
!
        name_eos = "adiabatic"
        eos      = "adi"

! initialize the indices of pressure and total energy
!
        ipr = 5
        ien = 5

! initialize the number of variables:
!
!   - 1 component of density
!   - 3 components of velocity
!   - 1 component of pressure
!
        nv  = 5

! allocate the variable indices
!
        allocate(ivars(NDIMS,nv), stat = status)

! fill up the variable indices
!
        ivars(1,:) = (/ idn, ivx, ivy, ivz, ipr /)
        ivars(2,:) = (/ idn, ivy, ivz, ivx, ipr /)
#if NDIMS == 3
        ivars(3,:) = (/ idn, ivz, ivx, ivy, ipr /)
#endif /* NDIMS == 3 */

! set pointers to subroutines
!
        prim2cons          => prim2cons_srhd_adi
        cons2prim          => cons2prim_srhd_adi
        fluxspeed          => fluxspeed_srhd_adi
        maxspeed           => maxspeed_srhd_adi
        get_maximum_speeds => get_maximum_speeds_srhd_adi

! warn about the unimplemented equation of state
!
      case default

        if (verbose) then
          write(*,*)
          write(*,"(1x,a)") "ERROR!"
          write(*,"(1x,a)") "The selected equation of state is not " //        &
                            "implemented for SR-HD: '" // trim(eos) // "'."
          write(*,"(1x,a)") "Available equations of state: 'adi'."
        end if
        status = 1

      end select

! choose the conserved to primitive variable conversion method
!
      select case(trim(c2p))

      case("1Dw", "1dw", "1DW", "1D(w)", "1D(W)")

! the type of equation of state
!
        name_c2p  = "1D(W)"

! set pointer to the conversion method
!
        nr_iterate => nr_iterate_srhd_adi_1dw

      case("2dwv", "2Dwv", "2D(w,v)", "2D(W,v)")

! the type of equation of state
!
        name_c2p  = "2D(W,v²)"

! set pointer to the conversion method
!
        nr_iterate => nr_iterate_srhd_adi_2dwv

      case("2dwu", "2Dwu", "2D(w,u)", "2D(W,u)")

! the type of equation of state
!
        name_c2p  = "2D(W,u²)"

! set pointer to the conversion method
!
        nr_iterate => nr_iterate_srhd_adi_2dwu

! warn about the unimplemented method
!
      case default

        if (verbose) then
          write(*,*)
          write(*,"(1x,a)") "ERROR!"
          write(*,"(1x,a)") "The selected conversion method is not " //        &
                            "implemented " // trim(c2p) // " for SR-HD."
          write(*,"(1x,a)") "Available conversions: '1Dw', '2Dwv', '2Dwu'."
        end if
        status = 1

      end select

! proceed if everything is fine
!
      if (status == 0) then

! initialize the number of fluxes
!
        nf = nv

! set passive scalars index limits and update the number of variables
!
        if (ns > 0) then
          isl = nv + 1
          isu = nv + ns
          nv  = isu
        end if

! allocate arrays for variable names
!
        allocate(pvars(nv), cvars(nv), stat = status)

        if (status == 0) then
! fill in the primitive variable names
!
          pvars(idn) = 'dens'
          pvars(ivx) = 'velx'
          pvars(ivy) = 'vely'
          pvars(ivz) = 'velz'
          if (ipr > 0) pvars(ipr) = 'pres'

! fill in the conservative variable names
!
          cvars(idn) = 'dens'
          cvars(imx) = 'momx'
          cvars(imy) = 'momy'
          cvars(imz) = 'momz'
          if (ien > 0) cvars(ien) = 'ener'

        end if ! status
      end if ! status

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

! the name of equation system
!
      name_eqsys = "Special Relativity MHD"
      eqsys      = "srmhd"

! set relativistic flag
!
      relativistic = .true.

! set magnetized flag
!
      magnetized = .true.

! initialize the indices of density, and velocity and momenta components
!
      idn = 1
      ivx = 2
      ivy = 3
      ivz = 4
      imx = 2
      imy = 3
      imz = 4

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

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

! the type of equation of state
!
        name_eos = "adiabatic"
        eos      = "adi"

! initialize the indices of pressure, total energy, magnetic field components
! and divergence potential
!
        ipr = 5
        ien = 5
        ibx = 6
        iby = 7
        ibz = 8
        ibp = 9

! initialize the number of variables:
!
!   - 1 component of density
!   - 3 components of velocity
!   - 1 component of pressure
!   - 3 components of magnetic field
!   - 1 component of magnetic divergence potential
!
        nv = 9

! allocate the variable indices
!
        allocate(ivars(NDIMS,nv), stat = status)

! fill up the variable indices
!
        ivars(1,:) = (/ idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp /)
        ivars(2,:) = (/ idn, ivy, ivz, ivx, ipr, iby, ibz, ibx, ibp /)
#if NDIMS == 3
        ivars(3,:) = (/ idn, ivz, ivx, ivy, ipr, ibz, ibx, iby, ibp /)
#endif /* NDIMS == 3 */

! set pointers to subroutines
!
        prim2cons          => prim2cons_srmhd_adi
        cons2prim          => cons2prim_srmhd_adi
        fluxspeed          => fluxspeed_srmhd_adi
        maxspeed           => maxspeed_srmhd_adi
        get_maximum_speeds => get_maximum_speeds_srmhd_adi

! warn about the unimplemented equation of state
!
      case default

        if (verbose) then
          write(*,*)
          write(*,"(1x,a)") "ERROR!"
          write(*,"(1x,a)") "The selected equation of state is not " //        &
                            "implemented for SR-MHD: '" // trim(eos) // "'."
          write(*,"(1x,a)") "Available equations of state: 'adi'."
        end if
        status = 1

      end select

! choose the conserved to primitive variable conversion method
!
      select case(trim(c2p))

      case("1Dw", "1dw", "1DW", "1D(w)", "1D(W)")

! the type of equation of state
!
        name_c2p  = "1D(W)"

! set pointer to the conversion method
!
        nr_iterate => nr_iterate_srmhd_adi_1dw

      case("2dwv", "2Dwv", "2D(w,v)", "2D(W,v)")

! the type of equation of state
!
        name_c2p  = "2D(W,v²)"

! set pointer to the conversion method
!
        nr_iterate => nr_iterate_srmhd_adi_2dwv

      case("2dwu", "2Dwu", "2D(w,u)", "2D(W,u)")

! the type of equation of state
!
        name_c2p  = "2D(W,u²)"

! set pointer to the conversion method
!
        nr_iterate => nr_iterate_srmhd_adi_2dwu

! warn about the unimplemented method
!
      case default

        if (verbose) then
          write(*,*)
          write(*,"(1x,a)") "ERROR!"
          write(*,"(1x,a)") "The selected conversion method is not " //        &
                            "implemented " // trim(c2p) // " for SR-MHD."
          write(*,"(1x,a)") "Available conversions: '1Dw', '2Dwv', '2Dwu'."
        end if
        status = 1

      end select

! proceed if everything is fine
!
      if (status == 0) then

! initialize the number of fluxes
!
        nf = nv

! set passive scalars index limits and update the number of variables
!
        if (ns > 0) then
          isl = nv + 1
          isu = nv + ns
          nv  = isu
        end if

! allocate arrays for variable names
!
        allocate(pvars(nv), cvars(nv), stat = status)

        if (status == 0) then

! fill in the primitive variable names
!
          pvars(idn) = 'dens'
          pvars(ivx) = 'velx'
          pvars(ivy) = 'vely'
          pvars(ivz) = 'velz'
          if (ipr > 0) pvars(ipr) = 'pres'
          pvars(ibx) = 'magx'
          pvars(iby) = 'magy'
          pvars(ibz) = 'magz'
          pvars(ibp) = 'bpsi'

! fill in the conservative variable names
!
          cvars(idn) = 'dens'
          cvars(imx) = 'momx'
          cvars(imy) = 'momy'
          cvars(imz) = 'momz'
          if (ien > 0) cvars(ien) = 'ener'
          cvars(ibx) = 'magx'
          cvars(iby) = 'magy'
          cvars(ibz) = 'magz'
          cvars(ibp) = 'bpsi'

        end if ! status
      end if ! status

!--- EQUATION SYSTEM NOT IMPLEMENTED ---
!
    case default

      if (verbose) then
        write(*,*)
        write(*,"(1x,a)") "ERROR!"
        write(*,"(1x,a)") "The selected equation system is not " //            &
                          "implemented: " // trim(eqsys) // "."
        write(*,"(1x,a)") "Available equation systems: 'hd' 'mhd'" //          &
                          ", 'srhd', 'srmhd'."
      end if
      status = 1

    end select

! allocate positivity indicators
!
    allocate(positive(nv), stat = status)

    if (status == 0) then

      positive( : ) = .false.
      positive(idn) = .true.
      if (ipr > 0) positive(ipr) = .true.

    end if

! allocate an array for errors
!
    allocate(errors(nf), stat = status)

    if (status == 0) errors(:) = 0.0d+00

! proceed if everything is fine
!
    if (status == 0) then

! generate the passive scalar names
!
      if (ns > 0) then
        do p = isl, isu
          write(pvars(p), '("ps",i2.2)') p - isl
          write(cvars(p), '("ps",i2.2)') p - isl
        end do
      end if

! obtain the adiabatic specific heat ratio
!
      call get_parameter("adiabatic_index", adiabatic_index)

! calculate additional parameters
!
      gammam1  = adiabatic_index - 1.0d+00
      gammam1i = 1.0d+00 / gammam1
      gammaxi  = gammam1 / adiabatic_index

! obtain the isothermal sound speed
!
      call get_parameter("sound_speed", csnd )

! calculate additional parameters
!
      csnd2 = csnd * csnd

! allocate array for the boundary values
!
      allocate(qpbnd(nv,3,2), stat = status)

      if (status == 0) then

! set the boundary values
!
        do p = 1, nv

! set the initial boundary values (1.0 for density and pressure, 0.0 otherwise)
!
          if (pvars(p) == "dens" .or. pvars(p) == "pres") then
            qpbnd(p,:,:) = 1.0d+00
          else
            qpbnd(p,:,:) = 0.0d+00
          end if

! read the boundary values from the parameter file
!
          call get_parameter(pvars(p) // "_bnd_xl", qpbnd(p,1,1))
          call get_parameter(pvars(p) // "_bnd_xr", qpbnd(p,1,2))
          call get_parameter(pvars(p) // "_bnd_yl", qpbnd(p,2,1))
          call get_parameter(pvars(p) // "_bnd_yr", qpbnd(p,2,2))
          call get_parameter(pvars(p) // "_bnd_zl", qpbnd(p,3,1))
          call get_parameter(pvars(p) // "_bnd_zr", qpbnd(p,3,2))

        end do ! over all variables

! get the minimum allowed density and pressure in the system, and the maximum
! Lorentz factor for special relativity
!
        call get_parameter("dmin", dmin)
        call get_parameter("pmin", pmin)
        call get_parameter("lmax", lmax)

! calculate the maximum speed corresponding to the maximum Lorentz factor
!
        vmax = 1.0d+00 - 1.0d+00 / lmax**2

! get the upper bound for the sonic Mach number
!
        call get_parameter("msmax", msmax)

! calculate the sonic Mach number factor
!
        msfac = 1.0d+00 / (adiabatic_index * msmax**2)

! get the tolerance
!
        call get_parameter("tolerance"       , tol  )

! get the maximum number of Newton-Raphson method iterations
!
        call get_parameter("nr_maxit"        , nrmax)
        call get_parameter("nr_extra"        , nrext)

! get the state correction flags
!
        call get_parameter("fix_unphysical_cells", unphysical_fix)

! check if the correction of unphysical cells is on
!
        select case(trim(unphysical_fix))
        case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
          fix_unphysical_cells = .true.
        case default
          fix_unphysical_cells = .false.
        end select

! get parameters for unphysical cells correction
!
        call get_parameter("ngavg", ngavg)
        call get_parameter("npavg", npavg)

! correct the above parameters to reasonable values
!
        ngavg = max(1, ngavg)
        npavg = max(2, npavg)

      end if ! status
    end if ! status

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

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    integer, intent(out) :: status
!
!-------------------------------------------------------------------------------
!
    status = 0

! release the procedure pointers
!
    nullify(prim2cons)
    nullify(cons2prim)
    nullify(fluxspeed)
    nullify(maxspeed)
    nullify(eigensystem_roe)

! deallocate variable indices
!
    if (allocated(ivars)) deallocate(ivars, stat = status)

! deallocate variable name arrays
!
    if (allocated(pvars)) deallocate(pvars, stat = status)
    if (allocated(cvars)) deallocate(cvars, stat = status)

! deallocate boundary values array
!
    if (allocated(qpbnd)) deallocate(qpbnd, stat = status)

! deallocate positivity indicator
!
    if (allocated(positive)) deallocate(positive, stat = status)

! deallocate the array for errors
!
    if (allocated(errors)) deallocate(errors, stat = status)

!-------------------------------------------------------------------------------
!
  end subroutine finalize_equations
!
!===============================================================================
!
! subroutine PRINT_EQUATIONS:
! --------------------------
!
!   Subroutine prints module parameters.
!
!   Arguments:
!
!     verbose - a logical flag turning the information printing;
!
!===============================================================================
!
  subroutine print_equations(verbose)

! include external procedures and variables
!
    use helpers, only : print_section, print_parameter

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    logical, intent(in) :: verbose

! local variables
!
    character(len= 80) :: sfmt
    character(len=255) :: msg
!
!-------------------------------------------------------------------------------
!
    if (verbose) then

      call print_section(verbose, "Physics")
      call print_parameter(verbose, "equation system"          , name_eqsys)
      call print_parameter(verbose, "equation of state"        , name_eos  )
      if (ipr > 0) then
        call print_parameter(verbose, "adiabatic index", adiabatic_index)
      else
        call print_parameter(verbose, "sound speed"    , csnd)
      end if
      call print_parameter(verbose, "number of variables"      , nv        )
      call print_parameter(verbose, "number of fluxes"         , nf        )
      call print_parameter(verbose, "number of passive scalars", ns        )
      write(sfmt,"(a,i0,a)") "(", nv, "(1x,a))"
      write(msg,sfmt) cvars
      call print_parameter(verbose, "conservative variables", msg       )
      write(msg,sfmt) pvars
      call print_parameter(verbose, "primitive variables"   , msg       )
      if (relativistic) then
        call print_parameter(verbose, "variable conversion"   , name_c2p  )
      end if
      if (fix_unphysical_cells) then
        call print_parameter(verbose, "fix unphysical cells"  , "on"      )
        call print_parameter(verbose, "ngavg"                 , ngavg     )
        call print_parameter(verbose, "npavg"                 , npavg     )
      else
        call print_parameter(verbose, "fix unphysical cells"  , "off"     )
      end if

    end if

!-------------------------------------------------------------------------------
!
  end subroutine print_equations
!
!===============================================================================
!
! subroutine RESET_MAXSPEED:
! -------------------------
!
!   Subroutine resets the maximum speed in the domain to zero.
!
!
!===============================================================================
!
  subroutine reset_maxspeed()

! local variables are not implicit by default
!
    implicit none
!
!-------------------------------------------------------------------------------
!
! reset the maximum speed
!
    cmax = 0.0d+00

!-------------------------------------------------------------------------------
!
  end subroutine reset_maxspeed
!
!===============================================================================
!
! function GET_MAXSPEED:
! -----------------
!
!   Function returns the maximum speed in the domain.
!
!
!===============================================================================
!
  real(kind=8) function get_maxspeed()

! local variables are not implicit by default
!
    implicit none
!
!-------------------------------------------------------------------------------
!
! return the maximum speed
!
    get_maxspeed = cmax

! return the value
!
    return

!-------------------------------------------------------------------------------
!
  end function get_maxspeed
!
!===============================================================================
!
! subroutine UPDATE_PRIMITIVE_VARIABLES:
! -------------------------------------
!
!   Subroutine updates primitive variables from their conservative
!   representation.  This process is done once after advance of the conserved
!   variables due to their evolution in time.
!
!   Arguments:
!
!     uu - the input array of conservative variables;
!     qq - the output array of primitive variables;
!
!===============================================================================
!
  subroutine update_primitive_variables(uu, qq, status)

    use coordinates, only : nn => bcells
    use coordinates, only : nb, ne, nbl, neu

    implicit none

    real(kind=8), dimension(:,:,:,:), intent(inout) :: uu
    real(kind=8), dimension(:,:,:,:), intent(inout) :: qq
    integer                         , intent(out)   :: status

    integer :: i, j, k

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

#if NDIMS == 2
    k = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
    do k = nb, ne
#endif /* NDIMS == 3 */
      do j = nb, ne

        call cons2prim(uu(1:nv,nb:ne,j,k), qq(1:nv,nb:ne,j,k), status)

        if (status /= 0) go to 100

      end do ! j = nb, ne
#if NDIMS == 3
    end do ! k = nb, ne
#endif /* NDIMS == 3 */

#if NDIMS == 3
    do i = nbl, 1, -1
      qq(1:nv, i,nb:ne,nb:ne) = qq(1:nv,nb,nb:ne,nb:ne)
    end do
    do i = neu, nn
      qq(1:nv, i,nb:ne,nb:ne) = qq(1:nv,ne,nb:ne,nb:ne)
    end do
    do j = nbl, 1, -1
      qq(1:nv,  :  , j,nb:ne) = qq(1:nv,  :  ,nb,nb:ne)
    end do
    do j = neu, nn
      qq(1:nv,  :  , j,nb:ne) = qq(1:nv,  :  ,ne,nb:ne)
    end do
    do k = nbl, 1, -1
      qq(1:nv,  :  ,  :  , k) = qq(1:nv,  :  ,  :  ,nb)
    end do
    do k = neu, nn
      qq(1:nv,  :  ,  :  , k) = qq(1:nv,  :  ,  :  ,ne)
    end do
#else /* NDIMS == 3 */
    do i = nbl, 1, -1
      qq(1:nv, i,nb:ne,  :  ) = qq(1:nv,nb,nb:ne,  :  )
    end do
    do i = neu, nn
      qq(1:nv, i,nb:ne,  :  ) = qq(1:nv,ne,nb:ne,  :  )
    end do
    do j = nbl, 1, -1
      qq(1:nv,  :  , j,  :  ) = qq(1:nv,  :  ,nb,  :  )
    end do
    do j = neu, nn
      qq(1:nv,  :  , j,  :  ) = qq(1:nv,  :  ,ne,  :  )
    end do
#endif /* NDIMS == 3 */

    100 continue

!-------------------------------------------------------------------------------
!
  end subroutine update_primitive_variables
!
!===============================================================================
!
! subroutine CORRECT_UNPHYSICAL_STATES:
! ------------------------------------
!
!   Subroutine seeks for unphysical states (cells with negative density or
!   pressure) and try to fix them by averaging their values from physical
!   neighbours.
!
!   Arguments:
!
!     it - the time step number;
!     id - the block id where the states are being checked;
!     qq - the output array of primitive variables;
!     uu - the input array of conservative variables;
!
!===============================================================================
!
  subroutine correct_unphysical_states(it, id, qq, uu)

    use coordinates, only : nn => bcells
    use helpers    , only : print_message

    implicit none

    integer(kind=4)                 , intent(in)    :: it, id
    real(kind=8), dimension(:,:,:,:), intent(inout) :: qq
    real(kind=8), dimension(:,:,:,:), intent(inout) :: uu

    character(len=255) :: msg, sfmt
    character(len=16)  :: sit, sid, snc
    integer            :: n, p, nc, np
    integer            :: i, il, iu
    integer            :: j, jl, ju
    integer            :: k
#if NDIMS == 3
    integer            :: kl, ku
#endif /* NDIMS == 3 */

#if NDIMS == 3
    logical, dimension(nn,nn,nn) :: physical
#else /* NDIMS == 3 */
    logical, dimension(nn,nn, 1) :: physical
#endif /* NDIMS == 3 */

    integer     , dimension(:,:), allocatable :: idx
    real(kind=8), dimension(:,:), allocatable :: q, u

    character(len=*), parameter :: loc = 'EQUATIONS::correct_unphysical_states()'

!-------------------------------------------------------------------------------
!
    np = 0
    i  = 1
    il = 1
    iu = 1
    j  = 1
    jl = 1
    ju = 1
#if NDIMS == 2
    k  = 1
    kl = 1
    ku = 1
#endif /* NDIMS == 2 */

! search for negative density or pressure
!
    physical(:,:,:) = qq(idn,:,:,:) > 0.0d+00
    if (ipr > 0) then
      physical(:,:,:) = physical(:,:,:) .and. qq(ipr,:,:,:) > 0.0d+00
    end if

! apply averaging for unphysical states
!
    if (.not. all(physical)) then

! count unphysical cells
!
      nc = count(.not. physical)

! inform about the encountered unphysical states
!
      write(sit,'(i15)') it
      write(sid,'(i15)') id
      write(snc,'(i15)') nc
      sfmt = '("Unphysical cells in block ID:",a," (",a," counted)'            &
             // ' at time step ",a,".")'
      write(msg,sfmt) trim(adjustl(sid)), trim(adjustl(snc)), trim(adjustl(sit))
      call print_message(loc, msg)

! allocate temporary vectors for unphysical states
!
      allocate(q(nv,nc), u(nv,nc), idx(3,nc))

! iterate over block cells
!
      n = 0
#if NDIMS == 3
      do k = 1, nn
#endif /* NDIMS == 3 */
        do j = 1, nn
          do i = 1, nn

! fix unphysical states
!
            if (.not. physical(i,j,k)) then

              n        = n + 1

              idx(:,n) = (/ i, j, k /)

! increase the region until we find at least three physical cells, but no more
! than 4 cells away
!
              np = 0
              p  = 1
              do while (np <= npavg .and. p <= ngavg)
                il = max( 1, i - p)
                iu = min(nn, i + p)
                jl = max( 1, j - p)
                ju = min(nn, j + p)
#if NDIMS == 3
                kl = max( 1, k - p)
                ku = min(nn, k + p)

                np = count(physical(il:iu,jl:ju,kl:ku))
#else /* NDIMS == 3 */

                np = count(physical(il:iu,jl:ju, 1   ))
#endif /* NDIMS == 3 */

                p  = p + 1
              end do

! average primitive variables
!
              if (np >= npavg) then

                do p = 1, nv
#if NDIMS == 3
                  q(p,n) = sum(qq(p,il:iu,jl:ju,kl:ku),                        &
                                            physical(il:iu,jl:ju,kl:ku))       &
                                    / np
#else /* NDIMS == 3 */
                  q(p,n) = sum(qq(p,il:iu,jl:ju, 1   ),                        &
                                            physical(il:iu,jl:ju, 1   ))       &
                                    / np
#endif /* NDIMS == 3 */
                end do

              else

! limit density or pressure to minimum value, since the averaging over
! neighbours failed
!
                msg = "Not sufficient number of physical neighbors!"
                call print_message(loc, msg)
                sfmt = '("Block ID:",a,", cell position = ( ",3(i4," ")," ).")'
                write(msg,sfmt) trim(adjustl(sid)), i, j, k
                call print_message(loc, msg)
                write(msg,"('Q = ',10(1x,1es24.16e3))") qq(:,i,j,k)
                call print_message(loc, msg)
                msg = "Applying lower bounds for positive variables."
                call print_message(loc, msg)

                q(:,n) = qq(:,i,j,k)
                q(idn ,n) = max(dmin, qq(idn,i,j,k))
                if (ipr > 0) q(ipr,n) = max(pmin, qq(ipr,i,j,k))

              end if ! not sufficient number of physical cells for averaging

            end if ! not physical

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

! convert the vector of primitive variables to conservative ones
!
      call prim2cons(q(:,1:nc), u(:,1:nc), .true.)

! update block variables
!
      do n = 1, nc
        i = idx(1,n)
        j = idx(2,n)
        k = idx(3,n)

        qq(:,i,j,k) = q(:,n)
        uu(:,i,j,k) = u(:,n)
      end do

! deallocate temporary vectors
!
      deallocate(q, u, idx)

    end if ! there are unphysical cells

!-------------------------------------------------------------------------------
!
  end subroutine correct_unphysical_states
!
!===============================================================================
!!
!!***  PRIVATE SUBROUTINES  ****************************************************
!!
!===============================================================================
!
!*******************************************************************************
!
! ISOTHERMAL HYDRODYNAMIC EQUATIONS
!
!*******************************************************************************
!
!===============================================================================
!
! subroutine PRIM2CONS_HD_ISO:
! ---------------------------
!
!   Subroutine converts primitive variables to their corresponding
!   conservative representation.
!
!   Arguments:
!
!     q - the input array of primitive variables;
!     u - the output array of conservative variables;
!     s - an optional flag indicating that passive scalars have
!         to be calculated too;
!
!===============================================================================
!
  subroutine prim2cons_hd_iso(q, u, s)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: q
    real(kind=8), dimension(:,:), intent(out) :: u
    logical     , optional      , intent(in)  :: s

! local variables
!
    integer :: i, p
!
!-------------------------------------------------------------------------------
!
! iterate over all positions
!
    do i = 1, size(q,2)

      u(idn,i) = q(idn,i)
      u(imx,i) = q(idn,i) * q(ivx,i)
      u(imy,i) = q(idn,i) * q(ivy,i)
      u(imz,i) = q(idn,i) * q(ivz,i)

    end do

! update primitive passive scalars
!
    if (ns > 0 .and. present(s)) then
      if (s) then
        do p = isl, isu
          u(p,:) = q(p,:) * u(idn,:)
        end do
      end if
    end if

!-------------------------------------------------------------------------------
!
  end subroutine prim2cons_hd_iso
!
!===============================================================================
!
! subroutine CONS2PRIM_HD_ISO:
! ---------------------------
!
!   Subroutine converts conservative variables to their corresponding
!   primitive representation.
!
!   Arguments:
!
!     u - the input array of conservative variables;
!     q - the output array of primitive variables;
!     s - the status flag;
!
!===============================================================================
!
  subroutine cons2prim_hd_iso(u, q, s)

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: u
    real(kind=8), dimension(:,:), intent(out) :: q
    integer                     , intent(out) :: s

    integer :: i, p

!-------------------------------------------------------------------------------
!
    s = 0

    do i = 1, size(u,2)

      if (u(idn,i) > 0.0d+00) then
        q(idn,i) = u(idn,i)
        q(ivx,i) = u(imx,i) / u(idn,i)
        q(ivy,i) = u(imy,i) / u(idn,i)
        q(ivz,i) = u(imz,i) / u(idn,i)
      else
        s = 1
        go to 100
      end if

    end do

    if (ns > 0) then
      do p = isl, isu
        q(p,:) = u(p,:) / u(idn,:)
      end do
    end if

    100 continue

!-------------------------------------------------------------------------------
!
  end subroutine cons2prim_hd_iso
!
!===============================================================================
!
! subroutine FLUXSPEED_HD_ISO:
! ---------------------------
!
!   Subroutine calculates physical fluxes and characteristic speeds from a
!   given equation system.
!
!   Arguments:
!
!     q - the input array of primitive variables;
!     u - the input array of conservative variables;
!     f - the output vector of fluxes;
!     c - the output vector of left- and right-going characteristic speeds;
!
!===============================================================================
!
  subroutine fluxspeed_hd_iso(q, u, f, c)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8), dimension(:,:)          , intent(in)  :: q, u
    real(kind=8), dimension(:,:)          , intent(out) :: f
    real(kind=8), dimension(:,:), optional, intent(out) :: c

! local variables
!
    integer :: i
!
!-------------------------------------------------------------------------------
!
! calculate the hydrodynamic fluxes
!
    do i = 1, size(q,2)

      f(idn,i) = u(imx,i)
      f(imx,i) = q(ivx,i) * u(imx,i)
      f(imy,i) = q(ivx,i) * u(imy,i)
      f(imz,i) = q(ivx,i) * u(imz,i)
      f(imx,i) = f(imx,i) + csnd2 * q(idn,i)

    end do

! calculate the characteristic speeds
!
    if (present(c)) then

      do i = 1, size(q,2)

        c(1,i) = q(ivx,i) - csnd
        c(2,i) = q(ivx,i) + csnd

      end do

    end if

!-------------------------------------------------------------------------------
!
  end subroutine fluxspeed_hd_iso
!
!===============================================================================
!
! function MAXSPEED_HD_ISO:
! ------------------------
!
!   Function scans the variable array and returns the maximum speed in within.
!
!   Arguments:
!
!     q - the array of primitive variables;
!
!
!===============================================================================
!
  function maxspeed_hd_iso(qq) result(maxspeed)

    use coordinates, only : nb, ne

    implicit none

    real(kind=8), dimension(:,:,:,:), intent(in) :: qq

    real(kind=8) :: maxspeed

    integer      :: i, j, k
    real(kind=8) :: vv, v

!-------------------------------------------------------------------------------
!
    maxspeed = 0.0d+00

#if NDIMS == 2
    k = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
    do k = nb, ne
#endif /* NDIMS == 3 */
      do j = nb, ne
        do i = nb, ne

! calculate the velocity amplitude
!
          vv = sum(qq(ivx:ivz,i,j,k) * qq(ivx:ivz,i,j,k))
          v  = sqrt(vv)

! calculate the maximum speed
!
          maxspeed = max(maxspeed, v + csnd)

        end do ! i = nb, ne
      end do ! j = nb, ne
#if NDIMS == 3
    end do ! k = nb, ne
#endif /* NDIMS == 3 */

    return

!-------------------------------------------------------------------------------
!
  end function maxspeed_hd_iso
!
!===============================================================================
!
! subroutine GET_MAXIMUM_SPEEDS_HD_ISO:
! ------------------------------------
!
!   Subroutine determines the maximum characteristic speed and eigenvalue
!   in the input array of primitive variables.
!
!   Arguments:
!
!     qq - the input array of primitive variables;
!     vm - the maximum physical speed;
!     cm - the maximum eigenvalue;
!
!===============================================================================
!
  subroutine get_maximum_speeds_hd_iso(qq, vm, cm)

    use coordinates, only : nb, ne

    implicit none

    real(kind=8), dimension(:,:,:,:), intent(in)  :: qq
    real(kind=8)                    , intent(out) :: vm, cm

    integer      :: i, j, k
    real(kind=8) :: vl, vu

!-------------------------------------------------------------------------------
!
    vm = 0.0d+00
    cm = 0.0d+00

#if NDIMS == 2
    k = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
    do k = nb, ne
#endif /* NDIMS == 3 */
      do j = nb, ne
        do i = nb, ne
          vl = minval(qq(ivx:ivz,i,j,k))
          vu = maxval(qq(ivx:ivz,i,j,k))
          vm = max(vm, abs(vl), abs(vu))
          cm = max(cm, abs(vl - csnd), abs(vu + csnd))
        end do
      end do
#if NDIMS == 3
    end do
#endif /* NDIMS == 3 */

!-------------------------------------------------------------------------------
!
  end subroutine get_maximum_speeds_hd_iso
!
!===============================================================================
!
! subroutine ESYSTEM_ROE_HD_ISO:
! -----------------------------
!
!   Subroutine computes eigenvalues and eigenvectors for a given set of
!   equations and input variables.
!
!   Arguments:
!
!     x - ratio of the perpendicular magnetic field component difference
!     y - ratio of the density
!     q - the intermediate Roe state vector;
!     c - the vector of eigenvalues;
!     r - the matrix of right eigenvectors;
!     l - the matrix of left eigenvectors;
!
!   References:
!
!     [1] Roe, P. L.
!         "Approximate Riemann Solvers, Parameter Vectors, and Difference
!          Schemes",
!         Journal of Computational Physics, 1981, 43, pp. 357-372
!     [2] Stone, J. M. & Gardiner, T. A.,
!         "ATHENA: A New Code for Astrophysical MHD",
!         The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177
!
!===============================================================================
!
  subroutine esystem_roe_hd_iso(x, y, q, c, r, l)

    implicit none

    real(kind=8)                , intent(in)    :: x, y
    real(kind=8), dimension(:)  , intent(in)    :: q
    real(kind=8), dimension(:)  , intent(inout) :: c
    real(kind=8), dimension(:,:), intent(inout) :: l, r

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

    real(kind=8) :: vc

!-------------------------------------------------------------------------------
!
    if (first) then
      lvec(:,:) = 0.0d+00
      rvec(:,:) = 0.0d+00

      lvec(ivx,1) = - 5.0d-01 / csnd
      lvec(ivy,2) =   1.0d+00
      lvec(ivz,3) =   1.0d+00
      lvec(ivx,4) = - lvec(ivx,1)

      rvec(1,idn) = 1.0d+00
      rvec(2,ivy) = 1.0d+00
      rvec(3,ivz) = 1.0d+00
      rvec(4,idn) = 1.0d+00

      first = .false.
    end if

! eigenvalues
!
    c(1) = q(ivx) - csnd
    c(2) = q(ivx)
    c(3) = q(ivx)
    c(4) = q(ivx) + csnd

! update the varying elements of the left eigenvectors matrix
!
    vc = 5.0d-01 * q(ivx) / csnd
    lvec(idn,1) = 5.0d-01 + vc
    lvec(idn,2) = - q(ivy)
    lvec(idn,3) = - q(ivz)
    lvec(idn,4) = 5.0d-01 - vc

! update the varying elements of the right eigenvectors matrix
!
    rvec(1,ivx) = c(1)
    rvec(1,ivy) = q(ivy)
    rvec(1,ivz) = q(ivz)

    rvec(4,ivx) = c(4)
    rvec(4,ivy) = q(ivy)
    rvec(4,ivz) = q(ivz)

! copy matrices of eigenvectors
!
    l(:,:) = lvec(:,:)
    r(:,:) = rvec(:,:)

!-------------------------------------------------------------------------------
!
  end subroutine esystem_roe_hd_iso
!
!*******************************************************************************
!
! ADIABATIC HYDRODYNAMIC EQUATIONS
!
!*******************************************************************************
!
!===============================================================================
!
! subroutine PRIM2CONS_HD_ADI:
! ---------------------------
!
!   Subroutine converts primitive variables to their corresponding
!   conservative representation.
!
!   Arguments:
!
!     q - the input array of primitive variables;
!     u - the output array of conservative variables;
!     s - an optional flag indicating that passive scalars have
!         to be calculated too;
!
!===============================================================================
!
  subroutine prim2cons_hd_adi(q, u, s)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: q
    real(kind=8), dimension(:,:), intent(out) :: u
    logical     , optional      , intent(in)  :: s

! local variables
!
    integer      :: i, p
    real(kind=8) :: ek, ei
!
!-------------------------------------------------------------------------------
!
! iterate over all positions
!
    do i = 1, size(q,2)

      u(idn,i) = q(idn,i)
      u(imx,i) = q(idn,i) * q(ivx,i)
      u(imy,i) = q(idn,i) * q(ivy,i)
      u(imz,i) = q(idn,i) * q(ivz,i)
      ek       = 0.5d+00 * (u(imx,i) * q(ivx,i) + u(imy,i) * q(ivy,i)          &
                                                + u(imz,i) * q(ivz,i))
      ei       = gammam1i * q(ipr,i)
      u(ien,i) = ei + ek

    end do

! update primitive passive scalars
!
    if (ns > 0 .and. present(s)) then
      if (s) then
        do p = isl, isu
          u(p,:) = q(p,:) * u(idn,:)
        end do
      end if
    end if

!-------------------------------------------------------------------------------
!
  end subroutine prim2cons_hd_adi
!
!===============================================================================
!
! subroutine CONS2PRIM_HD_ADI:
! ---------------------------
!
!   Subroutine converts conservative variables to their corresponding
!   primitive representation.
!
!   Arguments:
!
!     u - the input array of conservative variables;
!     q - the output array of primitive variables;
!     s - the status flag;
!
!===============================================================================
!
  subroutine cons2prim_hd_adi(u, q, s)

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: u
    real(kind=8), dimension(:,:), intent(out) :: q
    integer                     , intent(out) :: s

    integer      :: i, p
    real(kind=8) :: ek, ei

!-------------------------------------------------------------------------------
!
    s = 0

    do i = 1, size(u,2)

      if (u(idn,i) > 0.0d+00) then
        q(idn,i) = u(idn,i)
        q(ivx,i) = u(imx,i) / u(idn,i)
        q(ivy,i) = u(imy,i) / u(idn,i)
        q(ivz,i) = u(imz,i) / u(idn,i)
        ek       = 0.5d+00 * sum(u(imx:imz,i) * q(ivx:ivz,i))
        ei       = u(ien,i) - ek
        if (ei > 0.0d+00) then
          q(ipr,i) = gammam1 * ei
        else
          s = 1
          go to 100
        end if
      else
        s = 1
        go to 100
      end if
    end do

    if (ns > 0) then
      do p = isl, isu
        q(p,:) = u(p,:) / u(idn,:)
      end do
    end if

    100 continue

!-------------------------------------------------------------------------------
!
  end subroutine cons2prim_hd_adi
!
!===============================================================================
!
! subroutine FLUXSPEED_HD_ADI:
! ---------------------------
!
!   Subroutine calculates physical fluxes and characteristic speeds from a
!   given equation system.
!
!   Arguments:
!
!     q - the input array of primitive variables;
!     u - the input array of conservative variables;
!     f - the output vector of fluxes;
!     c - the output vector of left- and right-going characteristic speeds;
!
!===============================================================================
!
  subroutine fluxspeed_hd_adi(q, u, f, c)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8), dimension(:,:)          , intent(in)  :: q, u
    real(kind=8), dimension(:,:)          , intent(out) :: f
    real(kind=8), dimension(:,:), optional, intent(out) :: c

! local variables
!
    integer      :: i
    real(kind=8) :: cs
!
!-------------------------------------------------------------------------------
!
! calculate the hydrodynamic fluxes
!
    do i = 1, size(q,2)

      f(idn,i) = u(imx,i)
      f(imx,i) = q(ivx,i) * u(imx,i)
      f(imy,i) = q(ivx,i) * u(imy,i)
      f(imz,i) = q(ivx,i) * u(imz,i)
      f(imx,i) = f(imx,i) + q(ipr,i)
      f(ien,i) = q(ivx,i) * (u(ien,i) + q(ipr,i))

    end do

! calculate the characteristic speeds
!
    if (present(c)) then

      do i = 1, size(q,2)

        cs     = sqrt(adiabatic_index * q(ipr,i) / q(idn,i))

        c(1,i) = q(ivx,i) - cs
        c(2,i) = q(ivx,i) + cs

      end do

    end if

!-------------------------------------------------------------------------------
!
  end subroutine fluxspeed_hd_adi
!
!===============================================================================
!
! function MAXSPEED_HD_ADI:
! ------------------------
!
!   Function scans the variable array and returns the maximum speed in within.
!
!   Arguments:
!
!     q - the array of primitive variables;
!
!
!===============================================================================
!
  function maxspeed_hd_adi(qq) result(maxspeed)

    use coordinates, only : nb, ne

    implicit none

    real(kind=8), dimension(:,:,:,:), intent(in) :: qq

    real(kind=8) :: maxspeed

    integer      :: i, j, k
    real(kind=8) :: vv, v, c

!-------------------------------------------------------------------------------
!
    maxspeed = 0.0d+00

#if NDIMS == 2
    k = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
    do k = nb, ne
#endif /* NDIMS == 3 */
      do j = nb, ne
        do i = nb, ne

! calculate the velocity amplitude
!
          vv = sum(qq(ivx:ivz,i,j,k) * qq(ivx:ivz,i,j,k))
          v  = sqrt(vv)

! calculate the adiabatic speed of sound
!
          c = sqrt(adiabatic_index * qq(ipr,i,j,k) / qq(idn,i,j,k))

! calculate the maximum speed
!
          maxspeed = max(maxspeed, v + c)

        end do ! i = nb, ne
      end do ! j = nb, ne
#if NDIMS == 3
    end do ! k = nb, ne
#endif /* NDIMS == 3 */

    return

!-------------------------------------------------------------------------------
!
  end function maxspeed_hd_adi
!
!===============================================================================
!
! subroutine GET_MAXIMUM_SPEEDS_HD_ADI:
! ------------------------------------
!
!   Subroutine determines the maximum characteristic speed and eigenvalue
!   in the input array of primitive variables.
!
!   Arguments:
!
!     qq - the input array of primitive variables;
!     vm - the maximum physical speed;
!     cm - the maximum eigenvalue;
!
!===============================================================================
!
  subroutine get_maximum_speeds_hd_adi(qq, vm, cm)

    use coordinates, only : nb, ne

    implicit none

    real(kind=8), dimension(:,:,:,:), intent(in)  :: qq
    real(kind=8)                    , intent(out) :: vm, cm

    integer      :: i, j, k
    real(kind=8) :: vl, vu, cc

!-------------------------------------------------------------------------------
!
    vm = 0.0d+00
    cm = 0.0d+00

#if NDIMS == 2
    k = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
    do k = nb, ne
#endif /* NDIMS == 3 */
      do j = nb, ne
        do i = nb, ne
          vl = minval(qq(ivx:ivz,i,j,k))
          vu = maxval(qq(ivx:ivz,i,j,k))
          vm = max(vm, abs(vl), abs(vu))

          cc = sqrt(adiabatic_index * qq(ipr,i,j,k) / qq(idn,i,j,k))
          cm = max(cm, abs(vl - cc), abs(vu + cc))
        end do
      end do
#if NDIMS == 3
    end do
#endif /* NDIMS == 3 */

!-------------------------------------------------------------------------------
!
  end subroutine get_maximum_speeds_hd_adi
!
!===============================================================================
!
! subroutine ESYSTEM_ROE_HD_ADI:
! -----------------------------
!
!   Subroutine computes eigenvalues and eigenvectors for a given set of
!   equations and input variables.
!
!   Arguments:
!
!     x - ratio of the perpendicular magnetic field component difference
!     y - ratio of the density
!     q - the intermediate Roe state vector;
!     c - the vector of eigenvalues;
!     r - the matrix of right eigenvectors;
!     l - the matrix of left eigenvectors;
!
!   References:
!
!     [1] Roe, P. L.
!         "Approximate Riemann Solvers, Parameter Vectors, and Difference
!          Schemes",
!         Journal of Computational Physics, 1981, 43, pp. 357-372
!     [2] Stone, J. M. & Gardiner, T. A.,
!         "ATHENA: A New Code for Astrophysical MHD",
!         The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177
!
!===============================================================================
!
  subroutine esystem_roe_hd_adi(x, y, q, c, r, l)

    implicit none

    real(kind=8)                , intent(in)    :: x, y
    real(kind=8), dimension(:)  , intent(in)    :: q
    real(kind=8), dimension(:)  , intent(inout) :: c
    real(kind=8), dimension(:,:), intent(inout) :: l, r

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

    real(kind=8)  :: vv, vh, c2, cc, vc, fc, fh, f1, f2, fv, fx

!-------------------------------------------------------------------------------
!
    if (first) then
      lvec(:,:)   = 0.0d+00
      rvec(:,:)   = 0.0d+00

      lvec(ivy,2) = 1.0d+00
      lvec(ivz,3) = 1.0d+00

      rvec(1,idn) = 1.0d+00
      rvec(2,ivy) = 1.0d+00
      rvec(3,ivz) = 1.0d+00
      rvec(4,idn) = 1.0d+00
      rvec(5,idn) = 1.0d+00

      first = .false.
    end if

! calculate characteristic speeds and useful variables
!
    vv = sum(q(ivx:ivz)**2)
    vh = 5.0d-01 * vv
    c2 = gammam1 * (q(ien) - vh)
    cc = sqrt(c2)
    vc = q(ivx) * cc
    fc = gammam1 / c2
    fh = 5.0d-01 * fc
    f1 = 5.0d-01 * vc / c2
    f2 = 5.0d-01 * cc / c2
    fv = fh * vh
    fx = fh * q(ivx)

! eigenvalues
!
    c(1) = q(ivx) - cc
    c(2) = q(ivx)
    c(3) = q(ivx)
    c(4) = q(ivx)
    c(5) = q(ivx) + cc

! update the varying elements of the left eigenvectors matrix
!
    lvec(idn,1) =   fv + f1
    lvec(ivx,1) = - fx - f2
    lvec(ivy,1) = - fh * q(ivy)
    lvec(ivz,1) = - fh * q(ivz)
    lvec(ien,1) = fh

    lvec(idn,2) = - q(ivy)

    lvec(idn,3) = - q(ivz)

    lvec(idn,4) = 1.0d+00 - fh * vv
    lvec(ivx,4) =   fc * q(ivx)
    lvec(ivy,4) =   fc * q(ivy)
    lvec(ivz,4) =   fc * q(ivz)
    lvec(ien,4) = - fc

    lvec(idn,5) =   fv - f1
    lvec(ivx,5) = - fx + f2
    lvec(ivy,5) = lvec(ivy,1)
    lvec(ivz,5) = lvec(ivz,1)
    lvec(ien,5) = fh

! update the varying elements of the right eigenvectors matrix
!
    rvec(1,ivx) = q(ivx) - cc
    rvec(1,ivy) = q(ivy)
    rvec(1,ivz) = q(ivz)
    rvec(1,ien) = q(ien) - vc

    rvec(2,ien) = q(ivy)

    rvec(3,ien) = q(ivz)

    rvec(4,ivx) = q(ivx)
    rvec(4,ivy) = q(ivy)
    rvec(4,ivz) = q(ivz)
    rvec(4,ien) = vh

    rvec(5,ivx) = q(ivx) + cc
    rvec(5,ivy) = q(ivy)
    rvec(5,ivz) = q(ivz)
    rvec(5,ien) = q(ien) + vc

! copy matrices of eigenvectors
!
    l(:,:) = lvec(:,:)
    r(:,:) = rvec(:,:)

!-------------------------------------------------------------------------------
!
  end subroutine esystem_roe_hd_adi
!
!*******************************************************************************
!
! ISOTHERMAL MAGNETOHYDRODYNAMIC EQUATIONS
!
!*******************************************************************************
!
!===============================================================================
!
! subroutine PRIM2CONS_MHD_ISO:
! ----------------------------
!
!   Subroutine converts primitive variables to their corresponding
!   conservative representation.
!
!   Arguments:
!
!     q - the input array of primitive variables;
!     u - the output array of conservative variables;
!     s - an optional flag indicating that passive scalars have
!         to be calculated too;
!
!===============================================================================
!
  subroutine prim2cons_mhd_iso(q, u, s)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: q
    real(kind=8), dimension(:,:), intent(out) :: u
    logical     , optional      , intent(in)  :: s

! local variables
!
    integer :: i, p
!
!-------------------------------------------------------------------------------
!
! iterate over all positions
!
    do i = 1, size(q,2)

      u(idn,i) = q(idn,i)
      u(imx,i) = q(idn,i) * q(ivx,i)
      u(imy,i) = q(idn,i) * q(ivy,i)
      u(imz,i) = q(idn,i) * q(ivz,i)
      u(ibx,i) = q(ibx,i)
      u(iby,i) = q(iby,i)
      u(ibz,i) = q(ibz,i)
      u(ibp,i) = q(ibp,i)

    end do

! update primitive passive scalars
!
    if (ns > 0 .and. present(s)) then
      if (s) then
        do p = isl, isu
          u(p,:) = q(p,:) * u(idn,:)
        end do
      end if
    end if

!-------------------------------------------------------------------------------
!
  end subroutine prim2cons_mhd_iso
!
!===============================================================================
!
! subroutine CONS2PRIM_MHD_ISO:
! ----------------------------
!
!   Subroutine converts conservative variables to their corresponding
!   primitive representation.
!
!   Arguments:
!
!     u - the input array of conservative variables;
!     q - the output array of primitive variables;
!     s - the status flag;
!
!===============================================================================
!
  subroutine cons2prim_mhd_iso(u, q, s)

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: u
    real(kind=8), dimension(:,:), intent(out) :: q
    integer                     , intent(out) :: s

    integer :: i, p

!-------------------------------------------------------------------------------
!
    s = 0

    do i = 1, size(u,2)

      if (u(idn,i) > 0.0d+00) then
        q(idn,i) = u(idn,i)
        q(ivx,i) = u(imx,i) / u(idn,i)
        q(ivy,i) = u(imy,i) / u(idn,i)
        q(ivz,i) = u(imz,i) / u(idn,i)
        q(ibx,i) = u(ibx,i)
        q(iby,i) = u(iby,i)
        q(ibz,i) = u(ibz,i)
        q(ibp,i) = u(ibp,i)
      else
        s = 1
        go to 100
      end if

    end do

    if (ns > 0) then
      do p = isl, isu
        q(p,:) = u(p,:) / u(idn,:)
      end do
    end if

    100 continue

!-------------------------------------------------------------------------------
!
  end subroutine cons2prim_mhd_iso
!
!===============================================================================
!
! subroutine FLUXSPEED_MHD_ISO:
! ----------------------------
!
!   Subroutine calculates physical fluxes and characteristic speeds from a
!   given equation system.
!
!   Arguments:
!
!     q - the input array of primitive variables;
!     u - the input array of conservative variables;
!     f - the output vector of fluxes;
!     c - the output vector of left- and right-going characteristic speeds;
!
!===============================================================================
!
  subroutine fluxspeed_mhd_iso(q, u, f, c)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8), dimension(:,:)          , intent(in)  :: q, u
    real(kind=8), dimension(:,:)          , intent(out) :: f
    real(kind=8), dimension(:,:), optional, intent(out) :: c

! local variables
!
    integer      :: i
    real(kind=8) :: by2, bz2, pt
    real(kind=8) :: fa, fb, fc, cf

! local arrays
!
    real(kind=8), dimension(size(q,2)) :: bx2, bb
!
!-------------------------------------------------------------------------------
!
! calculate the magnetohydrodynamic fluxes
!
    do i = 1, size(q,2)

      bx2(i) = q(ibx,i) * q(ibx,i)
      by2    = q(iby,i) * q(iby,i)
      bz2    = q(ibz,i) * q(ibz,i)
      bb(i)  = bx2(i) + by2 + bz2
      pt     = csnd2 * q(idn,i) + 0.5d+00 * bb(i)

      f(idn,i) = u(imx,i)
      f(imx,i) = q(ivx,i) * u(imx,i) - bx2(i)
      f(imy,i) = q(ivx,i) * u(imy,i) - q(ibx,i) * q(iby,i)
      f(imz,i) = q(ivx,i) * u(imz,i) - q(ibx,i) * q(ibz,i)
      f(imx,i) = f(imx,i) + pt
      f(ibx,i) = q(ibp,i)
      f(iby,i) = q(ivx,i) * q(iby,i) - q(ibx,i) * q(ivy,i)
      f(ibz,i) = q(ivx,i) * q(ibz,i) - q(ibx,i) * q(ivz,i)
      f(ibp,i) = cmax2 * q(ibx,i)

    end do

! calculate the characteristic speeds
!
    if (present(c)) then

      do i = 1, size(q,2)

        fa = csnd2 * q(idn,i)
        fb = fa + bb(i)
        fc = max(0.0d+00, fb * fb - 4.0d+00 * fa * bx2(i))
        cf = sqrt(max(0.5d+00 * (fb + sqrt(fc)), bb(i)) / q(idn,i))

        c(1,i) = q(ivx,i) - cf
        c(2,i) = q(ivx,i) + cf

      end do

    end if

!-------------------------------------------------------------------------------
!
  end subroutine fluxspeed_mhd_iso
!
!===============================================================================
!
! function MAXSPEED_MHD_ISO:
! -------------------------
!
!   Function scans the variable array and returns the maximum speed in within.
!
!   Arguments:
!
!     q - the array of primitive variables;
!
!===============================================================================
!
  function maxspeed_mhd_iso(qq) result(maxspeed)

    use coordinates, only : nb, ne

    implicit none

    real(kind=8), dimension(:,:,:,:), intent(in) :: qq

    real(kind=8) :: maxspeed

    integer      :: i, j, k
    real(kind=8) :: vv, bb, v, c

!-------------------------------------------------------------------------------
!
    maxspeed = 0.0d+00

#if NDIMS == 2
    k = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
    do k = nb, ne
#endif /* NDIMS == 3 */
      do j = nb, ne
        do i = nb, ne

! calculate the velocity amplitude
!
          vv = sum(qq(ivx:ivz,i,j,k) * qq(ivx:ivz,i,j,k))
          v  = sqrt(vv)
          bb = sum(qq(ibx:ibz,i,j,k) * qq(ibx:ibz,i,j,k))

! calculate the fast magnetosonic speed
!
          c = sqrt(csnd2 + bb / qq(idn,i,j,k))

! calculate the maximum of speed
!
          maxspeed = max(maxspeed, v + c)

        end do ! i = nb, ne
      end do ! j = nb, ne
#if NDIMS == 3
    end do ! k = nb, ne
#endif /* NDIMS == 3 */

    return

!-------------------------------------------------------------------------------
!
  end function maxspeed_mhd_iso
!
!===============================================================================
!
! subroutine GET_MAXIMUM_SPEEDS_MHD_ISO:
! -------------------------------------
!
!   Subroutine determines the maximum characteristic speed and eigenvalue
!   in the input array of primitive variables.
!
!   Arguments:
!
!     qq - the input array of primitive variables;
!     vm - the maximum physical speed;
!     cm - the maximum eigenvalue;
!
!===============================================================================
!
  subroutine get_maximum_speeds_mhd_iso(qq, vm, cm)

    use coordinates, only : nb, ne

    implicit none

    real(kind=8), dimension(:,:,:,:), intent(in)  :: qq
    real(kind=8)                    , intent(out) :: vm, cm

    integer      :: i, j, k
    real(kind=8) :: vl, vu, cc, xx

    real(kind=8), dimension(3) :: bb

!-------------------------------------------------------------------------------
!
    vm = 0.0d+00
    cm = 0.0d+00

#if NDIMS == 2
    k = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
    do k = nb, ne
#endif /* NDIMS == 3 */
      do j = nb, ne
        do i = nb, ne
          vl = minval(qq(ivx:ivz,i,j,k))
          vu = maxval(qq(ivx:ivz,i,j,k))
          vm = max(vm, abs(vl), abs(vu))

          bb = qq(ibx:ibz,i,j,k) * qq(ibx:ibz,i,j,k) / qq(idn,i,j,k)
          xx = csnd2 + bb(1) + bb(2) + bb(3)
          cc = max(0.0d+00, xx**2 - 4.0d+00 * csnd2 * bb(1))
          cc = sqrt(5.0d-01 * (xx + sqrt(cc)))
          cm = max(cm, abs(vl - cc), abs(vu + cc))
        end do
      end do
#if NDIMS == 3
    end do
#endif /* NDIMS == 3 */

!-------------------------------------------------------------------------------
!
  end subroutine get_maximum_speeds_mhd_iso
!
!===============================================================================
!
! subroutine ESYSTEM_ROE_MHD_ISO:
! ------------------------------
!
!   Subroutine computes eigenvalues and eigenvectors for a given set of
!   equations and input variables.
!
!   Arguments:
!
!     x - ratio of the perpendicular magnetic field component difference
!     y - ratio of the density
!     q - the intermediate Roe state vector;
!     c - the vector of eigenvalues;
!     r - the matrix of right eigenvectors;
!     l - the matrix of left eigenvectors;
!
!   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] Balsara, D. S.
!         "Linearized Formulation of the Riemann Problem for Adiabatic and
!          Isothermal Magnetohydrodynamics",
!         The Astrophysical Journal Suplement Series, 1998, 116, pp. 119-131
!
!===============================================================================
!
  subroutine esystem_roe_mhd_iso(x, y, q, c, r, l)

    implicit none

    real(kind=8)                , intent(in)    :: x, y
    real(kind=8), dimension(:)  , intent(in)    :: q
    real(kind=8), dimension(:)  , intent(inout) :: c
    real(kind=8), dimension(:,:), intent(inout) :: l, r

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

    real(kind=8) :: ca, ca2, ct2, cf, cf2, cs, cs2, cf2_cs2
    real(kind=8) :: br, br2, brs, br2s, tsum, tdif, twid_c2, twid_c
    real(kind=8) :: bty, btz, btys, btzs, bt2s, alf, als, norm
    real(kind=8) :: sqrtd, sgn, qf, qs, af_prime, as_prime
    real(kind=8) :: cff, css, af, as, afpb, aspb, q2s, q3s, vqstr

!-------------------------------------------------------------------------------
!
    if (first) then
      lvec(:,:) = 0.0d+00
      rvec(:,:) = 0.0d+00

      first = .false.
    end if

! coefficients for eigenvalues
!
    ca2  = q(ibx) * q(ibx) / q(idn)
    ca   = sqrt(ca2)
    br2  = sum(q(iby:ibz)**2)
    br2s = br2 * y
    ct2  = br2s / q(idn)

    twid_c2 = csnd2 + x
    tsum    = ca2 + ct2 + twid_c2
    tdif    = ca2 + ct2 - twid_c2
    cf2_cs2 = sqrt(tdif * tdif + 4.0d+00 * twid_c2 * ct2)
    cf2     = 0.5d+00 * (tsum + cf2_cs2)
    cf      = sqrt(cf2)
    cs2     = twid_c2 * ca2 / cf2
    cs      = sqrt(cs2)

! eigenvalues
!
    c(1) = q(ivx) - cf
    c(2) = q(ivx) - ca
    c(3) = q(ivx) - cs
    c(4) = q(ivx)
    c(5) = q(ivx) + cs
    c(6) = q(ivx) + ca
    c(7) = q(ivx) + cf
    c(8) = c(7)

! calculate the eigenvectors only if the waves propagate in both direction
!
    if (c(1) >= 0.0d+00 .or. c(7) <= 0.0d+00) return

! remaining coefficients
!
    br    = sqrt(br2)
    brs   = sqrt(br2s)
    if (abs(br) > 0.0d+00) then
      bty = q(iby) / br
      btz = q(ibz) / br
    else
      bty = 1.0d+00
      btz = 0.0d+00
    end if
    btys = bty / sqrt(y)
    btzs = btz / sqrt(y)
    bt2s = btys * btys + btzs * btzs

    if (.not. abs(cf2 - cs2) > 0.0d+00) then
      alf = 1.0d+00
      als = 0.0d+00
    else if ((twid_c2 - cs2) <= 0.0d+00) then
      alf = 0.0d+00
      als = 1.0d+00
    else if ((cf2 - twid_c2) <= 0.0d+00) then
      alf = 1.0d+00
      als = 0.0d+00
    else
      alf = sqrt((twid_c2 - cs2) / (cf2 - cs2))
      als = sqrt((cf2 - twid_c2) / (cf2 - cs2))
    end if

    sqrtd    = sqrt(q(idn))
    sgn      = sign(1.0d+00, q(ibx))
    twid_c   = sqrt(twid_c2)
    qf       = cf * alf * sgn
    qs       = cs * als * sgn
    af_prime = twid_c * alf / sqrtd
    as_prime = twid_c * als / sqrtd

! === update the varying elements of the right eigenvectors matrix
!
! left-going fast wave
!
    rvec(1,idn) = alf
    rvec(1,ivx) = alf * c(1)
    rvec(1,ivy) = alf * q(ivy) + qs * btys
    rvec(1,ivz) = alf * q(ivz) + qs * btzs
    rvec(1,iby) = as_prime * btys
    rvec(1,ibz) = as_prime * btzs

! left-going Alfvèn wave
!
    rvec(2,ivy) = - btz
    rvec(2,ivz) =   bty
    rvec(2,iby) = - btz * sgn / sqrtd
    rvec(2,ibz) =   bty * sgn / sqrtd

! left-going slow wave
!
    rvec(3,idn) = als
    rvec(3,ivx) = als * c(3)
    rvec(3,ivy) = als * q(ivy) - qf * btys
    rvec(3,ivz) = als * q(ivz) - qf * btzs
    rvec(3,iby) = - af_prime * btys
    rvec(3,ibz) = - af_prime * btzs

! right-going slow wave
!
    rvec(5,idn) = als
    rvec(5,ivx) = als * c(5)
    rvec(5,ivy) = als * q(ivy) + qf * btys
    rvec(5,ivz) = als * q(ivz) + qf * btzs
    rvec(5,iby) = rvec(3,iby)
    rvec(5,ibz) = rvec(3,ibz)

! right-going Alfvèn wave
!
    rvec(6,ivy) =   btz
    rvec(6,ivz) = - bty
    rvec(6,iby) = rvec(2,iby)
    rvec(6,ibz) = rvec(2,ibz)

! right-going fast wave
!
    rvec(7,idn) = alf
    rvec(7,ivx) = alf * c(7)
    rvec(7,ivy) = alf * q(ivy) - qs * btys
    rvec(7,ivz) = alf * q(ivz) - qs * btzs
    rvec(7,iby) = rvec(1,iby)
    rvec(7,ibz) = rvec(1,ibz)

! === update the varying elements of the left eigenvectors matrix
!
    norm = 2.0d+00 * twid_c2
    cff  = alf * cf / norm
    css  = als * cs / norm
    qf   = qf / norm
    qs   = qs / norm
    af   = af_prime * q(idn) / norm
    as   = as_prime * q(idn) / norm
    afpb = af_prime * brs / norm
    aspb = as_prime * brs / norm

    q2s   = btys / bt2s
    q3s   = btzs / bt2s
    vqstr = q(ivy) * q2s + q(ivz) * q3s

! left-going fast wave
!
    lvec(idn,1) = cff * c(7) - qs * vqstr - aspb
    lvec(ivx,1) = - cff
    lvec(ivy,1) = qs * q2s
    lvec(ivz,1) = qs * q3s
    lvec(iby,1) = as * q2s
    lvec(ibz,1) = as * q3s

! left-going Alfvèn wave
!
    lvec(idn,2) =   0.5d+00 * (q(ivy) * btz - q(ivz) * bty)
    lvec(ivy,2) = - 0.5d+00 * btz
    lvec(ivz,2) =   0.5d+00 * bty
    lvec(iby,2) = - 0.5d+00 * sqrtd * btz * sgn
    lvec(ibz,2) =   0.5d+00 * sqrtd * bty * sgn

! left-going slow wave
!
    lvec(idn,3) =   css * c(5) + qf * vqstr + afpb
    lvec(ivx,3) = - css
    lvec(ivy,3) = - qf * q2s
    lvec(ivz,3) = - qf * q3s
    lvec(iby,3) = - af * q2s
    lvec(ibz,3) = - af * q3s

! right-going slow wave
!
    lvec(idn,5) = - css * c(3) - qf * vqstr + afpb
    lvec(ivx,5) =   css
    lvec(ivy,5) = - lvec(ivy,3)
    lvec(ivz,5) = - lvec(ivz,3)
    lvec(iby,5) =   lvec(iby,3)
    lvec(ibz,5) =   lvec(ibz,3)

! right-going Alfvèn wave
!
    lvec(idn,6) = - lvec(idn,2)
    lvec(ivy,6) = - lvec(ivy,2)
    lvec(ivz,6) = - lvec(ivz,2)
    lvec(iby,6) =   lvec(iby,2)
    lvec(ibz,6) =   lvec(ibz,2)

! right-going fast wave
!
    lvec(idn,7) = - cff * c(1) + qs * vqstr - aspb
    lvec(ivx,7) =   cff
    lvec(ivy,7) = - lvec(ivy,1)
    lvec(ivz,7) = - lvec(ivz,1)
    lvec(iby,7) =   lvec(iby,1)
    lvec(ibz,7) =   lvec(ibz,1)

! copy matrices of eigenvectors
!
    l(:,:) = lvec(:,:)
    r(:,:) = rvec(:,:)

!-------------------------------------------------------------------------------
!
  end subroutine esystem_roe_mhd_iso
!
!*******************************************************************************
!
! ADIABATIC MAGNETOHYDRODYNAMIC EQUATIONS
!
!*******************************************************************************
!
!===============================================================================
!
! subroutine PRIM2CONS_MHD_ADI:
! ----------------------------
!
!   Subroutine converts primitive variables to their corresponding
!   conservative representation.
!
!   Arguments:
!
!     q - the input array of primitive variables;
!     u - the output array of conservative variables;
!     s - an optional flag indicating that passive scalars have
!         to be calculated too;
!
!===============================================================================
!
  subroutine prim2cons_mhd_adi(q, u, s)

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: q
    real(kind=8), dimension(:,:), intent(out) :: u
    logical     , optional      , intent(in)  :: s

    integer      :: i, p
    real(kind=8) :: ei, ek, em, ep

!-------------------------------------------------------------------------------
!
    do i = 1, size(q,2)

      u(idn,i) = q(idn,i)
      u(imx,i) = q(idn,i) * q(ivx,i)
      u(imy,i) = q(idn,i) * q(ivy,i)
      u(imz,i) = q(idn,i) * q(ivz,i)
      u(ibx,i) = q(ibx,i)
      u(iby,i) = q(iby,i)
      u(ibz,i) = q(ibz,i)
      u(ibp,i) = q(ibp,i)
      ei       = gammam1i * q(ipr,i)
      ek       = sum(u(imx:imz,i) * q(ivx:ivz,i))
      em       = sum(q(ibx:ibz,i) * q(ibx:ibz,i))
      ep       =     q(ibp,i)     * q(ibp,i)
      u(ien,i) = ei + 5.0d-01 * (ek + em + ep)

    end do

    if (ns > 0 .and. present(s)) then
      if (s) then
        do p = isl, isu
          u(p,:) = q(p,:) * u(idn,:)
        end do
      end if
    end if

!-------------------------------------------------------------------------------
!
  end subroutine prim2cons_mhd_adi
!
!===============================================================================
!
! subroutine CONS2PRIM_MHD_ADI:
! ----------------------------
!
!   Subroutine converts conservative variables to their corresponding
!   primitive representation.
!
!   Arguments:
!
!     u - the input array of conservative variables;
!     q - the output array of primitive variables;
!     s - the status flag;
!
!===============================================================================
!
  subroutine cons2prim_mhd_adi(u, q, s)

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: u
    real(kind=8), dimension(:,:), intent(out) :: q
    integer                     , intent(out) :: s

    integer      :: i, p
    real(kind=8) :: ei, ek, em, ep

!-------------------------------------------------------------------------------
!
    s = 0

    do i = 1, size(u,2)

      if (u(idn,i) > 0.0d+00) then
        q(idn,i) = u(idn,i)
        q(ivx,i) = u(imx,i) / u(idn,i)
        q(ivy,i) = u(imy,i) / u(idn,i)
        q(ivz,i) = u(imz,i) / u(idn,i)
        q(ibx,i) = u(ibx,i)
        q(iby,i) = u(iby,i)
        q(ibz,i) = u(ibz,i)
        q(ibp,i) = u(ibp,i)
        ek       = sum(u(imx:imz,i) * q(ivx:ivz,i))
        em       = sum(q(ibx:ibz,i) * q(ibx:ibz,i))
        ep       =     q(ibp,i)     * q(ibp,i)
        ei       = u(ien,i) - 5.0d-01 * (ek + em + ep)
        if (ei > 0.0d+00) then
          q(ipr,i) = gammam1 * ei
        else
          s = 1
          go to 100
        end if
      else
        s = 1
        go to 100
      end if

    end do

    if (ns > 0) then
      do p = isl, isu
        q(p,:) = u(p,:) / u(idn,:)
      end do
    end if

    100 continue

!-------------------------------------------------------------------------------
!
  end subroutine cons2prim_mhd_adi
!
!===============================================================================
!
! subroutine FLUXSPEED_MHD_ADI:
! ----------------------------
!
!   Subroutine calculates physical fluxes and characteristic speeds from a
!   given equation system.
!
!   Arguments:
!
!     q - the input array of primitive variables;
!     u - the input array of conservative variables;
!     f - the output vector of fluxes;
!     c - the output vector of left- and right-going characteristic speeds;
!
!===============================================================================
!
  subroutine fluxspeed_mhd_adi(q, u, f, c)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8), dimension(:,:)          , intent(in)  :: q, u
    real(kind=8), dimension(:,:)          , intent(out) :: f
    real(kind=8), dimension(:,:), optional, intent(out) :: c

! local variables
!
    integer      :: i
    real(kind=8) :: by2, bz2, pt
    real(kind=8) :: vb
    real(kind=8) :: fa, fb, fc, cf

! local arrays
!
    real(kind=8), dimension(size(q,2)) :: bx2, bb
!
!-------------------------------------------------------------------------------
!
! calculate the magnetohydrodynamic fluxes
!
    do i = 1, size(q,2)

      bx2(i) = q(ibx,i) * q(ibx,i)
      by2    = q(iby,i) * q(iby,i)
      bz2    = q(ibz,i) * q(ibz,i)
      bb(i)  = bx2(i) + by2 + bz2
      vb  = sum(q(ivx:ivz,i) * q(ibx:ibz,i))
      pt  = q(ipr,i) + 0.5d+00 * bb(i)

      f(idn,i) = u(imx,i)
      f(imx,i) = q(ivx,i) * u(imx,i) - bx2(i)
      f(imy,i) = q(ivx,i) * u(imy,i) - q(ibx,i) * q(iby,i)
      f(imz,i) = q(ivx,i) * u(imz,i) - q(ibx,i) * q(ibz,i)
      f(imx,i) = f(imx,i) + pt
      f(ibx,i) = q(ibp,i)
      f(iby,i) = q(ivx,i) * q(iby,i) - q(ibx,i) * q(ivy,i)
      f(ibz,i) = q(ivx,i) * q(ibz,i) - q(ibx,i) * q(ivz,i)
      f(ibp,i) = cmax2 * q(ibx,i)
      f(ien,i) = q(ivx,i) * (u(ien,i) + pt) - q(ibx,i) * vb

    end do

! calculate the characteristic speeds
!
    if (present(c)) then

      do i = 1, size(q,2)

        fa = adiabatic_index * q(ipr,i)
        fb = fa + bb(i)
        fc = max(0.0d+00, fb * fb - 4.0d+00 * fa * bx2(i))
        cf = sqrt(max(0.5d+00 * (fb + sqrt(fc)), bb(i)) / q(idn,i))

        c(1,i) = q(ivx,i) - cf
        c(2,i) = q(ivx,i) + cf

      end do

    end if

!-------------------------------------------------------------------------------
!
  end subroutine fluxspeed_mhd_adi
!
!===============================================================================
!
! function MAXSPEED_MHD_ADI:
! -------------------------
!
!   Function scans the variable array and returns the maximum speed in within.
!
!   Arguments:
!
!     q - the array of primitive variables;
!
!===============================================================================
!
  function maxspeed_mhd_adi(qq) result(maxspeed)

    use coordinates, only : nb, ne

    implicit none

    real(kind=8), dimension(:,:,:,:), intent(in) :: qq

    real(kind=8) :: maxspeed

    integer      :: i, j, k
    real(kind=8) :: vv, bb, v, c

!-------------------------------------------------------------------------------
!
    maxspeed = 0.0d+00

#if NDIMS == 2
    k = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
    do k = nb, ne
#endif /* NDIMS == 3 */
      do j = nb, ne
        do i = nb, ne

! calculate the velocity amplitude
!
          vv = sum(qq(ivx:ivz,i,j,k) * qq(ivx:ivz,i,j,k))
          v  = sqrt(vv)
          bb = sum(qq(ibx:ibz,i,j,k) * qq(ibx:ibz,i,j,k))

! calculate the fast magnetosonic speed
!
          c = sqrt((adiabatic_index * qq(ipr,i,j,k) + bb) / qq(idn,i,j,k))

! calculate the maximum of speed
!
          maxspeed = max(maxspeed, v + c)

        end do ! i = nb, ne
      end do ! j = nb, ne
#if NDIMS == 3
    end do ! k = nb, ne
#endif /* NDIMS == 3 */

    return

!-------------------------------------------------------------------------------
!
  end function maxspeed_mhd_adi
!
!===============================================================================
!
! subroutine GET_MAXIMUM_SPEEDS_MHD_ADI:
! -------------------------------------
!
!   Subroutine determines the maximum characteristic speed and eigenvalue
!   in the input array of primitive variables.
!
!   Arguments:
!
!     qq - the input array of primitive variables;
!     vm - the maximum physical speed;
!     cm - the maximum eigenvalue;
!
!===============================================================================
!
  subroutine get_maximum_speeds_mhd_adi(qq, vm, cm)

    use coordinates, only : nb, ne

    implicit none

    real(kind=8), dimension(:,:,:,:), intent(in)  :: qq
    real(kind=8)                    , intent(out) :: vm, cm

    integer      :: i, j, k
    real(kind=8) :: vl, vu, aa, cc, xx

    real(kind=8), dimension(3) :: bb

!-------------------------------------------------------------------------------
!
    vm = 0.0d+00
    cm = 0.0d+00

#if NDIMS == 2
    k = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
    do k = nb, ne
#endif /* NDIMS == 3 */
      do j = nb, ne
        do i = nb, ne
          vl = minval(qq(ivx:ivz,i,j,k))
          vu = maxval(qq(ivx:ivz,i,j,k))
          vm = max(vm, abs(vl), abs(vu))

          aa = adiabatic_index * qq(ipr,i,j,k) / qq(idn,i,j,k)
          bb = qq(ibx:ibz,i,j,k) * qq(ibx:ibz,i,j,k) / qq(idn,i,j,k)
          xx = aa + bb(1) + bb(2) + bb(3)
          cc = max(0.0d+00, xx**2 - 4.0d+00 * aa * bb(1))
          cc = sqrt(5.0d-01 * (xx + sqrt(cc)))
          cm = max(cm, abs(vl - cc), abs(vu + cc))
        end do
      end do
#if NDIMS == 3
    end do
#endif /* NDIMS == 3 */

!-------------------------------------------------------------------------------
!
  end subroutine get_maximum_speeds_mhd_adi
!
!===============================================================================
!
! subroutine ESYSTEM_ROE_MHD_ADI:
! ------------------------------
!
!   Subroutine computes eigenvalues and eigenvectors for a given set of
!   equations and input variables.
!
!   Arguments:
!
!     x - ratio of the perpendicular magnetic field component difference
!     y - ratio of the density
!     q - the intermediate Roe state vector;
!     c - the vector of eigenvalues;
!     r - the matrix of right eigenvectors;
!     l - the matrix of left eigenvectors;
!
!   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] Balsara, D. S.
!         "Linearized Formulation of the Riemann Problem for Adiabatic and
!          Isothermal Magnetohydrodynamics",
!         The Astrophysical Journal Suplement Series, 1998, 116, pp. 119-131
!
!===============================================================================
!
  subroutine esystem_roe_mhd_adi(x, y, q, c, r, l)

    implicit none

    real(kind=8)                , intent(in)    :: x, y
    real(kind=8), dimension(:)  , intent(in)    :: q
    real(kind=8), dimension(:)  , intent(inout) :: c
    real(kind=8), dimension(:,:), intent(inout) :: l, r

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

    real(kind=8) :: ca, ca2, cf, cf2, cs, cs2, cf2_cs2, ct2
    real(kind=8) :: v2, br2, br, br2s, brs, hp
    real(kind=8) :: tsum, tdif, twid_a2, twid_a
    real(kind=8) :: bty, btz, btys, btzs, bt2s, vbet
    real(kind=8) :: alf, als, af_prime, as_prime
    real(kind=8) :: sqrtd, sgn, qf, qs, afpbb, aspbb
    real(kind=8) :: qa, qb, qc, qd, q2s, q3s
    real(kind=8) :: norm, cff, css, af, as, afpb, aspb, vqstr

    real(kind=8), parameter :: eps = epsilon(1.0d+00)

!-------------------------------------------------------------------------------
!
    if (first) then
      gammam2  = adiabatic_index - 2.0d+00

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

      first = .false.
    end if

! coefficients
!
    ca2     = q(ibx) * q(ibx) / q(idn)
    ca      = sqrt(ca2)
    v2      = sum(q(ivx:ivz)**2)
    br2     = sum(q(iby:ibz)**2)
    br2s    = (gammam1 - gammam2 * y) * br2
    hp      = q(ien) - (ca2 + br2 / q(idn))
    twid_a2 = max(eps, (gammam1 * (hp - 0.5d+00 * v2) - gammam2 * x))
    ct2     = br2s / q(idn)
    tsum    = ca2 + ct2 + twid_a2
    tdif    = ca2 + ct2 - twid_a2
    cf2_cs2 = sqrt((tdif * tdif + 4.0d+00 * twid_a2 * ct2))
    cf2     = 0.5d+00 * (tsum + cf2_cs2)
    cf      = sqrt(cf2)
    cs2     = twid_a2 * ca2 / cf2
    cs      = sqrt(cs2)

! eigenvalues
!
    c(1) = q(ivx) - cf
    c(2) = q(ivx) - ca
    c(3) = q(ivx) - cs
    c(4) = q(ivx)
    c(5) = q(ivx)
    c(6) = q(ivx) + cs
    c(7) = q(ivx) + ca
    c(8) = q(ivx) + cf
    c(9) = c(8)

! eigenvectors only for the case of waves propagating in both direction
!
    if (c(1) >= 0.0d+00 .or. c(8) <= 0.0d+00) return

! remaining coefficients
!
    br  = sqrt(br2)
    brs = sqrt(br2s)
    if (abs(br) > 0.0d+00) then
      bty = q(iby) / br
      btz = q(ibz) / br
    else
      bty = 1.0d+00
      btz = 0.0d+00
    end if
    btys = bty / sqrt(gammam1 - gammam2 * y)
    btzs = btz / sqrt(gammam1 - gammam2 * y)
    bt2s = btys * btys + btzs * btzs
    vbet = q(ivy) * btys + q(ivz) * btzs

    if ( .not. abs(cf2 - cs2) > 0.0d+00 ) then
      alf = 1.0d+00
      als = 0.0d+00
    else if ( (twid_a2 - cs2) <= 0.0d+00 ) then
      alf = 0.0d+00
      als = 1.0d+00
    else if ( (cf2 - twid_a2) <= 0.0d+00 ) then
      alf = 1.0d+00
      als = 0.0d+00
    else
      alf = sqrt((twid_a2 - cs2) / (cf2 - cs2))
      als = sqrt((cf2 - twid_a2) / (cf2 - cs2))
    end if

    sqrtd    = sqrt(q(idn))
    sgn      = sign(1.0d+00, q(ibx))
    twid_a   = sqrt(twid_a2)
    qf       = cf * alf * sgn
    qs       = cs * als * sgn
    af_prime = twid_a * alf / sqrtd
    as_prime = twid_a * als / sqrtd
    afpbb    = af_prime * brs * bt2s
    aspbb    = as_prime * brs * bt2s

! === update the varying elements of the right eigenvectors matrix
!
    rvec(1,idn) = alf
    rvec(3,idn) = als
    rvec(6,idn) = als
    rvec(8,idn) = alf

    rvec(1,ivx) = alf * c(1)
    rvec(3,ivx) = als * c(3)
    rvec(4,ivx) = q(ivx)
    rvec(6,ivx) = als * c(6)
    rvec(8,ivx) = alf * c(8)

    qa = alf * q(ivy)
    qb = als * q(ivy)
    qc = qs * btys
    qd = qf * btys

    rvec(1,ivy) = qa + qc
    rvec(2,ivy) = - btz
    rvec(3,ivy) = qb - qd
    rvec(4,ivy) = q(ivy)
    rvec(6,ivy) = qb + qd
    rvec(7,ivy) = btz
    rvec(8,ivy) = qa - qc

    qa = alf * q(ivz)
    qb = als * q(ivz)
    qc = qs * btzs
    qd = qf * btzs

    rvec(1,ivz) = qa + qc
    rvec(2,ivz) = bty
    rvec(3,ivz) = qb - qd
    rvec(4,ivz) = q(ivz)
    rvec(6,ivz) = qb + qd
    rvec(7,ivz) = - bty
    rvec(8,ivz) = qa - qc

    rvec(1,ipr) = alf * (hp - q(ivx) * cf) + qs * vbet + aspbb
    rvec(2,ipr) = -(q(ivy) * btz - q(ivz) * bty)
    rvec(3,ipr) = als * (hp - q(ivx) * cs) - qf * vbet - afpbb
    rvec(4,ipr) = 0.5d+00 * v2 + gammam2 * x / gammam1
    rvec(6,ipr) = als * (hp + q(ivx) * cs) + qf * vbet - afpbb
    rvec(7,ipr) = - rvec(2,ipr)
    rvec(8,ipr) = alf * (hp + q(ivx) * cf) - qs * vbet + aspbb

    rvec(1,iby) = as_prime * btys
    rvec(2,iby) = - btz * sgn / sqrtd
    rvec(3,iby) = - af_prime * btys
    rvec(6,iby) = rvec(3,iby)
    rvec(7,iby) = rvec(2,iby)
    rvec(8,iby) = rvec(1,iby)

    rvec(1,ibz) = as_prime * btzs
    rvec(2,ibz) = bty * sgn / sqrtd
    rvec(3,ibz) = - af_prime * btzs
    rvec(6,ibz) = rvec(3,ibz)
    rvec(7,ibz) = rvec(2,ibz)
    rvec(8,ibz) = rvec(1,ibz)

! === update the varying elements of the left eigenvectors matrix
!
    norm = 2.0d+00 * twid_a2
    cff  = alf * cf / norm
    css  = als * cs / norm
    qf   = qf / norm
    qs   = qs / norm
    af   = af_prime * q(idn) / norm
    as   = as_prime * q(idn) / norm
    afpb = af_prime * brs / norm
    aspb = as_prime * brs / norm

    norm  = norm / gammam1
    alf   = alf / norm
    als   = als / norm
    q2s   = btys / bt2s
    q3s   = btzs / bt2s
    vqstr = (q(ivy) * q2s + q(ivz) * q3s)

! left-going fast wave
!
    lvec(idn,1) =   alf * (v2 - hp) + cff * (cf + q(ivx)) - qs * vqstr - aspb
    lvec(ipr,1) =   alf
    lvec(ivx,1) = - alf * q(ivx) - cff
    lvec(ivy,1) = - alf * q(ivy) + qs * q2s
    lvec(ivz,1) = - alf * q(ivz) + qs * q3s
    lvec(iby,1) = as * q2s - alf * q(iby)
    lvec(ibz,1) = as * q3s - alf * q(ibz)

! left-going Alfvèn wave
!
    lvec(idn,2) =   5.0d-01 * (q(ivy) * btz - q(ivz) * bty)
    lvec(ivy,2) = - 5.0d-01 * btz
    lvec(ivz,2) =   5.0d-01 * bty
    lvec(iby,2) = - 5.0d-01 * sqrtd * btz * sgn
    lvec(ibz,2) =   5.0d-01 * sqrtd * bty * sgn

! left-going slow wave
!
    lvec(idn,3) = als * (v2 - hp) + css * (cs + q(ivx)) + qf * vqstr + afpb
    lvec(ipr,3) = als
    lvec(ivx,3) = - als * q(ivx) - css
    lvec(ivy,3) = - als * q(ivy) - qf * q2s
    lvec(ivz,3) = - als * q(ivz) - qf * q3s
    lvec(iby,3) = - af * q2s - als * q(iby)
    lvec(ibz,3) = - af * q3s - als * q(ibz)

! entropy wave
!
    lvec(idn,4) = 1.0d+00 - (5.0d-01 * v2 - gammam2 * x / gammam1) / twid_a2
    lvec(ipr,4) = - 1.0d+00 / twid_a2
    lvec(ivx,4) = q(ivx) / twid_a2
    lvec(ivy,4) = q(ivy) / twid_a2
    lvec(ivz,4) = q(ivz) / twid_a2
    lvec(iby,4) = q(iby) / twid_a2
    lvec(ibz,4) = q(ibz) / twid_a2

! right-going slow wave
!
    lvec(idn,6) = als * (v2 - hp) + css * (cs - q(ivx)) - qf * vqstr + afpb
    lvec(ipr,6) = als
    lvec(ivx,6) = - als * q(ivx) + css
    lvec(ivy,6) = - als * q(ivy) + qf * q2s
    lvec(ivz,6) = - als * q(ivz) + qf * q3s
    lvec(iby,6) = lvec(iby,3)
    lvec(ibz,6) = lvec(ibz,3)

! right-going Alfvèn wave
!
    lvec(idn,7) = - lvec(idn,2)
    lvec(ivy,7) = - lvec(ivy,2)
    lvec(ivz,7) = - lvec(ivz,2)
    lvec(iby,7) =   lvec(iby,2)
    lvec(ibz,7) =   lvec(ibz,2)

! right-going fast wave
!
    lvec(idn,8) = alf * (v2 - hp) + cff * (cf - q(ivx)) + qs * vqstr - aspb
    lvec(ipr,8) = alf
    lvec(ivx,8) = - alf * q(ivx) + cff
    lvec(ivy,8) = - alf * q(ivy) - qs * q2s
    lvec(ivz,8) = - alf * q(ivz) - qs * q3s
    lvec(iby,8) = lvec(iby,1)
    lvec(ibz,8) = lvec(ibz,1)

! copy matrices of eigenvectors
!
    l(:,:) = lvec(:,:)
    r(:,:) = rvec(:,:)

!-------------------------------------------------------------------------------
!
  end subroutine esystem_roe_mhd_adi
!
!*******************************************************************************
!
! ADIABATIC SPECIAL RELATIVITY HYDRODYNAMIC EQUATIONS
!
!*******************************************************************************
!
!===============================================================================
!
! subroutine PRIM2CONS_SRHD_ADI:
! -----------------------------
!
!   Subroutine converts primitive variables to their corresponding
!   conservative representation.
!
!   Arguments:
!
!     q - the input array of primitive variables;
!     u - the output array of conservative variables;
!     s - an optional flag indicating that passive scalars have
!         to be calculated too;
!
!===============================================================================
!
  subroutine prim2cons_srhd_adi(q, u, s)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: q
    real(kind=8), dimension(:,:), intent(out) :: u
    logical     , optional      , intent(in)  :: s

! local variables
!
    integer      :: i, p
    real(kind=8) :: vv, vm, vs, ww
!
!-------------------------------------------------------------------------------
!
! iterate over all positions
!
    do i = 1, size(q,2)

! calculate the square of velocity, the Lorentz factor and specific enthalpy
!
      vv  = sum(q(ivx:ivz,i) * q(ivx:ivz,i))
      vm  = 1.0d+00 - vv
      vs  = sqrt(vm)
      ww  = (q(idn,i) + q(ipr,i) / gammaxi) / vm

! calculate conservative variables
!
      u(idn,i) =      q(idn,i) / vs
      u(imx,i) = ww * q(ivx,i)
      u(imy,i) = ww * q(ivy,i)
      u(imz,i) = ww * q(ivz,i)
      u(ien,i) = ww - q(ipr,i) - u(idn,i)

    end do

! update primitive passive scalars
!
    if (ns > 0 .and. present(s)) then
      if (s) then
        do p = isl, isu
          u(p,:) = q(p,:) * u(idn,:)
        end do
      end if
    end if

!-------------------------------------------------------------------------------
!
  end subroutine prim2cons_srhd_adi
!
!===============================================================================
!
! subroutine CONS2PRIM_SRHD_ADI:
! -----------------------------
!
!   Subroutine converts conservative variables to their corresponding
!   primitive representation using an interative method.
!
!   Arguments:
!
!     u - the input array of conservative variables;
!     q - the output array of primitive variables;
!     s - the status flag;
!
!===============================================================================
!
  subroutine cons2prim_srhd_adi(u, q, s)

    use helpers, only : print_message

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: u
    real(kind=8), dimension(:,:), intent(out) :: q
    integer                     , intent(out) :: s

    logical      :: info
    integer      :: i, p
    real(kind=8) :: mm, bb, mb, en, dn
    real(kind=8) :: w , vv, vm, vs

    character(len=80) :: msg

    character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_srhd_adi()'

!-------------------------------------------------------------------------------
!
    s = 0

    do i = 1, size(u,2)

! prepare variables which do not change during the Newton-Ralphson iterations
!
      mm  = sum(u(imx:imz,i) * u(imx:imz,i))
      en  = u(ien,i) + u(idn,i)
      dn  = u(idn,i)

! find the exact W using the Newton-Ralphson interative method
!
      call nr_iterate(mm, bb, mb, en, dn, w, vv, info)

! if info is .true., the solution was found
!
      if (info) then

! prepare coefficients
!
        vm = 1.0d+00 - vv
        vs = sqrt(vm)

! calculate the primitive variables
!
        q(idn,i) = dn * vs
        q(ivx,i) = u(imx,i) / w
        q(ivy,i) = u(imy,i) / w
        q(ivz,i) = u(imz,i) / w
        q(ipr,i) = w - en

      else ! cannot find physical solution

        call print_message(loc, "Conversion to physical primitive state failed!")
        write(msg,"(a,5(1x,1es24.16e3))") "U = ", u(:,i)
        call print_message(loc, msg)
        write(msg,"(a,3(1x,1es24.16e3))") "D, |m|², E = ", dn, mm, en
        call print_message(loc, msg)

        s = 1
        go to 100

      end if

    end do

! update primitive passive scalars
!
    if (ns > 0) then
      do p = isl, isu
        q(p,:) = u(p,:) / u(idn,:)
      end do
    end if

    100 continue

!-------------------------------------------------------------------------------
!
  end subroutine cons2prim_srhd_adi
!
!===============================================================================
!
! subroutine FLUXSPEED_SRHD_ADI:
! -----------------------------
!
!   Subroutine calculates physical fluxes and characteristic speeds from a
!   given equation system.
!
!   Arguments:
!
!     q - the input array of primitive variables;
!     u - the input array of conservative variables;
!     f - the output vector of fluxes;
!     c - the output vector of left- and right-going characteristic speeds;
!
!   References:
!
!     [1] Mignone, A., Bodo, G.,
!         "An HLLC Riemann solver for relativistic flows - I. Hydrodynamics",
!         Monthly Notices of the Royal Astronomical Society, 2005, 364, 126-136
!
!===============================================================================
!
  subroutine fluxspeed_srhd_adi(q, u, f, c)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8), dimension(:,:)          , intent(in)  :: q, u
    real(kind=8), dimension(:,:)          , intent(out) :: f
    real(kind=8), dimension(:,:), optional, intent(out) :: c

! local variables
!
    integer      :: i
    real(kind=8) :: vv, ww, c2, ss, cc, fc
!
!-------------------------------------------------------------------------------
!
! calculate the relativistic hydrodynamic fluxes (eq. 2 in [1])
!
    do i = 1, size(q,2)

      f(idn,i) = u(idn,i) * q(ivx,i)
      f(imx,i) = u(imx,i) * q(ivx,i) + q(ipr,i)
      f(imy,i) = u(imy,i) * q(ivx,i)
      f(imz,i) = u(imz,i) * q(ivx,i)
      f(ien,i) = u(imx,i) - f(idn,i)

    end do

! calculate characteristic speeds (eq. 23 in [1])
!
    if (present(c)) then

      do i = 1, size(q,2)

        ww     = q(idn,i) + q(ipr,i) / gammaxi
        c2     = adiabatic_index * q(ipr,i) / ww
        vv     = sum(q(ivx:ivz,i) * q(ivx:ivz,i))
        ss     = c2 * (1.0d+00 - vv) / (1.0d+00 - c2)
        fc     = 1.0d+00 + ss
        cc     = sqrt(ss * (fc - q(ivx,i)**2))

        c(1,i) = (q(ivx,i) - cc) / fc
        c(2,i) = (q(ivx,i) + cc) / fc

      end do

    end if

!-------------------------------------------------------------------------------
!
  end subroutine fluxspeed_srhd_adi
!
!===============================================================================
!
! function MAXSPEED_SRHD_ADI:
! --------------------------
!
!   Function scans the variable array and returns the maximum speed in within.
!
!   Arguments:
!
!     q - the array of primitive variables;
!
!===============================================================================
!
  function maxspeed_srhd_adi(qq) result(maxspeed)

    use coordinates, only : nb, ne

    implicit none

    real(kind=8), dimension(:,:,:,:), intent(in) :: qq

    real(kind=8) :: maxspeed

    integer      :: i, j, k
    real(kind=8) :: vv, v, cc, ww, c2, ss, fc

!-------------------------------------------------------------------------------
!
    maxspeed = 0.0d+00

#if NDIMS == 2
    k = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
    do k = nb, ne
#endif /* NDIMS == 3 */
      do j = nb, ne
        do i = nb, ne

! calculate the velocity amplitude
!
          vv = sum(qq(ivx:ivz,i,j,k) * qq(ivx:ivz,i,j,k))
          v  = sqrt(vv)

! calculate the square of the sound speed
!
          ww = qq(idn,i,j,k) + qq(ipr,i,j,k) / gammaxi
          c2 = adiabatic_index * qq(ipr,i,j,k) / ww
          ss = c2 * (1.0d+00 - vv) / (1.0d+00 - c2)
          fc = 1.0d+00 + ss
          cc = sqrt(ss * (fc - vv))

! calculate the maximum of speed
!
          maxspeed = max(maxspeed, (v + cc) / fc)

        end do ! i = nb, ne
      end do ! j = nb, ne
#if NDIMS == 3
    end do ! k = nb, ne
#endif /* NDIMS == 3 */

    return

!-------------------------------------------------------------------------------
!
  end function maxspeed_srhd_adi
!
!===============================================================================
!
! subroutine GET_MAXIMUM_SPEEDS_SRHD_ADI:
! --------------------------------------
!
!   Subroutine determines the maximum characteristic speed and eigenvalue
!   in the input array of primitive variables.
!
!   Arguments:
!
!     qq - the input array of primitive variables;
!     vm - the maximum physical speed;
!     cm - the maximum eigenvalue;
!
!===============================================================================
!
  subroutine get_maximum_speeds_srhd_adi(qq, vm, cm)

    use coordinates, only : nb, ne

    implicit none

    real(kind=8), dimension(:,:,:,:), intent(in)  :: qq
    real(kind=8)                    , intent(out) :: vm, cm

    integer      :: i, j, k
    real(kind=8) :: vl, vu, vv, ww, aa, cc, ss, fc

!-------------------------------------------------------------------------------
!
    vm = 0.0d+00
    cm = 0.0d+00

#if NDIMS == 2
    k = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
    do k = nb, ne
#endif /* NDIMS == 3 */
      do j = nb, ne
        do i = nb, ne
          vl = minval(qq(ivx:ivz,i,j,k))
          vu = maxval(qq(ivx:ivz,i,j,k))
          vm = max(vm, abs(vl), abs(vu))
          vv = sum(qq(ivx:ivz,i,j,k) * qq(ivx:ivz,i,j,k))

          ww = qq(idn,i,j,k) + qq(ipr,i,j,k) / gammaxi
          aa = adiabatic_index * qq(ipr,i,j,k) / ww
          ss = aa * (1.0d+00 - vv) / (1.0d+00 - aa)
          fc = 1.0d+00 + ss
          cc = sqrt(ss * (fc - vv))
          cm = max(cm, abs(vl - cc), abs(vu + cc))
        end do
      end do
#if NDIMS == 3
    end do
#endif /* NDIMS == 3 */

!-------------------------------------------------------------------------------
!
  end subroutine get_maximum_speeds_srhd_adi
!
!===============================================================================
!
! subroutine NR_FUNCTION_SRHD_ADI_1D:
! ----------------------------------
!
!   Subroutine calculate the value of function
!
!     F(W) = W - P(W) - E
!
!   for a given enthalpy W. It is used to estimate the initial guess.
!
!   The pressure is
!
!     P(W) = (γ - 1)/γ (W - D / sqrt(1 - |v|²(W))) (1 - |v|²(W))
!
!   and the squared velocity is
!
!     |v|²(W) = |m|² / W²
!
!   Optional derivative is returned
!
!     dF(W)/dW = 1 - dP(W)/dW
!
!   Arguments:
!
!     mm, en, dn, w - input coefficients for |M|² E, D, and W, respectively;
!     f             - the value of function F(W);
!     df            - optional derivative F'(W);
!
!===============================================================================
!
  subroutine nr_function_srhd_adi_1d(mm, en, dn, w, f, df)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8)          , intent(in)    :: mm, en, dn, w
    real(kind=8)          , intent(out)   :: f
    real(kind=8), optional, intent(out)   :: df

! local variables
!
    real(kind=8) :: vv, vm, vs
!
!-------------------------------------------------------------------------------
!
    vv = mm / (w * w)
    vm = 1.0d+00 - vv
    vs = sqrt(vm)
    f  = (1.0d+00 - gammaxi * vm) * w + gammaxi * dn * vs - en
    if (present(df)) then
      df  = 1.0d+00 - gammaxi * (1.0d+00 + (1.0d+00 - dn / (vs * w)) * vv)
    end if

!-------------------------------------------------------------------------------
!
  end subroutine nr_function_srhd_adi_1d
!
!===============================================================================
!
! subroutine NR_ITERATE_SRHD_ADI_1DW:
! ----------------------------------
!
!   Subroutine finds a root W of equation
!
!     F(W) = W - P(W) - E = 0
!
!   using the Newton-Raphson 1Dw iterative method.
!
!   Arguments:
!
!     mm, en - input coefficients for |M|² and E, respectively;
!     bb, bm - input coefficients for |B|² and B.M, respectively;
!     w, vv  - input/output coefficients W and |v|²;
!     info   - the flag is .true. if the solution was found, otherwise
!              it is .false.;
!
!   References:
!
!     Noble, S. C., Gammie, C. F., McKinney, J. C, Del Zanna, L.,
!     "Primitive Variable Solvers for Conservative General Relativistic
!      Magnetohydrodynamics",
!     The Astrophysical Journal, 2006, vol. 641, pp. 626-637
!
!===============================================================================
!
  subroutine nr_iterate_srhd_adi_1dw(mm, bb, mb, en, dn, w, vv, info)

    use helpers, only : print_message

    implicit none

    real(kind=8), intent(in)    :: mm, bb, mb, en, dn
    real(kind=8), intent(inout) :: w, vv
    logical     , intent(out)   :: info

    logical      :: keep
    integer      :: it, cn
    real(kind=8) :: wl, wu, fl, fu
    real(kind=8) :: f, df, dw
    real(kind=8) :: err

    character(len=80) :: msg

    character(len=*), parameter :: loc = 'EQUATIONS::nr_iterate_srhd_adi_1dw()'

!-------------------------------------------------------------------------------
!
! check if the state is physical; this can save some time if unphysical state
! is considered
!
    wl = sqrt(mm + dn * dn)
    if (en > wl) then

! prepare the initial brackets
!
      wl = wl + gammaxi * pmin
      wu = en + pmin

! calculate the value of function for the lower bracket
!
      call nr_function_srhd_adi_1d(mm, en, dn, wl, fl)

! the lower bracket gives negative function, so there is chance it bounds
! the root
!
      if (fl < 0.0d+00) then

! make sure that the upper bracket is larger than the lower one and
! the function has positive value
!
        if (wu <= wl) wu = 2.0d+00 * wl

! check if the brackets bound the root region, if not proceed until
! opposite function signs are found for the brackets
!
        call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
        it   = nrmax
        keep = fl * fu > 0.0d+00
        do while (keep)
          it = it - 1
          wl = wu
          fl = fu
          wu = 2.0d+00 * wu
          call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
          keep = (fl * fu > 0.0d+00) .and. it > 0
        end do

! the upper bracket was found, so proceed with determining the root
!
        if (it > 0 .and. fu >= 0.0d+00) then

! estimate the value of enthalpy close to the root and corresponding v²
!
          w = wl - fl * (wu - wl) / (fu - fl)

! initialize iteration parameters
!
          info = .true.
          keep = .true.
          it   = nrmax
          cn   = nrext

! iterate using the Newton-Raphson method in order to find a root w of the
! function
!
          do while(keep)

! calculate F(W) and dF(W)/dW
!
            call nr_function_srhd_adi_1d(mm, en, dn, w, f, df)

! update brackets
!
            if (f > fl .and. f < 0.0d+00) then
              wl = w
              fl = f
            end if
            if (f < fu .and. f > 0.0d+00) then
              wu = w
              fu = f
            end if

! calculate the increment dW, update the solution, and estimate the error
!
            dw  = f / df
            w   = w - dw
            err = abs(dw / w)

! check the convergence, if the convergence is not reached, iterate until
! the maximum number of iteration is reached
!
            if (err < tol) then
              keep = cn > 0
              cn   = cn - 1
            else
              keep = it > 0
            end if

! if new W leaves the brackets, use the bisection method to estimate
! the new guess
!
            if (w < wl .or. w > wu) then
              w = 0.5d+00 * (wl + wu)
            end if

            it = it - 1
          end do ! NR iterations

! print information about failed convergence or unphysical variables
!
          if (err >= tol) then
            call print_message(loc, "Convergence not reached!")
            write(msg,"(a,1x,1es24.16e3)") "Error: ", err
            call print_message(loc, msg)
          end if

! calculate |V|² from W
!
          vv  = mm / (w * w)

        else ! the upper brack not found
          call print_message(loc, "Could not find the upper bracket!")
          info = .false.

        end if

      else if (fl > 0.0d+00) then ! the root cannot be found, since it is below the lower bracket
        call print_message(loc, "Positive function for lower bracket!")
        info = .false.
      else ! the lower bracket is a root, so return it
        w    = wl
        info = .true.
      end if

    else ! the state is unphysical
      call print_message(loc, "The state is not physical!")
      info = .false.
    end if

!-------------------------------------------------------------------------------
!
  end subroutine nr_iterate_srhd_adi_1dw
!
!===============================================================================
!
! subroutine NR_ITERATE_SRHD_ADI_2DWV:
! -----------------------------------
!
!   Subroutine finds a root (W,v²) of 2D equations
!
!     F(W,v²) = W - E - P(W,v²) = 0
!     G(W,v²) = W² v² - m²      = 0
!
!   using the Newton-Raphson 2D iterative method.
!
!   All evaluated equations incorporate already the pressure of the form
!
!     P(W,|v|²) = (γ - 1)/γ (W - Γ D) (1 - |v|²)
!
!   in order to optimize calculations.
!
!   Arguments:
!
!     mm, en - input coefficients for |M|² and E, respectively;
!     bb, bm - input coefficients for |B|² and B.M, respectively;
!     w, vv  - input/output coefficients W and |v|²;
!     info   - the flag is .true. if the solution was found, otherwise
!              it is .false.;
!
!   References:
!
!     Noble, S. C., Gammie, C. F., McKinney, J. C, Del Zanna, L.,
!     "Primitive Variable Solvers for Conservative General Relativistic
!      Magnetohydrodynamics",
!     The Astrophysical Journal, 2006, vol. 641, pp. 626-637
!
!===============================================================================
!
  subroutine nr_iterate_srhd_adi_2dwv(mm, bb, mb, en, dn, w, vv, info)

    use helpers, only : print_message

    implicit none

    real(kind=8), intent(in)    :: mm, bb, mb, en, dn
    real(kind=8), intent(inout) :: w, vv
    logical     , intent(out)   :: info

    logical      :: keep
    integer      :: it, cn
    real(kind=8) :: wl, wu, fl, fu
    real(kind=8) :: ww, vm, vs, gd, gv
    real(kind=8) :: f, dfw, dfv
    real(kind=8) :: g, dgw, dgv
    real(kind=8) :: det, jfw, jfv, jgw, jgv
    real(kind=8) :: dw, dv
    real(kind=8) :: err

    character(len=80) :: msg

    character(len=*), parameter :: loc = 'EQUATIONS::nr_iterate_srhd_adi_2dwv()'

!-------------------------------------------------------------------------------
!
! prepare the initial brackets
!
    wl   = sqrt(mm + dn * dn) + gammaxi * pmin
    wu   = en + pmin

! make sure that the upper bracket is larger than the lower one
!
    keep = wl >= wu
    it   = nrmax
    do while(keep)
      wu   = 2.0d+00 * wu
      it   = it - 1
      keep = (wl >= wu) .and. it > 0
    end do
    if (it <= 0) then
      call print_message(loc, "Could not find the upper limit for enthalpy!")
      info = .false.
      return
    end if

! check if the brackets bound the root region, if not proceed until
! opposite function signs are found for the brackets
!
    call nr_function_srhd_adi_1d(mm, en, dn, wl, fl)
    call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)

    keep = (fl * fu > 0.0d+00)
    it   = nrmax

    do while (keep)

      wl = wu
      fl = fu
      wu = 2.0d+00 * wu

      call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)

      it   = it - 1

      keep = (fl * fu > 0.0d+00) .and. it > 0

    end do
    if (it <= 0) then
      call print_message(loc, "No initial brackets found!")
      info = .false.
      return
    end if

! estimate the value of enthalpy close to the root and corresponding v²
!
    w  = wl - fl * (wu - wl) / (fu - fl)
    vv = mm / (w * w)

! initialize iteration parameters
!
    info = .true.
    keep = .true.
    it   = nrmax
    cn   = nrext

! iterate using the Newton-Raphson method in order to find the roots W and |v|²
! of functions
!
    do while(keep)

! calculate W², (1 - |v|²), and the Lorentz factor
!
      ww  = w * w
      vm  = 1.0d+00 - vv
      vs  = sqrt(vm)
      gd  = gammaxi * dn
      gv  = 1.0d+00 - gammaxi * vm

! calculate F(W,|v|²) and G(W,|v|²)
!
      f   = gv * w - en + gd * vs
      g   = vv * ww - mm

! calculate dF(W,|v|²)/dW and dF(W,|v|²)/d|v|²
!
      dfw = gv
      dfv = gammaxi * w - 0.5d+00 * gd / vs

! calculate dG(W,|v|²)/dW and dG(W,|v|²)/d|v|²
!
      dgw = 2.0d+00 * vv * w
      dgv = ww

! invert the Jacobian J = | dF/dW, dF/d|v|² |
!                         | dG/dW, dG/d|v|² |
!
      det = dfw * dgv - dfv * dgw

      jfw =   dgv / det
      jgw = - dfv / det
      jfv = - dgw / det
      jgv =   dfw / det

! calculate increments dW and d|v|²
!
      dw  = f * jfw + g * jgw
      dv  = f * jfv + g * jgv

! correct W and |v|²
!
      w   = w  - dw
      vv  = vv - dv

! check if the new enthalpy and velocity are physical
!
      if (w < wl) then
        write(msg,"(a,1x,2es24.16e3)") "Enthalpy smaller than the limit: ",    &
                                       w, wl
        call print_message(loc, msg)
        info = .false.
        return
      end if
      if (vv < 0.0d+00 .or. vv >= 1.0d+00) then
        write(msg,"(a,1x,1es24.16e3)") "Unphysical speed |v|²: ", vv
        call print_message(loc, msg)
        info = .false.
        return
      end if

! calculate the error
!
      err = max(abs(dw / w), abs(dv))

! check the convergence, if the convergence is not reached, iterate until
! the maximum number of iteration is reached
!
      if (err < tol) then
        keep = cn > 0
        cn   = cn - 1
      else
        keep = it > 0
      end if

! decrease the number of remaining iterations
!
      it = it - 1

    end do ! NR iterations

! print information about failed convergence or unphysical variables
!
    if (err >= tol) then
      write(msg,"(a,1x,1es24.16e3)") "Convergence not reached: ", err
      call print_message(loc, msg)
    end if

!-------------------------------------------------------------------------------
!
  end subroutine nr_iterate_srhd_adi_2dwv
!
!===============================================================================
!
! subroutine NR_ITERATE_SRHD_ADI_2DWU:
! -----------------------------------
!
!   Subroutine finds a root (W,u²) of 2D equations
!
!     F(W,u²) = (W - E - P(W,u²)) (u² + 1) = 0
!     G(W,u²) =  W² u² - (u² + 1) m²       = 0
!
!   using the Newton-Raphson 2D iterative method.
!
!   All evaluated equations incorporate already the pressure of the form
!
!     P(W,|u|²) = (γ - 1)/γ (W - Γ D) / (1 + |u|²)
!
!   in order to optimize calculations.
!
!   Arguments:
!
!     mm, en - input coefficients for |M|² and E, respectively;
!     bb, bm - input coefficients for |B|² and B.M, respectively;
!     w, vv  - input/output coefficients W and |v|²;
!     info   - the flag is .true. if the solution was found, otherwise
!              it is .false.;
!
!===============================================================================
!
  subroutine nr_iterate_srhd_adi_2dwu(mm, bb, mb, en, dn, w, vv, info)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8), intent(in)    :: mm, bb, mb, en, dn
    real(kind=8), intent(inout) :: w, vv
    logical     , intent(out)   :: info

! local variables
!
    logical      :: keep
    integer      :: it, cn
    real(kind=8) :: wl, wu, fl, fu
    real(kind=8) :: ww, uu, up, gm, gd
    real(kind=8) :: f, dfw, dfu
    real(kind=8) :: g, dgw, dgu
    real(kind=8) :: det, jfw, jfu, jgw, jgu
    real(kind=8) :: dw, du
    real(kind=8) :: err
!
!-------------------------------------------------------------------------------
!
! prepare the initial brackets
!
    wl   = sqrt(mm + dn * dn) + gammaxi * pmin
    wu   = en + pmin

! make sure that the upper bracket is larger than the lower one
!
    keep = wl >= wu
    it   = nrmax
    do while(keep)
      wu   = 2.0d+00 * wu
      it   = it - 1
      keep = (wl >= wu) .and. it > 0
    end do
    if (it <= 0) then
      write(*,*)
      write(*,"(a,1x,a)") "ERROR in"                                           &
                        , "EQUATIONS::nr_iterate_srhd_adi_1dw()"
      write(*,"(a)"     ) "Could not find the upper limit for enthalpy!"
      info = .false.
      return
    end if

! check if the brackets bound the root region, if not proceed until
! opposite function signs are found for the brackets
!
    call nr_function_srhd_adi_1d(mm, en, dn, wl, fl)
    call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)

    keep = (fl * fu > 0.0d+00)
    it   = nrmax

    do while (keep)

      wl = wu
      fl = fu
      wu = 2.0d+00 * wu

      call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)

      it   = it - 1

      keep = (fl * fu > 0.0d+00) .and. it > 0

    end do
    if (it <= 0) then
      write(*,*)
      write(*,"(a,1x,a)") "ERROR in"                                           &
                        , "EQUATIONS::nr_iterate_srhd_adi_2dwu()"
      write(*,"(a)"     ) "No initial brackets found!"
      info = .false.
      return
    end if

! estimate the value of enthalpy close to the root and corresponding u²
!
    w  = wl - fl * (wu - wl) / (fu - fl)
    uu = mm / (w * w - mm)

! initialize iteration parameters
!
    info = .true.
    keep = .true.
    it   = nrmax
    cn   = nrext

! iterate using the Newton-Raphson method in order to find the roots W and |u|²
! of functions
!
    do while(keep)

! calculate W², (1 + |u|²), and the Lorentz factor, and some repeated
! expressions
!
      ww  = w * w
      up  = 1.0d+00 + uu
      gm  = sqrt(up)
      gd  = gammaxi * dn

! calculate F(W,|u|²) and G(W,|u|²)
!
      f   = (up - gammaxi) * w - up * en + gm * gd
      g   = uu * ww - up * mm

! calculate dF(W,|u|²)/dW and dF(W,|u|²)/d|u|²
!
      dfw = up - gammaxi
      dfu = w - en + 0.5d+00 * gd / gm

! calculate dG(W,|u|²)/dW and dG(W,|u|²)/d|u|²
!
      dgw = 2.0d+00 * uu * w
      dgu = ww - mm

! invert the Jacobian J = | dF/dW, dF/d|u|² |
!                         | dG/dW, dG/d|u|² |
!
      det = dfw * dgu - dfu * dgw

      jfw =   dgu / det
      jgw = - dfu / det
      jfu = - dgw / det
      jgu =   dfw / det

! calculate increments dW and d|u|²
!
      dw  = f * jfw + g * jgw
      du  = f * jfu + g * jgu

! correct W and |u|²
!
      w   = w  - dw
      uu  = uu - du

! check if the new enthalpy gives physical pressure and velocity
!
      if (w < wl) then
        write(*,*)
        write(*,"(a,1x,a)"         ) "ERROR in"                                &
                                   , "EQUATIONS::nr_iterate_srhd_adi_2dwu()"
        write(*,"(a,1x,2es24.16e3)") "Enthalpy smaller than the limit: ", w, wl
        info = .false.
        return
      end if
      if (uu < 0.0d+00) then
        write(*,*)
        write(*,"(a,1x,a)"         ) "ERROR in"                                &
                                   , "EQUATIONS::nr_iterate_srhd_adi_2dwu()"
        write(*,"(a,1x,1es24.16e3)") "Unphysical speed |u|²: ", uu
        info = .false.
        return
      end if

! calculate the error
!
      err = max(abs(dw / w), abs(du))

! check the convergence, if the convergence is not reached, iterate until
! the maximum number of iteration is reached
!
      if (err < tol) then
        keep = cn > 0
        cn   = cn - 1
      else
        keep = it > 0
      end if

! decrease the number of remaining iterations
!
      it = it - 1

    end do ! NR iterations

! calculate |v|² from |u|²
!
    vv = uu / (1.0d+00 + uu)

! print information about failed convergence or unphysical variables
!
    if (err >= tol) then
      write(*,*)
      write(*,"(a,1x,a)"         ) "WARNING in"                                &
                                 , "EQUATIONS::nr_iterate_srhd_adi_2dwu()"
      write(*,"(a,1x,1es24.16e3)") "Convergence not reached: ", err
    end if

!-------------------------------------------------------------------------------
!
  end subroutine nr_iterate_srhd_adi_2dwu
!
!*******************************************************************************
!
! ADIABATIC SPECIAL RELATIVITY MAGNETOHYDRODYNAMIC EQUATIONS
!
!*******************************************************************************
!
!===============================================================================
!
! subroutine PRIM2CONS_SRMHD_ADI:
! ------------------------------
!
!   Subroutine converts primitive variables to their corresponding
!   conservative representation.
!
!   Arguments:
!
!     q - the input array of primitive variables;
!     u - the output array of conservative variables;
!     s - an optional flag indicating that passive scalars have
!         to be calculated too;
!
!===============================================================================
!
  subroutine prim2cons_srmhd_adi(q, u, s)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8), dimension(:,:), intent(in)  :: q
    real(kind=8), dimension(:,:), intent(out) :: u
    logical     , optional      , intent(in)  :: s

! local variables
!
    integer      :: i, p
    real(kind=8) :: vv, bb, vb, vm, vs, ww, wt
!
!-------------------------------------------------------------------------------
!
! iterate over all positions
!
    do i = 1, size(q,2)

! calculate the square of velocity, the quare of magnetic field, the scalar
! product of velocity and magnetic field, the Lorentz factor, specific and
! total enthalpies
!
      vv  = sum(q(ivx:ivz,i) * q(ivx:ivz,i))
      bb  = sum(q(ibx:ibz,i) * q(ibx:ibz,i))
      vb  = sum(q(ivx:ivz,i) * q(ibx:ibz,i))
      vm  = 1.0d+00 - vv
      vs  = sqrt(vm)
      ww  = (q(idn,i) + q(ipr,i) / gammaxi) / vm
      wt  = ww + bb

! calculate conservative variables
!
      u(idn,i) =      q(idn,i) / vs
      u(imx,i) = wt * q(ivx,i) - vb * q(ibx,i)
      u(imy,i) = wt * q(ivy,i) - vb * q(iby,i)
      u(imz,i) = wt * q(ivz,i) - vb * q(ibz,i)
      u(ibx,i) = q(ibx,i)
      u(iby,i) = q(iby,i)
      u(ibz,i) = q(ibz,i)
      u(ibp,i) = q(ibp,i)
      u(ien,i) = wt - q(ipr,i) - u(idn,i) - 0.5d+00 * (vm * bb + vb * vb)

    end do

! update primitive passive scalars
!
    if (ns > 0 .and. present(s)) then
      if (s) then
        do p = isl, isu
          u(p,:) = q(p,:) * u(idn,:)
        end do
      end if
    end if

!-------------------------------------------------------------------------------
!
  end subroutine prim2cons_srmhd_adi
!
!===============================================================================
!
! subroutine CONS2PRIM_SRMHD_ADI:
! ------------------------------
!
!   Subroutine converts conservative variables to their corresponding
!   primitive representation using an interative method.
!
!   Arguments:
!
!     u - the input array of conservative variables;
!     q - the output array of primitive variables;
!     s - the status flag;
!
!===============================================================================
!
  subroutine cons2prim_srmhd_adi(u, q, s)

    use helpers, only : print_message

    implicit none

    real(kind=8), dimension(:,:), intent(in)  :: u
    real(kind=8), dimension(:,:), intent(out) :: q
    integer                     , intent(out) :: s

    logical      :: info
    integer      :: i, p
    real(kind=8) :: mm, mb, bb, en, dn
    real(kind=8) :: w, wt, vv, vm, vs, fc

    character(len=80) :: msg

    character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_srmhd_adi()'

!-------------------------------------------------------------------------------
!
    s = 0

! iterate over all positions
!
    do i = 1, size(u,2)

! prepare variables which do not change during the Newton-Ralphson iterations
! (|B|², |M|² and B.M and their multiplications)
!
      mm  = sum(u(imx:imz,i) * u(imx:imz,i))
      mb  = sum(u(imx:imz,i) * u(ibx:ibz,i))
      bb  = sum(u(ibx:ibz,i) * u(ibx:ibz,i))
      en  = u(ien,i) + u(idn,i)
      dn  = u(idn,i)

! find the exact W using an Newton-Ralphson interative method
!
      call nr_iterate(mm, bb, mb, en, dn, w, vv, info)

! if info is .true., the solution was found
!
      if (info) then

! prepare coefficients
!
        vm = 1.0d+00 - vv
        vs = sqrt(vm)
        wt = w + bb
        fc = mb / w

! calculate the primitive variables
!
        q(idn,i) = dn * vs
        q(ivx,i) = (u(imx,i) + fc * u(ibx,i)) / wt
        q(ivy,i) = (u(imy,i) + fc * u(iby,i)) / wt
        q(ivz,i) = (u(imz,i) + fc * u(ibz,i)) / wt
        q(ibx,i) = u(ibx,i)
        q(iby,i) = u(iby,i)
        q(ibz,i) = u(ibz,i)
        q(ibp,i) = u(ibp,i)
        q(ipr,i) = w - en + 0.5d+00 * (bb + (bb * mm - mb * mb) / wt**2)

! check if the pressure is positive, if not, print a warning and replace it
! with the minimum allowed value pmin
!
        if (q(ipr,i) <= 0.0d+00) then

          call print_message(loc, "Conversion to physical primitive state " // &
                                  "resulted in negative pressure!")
          write(msg,"(a,9(1x,1es24.16e3))") "U = ", u(:,i)
          call print_message(loc, msg)
          write(msg,"(a,6(1x,1es24.16e3))") "D, |m|², m.B, |B|², E, W = ",     &
                                            dn, mm, mb, bb, en, w
          call print_message(loc, msg)

          s = 1
          go to 100

        end if ! p <= 0

      else ! unphysical state

        call print_message(loc, "Conversion to physical primitive state failed!")
        write(msg,"(a,9(1x,1es24.16e3))") "U = ", u(:,i)
        call print_message(loc, msg)
        write(msg,"(a,5(1x,1es24.16e3))") "D, |m|², m.B, |B|², E = ",          &
                                          dn, mm, mb, bb, en
        call print_message(loc, msg)

        s = 1
        go to 100

      end if ! unphysical state

    end do

    if (ns > 0) then
      do p = isl, isu
        q(p,:) = u(p,:) / u(idn,:)
      end do
    end if

    100 continue

!-------------------------------------------------------------------------------
!
  end subroutine cons2prim_srmhd_adi
!
!===============================================================================
!
! subroutine FLUXSPEED_SRMHD_ADI:
! ------------------------------
!
!   Subroutine calculates physical fluxes and characteristic speeds from a
!   given equation system.
!
!   Arguments:
!
!     q - the input array of primitive variables;
!     u - the input array of conservative variables;
!     f - the output vector of fluxes;
!     c - the output vector of left- and right-going characteristic speeds;
!
!   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
!     [2] van der Holst, B., Keppens, R., Meliani, Z.
!         "A multidimentional grid-adaptive relativistic magnetofluid code",
!         Computer Physics Communications, 2008, Volume 179, Pages 617-627
!
!===============================================================================
!
  subroutine fluxspeed_srmhd_adi(q, u, f, c)

    use algebra, only : quadratic, quartic
    use helpers, only : print_message

    implicit none

    real(kind=8), dimension(:,:)          , intent(in)  :: q, u
    real(kind=8), dimension(:,:)          , intent(out) :: f
    real(kind=8), dimension(:,:), optional, intent(out) :: c

    integer      :: i, nr
    real(kind=8) :: bb, vs
    real(kind=8) :: bx, by, bz, pm, pt
    real(kind=8) :: rh, v1, v2
    real(kind=8) :: ca, cc, c2, gn, rt, zm, zp
    real(kind=8) :: fa, fb, fc, fd, fe, ff, fg

    real(kind=8), dimension(size(q,2)) :: vv, vm, vb, b2

    real(kind=8), dimension(5) :: a
    real(kind=8), dimension(4) :: x

    character(len=80) :: msg

    character(len=*), parameter :: loc = 'EQUATIONS::fluxspeed_srmhd_adi()'
!
!-------------------------------------------------------------------------------
!
! iterate over all positions
!
    do i = 1, size(q,2)

! calculate the square of velocity, magnetic field and their scalar product
!
      vv(i) = sum(q(ivx:ivz,i) * q(ivx:ivz,i))
      bb    = sum(q(ibx:ibz,i) * q(ibx:ibz,i))
      vb(i) = sum(q(ivx:ivz,i) * q(ibx:ibz,i))

! calculate (1 - |V|²)
!
      vm(i) = 1.0d+00 - vv(i)

! calculate magnetic field components of the magnetic four-vector divided by
! the Lorentz factor (eq. 3 in [1])
!
      bx  = q(ibx,i) * vm(i) + vb(i) * q(ivx,i)
      by  = q(iby,i) * vm(i) + vb(i) * q(ivy,i)
      bz  = q(ibz,i) * vm(i) + vb(i) * q(ivz,i)

! calculate magnetic and total pressures (eq. 6 in [1])
!
      b2(i) = bb * vm(i) + vb(i) * vb(i)
      pm    = 0.5d+00 * b2(i)
      pt    = q(ipr,i) + pm

! calculate the relativistic hydrodynamic fluxes (eq. 13 in [1])
!
      f(idn,i) = u(idn,i) * q(ivx,i)
      f(imx,i) = u(imx,i) * q(ivx,i) - q(ibx,i) * bx + pt
      f(imy,i) = u(imy,i) * q(ivx,i) - q(ibx,i) * by
      f(imz,i) = u(imz,i) * q(ivx,i) - q(ibx,i) * bz
      f(ibx,i) = q(ibp,i)
      f(ibp,i) = cmax2 * q(ibx,i)
      f(iby,i) = q(ivx,i) * q(iby,i) - q(ibx,i) * q(ivy,i)
      f(ibz,i) = q(ivx,i) * q(ibz,i) - q(ibx,i) * q(ivz,i)
      f(ien,i) = u(imx,i) - f(idn,i)

    end do

    if (present(c)) then

! calculate the characteristic speeds
!
      do i = 1, size(q,2)

! check if the total velocity |V|² is larger than zero
!
        if (vv(i) > 0.0d+00) then

! calculate additional coefficients
!
          rh  = q(idn,i) + q(ipr,i) / gammaxi
          vs  = sqrt(vm(i))

! check if the normal component of magnetic field Bₓ is larger than zero
!
          if (abs(q(ibx,i)) > 0.0d+00) then ! Bₓ ≠ 0

! prepare parameters for this case
!
            c2 = adiabatic_index * q(ipr,i) / rh
            v1 = abs(q(ivx,i))
            v2 = v1 * v1

            fa = rh * (1.0d+00 - c2)
            fb = c2 * (rh - vb(i) * vb(i)) + b2(i)
            fc = sign(1.0d+00, q(ivx,i)) * q(ibx,i) * vs
            fd = c2 * fc * fc
            fe = 1.0d+00 - v2
            ff = c2 * vb(i) * fc
            fg = v1 * vs

! prepare polynomial coefficients
!
            a(5) = fa + fb * vm(i)
            a(4) = 2.0d+00 * (ff * vm(i) + fb * fg)
            a(3) = - fd * vm(i) + 4.0d+00 * ff * fg - fb * fe
            a(2) = - 2.0d+00 * (fd * fg + fe * ff)
            a(1) = fd * fe

! call the quartic solver
!
            nr = quartic(a(1:5), x(1:4))

! convert eigenvalues to charasteristic speeds
!
            x(1:nr) = sign(1.0d+00, q(ivx,i)) * (abs(v1) + x(1:nr) * vs)

          else ! Bₓ ≠ 0

! special case when Bₓ = 0, then the quartic equation reduces to quadratic one
!
! prepare parameters for this case
!
            c2 = adiabatic_index * q(ipr,i) / rh
            cc = (1.0d+00 - c2) / vm(i)
            gn = b2(i) - c2 * vb(i) * vb(i)

! prepare polynomial coefficients
!
            a(3) = rh * (c2 + cc) + gn
            a(2) = - 2.0d+00 * rh * cc * q(ivx,i)
            a(1) = rh * (cc * q(ivx,i)**2 - c2) - gn

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

          end if ! Bx ≠ 0

        else ! |V|² > 0

! special case when |V|² = 0 (Γ = 1), then the quartic equation reduces to
! bi-quartic one
!
! prepare parameters for this case
!
          rh = q(idn,i) + q(ipr,i) / gammaxi
          vs = sqrt(vm(i))
          rt = rh + b2(i)
          c2 = adiabatic_index * q(ipr,i) / rh
          ca = (q(ibx,i) * vs + vb(i) * q(ivx,i) / vs)**2

! prepare polynomial coefficients
!
          a(3) = 1.0d+00
          a(2) = - ((rh + ca) * c2 + b2(i)) / rt
          a(1) = c2 * ca / rt

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

! compute the roots
!
          if (nr > 0) then

            zm = min(x(1), x(2))
            zp = max(x(1), x(2))

            if (zm >= 0.0d+00) then

              zm   = sqrt(zm)
              zp   = sqrt(zp)

              x(1) =   zp
              x(2) =   zm
              x(3) = - zm
              x(4) = - zp

              nr   = 4

            else

              if (zp >= 0.0d+00) then

                zp   = sqrt(zp)

                x(1) =   zp
                x(2) = - zp

                nr = 2

              else

                x(:) = 0.0d+00

                nr   = 0

              end if
            end if
          end if

        end if ! |V|² > 0

! find the minimum and maximum characteristic speeds
!
        if (nr > 1) then

          c(1,i) = minval(x(1:nr))
          c(2,i) = maxval(x(1:nr))

          if (max(abs(c(1,i)), abs(c(2,i))) >= 1.0d+00) then
            call print_message(loc, "Estimation returned unphysical speeds!")
            write(msg,"(a,5(1es24.16))") "A = ", a(1:5)
            call print_message(loc, msg)
            write(msg,"(a,1i2)") "N = ", nr
            call print_message(loc, msg)
            write(msg,"(a,4(1es24.16))") "X = ", x(1:4)
            call print_message(loc, msg)
          end if

        else

! speed estimation failed, so we substitute the minimum and maximum physical
! speeds equal the speed of light
!
          c(1,i) = - 1.0d+00
          c(2,i) =   1.0d+00

        end if

      end do

    end if

!-------------------------------------------------------------------------------
!
  end subroutine fluxspeed_srmhd_adi
!
!===============================================================================
!
! function MAXSPEED_SRMHD_ADI:
! ---------------------------
!
!   Function scans the variable array and returns the maximum speed in within.
!
!   Arguments:
!
!     qq - the array of primitive variables;
!
!===============================================================================
!
  function maxspeed_srmhd_adi(qq) result(maxspeed)

! local variables are not implicit by default
!
    implicit none

! input arguments
!
    real(kind=8), dimension(:,:,:,:), intent(in) :: qq

! return value
!
    real(kind=8) :: maxspeed
!
!-------------------------------------------------------------------------------
!
    maxspeed = 1.0d+00

    return

!-------------------------------------------------------------------------------
!
  end function maxspeed_srmhd_adi
!
!===============================================================================
!
! subroutine GET_MAXIMUM_SPEEDS_SRMHD_ADI:
! ---------------------------------------
!
!   Subroutine determines the maximum characteristic speed and eigenvalue
!   in the input array of primitive variables.
!
!   Arguments:
!
!     qq - the input array of primitive variables;
!     vm - the maximum physical speed;
!     cm - the maximum eigenvalue;
!
!===============================================================================
!
  subroutine get_maximum_speeds_srmhd_adi(qq, vm, cm)

    use coordinates, only : nb, ne

    implicit none

    real(kind=8), dimension(:,:,:,:), intent(in)  :: qq
    real(kind=8)                    , intent(out) :: vm, cm

    integer      :: i, j, k
    real(kind=8) :: vl, vu

!-------------------------------------------------------------------------------
!
    vm = 0.0d+00
    cm = 1.0d+00

#if NDIMS == 2
    k = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
    do k = nb, ne
#endif /* NDIMS == 3 */
      do j = nb, ne
        do i = nb, ne
          vl = minval(qq(ivx:ivz,i,j,k))
          vu = maxval(qq(ivx:ivz,i,j,k))
          vm = max(vm, abs(vl), abs(vu))
        end do
      end do
#if NDIMS == 3
    end do
#endif /* NDIMS == 3 */

!-------------------------------------------------------------------------------
!
  end subroutine get_maximum_speeds_srmhd_adi
!
!===============================================================================
!
! subroutine NR_VELOCITY_SRMHD_ADI_1D:
! -----------------------------------
!
!   Subroutine calculates the squared velocity
!
!     |V|²(W)  = (|m|² W² + S² (2 W + |B|²)) / (W² (W² + |B|²)²)
!
!   and its derivative
!
!     |V|²'(W) = - 2 (|m|² W³ + 3 S² W² + 3 S² |B|² W + S² |B|⁴)
!                                                           / (W³ (W² + |B|²)³)
!
!   for a given enthalpy W.
!
!   Arguments:
!
!     mm, bb, mb, w - input coefficients for |M|², |B|², M.B, and W,
!                     respectively;
!     vv, dv        - the values of squared velocity |V|² and its derivative;
!
!===============================================================================
!
  subroutine nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv, dv)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8),           intent(in)  :: mm, bb, mb, w
    real(kind=8),           intent(out) :: vv
    real(kind=8), optional, intent(out) :: dv

! local variables
!
    real(kind=8) :: ss, ww, www, wt, wt2
!
!-------------------------------------------------------------------------------
!
! temporary variables
!
    ss  = mb * mb
    ww  = w * w
    www = ww * w
    wt  = w + bb
    wt2 = wt * wt

! the function and its derivative
!
    vv = (mm * ww + (w + wt) * ss) / (ww * wt2)
    if (present(dv)) then
      dv = - 2.0d+00 * (mm * www + (3.0d+00 * ww                               &
                           + (2.0d+00 * w + wt) * bb) * ss) / (www * wt2 * wt)
    end if

!-------------------------------------------------------------------------------
!
  end subroutine nr_velocity_srmhd_adi_1d
!
!===============================================================================
!
! subroutine NR_PRESSURE_SRMHD_ADI_1D:
! -----------------------------------
!
!   Subroutine calculates the pressure function
!
!     P(W)  = W  - E + ½ |B|² + ½ (|m|² |B|² - S²) / (W + |B|²)²
!
!   and its derivative
!
!     P'(W) = 1 - (|m|² |B|² - S²) / (W + |B|²)³
!
!   for a given enthalpy W.
!
!   This subroutine is used to find the minimum enthalpy for which the velocity
!   is physical and the pressure is positive.
!
!   Arguments:
!
!     mm, bb, mb, en, w - input coefficients for |M|², |B|², M.B, E,
!                         and W, respectively;
!     p, dp             - the values for the function P(W) and its
!                         derivative P'(W);
!
!===============================================================================
!
  subroutine nr_pressure_srmhd_adi_1d(mm, bb, mb, en, w, p, dp)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8)          , intent(in)  :: mm, bb, mb, en, w
    real(kind=8)          , intent(out) :: p
    real(kind=8), optional, intent(out) :: dp

! local variables
!
    real(kind=8) :: wt, wd, ss, fn
!
!-------------------------------------------------------------------------------
!
! temporary variables
!
    wt = w  + bb
    wd = wt * wt
    ss = mb * mb
    fn = (mm * bb - ss) / wd

! the pressure function and its derivative
!
    p  = w - en + 0.5d+00 * (bb + fn)
    if (present(dp)) dp = 1.0d+00 - fn / wt

!-------------------------------------------------------------------------------
!
  end subroutine nr_pressure_srmhd_adi_1d
!
!===============================================================================
!
! subroutine NR_FUNCTION_SRMHD_ADI_1D:
! -----------------------------------
!
!   Subroutine calculates the energy function
!
!     F(W)  = W - P(W) + ½ |B|² + ½ (|m|² |B|² - S²) / (W + |B|²)² - E
!
!   and its derivative
!
!     F'(W) = 1 - dP(W)/dW - (|m|² |B|² - S²) / (W + |B|²)³
!
!   for a given enthalpy W. It is used to estimate the initial guess.
!
!   Arguments:
!
!     mm, bb, mb, en, dn, w - input coefficients for |M|², |B|², M.B, E, D,
!                             and W, respectively;
!     f, df                 - the values of F(W) and its derivative;
!
!===============================================================================
!
  subroutine nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, w, f, df)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8),           intent(in)  :: mm, bb, mb, en, dn, w
    real(kind=8),           intent(out) :: f
    real(kind=8), optional, intent(out) :: df

! local variables
!
    real(kind=8) :: pr, dp, gm2, gm, dg
    real(kind=8) :: ww, wt, wd, ss, sw, ws, fn, dv, ds
!
!-------------------------------------------------------------------------------
!
! temporary variables
!
    ww  = w * w
    wt  = w + bb
    wd  = wt * wt
    ss  = mb * mb
    sw  = ss / ww
    ws  = (w + wt) * sw
    fn  = (mm * bb - ss) / wd
    dv  = wd - mm - ws
    ds  = sqrt(dv)

! calculate the Lorentz factor
!
    gm2 = wd / dv
    gm  = wt / ds

! calculate the pressure P(W) and energy function F(W)
!
    pr  = gammaxi * (w - gm * dn) / gm2
    f   = w - pr - en + 0.5d+00 * (bb + fn)

! if desired, calculate the derivatives dP(W)/dW and dF(W)/dW
!
    if (present(df)) then

      dg  = (1.0d+00 - wt * (wt - sw + ws / w) / dv) / ds
      dp  = gammaxi * (1.0d+00 - (2.0d+00 * w / gm - dn) * dg) / gm2
      df  = 1.0d+00 - dp - fn / wt

    end if ! df present

!-------------------------------------------------------------------------------
!
  end subroutine nr_function_srmhd_adi_1d
!
!===============================================================================
!
! subroutine NR_INITIAL_BRACKETS_SRMHD_ADI:
! ----------------------------------------
!
!   Subroutine finds the initial brackets and initial guess from
!   the positivity condition
!
!     W³ + (5/2 |B|² - E) W² + 2 (|B|² - E) |B|² W
!                          + 1/2 [(|B|⁴ - 2 |B|² E + |m|²) |B|² - S²] > 0
!
!   coming from the energy equation and
!
!     W⁴ + 2 |B|² W³ - (|m|² + D² - |B|⁴) W²
!                    - (2 S² + D² |B|²) W - (S² + D² |B|²) |B|² > 0
!
!   coming from the equation of state
!
!   using analytical. It takes the maximum estimated root as the lower bracket.
!   If the analytical estimation fails, the Newton-Raphson iterative method
!   is used.
!
!   Arguments:
!
!     mm, en - input coefficients for |M|² and E, respectively;
!     bb, mb - input coefficients for |B|² and m.B, respectively;
!     wl, wu - the lower and upper limits for the enthalpy;
!     wc     - the initial root guess;
!     info   - the flag is .true. if the initial brackets and guess were found,
!              otherwise it is .false.;
!
!===============================================================================
!
  subroutine nr_initial_brackets_srmhd_adi(mm, bb, mb, en, dn                  &
                                                     , wl, wu, wc, fl, fu, info)

    use algebra, only : cubic_normalized, quartic
    use helpers, only : print_message

    implicit none

    real(kind=8), intent(in)  :: mm, bb, mb, en, dn
    real(kind=8), intent(out) :: wl, wu, wc, fl, fu
    logical     , intent(out) :: info

    logical      :: keep
    integer      :: it, nr
    real(kind=8) :: dd, ss, ec
    real(kind=8) :: f , df
    real(kind=8) :: dw, err

    real(kind=8), dimension(5) :: a
    real(kind=8), dimension(4) :: x

    character(len=80) :: msg

    character(len=*), parameter :: loc = 'EQUATIONS::nr_initial_brackets_srmhd_adi()'

!-------------------------------------------------------------------------------
!
! calculate temporary variables
!
    dd   = dn * dn
    ss   = mb * mb
    ec   = en + pmin

! set the initial upper bracket
!
    wu   = en + pmin

! calculate the cubic equation coefficients for the positivity condition
! coming from the energy equation; the condition, in fact, finds the minimum
! enthalphy for which the pressure is equal to pmin
!
    a(3) = 2.5d+00 * bb - ec
    a(2) = 2.0d+00 * (bb - ec) * bb
    a(1) = 0.5d+00 * (((bb - 2.0d+00 * ec) * bb + mm) * bb - ss)

! solve the cubic equation
!
    nr = cubic_normalized(a(1:3), x(1:3))

! if solution was found, use the maximum root as the lower bracket
!
    if (nr > 0) then

      wl = x(nr)

! calculate the quartic equation coefficients for the positivity condition
! coming from the pressure equation
!
      a(5) = 1.0d+00
      a(4) = 2.0d+00 * bb
      a(3) = bb * bb - dd - mm
      a(2) = - 2.0d+00 * (ss + dd * bb)
      a(1) = - (ss + dd * bb) * bb

! solve the quartic equation
!
      nr = quartic(a(1:5), x(1:4))

! take the maximum ethalpy from both conditions to guarantee that the pressure
! obtains from any of those equations is positive
!
      if (nr > 0) wl = max(wl, x(nr))

    else ! nr = 0

! the root could not be found analytically, so use the iterative solver
! to find the lower bracket; as the initial guess use the initial upper bracket
!
      keep = .true.
      it   = nrmax
      wl   = wu
      do while(keep)
        call nr_pressure_srmhd_adi_1d(mm, bb, mb, ec, wl, f, df)
        dw   = f / df
        wl   = wl - dw
        err  = abs(dw / wl)
        it   = it - 1
        keep = (err > tol) .and. it > 0
      end do
      if (it <= 0) then
        call print_message(loc,                                                &
                         "Iterative solver failed to find the lower bracket!")
        write(msg,"(a,5(1x,1es24.16e3))") "D, |m|², m.B, |B|², E = ",          &
                                          dn, mm, mb, bb, en
        call print_message(loc, msg)
      end if

    end if ! nr > 0

! check if the energy function is negative for the lower limit
!
    call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, wl, fl)

    if (fl <= 0.0d+00) then

! make sure that the upper limit is larger than the lower one
!
      if (wu <= wl) wu = 2.0d+00 * wl

! check if the brackets bound the root region, if not proceed until
! opposite function signs are found for the brackets
!
      call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, wu, fu)
      it   = nrmax
      keep = fl * fu > 0.0d+00
      do while (keep)
        it = it - 1
        wl = wu
        fl = fu
        wu = 2.0d+00 * wu
        call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, wu, fu)
        keep = (fl * fu > 0.0d+00) .and. it > 0
      end do

! the upper bracket was found, so proceed with determining the root
!
      if (it > 0 .and. fu >= 0.0d+00) then

! estimate the enthalpy value close to the root
!
        wc = wl - fl * (wu - wl) / (fu - fl)

! we have good brackets and guess, so good to go
!
        info = .true.

      else ! the upper brack not found
        call print_message(loc, "Could not find the upper bracket!")
        write(msg,"(a,5(1x,1es24.16e3))") "D, |m|², m.B, |B|², E = ",          &
                                          dn, mm, mb, bb, en
        call print_message(loc, msg)
        info = .false.

      end if

    else ! the root cannot be found, since it is below the lower bracket
        call print_message(loc, "Positive function for lower bracket!")
        write(msg,"(a,6(1x,1es24.16e3))") "D, |m|², m.B, |B|², E, W = ",       &
                                          dn, mm, mb, bb, en, wl
        call print_message(loc, msg)
      info = .false.
    end if

!-------------------------------------------------------------------------------
!
  end subroutine nr_initial_brackets_srmhd_adi
!
!===============================================================================
!
! subroutine NR_ITERATE_SRMHD_ADI_1DW:
! -----------------------------------
!
!   Subroutine finds a root W of equation
!
!     F(W) = W - P(W) + ½ [(1 + |V|²) |B|² - S² / W²] - E = 0
!
!   using the Newton-Raphson 1Dw iterative method.
!
!   Arguments:
!
!     mm, en - input coefficients for |M|² and E, respectively;
!     bb, bm - input coefficients for |B|² and B.M, respectively;
!     w , vv - input/output coefficients W and |V|²;
!     info   - the flag is .true. if the solution was found, otherwise
!              it is .false.;
!
!   References:
!
!     Noble, S. C., Gammie, C. F., McKinney, J. C, Del Zanna, L.,
!     "Primitive Variable Solvers for Conservative General Relativistic
!      Magnetohydrodynamics",
!     The Astrophysical Journal, 2006, vol. 641, pp. 626-637
!
!===============================================================================
!
  subroutine nr_iterate_srmhd_adi_1dw(mm, bb, mb, en, dn, w, vv, info)

    use helpers, only : print_message

    implicit none

    real(kind=8), intent(in)    :: mm, bb, mb, en, dn
    real(kind=8), intent(inout) :: w, vv
    logical     , intent(out)   :: info

    logical      :: keep
    integer      :: it, cn
    real(kind=8) :: wl, wu, fl, fu
    real(kind=8) :: f , df, dw
    real(kind=8) :: err

    character(len=80) :: msg

    character(len=*), parameter :: loc = 'EQUATIONS::nr_iterate_srmhd_adi_1dw()'

!-------------------------------------------------------------------------------
!
! find the initial brackets and estimate the initial enthalpy
!
    call nr_initial_brackets_srmhd_adi(mm, bb, mb, en, dn                      &
                                                 , wl, wu, w, fl, fu, info)

! continue if brackets found
!
    if (info) then

! initialize iteration parameters
!
      info = .true.
      keep = .true.
      it   = nrmax
      cn   = nrext

! iterate using the Newton-Raphson method in order to find a root w of the
! function
!
      do while(keep)

! calculate F(W) and dF(W)/dW
!
        call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, w, f, df)

! update brackets
!
        if (f > fl .and. f < 0.0d+00) then
          wl = w
          fl = f
        end if
        if (f < fu .and. f > 0.0d+00) then
          wu = w
          fu = f
        end if

! calculate the increment dW, update the solution, and estimate the error
!
        dw  = f / df
        w   = w - dw
        err = abs(dw / w)

! check the convergence, if the convergence is not reached, iterate until
! the maximum number of iteration is reached
!
        if (err < tol) then
          keep = cn > 0
          cn   = cn - 1
        else
          keep = it > 0
        end if

! if new W lays out of the brackets, use the bisection method to estimate
! the new guess
!
        if (w < wl .or. w > wu) then
          w = 0.5d+00 * (wl + wu)
        end if

! decrease the number of remaining iterations
!
        it = it - 1

      end do ! NR iterations

! let know the user if the convergence failed
!
      if (err >= tol) then
        call print_message(loc, "Convergence not reached!")
        write(msg,"(a,1x,1es24.16e3)") "Error: ", err
        call print_message(loc, msg)
      end if

! calculate |V|² from W
!
      call nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv)

    end if ! correct brackets

!-------------------------------------------------------------------------------
!
  end subroutine nr_iterate_srmhd_adi_1dw
!
!===============================================================================
!
! subroutine NR_ITERATE_SRMHD_ADI_2DWV:
! ------------------------------------
!
!   Subroutine finds a root (W, |V|²) of equations
!
!     F(W,|V|²) = W - P + ½ [(1 + |V|²) |B|² - S² / W²] - E     = 0
!     G(W,|V|²) = |V|² (|B|² + W)² - S² (|B|² + 2W) / W² - |M|² = 0
!
!   using the Newton-Raphson 2D iterative method.
!
!   Arguments:
!
!     mm, en - input coefficients for |M|² and E, respectively;
!     bb, bm - input coefficients for |B|² and B.M, respectively;
!     w , vv - input/output coefficients W and |V|²;
!     info   - the flag is .true. if the solution was found, otherwise
!              it is .false.;
!
!   References:
!
!     Noble, S. C., Gammie, C. F., McKinney, J. C, Del Zanna, L.,
!     "Primitive Variable Solvers for Conservative General Relativistic
!      Magnetohydrodynamics",
!     The Astrophysical Journal, 2006, vol. 641, pp. 626-637
!
!===============================================================================
!
  subroutine nr_iterate_srmhd_adi_2dwv(mm, bb, mb, en, dn, w, vv, info)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8), intent(in)    :: mm, bb, mb, en, dn
    real(kind=8), intent(inout) :: w, vv
    logical     , intent(out)   :: info

! local variables
!
    logical      :: keep
    integer      :: it, cn
    real(kind=8) :: wl, wu, fl, fu
    real(kind=8) :: vm, vs, wt, mw, wt2
    real(kind=8) :: f, dfw, dfv
    real(kind=8) :: g, dgw, dgv
    real(kind=8) :: det, jfw, jfv, jgw, jgv
    real(kind=8) :: dv, dw
    real(kind=8) :: err
!
!-------------------------------------------------------------------------------
!
! find the initial brackets and estimate the initial enthalpy
!
    call nr_initial_brackets_srmhd_adi(mm, bb, mb, en, dn                      &
                                                     , wl, wu, w, fl, fu, info)

! if the brackets could not be found, return the lower bracket as the solution
!
    if (.not. info) then
      write(*,*)
      write(*,"(a,1x,a)"         ) "WARNING in"                                &
                                 , "EQUATIONS::nr_iterate_srmhd_adi_1dw()"
      write(*,"(a,1x)"           ) "The solution lays in unphysical regime."
      write(*,"(a,1x,1es24.16e3)") "Using the lower bracket as solution: ", wl

! use the lower bracket, since it guarantees the positive pressure
!
      w = wl

! calculate |V|² from W
!
      call nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv)

      info = .true.
      return

    end if

! and the corresponding |V|²
!
    call nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv)

! initialize iteration parameters
!
    info  = .true.
    keep  = .true.
    it    = nrmax
    cn    = nrext

! find root with the help of the Newton-Raphson method
!
    do while(keep)

! calculate (S/W)², Wt, Wt²
!
      mw  = (mb / w)**2
      wt  = w + bb
      wt2 = wt * wt

! prepare (1 - |V|²) and sqrt(1 - |V|²)
!
      vm  = 1.0d+00 - vv
      vs  = sqrt(vm)

! calculate functions F(W,|V|²) and G(W,|V|²)
!
      f   = w - en - gammaxi * (w * vm - dn * vs)                              &
                                        + 0.5d+00 * ((1.0d+00 + vv) * bb - mw)
      g   = vv * wt2 - (wt + w) * mw  - mm

! calculate derivatives dF(W,|V|²)/dW and dF(W,|V|²)/d|V|²
!
      dfw = 1.0d+00 - gammaxi * vm + mw / w
      dfv = - gammaxi * (0.5d+00 * dn / vs - w) + 0.5d+00 * bb

! calculate derivatives dG(W,|V|²)/dW and dG(W,|V|²)/d|V|²
!
      dgw = 2.0d+00 * wt * (vv + mw / w)
      dgv = wt2

! invert the Jacobian J = | dF/dW, dF/d|V|² |
!                         | dG/dW, dG/d|V|² |
!
      det = dfw * dgv - dfv * dgw

      jfw =   dgv / det
      jgw = - dfv / det
      jfv = - dgw / det
      jgv =   dfw / det

! calculate increments dW and d|V|²
!
      dw  = f * jfw + g * jgw
      dv  = f * jfv + g * jgv

! correct W and |V|²
!
      w   = w  - dw
      vv  = vv - dv

! check if the new enthalpy and velocity are physical
!
      if (w < wl) then
        write(*,*)
        write(*,"(a,1x,a)"         ) "ERROR in"                                &
                                   , "EQUATIONS::nr_iterate_srmhd_adi_2dwv()"
        write(*,"(a,1x,2es24.16e3)") "Enthalpy smaller than the limit: ", w, wl
        info = .false.
        return
      end if
      if (vv < 0.0d+00 .or. vv >= 1.0d+00) then
        write(*,*)
        write(*,"(a,1x,a)"         ) "ERROR in"                                &
                                   , "EQUATIONS::nr_iterate_srmhd_adi_2dwv()"
        write(*,"(a,1x,1es24.16e3)") "Unphysical speed |v|²: ", vv
        info = .false.
        return
      end if

! calculate the error
!
      err = max(abs(dw / w), abs(dv))

! check the convergence, if the convergence is not reached, iterate until
! the maximum number of iteration is reached
!
      if (err < tol) then
        keep = cn > 0
        cn   = cn - 1
      else
        keep = it > 0
      end if

! decrease the number of remaining iterations
!
      it = it - 1

    end do ! NR iterations

! let know the user if the convergence failed
!
    if (err >= tol) then
      write(*,*)
      write(*,"(a,1x,a)"         ) "WARNING in"                                &
                                 , "EQUATIONS::nr_iterate_srmhd_adi_2dwv()"
      write(*,"(a,1x,1es24.16e3)") "Convergence not reached: ", err
    end if

!-------------------------------------------------------------------------------
!
  end subroutine nr_iterate_srmhd_adi_2dwv
!
!===============================================================================
!
! subroutine NR_ITERATE_SRMHD_ADI_2DWU:
! ------------------------------------
!
!   Subroutine finds a root (W, |u|²) of equations
!
!     F(W,|u|²) = W - E - P + ½ [(1 + |u|² / (1 + |u|²)) |B|² - S² / W²]     = 0
!     G(W,|u|²) = (|B|² + W)² |u|² / (1 + |u|²) - (2W + |B|²) S² / W² - |M|² = 0
!
!   using the Newton-Raphson 2D iterative method.
!
!   Arguments:
!
!     mm, en - input coefficients for |M|² and E, respectively;
!     bb, bm - input coefficients for |B|² and B.M, respectively;
!     w , vv - input/output coefficients W and |v|²;
!     info   - the flag is .true. if the solution was found, otherwise
!              it is .false.;
!
!===============================================================================
!
  subroutine nr_iterate_srmhd_adi_2dwu(mm, bb, mb, en, dn, w, vv, info)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8), intent(in)    :: mm, bb, mb, en, dn
    real(kind=8), intent(inout) :: w, vv
    logical     , intent(out)   :: info

! local variables
!
    logical      :: keep
    integer      :: it, cn
    real(kind=8) :: wl, wu, fl, fu
    real(kind=8) :: uu, up, gm
    real(kind=8) :: ss, ww, wt, wt2, wd, wp
    real(kind=8) :: f, dfw, dfu
    real(kind=8) :: g, dgw, dgu
    real(kind=8) :: det, jfw, jfu, jgw, jgu
    real(kind=8) :: dw, du
    real(kind=8) :: err
!
!-------------------------------------------------------------------------------
!
! find the initial brackets and estimate the initial enthalpy
!
    call nr_initial_brackets_srmhd_adi(mm, bb, mb, en, dn                      &
                                                     , wl, wu, w, fl, fu, info)

! if the brackets could not be found, return the lower bracket as the solution
!
    if (.not. info) then
      write(*,*)
      write(*,"(a,1x,a)"         ) "WARNING in"                                &
                                 , "EQUATIONS::nr_iterate_srmhd_adi_1dw()"
      write(*,"(a,1x)"           ) "The solution lays in unphysical regime."
      write(*,"(a,1x,1es24.16e3)") "Using the lower bracket as solution: ", wl

! use the lower bracket, since it guarantees the positive pressure
!
      w = wl

! calculate |V|² from W
!
      call nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv)

      info = .true.
      return

    end if

! and the corresponding |u|²
!
    call nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv)
    uu = vv / (1.0d+00 - vv)

! initialize iteration parameters
!
    info  = .true.
    keep  = .true.
    it    = nrmax
    cn    = nrext

! find root with the help of the Newton-Raphson method
!
    do while(keep)

! prepare (1 + |u|²) and the Lorentz factor
!
      up  = 1.0d+00 + uu
      gm  = sqrt(up)

! calculate temporary variables
!
      ss  = mb * mb
      ww  = w  * w
      wt  = w + bb
      wt2 = wt * wt
      wd  = w - gm * dn
      wp  = wt / up

! calculate functions F(W,|u|²) and G(W,|u|²)
!
      f   = w - en - gammaxi * wd / up                                         &
                              + 0.5d+00 * (bb * (1.0d+00 + uu / up) - ss / ww)
      g   = wp * wt * uu - (w + wt) * ss / ww - mm

! calculate derivatives dF(W,|u|²)/dW and dF(W,|u|²)/d|u|²
!
      dfw = 1.0d+00 - gammaxi / up + ss / ww / w
      dfu = 0.5d+00 * (gammaxi * (w + wd) + bb) / up**2

! calculate derivatives dG(W,|u|²)/dW and dG(W,|u|²)/d|u|²
!
      dgw = 2.0d+00 * wt * (uu / up + ss / ww / w)
      dgu = wp * wp

! invert the Jacobian J = | dF/dW, dF/d|u|² |
!                         | dG/dW, dG/d|u|² |
!
      det = dfw * dgu - dfu * dgw

      jfw =   dgu / det
      jgw = - dfu / det
      jfu = - dgw / det
      jgu =   dfw / det

! calculate increments dW and d|u|²
!
      dw  = f * jfw + g * jgw
      du  = f * jfu + g * jgu

! correct W and |u|²
!
      w   = w  - dw
      uu  = uu - du

! check if the new enthalpy and velocity are physical
!
      if (w < wl) then
        write(*,*)
        write(*,"(a,1x,a)"         ) "ERROR in"                                &
                                   , "EQUATIONS::nr_iterate_srmhd_adi_2dwu()"
        write(*,"(a,1x,2es24.16e3)") "Enthalpy smaller than the limit: ", w, wl
        info = .false.
        return
      end if
      if (uu < 0.0d+00) then
        write(*,*)
        write(*,"(a,1x,a)"         ) "ERROR in"                                &
                                   , "EQUATIONS::nr_iterate_srmhd_adi_2dwu()"
        write(*,"(a,1x,1es24.16e3)") "Unphysical speed |u|²: ", uu
        info = .false.
        return
      end if

! calculate the error
!
      err = max(abs(dw / w), abs(du))

! check the convergence, if the convergence is not reached, iterate until
! the maximum number of iteration is reached
!
      if (err < tol) then
        keep = cn > 0
        cn   = cn - 1
      else
        keep = it > 0
      end if

! decrease the number of remaining iterations
!
      it = it - 1

    end do ! NR iterations

! calculate |v|² from |u|²
!
    vv = uu / (1.0d+00 + uu)

! let know the user if the convergence failed
!
    if (err >= tol) then
      write(*,*)
      write(*,"(a,1x,a)"         ) "WARNING in"                                &
                                 , "EQUATIONS::nr_iterate_srmhd_adi_2dwu()"
      write(*,"(a,1x,1es24.16e3)") "Convergence not reached: ", err
    end if

!-------------------------------------------------------------------------------
!
  end subroutine nr_iterate_srmhd_adi_2dwu

!===============================================================================
!
end module equations