!!******************************************************************************
!!
!!  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-2020 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: INTERPOLATIONS
!!
!!  This module provides subroutines to interpolate variables and reconstruct
!!  the Riemann states.
!!
!!
!!******************************************************************************
!
module interpolations

#ifdef PROFILE
! import external subroutines
!
  use timers, only : set_timer, start_timer, stop_timer
#endif /* PROFILE */

! module variables are not implicit by default
!
  implicit none

#ifdef PROFILE
! timer indices
!
  integer            , save :: imi, imr, imf, imc
#endif /* PROFILE */

! interfaces for procedure pointers
!
  abstract interface
    subroutine interfaces_iface(positive, h, q, qi)
      logical                           , intent(in)  :: positive
      real(kind=8), dimension(:)        , intent(in)  :: h
      real(kind=8), dimension(:,:,:)    , intent(in)  :: q
      real(kind=8), dimension(:,:,:,:,:), intent(out) :: qi
    end subroutine
    subroutine reconstruct_iface(h, fc, fl, fr)
      real(kind=8)              , intent(in)  :: h
      real(kind=8), dimension(:), intent(in)  :: fc
      real(kind=8), dimension(:), intent(out) :: fl, fr
    end subroutine
    function limiter_iface(x, a, b) result(c)
      real(kind=8), intent(in) :: x, a, b
      real(kind=8)             :: c
    end function
  end interface

! pointers to the reconstruction and limiter procedures
!
  procedure(interfaces_iface)  , pointer, save :: interfaces         => null()
  procedure(reconstruct_iface) , pointer, save :: reconstruct_states => null()
  procedure(limiter_iface)     , pointer, save :: limiter_tvd        => null()
  procedure(limiter_iface)     , pointer, save :: limiter_prol       => null()
  procedure(limiter_iface)     , pointer, save :: limiter_clip       => null()

! method names
!
  character(len=255), save :: name_rec  = ""
  character(len=255), save :: name_tlim = ""
  character(len=255), save :: name_plim = ""
  character(len=255), save :: name_clim = ""

! module parameters
!
  real(kind=8), save :: eps        = epsilon(1.0d+00)
  real(kind=8), save :: rad        = 0.5d+00

! monotonicity preserving reconstruction coefficients
!
  real(kind=8), save :: kappa      = 1.0d+00
  real(kind=8), save :: kbeta      = 1.0d+00

! number of ghost zones (required for compact schemes)
!
  integer     , save :: ng         = 2

! order of the reconstruction
!
  integer     , save :: order      = 1

! number of cells used in the Gaussian process reconstruction
!
  integer     , save :: ngp        = 5
  integer     , save :: mgp        = 1
  integer     , save :: dgp        = 9

! normal distribution width in the Gaussian process reconstruction
!
  real(kind=8), save :: sgp        = 9.0d+00

! PPM constant
!
  real(kind=8), save :: ppm_const  = 1.25d+00

! Gaussian process reconstruction coefficients vector
!
  real(kind=8), dimension(:)    , allocatable, save :: cgp
  real(kind=8), dimension(:,:,:), allocatable, save :: ugp

! flags for reconstruction corrections
!
  logical     , save :: mlp        = .false.
  logical     , save :: positivity = .false.
  logical     , save :: clip       = .false.

! by default everything is private
!
  private

! declare public subroutines
!
  public :: initialize_interpolations, finalize_interpolations
  public :: print_interpolations
  public :: interfaces, reconstruct, limiter_prol
  public :: fix_positivity
  public :: order

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
  contains
!
!===============================================================================
!
! subroutine INITIALIZE_INTERPOLATIONS:
! ------------------------------------
!
!   Subroutine initializes the interpolation module by reading the module
!   parameters.
!
!   Arguments:
!
!     nghosts - the number of ghosts cells at each side of a block;
!     verbose - a logical flag turning the information printing;
!     status  - an integer flag for error return value;
!
!===============================================================================
!
  subroutine initialize_interpolations(nghosts, verbose, status)

! include external procedures
!
    use iso_fortran_env, only : error_unit
    use parameters     , only : get_parameter

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    integer, intent(inout) :: nghosts
    logical, intent(in)    :: verbose
    integer, intent(out)   :: status

! local variables
!
    character(len=255) :: sreconstruction = "tvd"
    character(len=255) :: tlimiter        = "mm"
    character(len=255) :: plimiter        = "mm"
    character(len=255) :: climiter        = "mm"
    character(len=255) :: mlp_limiting    = "off"
    character(len=255) :: positivity_fix  = "off"
    character(len=255) :: clip_extrema    = "off"
    character(len= 16) :: stmp
    real(kind=8)       :: cfl             = 0.5d+00

! local parameters
!
    character(len=*), parameter :: loc = 'INTERPOLATIONS:initialize_interpolation()'
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! set timer descriptions
!
    call set_timer('interpolations:: initialization', imi)
    call set_timer('interpolations:: reconstruction', imr)
    call set_timer('interpolations:: fix positivity', imf)
    call set_timer('interpolations:: clip extrema'  , imc)

! start accounting time for module initialization/finalization
!
    call start_timer(imi)
#endif /* PROFILE */

! reset the status flag
!
    status = 0

! obtain the user defined interpolation methods and coefficients
!
    call get_parameter("reconstruction"      , sreconstruction)
    call get_parameter("limiter"             , tlimiter       )
    call get_parameter("fix_positivity"      , positivity_fix )
    call get_parameter("clip_extrema"        , clip_extrema   )
    call get_parameter("extrema_limiter"     , climiter       )
    call get_parameter("prolongation_limiter", plimiter       )
    call get_parameter("mlp_limiting"        , mlp_limiting   )
    call get_parameter("ngp"                 , ngp            )
    call get_parameter("sgp"                 , sgp            )
    call get_parameter("eps"                 , eps            )
    call get_parameter("limo3_rad"           , rad            )
    call get_parameter("kappa"               , kappa          )
    call get_parameter("kbeta"               , kbeta          )
    call get_parameter("ppm_const"           , ppm_const      )
    call get_parameter("cfl"                 , cfl            )

! calculate κ = (1 - ν) / ν
!
    kappa = min(kappa, (1.0d+00 - cfl) / cfl)

! check ngp
!
    if (mod(ngp,2) == 0 .or. ngp < 3) then
      if (verbose) then
        write(error_unit,"('[', a, ']: ', a)") trim(loc),                      &
                         "The parameter ngp has to be an odd integer >= 3. "// &
                         "Resetting to default value of ngp = 5."
      end if
      ngp = 5
    end if

! calculate mgp
!
    mgp   = (ngp - 1) / 2
    dgp   = ngp**NDIMS

! select the reconstruction method
!
    select case(trim(sreconstruction))
    case ("tvd", "TVD")
      name_rec           =  "2nd order TVD"
      interfaces         => interfaces_tvd
      reconstruct_states => reconstruct_tvd
    case ("weno3", "WENO3")
      name_rec           =  "3rd order WENO"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_weno3
      order              = 3
    case ("limo3", "LIMO3", "LimO3")
      name_rec           =  "3rd order logarithmic limited"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_limo3
      order              = 3
      eps                = max(1.0d-12, eps)
    case ("ppm", "PPM")
      name_rec           =  "3rd order PPM"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_ppm
      order              = 3
      nghosts            = max(nghosts, 4)
    case ("weno5z", "weno5-z", "WENO5Z", "WENO5-Z")
      name_rec           =  "5th order WENO-Z (Borges et al. 2008)"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_weno5z
      order              = 5
      nghosts            = max(nghosts, 4)
    case ("weno5yc", "weno5-yc", "WENO5YC", "WENO5-YC")
      name_rec           =  "5th order WENO-YC (Yamaleev & Carpenter 2009)"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_weno5yc
      order              = 5
      nghosts            = max(nghosts, 4)
    case ("weno5ns", "weno5-ns", "WENO5NS", "WENO5-NS")
      name_rec           =  "5th order WENO-NS (Ha et al. 2013)"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_weno5ns
      order              = 5
      nghosts            = max(nghosts, 4)
    case ("crweno5z", "crweno5-z", "CRWENO5Z", "CRWENO5-Z")
      name_rec           =  "5th order Compact WENO-Z"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_crweno5z
      order              = 5
      nghosts            = max(nghosts, 4)
    case ("crweno5yc", "crweno5-yc", "CRWENO5YC", "CRWENO5-YC")
      name_rec           =  "5th order Compact WENO-YC"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_crweno5yc
      order              = 5
      nghosts            = max(nghosts, 4)
    case ("crweno5ns", "crweno5-ns", "CRWENO5NS", "CRWENO5-NS")
      name_rec           =  "5th order Compact WENO-NS"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_crweno5ns
      order              = 5
      nghosts            = max(nghosts, 4)
    case ("mp5", "MP5")
      name_rec           =  "5th order Monotonicity Preserving"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_mp5
      order              = 5
      nghosts            = max(nghosts, 4)
    case ("mp7", "MP7")
      name_rec           =  "7th order Monotonicity Preserving"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_mp7
      order              = 7
      nghosts            = max(nghosts, 4)
    case ("mp9", "MP9")
      name_rec           =  "9th order Monotonicity Preserving"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_mp9
      order              = 9
      nghosts            = max(nghosts, 6)
    case ("mp5ld", "MP5LD")
      name_rec           =  "5th order Low-Dissipation Monotonicity Preserving"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_mp5ld
      order              = 5
      nghosts            = max(nghosts, 4)
    case ("mp7ld", "MP7LD")
      name_rec           =  "7th order Low-Dissipation Monotonicity Preserving"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_mp7ld
      order              = 7
      nghosts            = max(nghosts, 4)
    case ("mp9ld", "MP9LD")
      name_rec           =  "9th order Low-Dissipation Monotonicity Preserving"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_mp9ld
      order              = 9
      nghosts            = max(nghosts, 6)
    case ("crmp5", "CRMP5")
      name_rec           =  "5th order Compact Monotonicity Preserving"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_crmp5
      order              = 5
      nghosts            = max(nghosts, 4)
    case ("crmp5l", "crmp5ld", "CRMP5L", "CRMP5LD")
      name_rec           =  "5th order Low-Dissipation Compact Monotonicity Preserving"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_crmp5ld
      order              = 5
      nghosts            = max(nghosts, 4)
    case ("crmp7", "CRMP7")
      name_rec           =  "7th order Compact Monotonicity Preserving"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_crmp7
      order              = 7
      nghosts            = max(nghosts, 4)
    case ("crmp7l", "crmp7ld", "CRMP7L", "CRMP7LD")
      name_rec           =  "7th order Low-Dissipation Compact Monotonicity Preserving"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_crmp7ld
      order              = 7
      nghosts            = max(nghosts, 4)
    case ("ocmp5", "OCMP5")
      name_rec           =  "5th order Optimized Compact Monotonicity Preserving"
      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_ocmp5
      order              = 5
      nghosts            = max(nghosts, 4)
    case ("gp", "GP")
      write(stmp, '(f16.1)') sgp
      write(name_rec, '("Gaussian Process (",i1,"-point, δ=",a,")")') ngp      &
                                                         , trim(adjustl(stmp))

! allocate the Gaussian process reconstruction matrix and position vector
!
      allocate(cgp(ngp), stat = status)

! prepare matrix coefficients
!
      if (status == 0) call prepare_gp(status)

      interfaces         => interfaces_dir
      reconstruct_states => reconstruct_gp
      order              = ngp

      ng = 2
      do while(ng < (ngp + 1) / 2)
        ng = ng + 2
      end do
      nghosts = max(nghosts, ng)

    case ("mgp", "MGP")
      write(stmp, '(f16.1)') sgp
      write(name_rec, &
          '("Multidimensional Gaussian Process (",i1,"-point, δ=",a,")")')     &
                                                      ngp, trim(adjustl(stmp))

! allocate the Gaussian process reconstruction matrix and position vector
!
      allocate(ugp(dgp,2,NDIMS), stat = status)

! prepare matrix coefficients
!
      if (status == 0) call prepare_mgp(status)

      interfaces => interfaces_mgp

      ng = 2
      do while(ng < (ngp + 1) / 2)
        ng = ng + 2
      end do
      nghosts = max(nghosts, ng)

    case default

      if (verbose) then
        write(error_unit,"('[', a, ']: ', a)") trim(loc),                      &
                         "The selected reconstruction method is not "    //    &
                         "implemented: " // trim(sreconstruction) // "." //    &
                         "Available methods: 'tvd' 'limo3', 'ppm',"      //    &
                         " 'weno3', 'weno5z', 'weno5yc', 'weno5ns',"     //    &
                         " 'crweno5z', 'crweno5yc', 'crweno5ns','mp5', " //    &
                         " 'mp7', 'mp9', 'crmp5', 'crmp5ld', 'crmp7', 'gp', 'mgp'."
      end if
      status = 1

    end select

! select the TVD limiter
!
    select case(trim(tlimiter))
    case ("mm", "minmod")
      name_tlim          =  "minmod"
      limiter_tvd        => limiter_minmod
      order              = max(2, order)
    case ("mc", "monotonized_central")
      name_tlim          =  "monotonized central"
      limiter_tvd        => limiter_monotonized_central
      order              = max(2, order)
    case ("sb", "superbee")
      name_tlim          =  "superbee"
      limiter_tvd        => limiter_superbee
      order              = max(2, order)
    case ("vl", "vanleer")
      name_tlim          =  "van Leer"
      limiter_tvd        => limiter_vanleer
      order              = max(2, order)
    case ("va", "vanalbada")
      name_tlim          =  "van Albada"
      limiter_tvd        => limiter_vanalbada
      order              = max(2, order)
    case default
      name_tlim          =  "zero derivative"
      limiter_tvd        => limiter_zero
      order              = max(1, order)
    end select

! select the prolongation limiter
!
    select case(trim(plimiter))
    case ("mm", "minmod")
      name_plim          =  "minmod"
      limiter_prol       => limiter_minmod
    case ("mc", "monotonized_central")
      name_plim          =  "monotonized central"
      limiter_prol       => limiter_monotonized_central
    case ("sb", "superbee")
      name_plim          =  "superbee"
      limiter_prol       => limiter_superbee
    case ("vl", "vanleer")
      name_plim          =  "van Leer"
      limiter_prol       => limiter_vanleer
    case default
      name_plim          =  "zero derivative"
      limiter_prol       => limiter_zero
    end select

! select the clipping limiter
!
    select case(trim(climiter))
    case ("mm", "minmod")
      name_clim          =  "minmod"
      limiter_clip       => limiter_minmod
    case ("mc", "monotonized_central")
      name_clim          =  "monotonized central"
      limiter_clip       => limiter_monotonized_central
    case ("sb", "superbee")
      name_clim          =  "superbee"
      limiter_clip       => limiter_superbee
    case ("vl", "vanleer")
      name_clim          =  "van Leer"
      limiter_clip       => limiter_vanleer
    case default
      name_clim          =  "zero derivative"
      limiter_clip       => limiter_zero
    end select

! check additional reconstruction limiting
!
    select case(trim(mlp_limiting))
    case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
      mlp = .true.
    case default
      mlp = .false.
    end select
    select case(trim(positivity_fix))
    case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
      positivity = .true.
    case default
      positivity = .false.
    end select
    select case(trim(clip_extrema))
    case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
      clip = .true.
    case default
      clip = .false.
    end select

! set the internal number of ghost cells
!
    ng = nghosts

#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
    call stop_timer(imi)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine initialize_interpolations
!
!===============================================================================
!
! subroutine FINALIZE_INTERPOLATIONS:
! ----------------------------------
!
!   Subroutine finalizes the interpolation module by releasing all memory used
!   by its module variables.
!
!   Arguments:
!
!     status  - an integer flag for error return value;
!
!===============================================================================
!
  subroutine finalize_interpolations(status)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    integer, intent(out) :: status
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for module initialization/finalization
!
    call start_timer(imi)
#endif /* PROFILE */

! reset the status flag
!
    status = 0

! deallocate Gaussian process reconstruction coefficient vector if used
!
    if (allocated(cgp)) deallocate(cgp, stat = status)
    if (allocated(ugp)) deallocate(ugp, stat = status)

! release the procedure pointers
!
    nullify(reconstruct_states)
    nullify(limiter_tvd)
    nullify(limiter_prol)
    nullify(limiter_clip)

#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
    call stop_timer(imi)
#endif /* PROFILE */

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

! import 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
!
!-------------------------------------------------------------------------------
!
    if (verbose) then
      call print_section(verbose, "Interpolations")
      call print_parameter(verbose, "reconstruction"      , name_rec )
      call print_parameter(verbose, "TVD limiter"         , name_tlim)
      call print_parameter(verbose, "prolongation limiter", name_plim)
      if (mlp) then
        call print_parameter(verbose, "MLP limiting"        , "on"     )
      else
        call print_parameter(verbose, "MLP limiting"        , "off"    )
      end if
      if (positivity) then
        call print_parameter(verbose, "fix positivity"      , "on"     )
      else
        call print_parameter(verbose, "fix positivity"      , "off"    )
      end if
      if (clip) then
        call print_parameter(verbose, "clip extrema"        , "on"     )
        call print_parameter(verbose, "extrema limiter"     , name_clim)
      else
        call print_parameter(verbose, "clip extrema"        , "off"    )
      end if

    end if

!-------------------------------------------------------------------------------
!
  end subroutine print_interpolations
!
!===============================================================================
!
! subroutine PREPARE_MGP:
! ---------------------
!
!   Subroutine prepares matrixes for the multidimensional
!   Gaussian Process (GP) method.
!
!===============================================================================
!
  subroutine prepare_mgp(iret)

! include external procedures
!
    use algebra        , only : max_prec, invert
    use constants      , only : pi
    use iso_fortran_env, only : error_unit

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    integer, intent(inout) :: iret

! local variables
!
    logical             :: flag
    integer             :: i, j, i1, j1, i2, j2
#if NDIMS == 3
    integer             :: k1, k2
#endif /* NDIMS == 3 */
    real(kind=max_prec) :: sig, fc, fx, fy, xl, xr, yl, yr
#if NDIMS == 3
    real(kind=max_prec) :: fz, zl, zr
#endif /* NDIMS == 3 */

! local arrays for derivatives
!
    real(kind=max_prec), dimension(:,:)  , allocatable :: cov, inv
    real(kind=max_prec), dimension(:,:,:), allocatable :: xgp

! local parameters
!
    character(len=*), parameter :: loc = 'INTERPOLATIONS::prepare_mgp()'
!
!-------------------------------------------------------------------------------
!
! calculate normal distribution sigma
!
    sig = sqrt(2.0d+00) * sgp

! allocate the convariance matrix and interpolation position vector
!
    allocate(cov(dgp,dgp))
    allocate(inv(dgp,dgp))
    allocate(xgp(dgp,2,NDIMS))

! prepare the covariance matrix
!
    fc = 0.5d+00 * sqrt(pi) * sig
    i = 0
#if NDIMS == 3
    do k1 = - mgp, mgp
#endif /* NDIMS == 3 */
      do j1 = - mgp, mgp
        do i1 = - mgp, mgp
          i = i + 1
          j = 0
#if NDIMS == 3
          do k2 = - mgp, mgp
#endif /* NDIMS == 3 */
            do j2 = - mgp, mgp
              do i2 = - mgp, mgp
                j = j + 1

                xl = (1.0d+00 * (i1 - i2) - 0.5d+00) / sig
                xr = (1.0d+00 * (i1 - i2) + 0.5d+00) / sig
                yl = (1.0d+00 * (j1 - j2) - 0.5d+00) / sig
                yr = (1.0d+00 * (j1 - j2) + 0.5d+00) / sig
#if NDIMS == 3
                zl = (1.0d+00 * (k1 - k2) - 0.5d+00) / sig
                zr = (1.0d+00 * (k1 - k2) + 0.5d+00) / sig
#endif /* NDIMS == 3 */

                cov(i,j) = fc * (erf(xr) - erf(xl)) * (erf(yr) - erf(yl))
#if NDIMS == 3
                cov(i,j) = cov(i,j) * (erf(zr) - erf(zl))
#endif /* NDIMS == 3 */
              end do
            end do
          end do
        end do
#if NDIMS == 3
      end do
    end do
#endif /* NDIMS == 3 */

! prepare the interpolation position vector
!
    i = 0
#if NDIMS == 3
    do k1 = - mgp, mgp
#endif /* NDIMS == 3 */
      do j1 = - mgp, mgp
        do i1 = - mgp, mgp
          i = i + 1

          xl = (1.0d+00 * i1 - 0.5d+00) / sig
          xr = (1.0d+00 * i1 + 0.5d+00) / sig
          yl = (1.0d+00 * j1 - 0.5d+00) / sig
          yr = (1.0d+00 * j1 + 0.5d+00) / sig
#if NDIMS == 3
          zl = (1.0d+00 * k1 - 0.5d+00) / sig
          zr = (1.0d+00 * k1 + 0.5d+00) / sig
#endif /* NDIMS == 3 */

          fx = erf(xr) - erf(xl)
          fy = erf(yr) - erf(yl)
#if NDIMS == 3
          fz = erf(zr) - erf(zl)

          xgp(i,1,1) = exp(- xl**2) * fy * fz
          xgp(i,2,1) = exp(- xr**2) * fy * fz
          xgp(i,1,2) = exp(- yl**2) * fx * fz
          xgp(i,2,2) = exp(- yr**2) * fx * fz
          xgp(i,1,3) = exp(- zl**2) * fx * fy
          xgp(i,2,3) = exp(- zr**2) * fx * fy
#else /* NDIMS == 3 */
          xgp(i,1,1) = exp(- xl**2) * fy
          xgp(i,2,1) = exp(- xr**2) * fy
          xgp(i,1,2) = exp(- yl**2) * fx
          xgp(i,2,2) = exp(- yr**2) * fx
#endif /* NDIMS == 3 */
        end do
      end do
#if NDIMS == 3
    end do
#endif /* NDIMS == 3 */

! invert the matrix
!
    call invert(dgp, cov(1:dgp,1:dgp), inv(1:dgp,1:dgp), flag)

! prepare the interpolation coefficients vector
!
    do j = 1, NDIMS
      do i = 1, 2
        ugp(1:dgp,i,j) = real(matmul(xgp(1:dgp,i,j), inv(1:dgp,1:dgp)), kind=8)
      end do
    end do

! deallocate the convariance matrix and interpolation position vector
!
    deallocate(cov)
    deallocate(inv)
    deallocate(xgp)

! check if the matrix was inverted successfully
!
    if (.not. flag) then
      write(error_unit,"('[',a,']: ',a)") trim(loc)                            &
                      , "Could not invert the covariance matrix!"
      iret = 201
    end if

!-------------------------------------------------------------------------------
!
  end subroutine prepare_mgp
!
!===============================================================================
!
! subroutine INTERFACES_TVD:
! -------------------------
!
!   Subroutine reconstructs both side interfaces of variable using TVD methods.
!
!   Arguments:
!
!     positive - the variable positivity flag;
!     h        - the spatial step;
!     q        - the variable array;
!     qi       - the array of reconstructed interfaces (2 in each direction);
!
!===============================================================================
!
  subroutine interfaces_tvd(positive, h, q, qi)

! include external procedures
!
    use coordinates    , only : nb, ne, nbl, neu
    use iso_fortran_env, only : error_unit

! local variables are not implicit by default
!
    implicit none

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

! local variables
!
    integer :: i  , j  , k = 1
    integer :: im1, jm1, ip1, jp1
#if NDIMS == 3
    integer :: km1, kp1
#endif /* NDIMS == 3 */

! local vectors
!
    real(kind=8), dimension(NDIMS) :: dql, dqr, dq

! local parameters
!
    character(len=*), parameter :: loc = 'INTERPOLATIONS::interfaces_tvd()'
!
!-------------------------------------------------------------------------------
!
! copy ghost zones
!
    do j = 1, NDIMS
      do i = 1, 2
        qi(  :nb,  :  ,  :  ,i,j) = q(  :nb,  :  ,  :  )
        qi(ne:  ,  :  ,  :  ,i,j) = q(ne:  ,  :  ,  :  )
        qi(nb:ne,  :nb,  :  ,i,j) = q(nb:ne,  :nb,  :  )
        qi(nb:ne,ne:  ,  :  ,i,j) = q(nb:ne,ne:  ,  :  )
#if NDIMS == 3
        qi(nb:ne,nb:ne,  :nb,i,j) = q(nb:ne,nb:ne,  :nb)
        qi(nb:ne,nb:ne,ne:  ,i,j) = q(nb:ne,nb:ne,ne:  )
#endif /* NDIMS == 3 */
      end do
    end do

! interpolate interfaces
!
#if NDIMS == 3
    do k = nbl, neu
      km1 = k - 1
      kp1 = k + 1
#endif /* NDIMS == 3 */
      do j = nbl, neu
        jm1 = j - 1
        jp1 = j + 1
        do i = nbl, neu
          im1 = i - 1
          ip1 = i + 1

! calculate the TVD derivatives
!
          dql(1) = q(i  ,j,k) - q(im1,j,k)
          dqr(1) = q(ip1,j,k) - q(i  ,j,k)
          dq (1) = limiter_tvd(0.5d+00, dql(1), dqr(1))

          dql(2) = q(i,j  ,k) - q(i,jm1,k)
          dqr(2) = q(i,jp1,k) - q(i,j  ,k)
          dq (2) = limiter_tvd(0.5d+00, dql(2), dqr(2))

#if NDIMS == 3
          dql(3) = q(i,j,k  ) - q(i,j,km1)
          dqr(3) = q(i,j,kp1) - q(i,j,k  )
          dq (3) = limiter_tvd(0.5d+00, dql(3), dqr(3))
#endif /* NDIMS == 3 */

! limit the derivatives if they produce negative interpolation for positive
! variables
!
          if (positive) then
            if (q(i,j,k) > 0.0d+00) then
              do while (q(i,j,k) <= sum(abs(dq(1:NDIMS))))
                dq(:) = 0.5d+00 * dq(:)
              end do
            else
              write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc)              &
                    , "Positive variable is not positive at ( ", i, j, k, " )"
              dq(:) = 0.0d+00
            end if
          end if

! interpolate states
!
          qi(i  ,j,k,1,1) = q(i,j,k) + dq(1)
          qi(im1,j,k,2,1) = q(i,j,k) - dq(1)

          qi(i,j  ,k,1,2) = q(i,j,k) + dq(2)
          qi(i,jm1,k,2,2) = q(i,j,k) - dq(2)

#if NDIMS == 3
          qi(i,j,k  ,1,3) = q(i,j,k) + dq(3)
          qi(i,j,km1,2,3) = q(i,j,k) - dq(3)
#endif /* NDIMS == 3 */

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

! apply Multi-dimensional Limiting Process
!
    if (mlp) call mlp_limiting(q(:,:,:), qi(:,:,:,:,:))

!-------------------------------------------------------------------------------
!
  end subroutine interfaces_tvd
!
!===============================================================================
!
! subroutine INTERFACES_DIR:
! -------------------------
!
!   Subroutine reconstructs both side interfaces of variable separately
!   along each direction.
!
!   Arguments:
!
!     positive - the variable positivity flag;
!     h        - the spatial step;
!     q        - the variable array;
!     qi       - the array of reconstructed interfaces (2 in each direction);
!
!===============================================================================
!
  subroutine interfaces_dir(positive, h, q, qi)

! include external procedures
!
    use coordinates, only : nb, ne, nbl, neu

! local variables are not implicit by default
!
    implicit none

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

! local variables
!
    integer :: i, j, k = 1
!
!-------------------------------------------------------------------------------
!
! copy ghost zones
!
    do j = 1, NDIMS
      do i = 1, 2
        qi(  :nb,  :  ,  :  ,i,j) = q(  :nb,  :  ,  :  )
        qi(ne:  ,  :  ,  :  ,i,j) = q(ne:  ,  :  ,  :  )
        qi(nb:ne,  :nb,  :  ,i,j) = q(nb:ne,  :nb,  :  )
        qi(nb:ne,ne:  ,  :  ,i,j) = q(nb:ne,ne:  ,  :  )
#if NDIMS == 3
        qi(nb:ne,nb:ne,  :nb,i,j) = q(nb:ne,nb:ne,  :nb)
        qi(nb:ne,nb:ne,ne:  ,i,j) = q(nb:ne,nb:ne,ne:  )
#endif /* NDIMS == 3 */
      end do
    end do

! interpolate interfaces
!
#if NDIMS == 3
    do k = nbl, neu
#endif /* NDIMS == 3 */
      do j = nbl, neu
        call reconstruct(h(1), q(:,j,k), qi(:,j,k,1,1), qi(:,j,k,2,1))
      end do ! j = nbl, neu
      do i = nbl, neu
        call reconstruct(h(2), q(i,:,k), qi(i,:,k,1,2), qi(i,:,k,2,2))
      end do ! i = nbl, neu
#if NDIMS == 3
    end do ! k = nbl, neu
    do j = nbl, neu
      do i = nbl, neu
        call reconstruct(h(3), q(i,j,:), qi(i,j,:,1,3), qi(i,j,:,2,3))
      end do ! i = nbl, neu
    end do ! j = nbl, neu
#endif /* NDIMS == 3 */

! make sure the interface states are positive for positive variables
!
    if (positive) then

#if NDIMS == 3
      do k = nbl, neu
#endif /* NDIMS == 3 */
        do j = nbl, neu
          call fix_positivity(q(:,j,k), qi(:,j,k,1,1), qi(:,j,k,2,1))
        end do ! j = nbl, neu
        do i = nbl, neu
          call fix_positivity(q(i,:,k), qi(i,:,k,1,2), qi(i,:,k,2,2))
        end do ! i = nbl, neu
#if NDIMS == 3
      end do ! k = nbl, neu
      do j = nbl, neu
        do i = nbl, neu
          call fix_positivity(q(i,j,:), qi(i,j,:,1,3), qi(i,j,:,2,3))
        end do ! i = nbl, neu
      end do ! j = nbl, neu
#endif /* NDIMS == 3 */

    end if

! apply Multi-dimensional Limiting Process
!
    if (mlp) call mlp_limiting(q(:,:,:), qi(:,:,:,:,:))

!-------------------------------------------------------------------------------
!
  end subroutine interfaces_dir
!
!===============================================================================
!
! subroutine INTERFACES_MGP:
! -------------------------
!
!   Subroutine reconstructs both side interfaces of variable using
!   multidimensional Gaussian Process method.
!
!   Arguments:
!
!     positive - the variable positivity flag;
!     h        - the spatial step;
!     q        - the variable array;
!     qi       - the array of reconstructed interfaces (2 in each direction);
!
!===============================================================================
!
  subroutine interfaces_mgp(positive, h, q, qi)

! include external procedures
!
    use coordinates    , only : nb, ne, nbl, neu
    use iso_fortran_env, only : error_unit

! local variables are not implicit by default
!
    implicit none

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

! local variables
!
    logical :: flag
    integer :: i, il, iu, im1, ip1
    integer :: j, jl, ju, jm1, jp1
    integer :: k = 1
#if NDIMS == 3
    integer :: kl, ku, km1, kp1
#endif /* NDIMS == 3 */

! local arrays for derivatives
!
    real(kind=8), dimension(NDIMS) :: dql, dqr, dq
    real(kind=8), dimension(dgp)   :: u

! local parameters
!
    character(len=*), parameter :: loc = 'INTERPOLATIONS::interfaces_mgp()'
!
!-------------------------------------------------------------------------------
!
! copy ghost zones
!
    do j = 1, NDIMS
      do i = 1, 2
        qi(  :nb,  :  ,  :  ,i,j) = q(  :nb,  :  ,  :  )
        qi(ne:  ,  :  ,  :  ,i,j) = q(ne:  ,  :  ,  :  )
        qi(nb:ne,  :nb,  :  ,i,j) = q(nb:ne,  :nb,  :  )
        qi(nb:ne,ne:  ,  :  ,i,j) = q(nb:ne,ne:  ,  :  )
#if NDIMS == 3
        qi(nb:ne,nb:ne,  :nb,i,j) = q(nb:ne,nb:ne,  :nb)
        qi(nb:ne,nb:ne,ne:  ,i,j) = q(nb:ne,nb:ne,ne:  )
#endif /* NDIMS == 3 */
      end do
    end do

! interpolate interfaces using precomputed interpolation vectors
!
#if NDIMS == 3
    do k = nbl, neu
      kl  = k - mgp
      ku  = k + mgp
      km1 = k - 1
      kp1 = k + 1
#endif /* NDIMS == 3 */
      do j = nbl, neu
        jl  = j - mgp
        ju  = j + mgp
        jm1 = j - 1
        jp1 = j + 1
        do i = nbl, neu
          il  = i - mgp
          iu  = i + mgp
          im1 = i - 1
          ip1 = i + 1

#if NDIMS == 3
          u(:) = reshape(q(il:iu,jl:ju,kl:ku), (/ dgp /)) - q(i,j,k)

          qi(i  ,j,k,1,1) = sum(ugp(1:dgp,1,1) * u(1:dgp)) + q(i,j,k)
          qi(im1,j,k,2,1) = sum(ugp(1:dgp,2,1) * u(1:dgp)) + q(i,j,k)
          qi(i,j  ,k,1,2) = sum(ugp(1:dgp,1,2) * u(1:dgp)) + q(i,j,k)
          qi(i,jm1,k,2,2) = sum(ugp(1:dgp,2,2) * u(1:dgp)) + q(i,j,k)
          qi(i,j,k  ,1,3) = sum(ugp(1:dgp,1,3) * u(1:dgp)) + q(i,j,k)
          qi(i,j,km1,2,3) = sum(ugp(1:dgp,2,3) * u(1:dgp)) + q(i,j,k)
#else /* NDIMS == 3 */
          u(:) = reshape(q(il:iu,jl:ju,k    ), (/ dgp /)) - q(i,j,k)

          qi(i  ,j,k,1,1) = sum(ugp(1:dgp,1,1) * u(1:dgp)) + q(i,j,k)
          qi(im1,j,k,2,1) = sum(ugp(1:dgp,2,1) * u(1:dgp)) + q(i,j,k)
          qi(i,j  ,k,1,2) = sum(ugp(1:dgp,1,2) * u(1:dgp)) + q(i,j,k)
          qi(i,jm1,k,2,2) = sum(ugp(1:dgp,2,2) * u(1:dgp)) + q(i,j,k)
#endif /* NDIMS == 3 */

! if the interpolation is not monotonic, apply a TVD slope
!
          flag =           ((qi(i  ,j,k,1,1) - q(ip1,j,k))                     &
                          * (qi(i  ,j,k,1,1) - q(i  ,j,k)) > eps)
          flag = flag .or. ((qi(im1,j,k,2,1) - q(im1,j,k))                     &
                          * (qi(im1,j,k,2,1) - q(i  ,j,k)) > eps)
          flag = flag .or. ((qi(i,j  ,k,1,2) - q(i,jp1,k))                     &
                          * (qi(i,j  ,k,1,2) - q(i,j  ,k)) > eps)
          flag = flag .or. ((qi(i,jm1,k,2,2) - q(i,jm1,k))                     &
                          * (qi(i,jm1,k,2,2) - q(i,j  ,k)) > eps)
#if NDIMS == 3
          flag = flag .or. ((qi(i,j,k  ,1,3) - q(i,j,kp1))                     &
                          * (qi(i,j,k  ,1,3) - q(i,j,k  )) > eps)
          flag = flag .or. ((qi(i,j,km1,2,3) - q(i,j,km1))                     &
                          * (qi(i,j,km1,2,3) - q(i,j,k  )) > eps)
#endif /* NDIMS == 3 */

          if (flag) then

! calculate the TVD derivatives
!
            dql(1) = q(i  ,j,k) - q(im1,j,k)
            dqr(1) = q(ip1,j,k) - q(i  ,j,k)
            dq (1) = limiter_tvd(0.5d+00, dql(1), dqr(1))

            dql(2) = q(i,j  ,k) - q(i,jm1,k)
            dqr(2) = q(i,jp1,k) - q(i,j  ,k)
            dq (2) = limiter_tvd(0.5d+00, dql(2), dqr(2))

#if NDIMS == 3
            dql(3) = q(i,j,k  ) - q(i,j,km1)
            dqr(3) = q(i,j,kp1) - q(i,j,k  )
            dq (3) = limiter_tvd(0.5d+00, dql(3), dqr(3))
#endif /* NDIMS == 3 */

! limit the derivatives if they produce negative interpolation for positive
! variables
!
            if (positive) then
              if (q(i,j,k) > 0.0d+00) then
                do while (q(i,j,k) <= sum(abs(dq(1:NDIMS))))
                  dq(:) = 0.5d+00 * dq(:)
                end do
              else
                write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc)            &
                    , "Positive variable is not positive at ( ", i, j, k, " )"
                dq(:) = 0.0d+00
              end if
            end if

! interpolate states
!
            qi(i  ,j,k,1,1) = q(i,j,k) + dq(1)
            qi(im1,j,k,2,1) = q(i,j,k) - dq(1)

            qi(i,j  ,k,1,2) = q(i,j,k) + dq(2)
            qi(i,jm1,k,2,2) = q(i,j,k) - dq(2)

#if NDIMS == 3
            qi(i,j,k  ,1,3) = q(i,j,k) + dq(3)
            qi(i,j,km1,2,3) = q(i,j,k) - dq(3)
#endif /* NDIMS == 3 */

          end if

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

! apply Multi-dimensional Limiting Process
!
    if (mlp) call mlp_limiting(q(:,:,:), qi(:,:,:,:,:))

!-------------------------------------------------------------------------------
!
  end subroutine interfaces_mgp
!
!===============================================================================
!
! subroutine MLP_LIMITING:
! -----------------------
!
!   Subroutine applies the multi-dimensional limiting process to
!   the reconstructed states in order to control oscillations due to
!   the Gibbs phenomena near discontinuities.
!
!   Arguments:
!
!     q        - the variable array;
!     qi       - the array of reconstructed interfaces (2 in each direction);
!
!   References:
!
!     [1] Gerlinger, P.,
!         "Multi-dimensional limiting for high-order schemes including
!          turbulence and combustion",
!         Journal of Computational Physics,
!         2012, vol. 231, pp. 2199-2228,
!         http://dx.doi.org/10.1016/j.jcp.2011.10.024
!
!===============================================================================
!
  subroutine mlp_limiting(q, qi)

! include external procedures
!
    use coordinates, only : nn => bcells

! local variables are not implicit by default
!
    implicit none

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

! local variables
!
    integer      :: i, im1, ip1
    integer      :: j, jm1, jp1
    integer      :: k = 1
#if NDIMS == 3
    integer      :: km1, kp1
#endif /* NDIMS == 3 */
    integer      :: m
#if NDIMS == 3
    integer      :: n, np1, np2
#endif /* NDIMS == 3 */
    real(kind=8) :: qmn, qmx, dqc
    real(kind=8) :: fc, fl

! local vectors
!
    real(kind=8), dimension(NDIMS) :: dql, dqr, dq
    real(kind=8), dimension(NDIMS) :: dqm, ap
#if NDIMS == 3
    real(kind=8), dimension(NDIMS) :: hh, uu
#endif /* NDIMS == 3 */
!
!-------------------------------------------------------------------------------
!
#if NDIMS == 3
    do k = 1, nn
      km1 = max( 1, k - 1)
      kp1 = min(nn, k + 1)
#endif /* NDIMS == 3 */
      do j = 1, nn
        jm1 = max( 1, j - 1)
        jp1 = min(nn, j + 1)
        do i = 1, nn
          im1 = max( 1, i - 1)
          ip1 = min(nn, i + 1)

! calculate the minmod TVD derivatives
!
          dql(1) = q(i  ,j,k) - q(im1,j,k)
          dqr(1) = q(ip1,j,k) - q(i  ,j,k)
          dq (1) = limiter_minmod(0.5d+00, dql(1), dqr(1))

          dql(2) = q(i,j  ,k) - q(i,jm1,k)
          dqr(2) = q(i,jp1,k) - q(i,j  ,k)
          dq (2) = limiter_minmod(0.5d+00, dql(2), dqr(2))

#if NDIMS == 3
          dql(3) = q(i,j,k  ) - q(i,j,km1)
          dqr(3) = q(i,j,kp1) - q(i,j,k  )
          dq (3) = limiter_minmod(0.5d+00, dql(3), dqr(3))
#endif /* NDIMS == 3 */

! calculate dqc
!
#if NDIMS == 3
          qmn = minval(q(im1:ip1,jm1:jp1,km1:kp1))
          qmx = maxval(q(im1:ip1,jm1:jp1,km1:kp1))
#else /* NDIMS == 3 */
          qmn = minval(q(im1:ip1,jm1:jp1,k))
          qmx = maxval(q(im1:ip1,jm1:jp1,k))
#endif /* NDIMS == 3 */
          dqc = min(qmx - q(i,j,k), q(i,j,k) - qmn)

! if needed, apply the multi-dimensional limiting process
!
          if (sum(abs(dq(1:NDIMS))) > dqc) then

            dqm(1) = abs(q(ip1,j,k) - q(im1,j,k))
            dqm(2) = abs(q(i,jp1,k) - q(i,jm1,k))
#if NDIMS == 3
            dqm(3) = abs(q(i,j,kp1) - q(i,j,km1))
#endif /* NDIMS == 3 */

            fc = dqc / max(1.0d-16, sum(abs(dqm(1:NDIMS))))
            do m = 1, NDIMS
              ap(m) = fc * abs(dqm(m))
            end do

#if NDIMS == 3
            do n = 1, NDIMS
              hh(n)   = max(ap(n) - abs(dq(n)), 0.0d+00)
              np1     = mod(n    , NDIMS) + 1
              np2     = mod(n + 1, NDIMS) + 1
              uu(n)   = ap(n)   -           hh(n)
              uu(np1) = ap(np1) + 0.5d+00 * hh(n)
              uu(np2) = ap(np2) + 0.5d+00 * hh(n)
              fc      = hh(n) / (hh(n) + 1.0d-16)
              fl      = fc * (max(ap(np1) - abs(dq(np1)), 0.0d+00)             &
                            - max(ap(np2) - abs(dq(np2)), 0.0d+00))
              ap(n  ) = uu(n)
              ap(np1) = uu(np1) - fl
              ap(np2) = uu(np2) + fl
            end do
#else /* NDIMS == 3 */
            fl     = max(ap(1) - abs(dq(1)), 0.0d+00)                          &
                   - max(ap(2) - abs(dq(2)), 0.0d+00)
            ap(1)  = ap(1) - fl
            ap(2)  = ap(2) + fl
#endif /* NDIMS == 3 */

            do m = 1, NDIMS
              dq(m) = sign(ap(m), dq(m))
            end do

          end if

! update the interpolated states
!
          dql(1) = qi(i  ,j,k,1,1) - q(i,j,k)
          dqr(1) = qi(im1,j,k,2,1) - q(i,j,k)
          if (max(abs(dql(1)), abs(dqr(1))) > abs(dq(1))) then
            qi(i  ,j,k,1,1) = q(i,j,k) + dq(1)
            qi(im1,j,k,2,1) = q(i,j,k) - dq(1)
          end if

          dql(2) = qi(i,j  ,k,1,2) - q(i,j,k)
          dqr(2) = qi(i,jm1,k,2,2) - q(i,j,k)
          if (max(abs(dql(2)), abs(dqr(2))) > abs(dq(2))) then
            qi(i,j  ,k,1,2) = q(i,j,k) + dq(2)
            qi(i,jm1,k,2,2) = q(i,j,k) - dq(2)
          end if

#if NDIMS == 3
          dql(3) = qi(i,j,k  ,1,3) - q(i,j,k)
          dqr(3) = qi(i,j,km1,2,3) - q(i,j,k)
          if (max(abs(dql(3)), abs(dqr(3))) > abs(dq(3))) then
            qi(i,j,k  ,1,3) = q(i,j,k) + dq(3)
            qi(i,j,km1,2,3) = q(i,j,k) - dq(3)
          end if
#endif /* NDIMS == 3 */

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

!-------------------------------------------------------------------------------
!
  end subroutine mlp_limiting
!
!===============================================================================
!
! subroutine RECONSTRUCT:
! ----------------------
!
!   Subroutine calls a reconstruction procedure, depending on the compilation
!   flag SPACE, in order to interpolate the left and right states from their
!   cell integrals.  These states are required by any approximate Riemann
!   solver.
!
!   Arguments:
!
!     h  - the spatial step; this is required for some reconstruction methods;
!     f  - the input vector of cell averaged values;
!     fl - the left side state reconstructed for location (i+1/2);
!     fr - the right side state reconstructed for location (i+1/2);
!
!===============================================================================
!
  subroutine reconstruct(h, f, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: f
    real(kind=8), dimension(:), intent(out) :: fl, fr
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for reconstruction
!
    call start_timer(imr)
#endif /* PROFILE */

! reconstruct the states using the selected subroutine
!
    call reconstruct_states(h, f(:), fl(:), fr(:))

! correct the reconstruction near extrema by clipping them in order to improve
! the stability of scheme
!
    if (clip) call clip_extrema(f(:), fl(:), fr(:))

#ifdef PROFILE
! stop accounting time for reconstruction
!
    call stop_timer(imr)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct
!
!===============================================================================
!
! subroutine RECONSTRUCT_TVD:
! --------------------------
!
!   Subroutine reconstructs the interface states using the second order TVD
!   method with a selected limiter.
!
!   Arguments are described in subroutine reconstruct().
!
!
!===============================================================================
!
  subroutine reconstruct_tvd(h, fc, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer      :: n, i, im1, ip1
    real(kind=8) :: df, dfl, dfr
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

! calculate the left- and right-side interface interpolations
!
    do i = 2, n - 1

! calculate left and right indices
!
      im1     = i - 1
      ip1     = i + 1

! calculate left and right side derivatives
!
      dfl     = fc(i  ) - fc(im1)
      dfr     = fc(ip1) - fc(i  )

! obtain the TVD limited derivative
!
      df      = limiter_tvd(0.5d+00, dfl, dfr)

! update the left and right-side interpolation states
!
      fl(i  ) = fc(i) + df
      fr(im1) = fc(i) - df

    end do ! i = 2, n - 1

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_tvd
!
!===============================================================================
!
! subroutine RECONSTRUCT_WENO3:
! ----------------------------
!
!   Subroutine reconstructs the interface states using the third order
!   Weighted Essentially Non-Oscillatory (WENO) method.
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Yamaleev & Carpenter, 2009, J. Comput. Phys., 228, 3025
!
!===============================================================================
!
  subroutine reconstruct_weno3(h, fc, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer      :: n, i, im1, ip1
    real(kind=8) :: bp, bm, ap, am, wp, wm, ww
    real(kind=8) :: dfl, dfr, fp, fm, h2

! selection weights
!
    real(kind=8), parameter :: dp = 2.0d+00 / 3.0d+00, dm = 1.0d+00 / 3.0d+00
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

! prepare common parameters
!
    h2 = h * h

! iterate along the vector
!
    do i = 2, n - 1

! prepare neighbour indices
!
      im1     = i - 1
      ip1     = i + 1

! calculate the left and right derivatives
!
      dfl     = fc(i  ) - fc(im1)
      dfr     = fc(ip1) - fc(i  )

! calculate coefficient omega
!
      ww      = (dfr - dfl)**2

! calculate corresponding betas
!
      bp      = dfr * dfr
      bm      = dfl * dfl

! calculate improved alphas
!
      ap      = 1.0d+00 + ww / (bp + h2)
      am      = 1.0d+00 + ww / (bm + h2)

! calculate weights
!
      wp      = dp * ap
      wm      = dm * am
      ww      = 2.0d+00 * (wp + wm)

! calculate central interpolation
!
      fp      =   fc(i  ) +           fc(ip1)

! calculate left side interpolation
!
      fm      = - fc(im1) + 3.0d+00 * fc(i  )

! calculate the left state
!
      fl( i ) = (wp * fp + wm * fm) / ww

! calculate weights
!
      wp      = dp * am
      wm      = dm * ap
      ww      = 2.0d+00 * (wp + wm)

! calculate central interpolation
!
      fp      =   fc(i  ) +           fc(im1)

! calculate right side interpolation
!
      fm      = - fc(ip1) + 3.0d+00 * fc(i  )

! calculate the right state
!
      fr(im1) = (wp * fp + wm * fm) / ww

    end do ! i = 2, n - 1

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_weno3
!
!===============================================================================
!
! subroutine RECONSTRUCT_LIMO3:
! ----------------------------
!
!   Subroutine reconstructs the interface states using the third order method
!   with a limiter function LimO3.
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Cada, M. & Torrilhon, M.,
!         "Compact third-order limiter functions for finite volume methods",
!         Journal of Computational Physics, 2009, 228, 4118-4145
!     [2] Mignone, A., Tzeferacos, P., & Bodo, G.,
!         "High-order conservative finite divergence GLM-MHD schemes for
!          cell-centered MHD",
!         Journal of Computational Physics, 2010, 229, 5896-5920
!
!===============================================================================
!
  subroutine reconstruct_limo3(h, fc, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer      :: n, i, im1, ip1
    real(kind=8) :: dfl, dfr
    real(kind=8) :: th, et, f1, f2, xl, xi, rdx, rdx2
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

! prepare parameters
!
    rdx   = rad * h
    rdx2  = rdx * rdx

! iterate over positions and interpolate states
!
    do i = 2, n - 1

! prepare neighbour indices
!
      im1 = i - 1
      ip1 = i + 1

! prepare left and right differences
!
      dfl = fc(i  ) - fc(im1)
      dfr = fc(ip1) - fc(i  )

! calculate the indicator function (eq. 3.17 in [1])
!
      et = (dfl * dfl + dfr * dfr) / rdx2

! the switching function (embedded in eq. 3.22 in [1], eq. 32 in [2])
!
      xi = max(0.0d+00, 0.5d+00 * min(2.0d+00, 1.0d+00 + (et - 1.0d+00) / eps))
      xl = 1.0d+00 - xi

! calculate values at i + ½
!
      if (abs(dfr) > eps) then

! calculate the slope ratio (eq. 2.8 in [1])
!
        th = dfl / dfr

! calculate the quadratic reconstruction (eq. 3.8 in [1], divided by 2)
!
        f1 = (2.0d+00 + th) / 6.0d+00

! calculate the third order limiter (eq. 3.13 in [1], cofficients divided by 2)
!
        if (th >= 0.0d+00) then
          f2 = max(0.0d+00, min(f1, th, 0.8d+00))
        else
          f2 = max(0.0d+00, min(f1, - 0.25d+00 * th))
        end if

! interpolate the left state (eq. 3.5 in [1], eq. 30 in [2])
!
        fl(i) = fc(i) + dfr * (xl * f1 + xi * f2)

      else

        fl(i) = fc(i)

      end if

! calculate values at i - ½
!
      if (abs(dfl) > eps) then

! calculate the slope ratio (eq. 2.8 in [1])
!
        th = dfr / dfl

! calculate the quadratic reconstruction (eq. 3.8 in [1], divided by 2)
!
        f1 = (2.0d+00 + th) / 6.0d+00

! calculate the third order limiter (eq. 3.13 in [1], cofficients divided by 2)
!
        if (th >= 0.0d+00) then
          f2 = max(0.0d+00, min(f1, th, 0.8d+00))
        else
          f2 = max(0.0d+00, min(f1, - 0.25d+00 * th))
        end if

! interpolate the right state (eq. 3.5 in [1], eq. 30 in [2])
!
        fr(im1) = fc(i) - dfl * (xl * f1 + xi * f2)

      else

        fr(im1) = fc(i)

      end if

    end do ! i = 2, n - 1

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_limo3
!
!===============================================================================
!
! subroutine RECONSTRUCT_PPM:
! --------------------------
!
!   Subroutine reconstructs the interface states using the third order
!   Piecewise Parabolic Method (PPM) with new limiters. This version lacks
!   the support for flattening important for keeping the spurious oscillations
!   near strong shocks/discontinuties under control.
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Colella, P., Woodward, P. R.,
!         "The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations",
!         Journal of Computational Physics, 1984, vol. 54, pp. 174-201,
!         https://dx.doi.org/10.1016/0021-9991(84)90143-8
!     [2] Colella, P., Sekora, M. D.,
!         "A limiter for PPM that preserves accuracy at smooth extrema",
!         Journal of Computational Physics, 2008, vol. 227, pp. 7069-7076,
!         https://doi.org/10.1016/j.jcp.2008.03.034
!
!===============================================================================
!
  subroutine reconstruct_ppm(h, fc, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    logical      :: cm, cp, ext
    integer      :: n, i, im1, ip1, im2
    real(kind=8) :: df2c, df2m, df2p, df2s, df2l
    real(kind=8) :: dfm, dfp, dcm, dcp
    real(kind=8) :: alm, alp, amx, sgn, del

! local arrays
!
    real(kind=8), dimension(size(fc)) :: fi, df2
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

! calculate the high-order interface interpolation; eq. (16) in [2]
!
    fi(1  ) = 5.0d-01 * (fc(1) + fc(2))
    do i = 2, n - 2
      fi(i) = (7.0d+00 * (fc(i) + fc(i+1)) - (fc(i-1) + fc(i+2))) / 1.2d+01
    end do
    fi(n-1) = 5.0d-01 * (fc(n-1) + fc(n))
    fi(n  ) = fc(n)

! calculate second-order central derivative
!
    df2(1) = 0.0d+00
    do i = 2, n - 1
      df2(i) = (fc(i+1) + fc(i-1)) - 2.0d+00 * fc(i)
    end do
    df2(n) = 0.0d+00

! limit the interpolation; eqs. (18) and (19) in [2]
!
    do i = 2, n - 2
      im1 = i - 1
      ip1 = i + 1

      if ((fc(ip1) - fi(i)) * (fi(i) - fc(i)) < 0.0d+00) then
        df2c = 3.0d+00 * ((fc(ip1) + fc(i  )) - 2.0d+00 * fi(i))
        df2m = df2(i  )
        df2p = df2(ip1)
        if (min(df2c, df2m, df2p) * max(df2c, df2m, df2p) > 0.0d+00) then
          df2l = sign(min(ppm_const * min(abs(df2m), abs(df2p)), abs(df2c)), df2c)
        else
          df2l = 0.0d+00
        end if
        fi(i) = 5.0d-01 * (fc(i) + fc(ip1)) - df2l / 6.0d+00
      end if
    end do

! iterate along the vector
!
    do i = 2, n - 1
      im1 = i - 1
      ip1 = i + 1

! limit states if local extremum is detected or the interpolation is not
! monotone
!
      alm = fi(im1) - fc(i)
      alp = fi(i  ) - fc(i)
      cm  = abs(alm) >= 2.0d+00 * abs(alp)
      cp  = abs(alp) >= 2.0d+00 * abs(alm)
      ext = .false.

! check if we have local extremum here
!
      if (alm * alp > 0.0d+00) then
        ext = .true.
      else if (cm .or. cp) then
        im2 = max(1, i - 1)

        dfm = fi(im1) - fi(im2)
        dfp = fi(ip1) - fi(i  )
        dcm = fc(i  ) - fc(im1)
        dcp = fc(ip1) - fc(i  )
        if (min(abs(dfm),abs(dfp)) >= min(abs(dcm),abs(dcp))) then
          ext = dfm * dfp < 0.0d+00
        else
          ext = dcm * dcp < 0.0d+00
        end if
      end if

! limit states if the local extremum is detected
!
      if (ext) then
        df2s = 6.0d+00 * (alm + alp)
        df2m = df2(im1)
        df2c = df2(i  )
        df2p = df2(ip1)
        if (min(df2s, df2c, df2m, df2p)                                        &
                      * max(df2s, df2c, df2m, df2p) > 0.0d+00) then
          df2l = sign(min(ppm_const * min(abs(df2m), abs(df2c), abs(df2p))     &
                                                   , abs(df2s)), df2s)
          if (abs(df2l) > 0.0d+00) then
            alm = alm * df2l / df2s
            alp = alp * df2l / df2s
          else
            alm = 0.0d+00
            alp = 0.0d+00
          end if
        else
          alm = 0.0d+00
          alp = 0.0d+00
        end if
      else

! the interpolation is not monotonic so apply additional limiter
!
        if (cp) then
          sgn = sign(1.0d+00, alp)
          amx = - 2.5d-01 * alp**2 / (alp + alm)
          dfp = fc(ip1) - fc(i)
          if (sgn * amx >= sgn * dfp) then
            del = dfp * (dfp - alm)
            if (del >= 0.0d+00) then
              alp = - 2.0d+00 * (dfp + sgn * sqrt(del))
            else
              alp = - 2.0d+00 * alm
            end if
          end if
        end if

        if (cm) then
          sgn = sign(1.0d+00, alm)
          amx = - 2.5d-01 * alm**2 / (alp + alm)
          dfm = fc(im1) - fc(i)
          if (sgn * amx >= sgn * dfm) then
            del = dfm * (dfm - alp)
            if (del >= 0.0d+00) then
              alm = - 2.0d+00 * (dfm + sgn * sqrt(del))
            else
              alm = - 2.0d+00 * alp
            end if
          end if
        end if
      end if

! update the states
!
      fr(im1) = fc(i) + alm
      fl(i  ) = fc(i) + alp

    end do ! i = 2, n

! update the interpolation of the first and last points
!
    fl(1  ) = fi(1)
    fr(n-1) = fi(n)
    fl(n  ) = fi(n)
    fr(n  ) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_ppm
!
!===============================================================================
!
! subroutine RECONSTRUCT_WENO5Z:
! -----------------------------
!
!   Subroutine reconstructs the interface states using the fifth order
!   Explicit Weighted Essentially Non-Oscillatory (WENO5) method with
!   stencil weights by Borges et al. (2008).
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Borges, R., Carmona, M., Costa, B., & Don, W.-S.,
!         "An improved weighted essentially non-oscillatory scheme for
!          hyperbolic conservation laws"
!         Journal of Computational Physics,
!         2008, vol. 227, pp. 3191-3211,
!         http://dx.doi.org/10.1016/j.jcp.2007.11.038
!
!===============================================================================
!
  subroutine reconstruct_weno5z(h, fc, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer      :: n, i, im1, ip1, im2, ip2
    real(kind=8) :: bl, bc, br, tt, df
    real(kind=8) :: al, ac, ar
    real(kind=8) :: wl, wc, wr, ww
    real(kind=8) :: ql, qc, qr

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc)) :: dfm, dfp, df2

! smoothness indicator coefficients
!
    real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 2.5d-01

! weight coefficients
!
    real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01

! interpolation coefficients
!
    real(kind=8), parameter :: a11 =   2.0d+00 / 6.0d+00                       &
                             , a12 = - 7.0d+00 / 6.0d+00                       &
                             , a13 =   1.1d+01 / 6.0d+00
    real(kind=8), parameter :: a21 = - 1.0d+00 / 6.0d+00                       &
                             , a22 =   5.0d+00 / 6.0d+00                       &
                             , a23 =   2.0d+00 / 6.0d+00
    real(kind=8), parameter :: a31 =   2.0d+00 / 6.0d+00                       &
                             , a32 =   5.0d+00 / 6.0d+00                       &
                             , a33 = - 1.0d+00 / 6.0d+00
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

! calculate the left and right derivatives
!
    do i = 1, n - 1
      ip1      = i + 1
      dfp(i  ) = fc(ip1) - fc(i)
      dfm(ip1) = dfp(i)
    end do
    dfm(1) = dfp(1)
    dfp(n) = dfm(n)

! calculate the absolute value of the second derivative
!
    df2(:) = c1 * (dfp(:) - dfm(:))**2

! iterate along the vector
!
    do i = 3, n - 2

! prepare neighbour indices
!
      im2 = i - 2
      im1 = i - 1
      ip1 = i + 1
      ip2 = i + 2

! calculate βₖ (eqs. 9-11 in [1])
!
      bl  = df2(im1) + c2 * (3.0d+00 * dfm(i  ) - dfm(im1))**2
      bc  = df2(i  ) + c2 * (          dfp(i  ) + dfm(i  ))**2
      br  = df2(ip1) + c2 * (3.0d+00 * dfp(i  ) - dfp(ip1))**2

! calculate τ (below eq. 25 in [1])
!
      tt  = abs(br - bl)

! calculate αₖ (eq. 28 in [1])
!
      al  = 1.0d+00 + tt / (bl + eps)
      ac  = 1.0d+00 + tt / (bc + eps)
      ar  = 1.0d+00 + tt / (br + eps)

! calculate weights
!
      wl  = dl * al
      wc  = dc * ac
      wr  = dr * ar
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the left state
!
      ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i  )
      qc = a21 * fc(im1) + a22 * fc(i  ) + a23 * fc(ip1)
      qr = a31 * fc(i  ) + a32 * fc(ip1) + a33 * fc(ip2)

! calculate the left state
!
      fl(i  ) = (wl * ql + wr * qr) + wc * qc

! normalize weights
!
      wl  = dl * ar
      wc  = dc * ac
      wr  = dr * al
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the right state
!
      ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i  )
      qc = a21 * fc(ip1) + a22 * fc(i  ) + a23 * fc(im1)
      qr = a31 * fc(i  ) + a32 * fc(im1) + a33 * fc(im2)

! calculate the right state
!
      fr(im1) = (wl * ql + wr * qr) + wc * qc

    end do ! i = 3, n - 2

! update the interpolation of the first and last two points
!
    fl(1)   = 0.5d+00 * (fc(1) + fc(2))
    df      = limiter_tvd(0.5d+00, dfm(2), dfp(2))
    fr(1)   = fc(2) - df
    fl(2)   = fc(2) + df
    i       = n - 1
    df      = limiter_tvd(0.5d+00, dfm(i), dfp(i))
    fr(i-1) = fc(i) - df
    fl(i)   = fc(i) + df
    fr(i)   = 0.5d+00 * (fc(i) + fc(n))
    fl(n)   = fc(n)
    fr(n)   = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_weno5z
!
!===============================================================================
!
! subroutine RECONSTRUCT_WENO5YC:
! ------------------------------
!
!   Subroutine reconstructs the interface states using the fifth order
!   Explicit Weighted Essentially Non-Oscillatory (WENO5) method with
!   stencil weights by Yamaleev & Carpenter (2009).
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Yamaleev, N. K. & Carpenter, H. C.,
!         "A Systematic Methodology for Constructing High-Order Energy Stable
!          WENO Schemes"
!         Journal of Computational Physics,
!         2009, vol. 228, pp. 4248-4272,
!         http://dx.doi.org/10.1016/j.jcp.2009.03.002
!
!===============================================================================
!
  subroutine reconstruct_weno5yc(h, fc, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer      :: n, i, im1, ip1, im2, ip2
    real(kind=8) :: bl, bc, br, tt, df
    real(kind=8) :: al, ac, ar
    real(kind=8) :: wl, wc, wr, ww
    real(kind=8) :: ql, qc, qr

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc)) :: dfm, dfp, df2

! smoothness indicator coefficients
!
    real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 2.5d-01

! weight coefficients
!
    real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01

! interpolation coefficients
!
    real(kind=8), parameter :: a11 =   2.0d+00 / 6.0d+00                       &
                             , a12 = - 7.0d+00 / 6.0d+00                       &
                             , a13 =   1.1d+01 / 6.0d+00
    real(kind=8), parameter :: a21 = - 1.0d+00 / 6.0d+00                       &
                             , a22 =   5.0d+00 / 6.0d+00                       &
                             , a23 =   2.0d+00 / 6.0d+00
    real(kind=8), parameter :: a31 =   2.0d+00 / 6.0d+00                       &
                             , a32 =   5.0d+00 / 6.0d+00                       &
                             , a33 = - 1.0d+00 / 6.0d+00
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

! calculate the left and right derivatives
!
    do i = 1, n - 1
      ip1      = i + 1
      dfp(i  ) = fc(ip1) - fc(i)
      dfm(ip1) = dfp(i)
    end do
    dfm(1) = dfp(1)
    dfp(n) = dfm(n)

! calculate the absolute value of the second derivative
!
    df2(:) = c1 * (dfp(:) - dfm(:))**2

! iterate along the vector
!
    do i = 3, n - 2

! prepare neighbour indices
!
      im2 = i - 2
      im1 = i - 1
      ip1 = i + 1
      ip2 = i + 2

! calculate βₖ (eq. 19 in [1])
!
      bl  = df2(im1) + c2 * (3.0d+00 * dfm(i  ) - dfm(im1))**2
      bc  = df2(i  ) + c2 * (          dfp(i  ) + dfm(i  ))**2
      br  = df2(ip1) + c2 * (3.0d+00 * dfp(i  ) - dfp(ip1))**2

! calculate τ (below eq. 20 in [1])
!
      tt  = (6.0d+00 * fc(i) - 4.0d+00 * (fc(im1) + fc(ip1))                   &
                                       + (fc(im2) + fc(ip2)))**2

! calculate αₖ (eqs. 18 or 58 in [1])
!
      al  = 1.0d+00 + tt / (bl + eps)
      ac  = 1.0d+00 + tt / (bc + eps)
      ar  = 1.0d+00 + tt / (br + eps)

! calculate weights
!
      wl  = dl * al
      wc  = dc * ac
      wr  = dr * ar
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the left state
!
      ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i  )
      qc = a21 * fc(im1) + a22 * fc(i  ) + a23 * fc(ip1)
      qr = a31 * fc(i  ) + a32 * fc(ip1) + a33 * fc(ip2)

! calculate the left state
!
      fl(i  ) = (wl * ql + wr * qr) + wc * qc

! normalize weights
!
      wl  = dl * ar
      wc  = dc * ac
      wr  = dr * al
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the right state
!
      ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i  )
      qc = a21 * fc(ip1) + a22 * fc(i  ) + a23 * fc(im1)
      qr = a31 * fc(i  ) + a32 * fc(im1) + a33 * fc(im2)

! calculate the right state
!
      fr(im1) = (wl * ql + wr * qr) + wc * qc

    end do ! i = 3, n - 2

! update the interpolation of the first and last two points
!
    fl(1)   = 0.5d+00 * (fc(1) + fc(2))
    df      = limiter_tvd(0.5d+00, dfm(2), dfp(2))
    fr(1)   = fc(2) - df
    fl(2)   = fc(2) + df
    i       = n - 1
    df      = limiter_tvd(0.5d+00, dfm(i), dfp(i))
    fr(i-1) = fc(i) - df
    fl(i)   = fc(i) + df
    fr(i)   = 0.5d+00 * (fc(i) + fc(n))
    fl(n)   = fc(n)
    fr(n)   = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_weno5yc
!
!===============================================================================
!
! subroutine RECONSTRUCT_WENO5NS:
! ------------------------------
!
!   Subroutine reconstructs the interface states using the fifth order
!   Explicit Weighted Essentially Non-Oscillatory (WENO5) method with new
!   smoothness indicators and stencil weights by Ha et al. (2013).
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Ha, Y., Kim, C. H., Lee, Y. J., & Yoon, J.,
!         "An improved weighted essentially non-oscillatory scheme with a new
!          smoothness indicator",
!         Journal of Computational Physics,
!         2013, vol. 232, pp. 68-86
!         http://dx.doi.org/10.1016/j.jcp.2012.06.016
!
!===============================================================================
!
  subroutine reconstruct_weno5ns(h, fc, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer      :: n, i, im1, ip1, im2, ip2
    real(kind=8) :: bl, bc, br
    real(kind=8) :: al, ac, ar, aa
    real(kind=8) :: wl, wc, wr
    real(kind=8) :: df, lq, l3, zt
    real(kind=8) :: ql, qc, qr

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc)) :: dfm, dfp, df2

! weight coefficients
!
    real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01

! interpolation coefficients
!
    real(kind=8), parameter :: a11 =   2.0d+00 / 6.0d+00                       &
                             , a12 = - 7.0d+00 / 6.0d+00                       &
                             , a13 =   1.1d+01 / 6.0d+00
    real(kind=8), parameter :: a21 = - 1.0d+00 / 6.0d+00                       &
                             , a22 =   5.0d+00 / 6.0d+00                       &
                             , a23 =   2.0d+00 / 6.0d+00
    real(kind=8), parameter :: a31 =   2.0d+00 / 6.0d+00                       &
                             , a32 =   5.0d+00 / 6.0d+00                       &
                             , a33 = - 1.0d+00 / 6.0d+00

! the free parameter for smoothness indicators (see Eq. 3.6 in [1])
!
    real(kind=8), parameter :: xi  =   4.0d-01
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

! calculate the left and right derivatives
!
    do i = 1, n - 1
      ip1      = i + 1
      dfp(i  ) = fc(ip1) - fc(i)
      dfm(ip1) = dfp(i)
    end do
    dfm(1) = dfp(1)
    dfp(n) = dfm(n)

! calculate the absolute value of the second derivative
!
    df2(:) = 0.5d+00 * abs(dfp(:) - dfm(:))

! iterate along the vector
!
    do i = 3, n - 2

! prepare neighbour indices
!
      im2 = i - 2
      im1 = i - 1
      ip1 = i + 1
      ip2 = i + 2

! calculate βₖ (eq. 3.6 in [1])
!
      df  = abs(dfp(i))
      lq  = xi * df
      bl  = df2(im1) + xi * abs(2.0d+00 * dfm(i) - dfm(im1))
      bc  = df2(i  ) + lq
      br  = df2(ip1) + lq

! calculate ζ (below eq. 3.6 in [1])
!
      l3  = df**3
      zt  = 0.5d+00 * ((bl - br)**2 + (l3 / (1.0d+00 + l3))**2)

! calculate αₖ (eq. 3.9 in [4])
!
      al  = dl * (1.0d+00 + zt / (bl + eps)**2)
      ac  = dc * (1.0d+00 + zt / (bc + eps)**2)
      ar  = dr * (1.0d+00 + zt / (br + eps)**2)

! calculate weights
!
      aa  = (al + ar) + ac
      wl  = al / aa
      wr  = ar / aa
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the left state
!
      ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i  )
      qc = a21 * fc(im1) + a22 * fc(i  ) + a23 * fc(ip1)
      qr = a31 * fc(i  ) + a32 * fc(ip1) + a33 * fc(ip2)

! calculate the left state
!
      fl(i  ) = (wl * ql + wr * qr) + wc * qc

! calculate βₖ (eq. 3.6 in [1])
!
      df  = abs(dfm(i))
      lq  = xi * df
      bl  = df2(ip1) + xi * abs(2.0d+00 * dfp(i) - dfp(ip1))
      bc  = df2(i  ) + lq
      br  = df2(im1) + lq

! calculate ζ (below eq. 3.6 in [1])

      l3  = df**3
      zt  = 0.5d+00 * ((bl - br)**2 + (l3 / (1.0d+00 + l3))**2)

! calculate αₖ (eq. 3.9 in [4])
!
      al  = dl * (1.0d+00 + zt / (bl + eps)**2)
      ac  = dc * (1.0d+00 + zt / (bc + eps)**2)
      ar  = dr * (1.0d+00 + zt / (br + eps)**2)

! normalize weights
!
      aa  = (al + ar) + ac
      wl  = al / aa
      wr  = ar / aa
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the right state
!
      ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i  )
      qc = a21 * fc(ip1) + a22 * fc(i  ) + a23 * fc(im1)
      qr = a31 * fc(i  ) + a32 * fc(im1) + a33 * fc(im2)

! calculate the right state
!
      fr(im1) = (wl * ql + wr * qr) + wc * qc

    end do ! i = 3, n - 2

! update the interpolation of the first and last two points
!
    fl(1)   = 0.5d+00 * (fc(1) + fc(2))
    df      = limiter_tvd(0.5d+00, dfm(2), dfp(2))
    fr(1)   = fc(2) - df
    fl(2)   = fc(2) + df
    i       = n - 1
    df      = limiter_tvd(0.5d+00, dfm(i), dfp(i))
    fr(i-1) = fc(i) - df
    fl(i)   = fc(i) + df
    fr(i)   = 0.5d+00 * (fc(i) + fc(n))
    fl(n)   = fc(n)
    fr(n)   = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_weno5ns
!
!===============================================================================
!
! subroutine RECONSTRUCT_CRWENO5Z:
! -------------------------------
!
!   Subroutine reconstructs the interface states using the fifth order
!   Compact-Reconstruction Weighted Essentially Non-Oscillatory (CRWENO)
!   method and smoothness indicators by Borges et al. (2008).
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Ghosh, D. & Baeder, J. D.,
!         "Compact Reconstruction Schemes with Weighted ENO Limiting for
!          Hyperbolic Conservation Laws"
!         SIAM Journal on Scientific Computing,
!         2012, vol. 34, no. 3, pp. A1678-A1706,
!         http://dx.doi.org/10.1137/110857659
!     [2] Ghosh, D. & Baeder, J. D.,
!         "Weighted Non-linear Compact Schemes for the Direct Numerical
!          Simulation of Compressible, Turbulent Flows"
!         Journal on Scientific Computing,
!         2014,
!         http://dx.doi.org/10.1007/s10915-014-9818-0
!     [3] Borges, R., Carmona, M., Costa, B., & Don, W.-S.,
!         "An improved weighted essentially non-oscillatory scheme for
!          hyperbolic conservation laws"
!         Journal of Computational Physics,
!         2008, vol. 227, pp. 3191-3211,
!         http://dx.doi.org/10.1016/j.jcp.2007.11.038
!
!===============================================================================
!
  subroutine reconstruct_crweno5z(h, fc, fl, fr)

! include external procedures
!
    use algebra   , only : tridiag

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer      :: n, i, im1, ip1, im2, ip2
    real(kind=8) :: bl, bc, br, tt
    real(kind=8) :: wl, wc, wr, ww
    real(kind=8) :: ql, qc, qr

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc))   :: dfm, dfp, df2
    real(kind=8), dimension(size(fc))   :: al, ac, ar
    real(kind=8), dimension(size(fc))   :: u
    real(kind=8), dimension(size(fc),2) :: a, b, c, r

! smoothness indicator coefficients
!
    real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 2.5d-01

! weight coefficients for implicit (c) and explicit (d) interpolations
!
    real(kind=8), parameter :: cl = 1.0d+00 / 9.0d+00
    real(kind=8), parameter :: cc = 5.0d+00 / 9.0d+00
    real(kind=8), parameter :: cr = 1.0d+00 / 3.0d+00
    real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01

! implicit method coefficients
!
    real(kind=8), parameter :: dq = 5.0d-01

! interpolation coefficients
!
    real(kind=8), parameter :: a11 =   2.0d+00 / 6.0d+00                       &
                             , a12 = - 7.0d+00 / 6.0d+00                       &
                             , a13 =   1.1d+01 / 6.0d+00
    real(kind=8), parameter :: a21 = - 1.0d+00 / 6.0d+00                       &
                             , a22 =   5.0d+00 / 6.0d+00                       &
                             , a23 =   2.0d+00 / 6.0d+00
    real(kind=8), parameter :: a31 =   2.0d+00 / 6.0d+00                       &
                             , a32 =   5.0d+00 / 6.0d+00                       &
                             , a33 = - 1.0d+00 / 6.0d+00
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

! calculate the left and right derivatives
!
    do i = 1, n - 1
      ip1      = i + 1
      dfp(i  ) = fc(ip1) - fc(i)
      dfm(ip1) = dfp(i)
    end do
    dfm(1) = dfp(1)
    dfp(n) = dfm(n)

! calculate the absolute value of the second derivative
!
    df2(:) = c1 * (dfp(:) - dfm(:))**2

! prepare smoothness indicators
!
    do i = 2, n - 1

! prepare neighbour indices
!
      im1   = i - 1
      ip1   = i + 1

! calculate βₖ (eqs. 9-11 in [1])
!
      bl    = df2(im1) + c2 * (3.0d+00 * dfm(i  ) - dfm(im1))**2
      bc    = df2(i  ) + c2 * (          dfp(i  ) + dfm(i  ))**2
      br    = df2(ip1) + c2 * (3.0d+00 * dfp(i  ) - dfp(ip1))**2

! calculate τ (below eq. 25 in [1])
!
      tt    = abs(br - bl)

! calculate αₖ (eq. 28 in [1])
!
      al(i) = 1.0d+00 + tt / (bl + eps)
      ac(i) = 1.0d+00 + tt / (bc + eps)
      ar(i) = 1.0d+00 + tt / (br + eps)

    end do ! i = 2, n - 1

! prepare tridiagonal system coefficients
!
    do i = ng, n - ng + 1

! prepare neighbour indices
!
      im1 = i - 1
      ip1 = i + 1

! calculate weights
!
      wl  = cl * al(i)
      wc  = cc * ac(i)
      wr  = cr * ar(i)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate tridiagonal matrix coefficients
!
      a(i,1) = 2.0d+00 * wl +            wc
      b(i,1) =           wl + 2.0d+00 * (wc + wr)
      c(i,1) =           wr

! prepare right hand side of tridiagonal equation
!
      r(i,1) = (wl * fc(im1) + (5.0d+00 * (wl + wc) + wr) * fc(i  )            &
                                    + (wc + 5.0d+00 * wr) * fc(ip1)) * dq

! calculate weights
!
      wl  = cl * ar(i)
      wc  = cc * ac(i)
      wr  = cr * al(i)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate tridiagonal matrix coefficients
!
      a(i,2) =           wr
      b(i,2) =           wl + 2.0d+00 * (wc + wr)
      c(i,2) = 2.0d+00 * wl +            wc

! prepare right hand side of tridiagonal equation
!
      r(i,2) = (wl * fc(ip1) + (5.0d+00 * (wl + wc) + wr) * fc(i  )             &
                                    + (wc + 5.0d+00 * wr) * fc(im1)) * dq

    end do ! i = ng, n - ng + 1

! interpolate ghost zones using explicit solver (left-side reconstruction)
!
    do i = 2, ng

! prepare neighbour indices
!
      im2 = max(1, i - 2)
      im1 = i - 1
      ip1 = i + 1
      ip2 = i + 2

! calculate weights
!
      wl  = dl * al(i)
      wc  = dc * ac(i)
      wr  = dr * ar(i)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the left state
!
      ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i  )
      qc = a21 * fc(im1) + a22 * fc(i  ) + a23 * fc(ip1)
      qr = a31 * fc(i  ) + a32 * fc(ip1) + a33 * fc(ip2)

! calculate the left state
!
      fl(i) = (wl * ql + wr * qr) + wc * qc

! prepare coefficients of the tridiagonal system
!
      a(i,1) = 0.0d+00
      b(i,1) = 1.0d+00
      c(i,1) = 0.0d+00
      r(i,1) = fl(i)

    end do ! i = 2, ng
    a(1,1) = 0.0d+00
    b(1,1) = 1.0d+00
    c(1,1) = 0.0d+00
    r(1,1) = 0.5d+00 * (fc(1) + fc(2))

! interpolate ghost zones using explicit solver (left-side reconstruction)
!
    do i = n - ng, n - 1

! prepare neighbour indices
!
      im2 = i - 2
      im1 = i - 1
      ip1 = i + 1
      ip2 = min(n, i + 2)

! calculate weights
!
      wl  = dl * al(i)
      wc  = dc * ac(i)
      wr  = dr * ar(i)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the left state
!
      ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i  )
      qc = a21 * fc(im1) + a22 * fc(i  ) + a23 * fc(ip1)
      qr = a31 * fc(i  ) + a32 * fc(ip1) + a33 * fc(ip2)

! calculate the left state
!
      fl(i) = (wl * ql + wr * qr) + wc * qc

! prepare coefficients of the tridiagonal system
!
      a(i,1) = 0.0d+00
      b(i,1) = 1.0d+00
      c(i,1) = 0.0d+00
      r(i,1) = fl(i)

    end do ! i = n - ng, n - 1
    a(n,1) = 0.0d+00
    b(n,1) = 1.0d+00
    c(n,1) = 0.0d+00
    r(n,1) = fc(n)

! interpolate ghost zones using explicit solver (right-side reconstruction)
!
    do i = 2, ng + 1

! prepare neighbour indices
!
      im2 = max(1, i - 2)
      im1 = i - 1
      ip1 = i + 1
      ip2 = i + 2

! normalize weights
!
      wl  = dl * ar(i)
      wc  = dc * ac(i)
      wr  = dr * al(i)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the right state
!
      ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i  )
      qc = a21 * fc(ip1) + a22 * fc(i  ) + a23 * fc(im1)
      qr = a31 * fc(i  ) + a32 * fc(im1) + a33 * fc(im2)

! calculate the right state
!
      fr(i) = (wl * ql + wr * qr) + wc * qc

! prepare coefficients of the tridiagonal system
!
      a(i,2) = 0.0d+00
      b(i,2) = 1.0d+00
      c(i,2) = 0.0d+00
      r(i,2) = fr(i)

    end do ! i = 2, ng + 1
    a(1,2) = 0.0d+00
    b(1,2) = 1.0d+00
    c(1,2) = 0.0d+00
    r(1,2) = fc(1)

! interpolate ghost zones using explicit solver (right-side reconstruction)
!
    do i = n - ng + 1, n - 1

! prepare neighbour indices
!
      im2 = max(1, i - 2)
      im1 = max(1, i - 1)
      ip1 = min(n, i + 1)
      ip2 = min(n, i + 2)

! normalize weights
!
      wl  = dl * ar(i)
      wc  = dc * ac(i)
      wr  = dr * al(i)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the right state
!
      ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i  )
      qc = a21 * fc(ip1) + a22 * fc(i  ) + a23 * fc(im1)
      qr = a31 * fc(i  ) + a32 * fc(im1) + a33 * fc(im2)

! calculate the right state
!
      fr(i) = (wl * ql + wr * qr) + wc * qc

! prepare coefficients of the tridiagonal system
!
      a(i,2) = 0.0d+00
      b(i,2) = 1.0d+00
      c(i,2) = 0.0d+00
      r(i,2) = fr(i)

    end do ! i = n - ng + 1, n - 1
    a(n,2) = 0.0d+00
    b(n,2) = 1.0d+00
    c(n,2) = 0.0d+00
    r(n,2) = 0.5d+00 * (fc(n-1) + fc(n))

! solve the tridiagonal system of equations for the left-side interpolation
!
    call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n))

! substitute the left-side values
!
    fl(1:n  ) = u(1:n)

! solve the tridiagonal system of equations for the left-side interpolation
!
    call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n))

! substitute the right-side values
!
    fr(1:n-1) = u(2:n)

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_crweno5z
!
!===============================================================================
!
! subroutine RECONSTRUCT_CRWENO5YC:
! --------------------------------
!
!   Subroutine reconstructs the interface states using the fifth order
!   Compact-Reconstruction Weighted Essentially Non-Oscillatory (CRWENO)
!   method and smoothness indicators by Yamaleev & Carpenter (2009).
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Ghosh, D. & Baeder, J. D.,
!         "Compact Reconstruction Schemes with Weighted ENO Limiting for
!          Hyperbolic Conservation Laws"
!         SIAM Journal on Scientific Computing,
!         2012, vol. 34, no. 3, pp. A1678-A1706,
!         http://dx.doi.org/10.1137/110857659
!     [2] Ghosh, D. & Baeder, J. D.,
!         "Weighted Non-linear Compact Schemes for the Direct Numerical
!          Simulation of Compressible, Turbulent Flows"
!         Journal on Scientific Computing,
!         2014,
!         http://dx.doi.org/10.1007/s10915-014-9818-0
!     [3] Yamaleev, N. K. & Carpenter, H. C.,
!         "A Systematic Methodology for Constructing High-Order Energy Stable
!          WENO Schemes"
!         Journal of Computational Physics,
!         2009, vol. 228, pp. 4248-4272,
!         http://dx.doi.org/10.1016/j.jcp.2009.03.002
!
!===============================================================================
!
  subroutine reconstruct_crweno5yc(h, fc, fl, fr)

! include external procedures
!
    use algebra   , only : tridiag

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer      :: n, i, im1, ip1, im2, ip2
    real(kind=8) :: bl, bc, br, tt
    real(kind=8) :: wl, wc, wr, ww
    real(kind=8) :: ql, qc, qr

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc))   :: dfm, dfp, df2
    real(kind=8), dimension(size(fc))   :: al, ac, ar
    real(kind=8), dimension(size(fc))   :: u
    real(kind=8), dimension(size(fc),2) :: a, b, c, r

! smoothness indicator coefficients
!
    real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 2.5d-01

! weight coefficients for implicit (c) and explicit (d) interpolations
!
    real(kind=8), parameter :: cl = 1.0d+00 / 9.0d+00
    real(kind=8), parameter :: cc = 5.0d+00 / 9.0d+00
    real(kind=8), parameter :: cr = 1.0d+00 / 3.0d+00
    real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01

! implicit method coefficients
!
    real(kind=8), parameter :: dq = 5.0d-01

! interpolation coefficients
!
    real(kind=8), parameter :: a11 =   2.0d+00 / 6.0d+00                       &
                             , a12 = - 7.0d+00 / 6.0d+00                       &
                             , a13 =   1.1d+01 / 6.0d+00
    real(kind=8), parameter :: a21 = - 1.0d+00 / 6.0d+00                       &
                             , a22 =   5.0d+00 / 6.0d+00                       &
                             , a23 =   2.0d+00 / 6.0d+00
    real(kind=8), parameter :: a31 =   2.0d+00 / 6.0d+00                       &
                             , a32 =   5.0d+00 / 6.0d+00                       &
                             , a33 = - 1.0d+00 / 6.0d+00
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

! calculate the left and right derivatives
!
    do i = 1, n - 1
      ip1      = i + 1
      dfp(i  ) = fc(ip1) - fc(i)
      dfm(ip1) = dfp(i)
    end do
    dfm(1) = dfp(1)
    dfp(n) = dfm(n)

! calculate the absolute value of the second derivative
!
    df2(:) = c1 * (dfp(:) - dfm(:))**2

! prepare smoothness indicators
!
    do i = 2, n - 1

! prepare neighbour indices
!
      im2   = max(1, i - 2)
      im1   = i - 1
      ip1   = i + 1
      ip2   = min(n, i + 2)

! calculate βₖ (eqs. 9-11 in [1])
!
      bl    = df2(im1) + c2 * (3.0d+00 * dfm(i  ) - dfm(im1))**2
      bc    = df2(i  ) + c2 * (          dfp(i  ) + dfm(i  ))**2
      br    = df2(ip1) + c2 * (3.0d+00 * dfp(i  ) - dfp(ip1))**2

! calculate τ (below eq. 64 in [3])
!
      tt  = (6.0d+00 * fc(i) + (fc(im2) + fc(ip2))                             &
                   - 4.0d+00 * (fc(im1) + fc(ip1)))**2

! calculate αₖ (eq. 28 in [1])
!
      al(i) = 1.0d+00 + tt / (bl + eps)
      ac(i) = 1.0d+00 + tt / (bc + eps)
      ar(i) = 1.0d+00 + tt / (br + eps)

    end do ! i = 2, n - 1

! prepare tridiagonal system coefficients
!
    do i = ng, n - ng + 1

! prepare neighbour indices
!
      im1 = i - 1
      ip1 = i + 1

! calculate weights
!
      wl  = cl * al(i)
      wc  = cc * ac(i)
      wr  = cr * ar(i)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate tridiagonal matrix coefficients
!
      a(i,1) = 2.0d+00 * wl +            wc
      b(i,1) =           wl + 2.0d+00 * (wc + wr)
      c(i,1) =           wr

! prepare right hand side of tridiagonal equation
!
      r(i,1) = (wl * fc(im1) + (5.0d+00 * (wl + wc) + wr) * fc(i  )            &
                                    + (wc + 5.0d+00 * wr) * fc(ip1)) * dq

! calculate weights
!
      wl  = cl * ar(i)
      wc  = cc * ac(i)
      wr  = cr * al(i)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate tridiagonal matrix coefficients
!
      a(i,2) =           wr
      b(i,2) =           wl + 2.0d+00 * (wc + wr)
      c(i,2) = 2.0d+00 * wl +            wc

! prepare right hand side of tridiagonal equation
!
      r(i,2) = (wl * fc(ip1) + (5.0d+00 * (wl + wc) + wr) * fc(i  )            &
                                    + (wc + 5.0d+00 * wr) * fc(im1)) * dq

    end do ! i = ng, n - ng + 1

! interpolate ghost zones using explicit solver (left-side reconstruction)
!
    do i = 2, ng

! prepare neighbour indices
!
      im2 = max(1, i - 2)
      im1 = i - 1
      ip1 = i + 1
      ip2 = i + 2

! calculate weights
!
      wl  = dl * al(i)
      wc  = dc * ac(i)
      wr  = dr * ar(i)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the left state
!
      ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i  )
      qc = a21 * fc(im1) + a22 * fc(i  ) + a23 * fc(ip1)
      qr = a31 * fc(i  ) + a32 * fc(ip1) + a33 * fc(ip2)

! calculate the left state
!
      fl(i) = (wl * ql + wr * qr) + wc * qc

! prepare coefficients of the tridiagonal system
!
      a(i,1) = 0.0d+00
      b(i,1) = 1.0d+00
      c(i,1) = 0.0d+00
      r(i,1) = fl(i)

    end do ! i = 2, ng
    a(1,1) = 0.0d+00
    b(1,1) = 1.0d+00
    c(1,1) = 0.0d+00
    r(1,1) = 0.5d+00 * (fc(1) + fc(2))

! interpolate ghost zones using explicit solver (left-side reconstruction)
!
    do i = n - ng, n - 1

! prepare neighbour indices
!
      im2 = i - 2
      im1 = i - 1
      ip1 = i + 1
      ip2 = min(n, i + 2)

! calculate weights
!
      wl  = dl * al(i)
      wc  = dc * ac(i)
      wr  = dr * ar(i)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the left state
!
      ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i  )
      qc = a21 * fc(im1) + a22 * fc(i  ) + a23 * fc(ip1)
      qr = a31 * fc(i  ) + a32 * fc(ip1) + a33 * fc(ip2)

! calculate the left state
!
      fl(i) = (wl * ql + wr * qr) + wc * qc

! prepare coefficients of the tridiagonal system
!
      a(i,1) = 0.0d+00
      b(i,1) = 1.0d+00
      c(i,1) = 0.0d+00
      r(i,1) = fl(i)

    end do ! i = n - ng, n - 1
    a(n,1) = 0.0d+00
    b(n,1) = 1.0d+00
    c(n,1) = 0.0d+00
    r(n,1) = fc(n)

! interpolate ghost zones using explicit solver (right-side reconstruction)
!
    do i = 2, ng + 1

! prepare neighbour indices
!
      im2 = max(1, i - 2)
      im1 = i - 1
      ip1 = i + 1
      ip2 = i + 2

! normalize weights
!
      wl  = dl * ar(i)
      wc  = dc * ac(i)
      wr  = dr * al(i)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the right state
!
      ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i  )
      qc = a21 * fc(ip1) + a22 * fc(i  ) + a23 * fc(im1)
      qr = a31 * fc(i  ) + a32 * fc(im1) + a33 * fc(im2)

! calculate the right state
!
      fr(i) = (wl * ql + wr * qr) + wc * qc

! prepare coefficients of the tridiagonal system
!
      a(i,2) = 0.0d+00
      b(i,2) = 1.0d+00
      c(i,2) = 0.0d+00
      r(i,2) = fr(i)

    end do ! i = 2, ng + 1
    a(1,2) = 0.0d+00
    b(1,2) = 1.0d+00
    c(1,2) = 0.0d+00
    r(1,2) = fc(1)

! interpolate ghost zones using explicit solver (right-side reconstruction)
!
    do i = n - ng + 1, n - 1

! prepare neighbour indices
!
      im2 = i - 2
      im1 = i - 1
      ip1 = i + 1
      ip2 = min(n, i + 2)

! normalize weights
!
      wl  = dl * ar(i)
      wc  = dc * ac(i)
      wr  = dr * al(i)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the right state
!
      ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i  )
      qc = a21 * fc(ip1) + a22 * fc(i  ) + a23 * fc(im1)
      qr = a31 * fc(i  ) + a32 * fc(im1) + a33 * fc(im2)

! calculate the right state
!
      fr(i) = (wl * ql + wr * qr) + wc * qc

! prepare coefficients of the tridiagonal system
!
      a(i,2) = 0.0d+00
      b(i,2) = 1.0d+00
      c(i,2) = 0.0d+00
      r(i,2) = fr(i)

    end do ! i = n - ng + 1, n - 1
    a(n,2) = 0.0d+00
    b(n,2) = 1.0d+00
    c(n,2) = 0.0d+00
    r(n,2) = 0.5d+00 * (fc(n-1) + fc(n))

! solve the tridiagonal system of equations for the left-side interpolation
!
    call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n))

! substitute the left-side values
!
    fl(1:n  ) = u(1:n)

! solve the tridiagonal system of equations for the left-side interpolation
!
    call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n))

! substitute the right-side values
!
    fr(1:n-1) = u(2:n)

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_crweno5yc
!
!===============================================================================
!
! subroutine RECONSTRUCT_CRWENO5NS:
! --------------------------------
!
!   Subroutine reconstructs the interface states using the fifth order
!   Compact-Reconstruction Weighted Essentially Non-Oscillatory (CRWENO)
!   method combined with the smoothness indicators by Ha et al. (2013).
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Ghosh, D. & Baeder, J. D.,
!         "Compact Reconstruction Schemes with Weighted ENO Limiting for
!          Hyperbolic Conservation Laws"
!         SIAM Journal on Scientific Computing,
!         2012, vol. 34, no. 3, pp. A1678-A1706,
!         http://dx.doi.org/10.1137/110857659
!     [2] Ghosh, D. & Baeder, J. D.,
!         "Weighted Non-linear Compact Schemes for the Direct Numerical
!          Simulation of Compressible, Turbulent Flows"
!         Journal on Scientific Computing,
!         2014,
!         http://dx.doi.org/10.1007/s10915-014-9818-0
!     [3] Ha, Y., Kim, C. H., Lee, Y. J., & Yoon, J.,
!         "An improved weighted essentially non-oscillatory scheme with a new
!          smoothness indicator",
!         Journal of Computational Physics,
!         2013, vol. 232, pp. 68-86
!         http://dx.doi.org/10.1016/j.jcp.2012.06.016
!
!===============================================================================
!
  subroutine reconstruct_crweno5ns(h, fc, fl, fr)

! include external procedures
!
    use algebra   , only : tridiag

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer      :: n, i, im1, ip1, im2, ip2
    real(kind=8) :: bl, bc, br
    real(kind=8) :: wl, wc, wr, ww
    real(kind=8) :: df, lq, l3, zt
    real(kind=8) :: ql, qc, qr

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc))   :: dfm, dfp, df2
    real(kind=8), dimension(size(fc),2) :: al, ac, ar
    real(kind=8), dimension(size(fc))   :: u
    real(kind=8), dimension(size(fc),2) :: a, b, c, r

! the free parameter for smoothness indicators (see eq. 3.6 in [3])
!
    real(kind=8), parameter :: xi = 4.0d-01

! weight coefficients for implicit (c) and explicit (d) interpolations
!
    real(kind=8), parameter :: cl = 1.0d+00 / 9.0d+00
    real(kind=8), parameter :: cc = 5.0d+00 / 9.0d+00
    real(kind=8), parameter :: cr = 1.0d+00 / 3.0d+00
    real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01

! implicit method coefficients
!
    real(kind=8), parameter :: dq = 5.0d-01

! 3rd order interpolation coefficients for three stencils
!
    real(kind=8), parameter :: a11 =   2.0d+00 / 6.0d+00                       &
                             , a12 = - 7.0d+00 / 6.0d+00                       &
                             , a13 =   1.1d+01 / 6.0d+00
    real(kind=8), parameter :: a21 = - 1.0d+00 / 6.0d+00                       &
                             , a22 =   5.0d+00 / 6.0d+00                       &
                             , a23 =   2.0d+00 / 6.0d+00
    real(kind=8), parameter :: a31 =   2.0d+00 / 6.0d+00                       &
                             , a32 =   5.0d+00 / 6.0d+00                       &
                             , a33 = - 1.0d+00 / 6.0d+00
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

! calculate the left and right derivatives
!
    do i = 1, n - 1
      ip1      = i + 1
      dfp(i  ) = fc(ip1) - fc(i)
      dfm(ip1) = dfp(i)
    end do
    dfm(1) = dfp(1)
    dfp(n) = dfm(n)

! calculate the absolute value of the second derivative
!
    df2(:) = 0.5d+00 * abs(dfp(:) - dfm(:))

! prepare smoothness indicators
!
    do i = 2, n - 1

! prepare neighbour indices
!
      im1  = i - 1
      ip1  = i + 1

! calculate βₖ
!
      df  = abs(dfp(i))
      lq  = xi * df
      bl  = df2(im1) + xi * abs(2.0d+00 * dfm(i) - dfm(im1))
      bc  = df2(i  ) + lq
      br  = df2(ip1) + lq

! calculate ζ
!
      l3  = df**3
      zt  = 0.5d+00 * ((bl - br)**2 + (l3 / (1.0d+00 + l3))**2)

! calculate αₖ
!
      al(i,1) = 1.0d+00 + zt / (bl + eps)**2
      ac(i,1) = 1.0d+00 + zt / (bc + eps)**2
      ar(i,1) = 1.0d+00 + zt / (br + eps)**2

! calculate βₖ
!
      df  = abs(dfm(i))
      lq  = xi * df
      bl  = df2(im1) + lq
      bc  = df2(i  ) + lq
      br  = df2(ip1) + xi * abs(2.0d+00 * dfp(i) - dfp(ip1))

! calculate ζ

      l3  = df**3
      zt  = 0.5d+00 * ((bl - br)**2 + (l3 / (1.0d+00 + l3))**2)

! calculate αₖ
!
      al(i,2) = 1.0d+00 + zt / (bl + eps)**2
      ac(i,2) = 1.0d+00 + zt / (bc + eps)**2
      ar(i,2) = 1.0d+00 + zt / (br + eps)**2

    end do ! i = 2, n - 1

! prepare tridiagonal system coefficients
!
    do i = ng, n - ng + 1

! prepare neighbour indices
!
      im1 = i - 1
      ip1 = i + 1

! calculate weights
!
      wl  = cl * al(i,1)
      wc  = cc * ac(i,1)
      wr  = cr * ar(i,1)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate tridiagonal matrix coefficients
!
      a(i,1) = 2.0d+00 * wl +            wc
      b(i,1) =           wl + 2.0d+00 * (wc + wr)
      c(i,1) =           wr

! prepare right hand side of tridiagonal equation
!
      r(i,1) = (wl * fc(im1) + (5.0d+00 * (wl + wc) + wr) * fc(i  )            &
                                    + (wc + 5.0d+00 * wr) * fc(ip1)) * dq

! calculate weights
!
      wl  = cl * ar(i,2)
      wc  = cc * ac(i,2)
      wr  = cr * al(i,2)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate tridiagonal matrix coefficients
!
      a(i,2) =           wr
      b(i,2) =           wl + 2.0d+00 * (wc + wr)
      c(i,2) = 2.0d+00 * wl +            wc

! prepare right hand side of tridiagonal equation
!
      r(i,2) = (wl * fc(ip1) + (5.0d+00 * (wl + wc) + wr) * fc(i  )            &
                                    + (wc + 5.0d+00 * wr) * fc(im1)) * dq

    end do ! i = ng, n - ng + 1

! interpolate ghost zones using explicit solver (left-side reconstruction)
!
    do i = 2, ng

! prepare neighbour indices
!
      im2 = max(1, i - 2)
      im1 = i - 1
      ip1 = i + 1
      ip2 = i + 2

! calculate weights
!
      wl  = dl * al(i,1)
      wc  = dc * ac(i,1)
      wr  = dr * ar(i,1)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the left state
!
      ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i  )
      qc = a21 * fc(im1) + a22 * fc(i  ) + a23 * fc(ip1)
      qr = a31 * fc(i  ) + a32 * fc(ip1) + a33 * fc(ip2)

! calculate the left state
!
      fl(i) = (wl * ql + wr * qr) + wc * qc

! prepare coefficients of the tridiagonal system
!
      a(i,1) = 0.0d+00
      b(i,1) = 1.0d+00
      c(i,1) = 0.0d+00
      r(i,1) = fl(i)

    end do ! i = 2, ng
    a(1,1) = 0.0d+00
    b(1,1) = 1.0d+00
    c(1,1) = 0.0d+00
    r(1,1) = 0.5d+00 * (fc(1) + fc(2))

! interpolate ghost zones using explicit solver (left-side reconstruction)
!
    do i = n - ng, n - 1

! prepare neighbour indices
!
      im2 = i - 2
      im1 = i - 1
      ip1 = i + 1
      ip2 = min(n, i + 2)

! calculate weights
!
      wl  = dl * al(i,1)
      wc  = dc * ac(i,1)
      wr  = dr * ar(i,1)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the left state
!
      ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i  )
      qc = a21 * fc(im1) + a22 * fc(i  ) + a23 * fc(ip1)
      qr = a31 * fc(i  ) + a32 * fc(ip1) + a33 * fc(ip2)

! calculate the left state
!
      fl(i) = (wl * ql + wr * qr) + wc * qc

! prepare coefficients of the tridiagonal system
!
      a(i,1) = 0.0d+00
      b(i,1) = 1.0d+00
      c(i,1) = 0.0d+00
      r(i,1) = fl(i)

    end do ! i = n - ng, n - 1
    a(n,1) = 0.0d+00
    b(n,1) = 1.0d+00
    c(n,1) = 0.0d+00
    r(n,1) = fc(n)

! interpolate ghost zones using explicit solver (right-side reconstruction)
!
    do i = 2, ng + 1

! prepare neighbour indices
!
      im2 = max(1, i - 2)
      im1 = i - 1
      ip1 = i + 1
      ip2 = i + 2

! normalize weights
!
      wl  = dl * ar(i,2)
      wc  = dc * ac(i,2)
      wr  = dr * al(i,2)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the right state
!
      ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i  )
      qc = a21 * fc(ip1) + a22 * fc(i  ) + a23 * fc(im1)
      qr = a31 * fc(i  ) + a32 * fc(im1) + a33 * fc(im2)

! calculate the right state
!
      fr(i) = (wl * ql + wr * qr) + wc * qc

! prepare coefficients of the tridiagonal system
!
      a(i,2) = 0.0d+00
      b(i,2) = 1.0d+00
      c(i,2) = 0.0d+00
      r(i,2) = fr(i)

    end do ! i = 2, ng + 1
    a(1,2) = 0.0d+00
    b(1,2) = 1.0d+00
    c(1,2) = 0.0d+00
    r(1,2) = fc(1)

! interpolate ghost zones using explicit solver (right-side reconstruction)
!
    do i = n - ng + 1, n - 1

! prepare neighbour indices
!
      im2 = i - 2
      im1 = i - 1
      ip1 = i + 1
      ip2 = min(n, i + 2)

! normalize weights
!
      wl  = dl * ar(i,2)
      wc  = dc * ac(i,2)
      wr  = dr * al(i,2)
      ww  = (wl + wr) + wc
      wl  = wl / ww
      wr  = wr / ww
      wc  = 1.0d+00 - (wl + wr)

! calculate the interpolations of the right state
!
      ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i  )
      qc = a21 * fc(ip1) + a22 * fc(i  ) + a23 * fc(im1)
      qr = a31 * fc(i  ) + a32 * fc(im1) + a33 * fc(im2)

! calculate the right state
!
      fr(i) = (wl * ql + wr * qr) + wc * qc

! prepare coefficients of the tridiagonal system
!
      a(i,2) = 0.0d+00
      b(i,2) = 1.0d+00
      c(i,2) = 0.0d+00
      r(i,2) = fr(i)

    end do ! i = n - ng + 1, n - 1
    a(n,2) = 0.0d+00
    b(n,2) = 1.0d+00
    c(n,2) = 0.0d+00
    r(n,2) = 0.5d+00 * (fc(n-1) + fc(n))

! solve the tridiagonal system of equations for the left-side interpolation
!
    call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n))

! substitute the left-side values
!
    fl(1:n  ) = u(1:n)

! solve the tridiagonal system of equations for the left-side interpolation
!
    call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n))

! substitute the right-side values
!
    fr(1:n-1) = u(2:n)

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_crweno5ns
!
!===============================================================================
!
! subroutine RECONSTRUCT_MP5:
! --------------------------
!
!   Subroutine reconstructs the interface states using the fifth order
!   Monotonicity Preserving (MP) method.
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Suresh, A. & Huynh, H. T.,
!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
!          Time Stepping"
!         Journal on Computational Physics,
!         1997, vol. 136, pp. 83-99,
!         http://dx.doi.org/10.1006/jcph.1997.5745
!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
!         "A 5th order monotonicity-preserving upwind compact difference
!          scheme",
!         Science China Physics, Mechanics and Astronomy,
!         Volume 54, Issue 3, pp. 511-522,
!         http://dx.doi.org/10.1007/s11433-010-4220-x
!
!===============================================================================
!
  subroutine reconstruct_mp5(h, fc, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer :: n, i

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc)) :: fi
    real(kind=8), dimension(size(fc)) :: u

! local parameters
!
    real(kind=8), dimension(5), parameter ::                                   &
                                ce5 = (/ 2.0d+00,-1.3d+01, 4.7d+01             &
                                                , 2.7d+01,-3.0d+00 /) / 6.0d+01
    real(kind=8), dimension(3), parameter :: &
                                ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00 /) / 6.0d+00
    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

!! === left-side interpolation ===
!!
! reconstruct the interface state using the 5th order interpolation
!
    do i = 3, n - 2
      u(i) = sum(ce5(:) * fc(i-2:i+2))
    end do

! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
    u(  1) = sum(ce2(:) * fc(  1:  2))
    u(  2) = sum(ce3(:) * fc(  1:  3))
    u(n-1) = sum(ce3(:) * fc(n-2:  n))
    u(n  ) =              fc(n      )

! apply the monotonicity preserving limiting
!
    call mp_limiting(fc(:), u(:))

! copy the interpolation to the respective vector
!
    fl(1:n) = u(1:n)

!! === right-side interpolation ===
!!
! invert the cell-centered value vector
!
    fi(1:n) = fc(n:1:-1)

! reconstruct the interface state using the 5th order interpolation
!
    do i = 3, n - 2
      u(i) = sum(ce5(:) * fi(i-2:i+2))
    end do

! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
    u(  1) = sum(ce2(:) * fi(  1:  2))
    u(  2) = sum(ce3(:) * fi(  1:  3))
    u(n-1) = sum(ce3(:) * fi(n-2:  n))
    u(n  ) =              fi(n      )

! apply the monotonicity preserving limiting
!
    call mp_limiting(fi(:), u(:))

! copy the interpolation to the respective vector
!
    fr(1:n-1) = u(n-1:1:-1)

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_mp5
!
!===============================================================================
!
! subroutine RECONSTRUCT_MP7:
! --------------------------
!
!   Subroutine reconstructs the interface states using the seventh order
!   Monotonicity Preserving (MP) method.
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Suresh, A. & Huynh, H. T.,
!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
!          Time Stepping"
!         Journal on Computational Physics,
!         1997, vol. 136, pp. 83-99,
!         http://dx.doi.org/10.1006/jcph.1997.5745
!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
!         "A 5th order monotonicity-preserving upwind compact difference
!          scheme",
!         Science China Physics, Mechanics and Astronomy,
!         Volume 54, Issue 3, pp. 511-522,
!         http://dx.doi.org/10.1007/s11433-010-4220-x
!
!===============================================================================
!
  subroutine reconstruct_mp7(h, fc, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer :: n, i

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc)) :: fi
    real(kind=8), dimension(size(fc)) :: u

! local parameters
!
    real(kind=8), dimension(7), parameter ::                                   &
                                ce7 = (/-3.0d+00, 2.5d+01,-1.01d+02, 3.19d+02  &
                                      , 2.14d+02,-3.8d+01, 4.0d+00 /) / 4.2d+02
    real(kind=8), dimension(5), parameter ::                                   &
                                ce5 = (/ 2.0d+00,-1.3d+01, 4.7d+01             &
                                                , 2.7d+01,-3.0d+00 /) / 6.0d+01
    real(kind=8), dimension(3), parameter :: &
                                ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00 /) / 6.0d+00
    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

!! === left-side interpolation ===
!!
! reconstruct the interface state using the 5th order interpolation
!
    do i = 4, n - 3
      u(i) = sum(ce7(:) * fc(i-3:i+3))
    end do

! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
    u(  1) = sum(ce2(:) * fc(  1:  2))
    u(  2) = sum(ce3(:) * fc(  1:  3))
    u(  3) = sum(ce5(:) * fc(  1:  5))
    u(n-2) = sum(ce5(:) * fc(n-4:  n))
    u(n-1) = sum(ce3(:) * fc(n-2:  n))
    u(n  ) =              fc(n      )

! apply the monotonicity preserving limiting
!
    call mp_limiting(fc(:), u(:))

! copy the interpolation to the respective vector
!
    fl(1:n) = u(1:n)

!! === right-side interpolation ===
!!
! invert the cell-centered value vector
!
    fi(1:n) = fc(n:1:-1)

! reconstruct the interface state using the 5th order interpolation
!
    do i = 4, n - 3
      u(i) = sum(ce7(:) * fi(i-3:i+3))
    end do

! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
    u(  1) = sum(ce2(:) * fi(  1:  2))
    u(  2) = sum(ce3(:) * fi(  1:  3))
    u(  3) = sum(ce5(:) * fi(  1:  5))
    u(n-2) = sum(ce5(:) * fi(n-4:  n))
    u(n-1) = sum(ce3(:) * fi(n-2:  n))
    u(n  ) =              fi(n      )

! apply the monotonicity preserving limiting
!
    call mp_limiting(fi(:), u(:))

! copy the interpolation to the respective vector
!
    fr(1:n-1) = u(n-1:1:-1)

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_mp7
!
!===============================================================================
!
! subroutine RECONSTRUCT_MP9:
! --------------------------
!
!   Subroutine reconstructs the interface states using the ninth order
!   Monotonicity Preserving (MP) method.
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Suresh, A. & Huynh, H. T.,
!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
!          Time Stepping"
!         Journal on Computational Physics,
!         1997, vol. 136, pp. 83-99,
!         http://dx.doi.org/10.1006/jcph.1997.5745
!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
!         "A 5th order monotonicity-preserving upwind compact difference
!          scheme",
!         Science China Physics, Mechanics and Astronomy,
!         Volume 54, Issue 3, pp. 511-522,
!         http://dx.doi.org/10.1007/s11433-010-4220-x
!
!===============================================================================
!
  subroutine reconstruct_mp9(h, fc, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer :: n, i

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc)) :: fi
    real(kind=8), dimension(size(fc)) :: u

! local parameters
!
    real(kind=8), dimension(9), parameter ::                                   &
             ce9 = (/ 4.000d+00,-4.100d+01, 1.990d+02,-6.410d+02, 1.879d+03,   &
                      1.375d+03,-3.050d+02, 5.500d+01,-5.000d+00 /) / 2.520d+03
    real(kind=8), dimension(7), parameter ::                                   &
             ce7 = (/-3.000d+00, 2.500d+01,-1.010d+02, 3.190d+02, 2.140d+02,   &
                     -3.800d+01, 4.000d+00 /) / 4.200d+02
    real(kind=8), dimension(5), parameter ::                                   &
             ce5 = (/ 2.0d+00,-1.3d+01, 4.7d+01, 2.7d+01,-3.0d+00 /) / 6.0d+01
    real(kind=8), dimension(3), parameter ::                                   &
             ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00 /) / 6.0d+00
    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

!! === left-side interpolation ===
!!
! reconstruct the interface state using the 9th order interpolation
!
    do i = 5, n - 4
      u(i) = sum(ce9(:) * fc(i-4:i+4))
    end do

! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
    u(  1) = sum(ce2(:) * fc(  1:  2))
    u(  2) = sum(ce3(:) * fc(  1:  3))
    u(  3) = sum(ce5(:) * fc(  1:  5))
    u(  4) = sum(ce7(:) * fc(  1:  7))
    u(n-3) = sum(ce7(:) * fc(n-6:  n))
    u(n-2) = sum(ce5(:) * fc(n-4:  n))
    u(n-1) = sum(ce3(:) * fc(n-2:  n))
    u(n  ) =              fc(n      )

! apply the monotonicity preserving limiting
!
    call mp_limiting(fc(:), u(:))

! copy the interpolation to the respective vector
!
    fl(1:n) = u(1:n)

!! === right-side interpolation ===
!!
! invert the cell-centered value vector
!
    fi(1:n) = fc(n:1:-1)

! reconstruct the interface state using the 9th order interpolation
!
    do i = 5, n - 4
      u(i) = sum(ce9(:) * fi(i-4:i+4))
    end do

! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
    u(  1) = sum(ce2(:) * fi(  1:  2))
    u(  2) = sum(ce3(:) * fi(  1:  3))
    u(  3) = sum(ce5(:) * fi(  1:  5))
    u(  4) = sum(ce7(:) * fi(  1:  7))
    u(n-3) = sum(ce7(:) * fi(n-6:  n))
    u(n-2) = sum(ce5(:) * fi(n-4:  n))
    u(n-1) = sum(ce3(:) * fi(n-2:  n))
    u(n  ) =              fi(n      )

! apply the monotonicity preserving limiting
!
    call mp_limiting(fi(:), u(:))

! copy the interpolation to the respective vector
!
    fr(1:n-1) = u(n-1:1:-1)

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_mp9
!
!===============================================================================
!
! subroutine RECONSTRUCT_MP5LD:
! ----------------------------
!
!   Subroutine reconstructs the interface states using the fifth order
!   low dissipation Monotonicity Preserving (MP) method.
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Suresh, A. & Huynh, H. T.,
!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
!          Time Stepping"
!         Journal on Computational Physics,
!         1997, vol. 136, pp. 83-99,
!         http://dx.doi.org/10.1006/jcph.1997.5745
!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
!         "A 5th order monotonicity-preserving upwind compact difference
!          scheme",
!         Science China Physics, Mechanics and Astronomy,
!         Volume 54, Issue 3, pp. 511-522,
!         http://dx.doi.org/10.1007/s11433-010-4220-x
!
!===============================================================================
!
  subroutine reconstruct_mp5ld(h, fc, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer :: n, i

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc)) :: fi
    real(kind=8), dimension(size(fc)) :: u

! local parameters
!
    real(kind=8), dimension(6), parameter ::                                   &
                              ce5 = [ 1.20d+01,-8.10d+01, 3.09d+02,            &
                                      2.09d+02,-3.10d+01, 2.00d+00 ] / 4.2d+02
    real(kind=8), dimension(3), parameter ::                                   &
                              ce3 = [-1.00d+00, 5.00d+00, 2.00d+00 ] / 6.0d+00
    real(kind=8), dimension(2), parameter ::                                   &
                              ce2 = [ 5.0d-01, 5.0d-01 ]
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

!! === left-side interpolation ===
!!
! reconstruct the interface state using the 5th order interpolation
!
    do i = 3, n - 3
      u(i) = sum(ce5(:) * fc(i-2:i+3))
    end do

! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
    u(  1) = sum(ce2(:) * fc(  1:  2))
    u(  2) = sum(ce3(:) * fc(  1:  3))
    u(n-2) = sum(ce3(:) * fc(n-3:n-1))
    u(n-1) = sum(ce3(:) * fc(n-2:n  ))
    u(n  ) =              fc(n      )

! apply the monotonicity preserving limiting
!
    call mp_limiting(fc(:), u(:))

! copy the interpolation to the respective vector
!
    fl(1:n) = u(1:n)

!! === right-side interpolation ===
!!
! invert the cell-centered value vector
!
    fi(1:n) = fc(n:1:-1)

! reconstruct the interface state using the 5th order interpolation
!
    do i = 3, n - 3
      u(i) = sum(ce5(:) * fi(i-2:i+3))
    end do

! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
    u(  1) = sum(ce2(:) * fi(  1:  2))
    u(  2) = sum(ce3(:) * fi(  1:  3))
    u(n-2) = sum(ce3(:) * fi(n-3:n-1))
    u(n-1) = sum(ce3(:) * fi(n-2:n  ))
    u(n  ) =              fi(n      )

! apply the monotonicity preserving limiting
!
    call mp_limiting(fi(:), u(:))

! copy the interpolation to the respective vector
!
    fr(1:n-1) = u(n-1:1:-1)

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_mp5ld
!
!===============================================================================
!
! subroutine RECONSTRUCT_MP7LD:
! ----------------------------
!
!   Subroutine reconstructs the interface states using the seventh order
!   low dissipation Monotonicity Preserving (MP) method.
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Suresh, A. & Huynh, H. T.,
!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
!          Time Stepping"
!         Journal on Computational Physics,
!         1997, vol. 136, pp. 83-99,
!         http://dx.doi.org/10.1006/jcph.1997.5745
!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
!         "A 5th order monotonicity-preserving upwind compact difference
!          scheme",
!         Science China Physics, Mechanics and Astronomy,
!         Volume 54, Issue 3, pp. 511-522,
!         http://dx.doi.org/10.1007/s11433-010-4220-x
!
!===============================================================================
!
  subroutine reconstruct_mp7ld(h, fc, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer :: n, i

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc)) :: fi
    real(kind=8), dimension(size(fc)) :: u

! local parameters
!
    real(kind=8), dimension(8), parameter ::                                   &
                    ce7 = [-8.00d+00, 6.80d+01,-2.82d+02, 9.22d+02,            &
                            6.77d+02,-1.35d+02, 1.9d+01, -1.0d+00 ] / 1.26d+03
    real(kind=8), dimension(5), parameter ::                                   &
                    ce5 = [ 2.0d+00,-1.3d+01, 4.7d+01                          &
                                            , 2.7d+01,-3.0d+00 ] / 6.0d+01
    real(kind=8), dimension(3), parameter ::                                   &
                    ce3 = [-1.0d+00, 5.0d+00, 2.0d+00 ] / 6.0d+00
    real(kind=8), dimension(2), parameter ::                                   &
                    ce2 = [ 5.0d-01, 5.0d-01 ]
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

!! === left-side interpolation ===
!!
! reconstruct the interface state using the 5th order interpolation
!
    do i = 4, n - 4
      u(i) = sum(ce7(:) * fc(i-3:i+4))
    end do

! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
    u(  1) = sum(ce2(:) * fc(  1:  2))
    u(  2) = sum(ce3(:) * fc(  1:  3))
    u(  3) = sum(ce5(:) * fc(  1:  5))
    u(n-3) = sum(ce5(:) * fc(n-5:n-1))
    u(n-2) = sum(ce5(:) * fc(n-4:n  ))
    u(n-1) = sum(ce3(:) * fc(n-2:n  ))
    u(n  ) =              fc(n      )

! apply the monotonicity preserving limiting
!
    call mp_limiting(fc(:), u(:))

! copy the interpolation to the respective vector
!
    fl(1:n) = u(1:n)

!! === right-side interpolation ===
!!
! invert the cell-centered value vector
!
    fi(1:n) = fc(n:1:-1)

! reconstruct the interface state using the 5th order interpolation
!
    do i = 4, n - 4
      u(i) = sum(ce7(:) * fi(i-3:i+4))
    end do

! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
    u(  1) = sum(ce2(:) * fi(  1:  2))
    u(  2) = sum(ce3(:) * fi(  1:  3))
    u(  3) = sum(ce5(:) * fi(  1:  5))
    u(n-3) = sum(ce5(:) * fi(n-5:n-1))
    u(n-2) = sum(ce5(:) * fi(n-4:n  ))
    u(n-1) = sum(ce3(:) * fi(n-2:n  ))
    u(n  ) =              fi(n      )

! apply the monotonicity preserving limiting
!
    call mp_limiting(fi(:), u(:))

! copy the interpolation to the respective vector
!
    fr(1:n-1) = u(n-1:1:-1)

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_mp7ld
!
!===============================================================================
!
! subroutine RECONSTRUCT_MP9LD:
! ----------------------------
!
!   Subroutine reconstructs the interface states using the ninth order
!   low dissipation Monotonicity Preserving (MP) method.
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Suresh, A. & Huynh, H. T.,
!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
!          Time Stepping"
!         Journal on Computational Physics,
!         1997, vol. 136, pp. 83-99,
!         http://dx.doi.org/10.1006/jcph.1997.5745
!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
!         "A 5th order monotonicity-preserving upwind compact difference
!          scheme",
!         Science China Physics, Mechanics and Astronomy,
!         Volume 54, Issue 3, pp. 511-522,
!         http://dx.doi.org/10.1007/s11433-010-4220-x
!
!===============================================================================
!
  subroutine reconstruct_mp9ld(h, fc, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer :: n, i

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc)) :: fi
    real(kind=8), dimension(size(fc)) :: u

! local parameters
!
    real(kind=8), dimension(10), parameter ::                                  &
           ce9 = [ 4.0000d+01,-4.1500d+02, 2.0450d+03,-6.7150d+03,             &
                   2.0165d+04, 1.5629d+04,-3.6910d+03, 7.4900d+02,             &
                  -9.1000d+01, 4.0000d+00 ] / 2.7720d+04
    real(kind=8), dimension(7), parameter ::                                   &
           ce7 = [-3.000d+00, 2.500d+01,-1.010d+02, 3.190d+02, 2.140d+02,      &
                  -3.800d+01, 4.000d+00 ] / 4.200d+02
    real(kind=8), dimension(5), parameter ::                                   &
           ce5 = [ 2.0d+00,-1.3d+01, 4.7d+01, 2.7d+01,-3.0d+00 ] / 6.0d+01
    real(kind=8), dimension(3), parameter ::                                   &
           ce3 = [-1.0d+00, 5.0d+00, 2.0d+00 ] / 6.0d+00
    real(kind=8), dimension(2), parameter ::                                   &
           ce2 = [ 5.0d-01, 5.0d-01 ]
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

!! === left-side interpolation ===
!!
! reconstruct the interface state using the 9th order interpolation
!
    do i = 5, n - 5
      u(i) = sum(ce9(:) * fc(i-4:i+5))
    end do

! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
    u(  1) = sum(ce2(:) * fc(  1:  2))
    u(  2) = sum(ce3(:) * fc(  1:  3))
    u(  3) = sum(ce5(:) * fc(  1:  5))
    u(  4) = sum(ce7(:) * fc(  1:  7))
    u(n-4) = sum(ce7(:) * fc(n-7:n-1))
    u(n-3) = sum(ce7(:) * fc(n-6:n  ))
    u(n-2) = sum(ce5(:) * fc(n-4:n  ))
    u(n-1) = sum(ce3(:) * fc(n-2:n  ))
    u(n  ) =              fc(n      )

! apply the monotonicity preserving limiting
!
    call mp_limiting(fc(:), u(:))

! copy the interpolation to the respective vector
!
    fl(1:n) = u(1:n)

!! === right-side interpolation ===
!!
! invert the cell-centered value vector
!
    fi(1:n) = fc(n:1:-1)

! reconstruct the interface state using the 9th order interpolation
!
    do i = 5, n - 5
      u(i) = sum(ce9(:) * fi(i-4:i+5))
    end do

! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
    u(  1) = sum(ce2(:) * fi(  1:  2))
    u(  2) = sum(ce3(:) * fi(  1:  3))
    u(  3) = sum(ce5(:) * fi(  1:  5))
    u(  4) = sum(ce7(:) * fi(  1:  7))
    u(n-4) = sum(ce7(:) * fi(n-7:n-1))
    u(n-3) = sum(ce7(:) * fi(n-6:n  ))
    u(n-2) = sum(ce5(:) * fi(n-4:n  ))
    u(n-1) = sum(ce3(:) * fi(n-2:n  ))
    u(n  ) =              fi(n      )

! apply the monotonicity preserving limiting
!
    call mp_limiting(fi(:), u(:))

! copy the interpolation to the respective vector
!
    fr(1:n-1) = u(n-1:1:-1)

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_mp9ld
!
!===============================================================================
!
! subroutine RECONSTRUCT_CRMP5:
! ----------------------------
!
!   Subroutine reconstructs the interface states using the fifth order
!   Compact Reconstruction Monotonicity Preserving (CRMP) method.
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Suresh, A. & Huynh, H. T.,
!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
!          Time Stepping"
!         Journal on Computational Physics,
!         1997, vol. 136, pp. 83-99,
!         http://dx.doi.org/10.1006/jcph.1997.5745
!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
!         "A 5th order monotonicity-preserving upwind compact difference
!          scheme",
!         Science China Physics, Mechanics and Astronomy,
!         Volume 54, Issue 3, pp. 511-522,
!         http://dx.doi.org/10.1007/s11433-010-4220-x
!
!===============================================================================
!
  subroutine reconstruct_crmp5(h, fc, fl, fr)

! include external procedures
!
    use algebra   , only : tridiag

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer :: n, i

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc)) :: fi
    real(kind=8), dimension(size(fc)) :: a, b, c
    real(kind=8), dimension(size(fc)) :: r
    real(kind=8), dimension(size(fc)) :: u

! local parameters
!
    real(kind=8), dimension(3), parameter ::                                   &
                                di  = (/ 3.0d-01, 6.0d-01, 1.0d-01 /)
    real(kind=8), dimension(3), parameter ::                                   &
                                ci5 = (/ 1.0d+00, 1.9d+01, 1.0d+01 /) / 3.0d+01
    real(kind=8), dimension(5), parameter ::                                   &
                                ce5 = (/ 2.0d+00,-1.3d+01, 4.7d+01             &
                                                , 2.7d+01,-3.0d+00 /) / 6.0d+01
    real(kind=8), dimension(3), parameter :: &
                                ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00 /) / 6.0d+00
    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

! prepare the diagonals of the tridiagonal matrix
!
    do i = 1, ng
      a(i) = 0.0d+00
      b(i) = 1.0d+00
      c(i) = 0.0d+00
    end do
    do i = ng + 1, n - ng - 1
      a(i) = di(1)
      b(i) = di(2)
      c(i) = di(3)
    end do
    do i = n - ng, n
      a(i) = 0.0d+00
      b(i) = 1.0d+00
      c(i) = 0.0d+00
    end do

!! === left-side interpolation ===
!!
! prepare the tridiagonal system coefficients for the interior part
!
    do i = ng, n - ng + 1
      r(i) = sum(ci5(:) * fc(i-1:i+1))
    end do

! interpolate ghost zones using the explicit method
!
    r(  1) = sum(ce2(:) * fc(  1:  2))
    r(  2) = sum(ce3(:) * fc(  1:  3))
    do i = 3, ng
      r(i) = sum(ce5(:) * fc(i-2:i+2))
    end do
    do i = n - ng, n - 2
      r(i) = sum(ce5(:) * fc(i-2:i+2))
    end do
    r(n-1) = sum(ce3(:) * fc(n-2:  n))
    r(n  ) =              fc(n      )

! solve the tridiagonal system of equations
!
    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))

! apply the monotonicity preserving limiting
!
    call mp_limiting(fc(:), u(:))

! copy the interpolation to the respective vector
!
    fl(1:n) = u(1:n)

!! === right-side interpolation ===
!!
! invert the cell-centered value vector
!
    fi(1:n) = fc(n:1:-1)

! prepare the tridiagonal system coefficients for the interior part
!
    do i = ng, n - ng + 1
      r(i) = sum(ci5(:) * fi(i-1:i+1))
    end do ! i = ng, n - ng + 1

! interpolate ghost zones using the explicit method
!
    r(  1) = sum(ce2(:) * fi(  1:  2))
    r(  2) = sum(ce3(:) * fi(  1:  3))
    do i = 3, ng
      r(i) = sum(ce5(:) * fi(i-2:i+2))
    end do
    do i = n - ng, n - 2
      r(i) = sum(ce5(:) * fi(i-2:i+2))
    end do
    r(n-1) = sum(ce3(:) * fi(n-2:  n))
    r(n  ) =              fi(n      )

! solve the tridiagonal system of equations
!
    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))

! apply the monotonicity preserving limiting
!
    call mp_limiting(fi(:), u(:))

! copy the interpolation to the respective vector
!
    fr(1:n-1) = u(n-1:1:-1)

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_crmp5
!
!===============================================================================
!
! subroutine RECONSTRUCT_CRMP5LD:
! ------------------------------
!
!   Subroutine reconstructs the interface states using the fifth order
!   Low-Dissipation Compact Reconstruction Monotonicity Preserving (CRMP)
!   method.
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Suresh, A. & Huynh, H. T.,
!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
!          Time Stepping"
!         Journal on Computational Physics,
!         1997, vol. 136, pp. 83-99,
!         http://dx.doi.org/10.1006/jcph.1997.5745
!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
!         "A 5th order monotonicity-preserving upwind compact difference
!          scheme",
!         Science China Physics, Mechanics and Astronomy,
!         Volume 54, Issue 3, pp. 511-522,
!         http://dx.doi.org/10.1007/s11433-010-4220-x
!     [3] Ghosh, D. & Baeder, J.,
!         "Compact Reconstruction Schemes With Weighted ENO Limiting For
!          Hyperbolic Conservation Laws",
!         SIAM Journal on Scientific Computing,
!         2012, vol. 34, no. 3, pp. A1678-A1705,
!         http://dx.doi.org/10.1137/110857659
!
!===============================================================================
!
  subroutine reconstruct_crmp5ld(h, fc, fl, fr)

! include external procedures
!
    use algebra   , only : tridiag

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer :: n, i

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc)) :: fi
    real(kind=8), dimension(size(fc)) :: a, b, c
    real(kind=8), dimension(size(fc)) :: r
    real(kind=8), dimension(size(fc)) :: u

! local parameters
!
    real(kind=8), dimension(3), parameter ::                                   &
                                di  = (/ 2.5d-01, 6.0d-01, 1.5d-01 /)
    real(kind=8), dimension(4), parameter ::                                   &
                                ci5 = (/ 3.0d+00, 6.7d+01, 4.9d+01             &
                                                         , 1.0d+00 /) / 1.2d+02
    real(kind=8), dimension(5), parameter ::                                   &
                                ce5 = (/ 2.0d+00,-1.3d+01, 4.7d+01             &
                                                , 2.7d+01,-3.0d+00 /) / 6.0d+01
    real(kind=8), dimension(3), parameter :: &
                                ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00 /) / 6.0d+00
    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

! prepare the diagonals of the tridiagonal matrix
!
    do i = 1, ng
      a(i) = 0.0d+00
      b(i) = 1.0d+00
      c(i) = 0.0d+00
    end do
    do i = ng + 1, n - ng - 1
      a(i) = di(1)
      b(i) = di(2)
      c(i) = di(3)
    end do
    do i = n - ng, n
      a(i) = 0.0d+00
      b(i) = 1.0d+00
      c(i) = 0.0d+00
    end do

!! === left-side interpolation ===
!!
! prepare the tridiagonal system coefficients for the interior part
!
    do i = ng, n - ng + 1
      r(i) = sum(ci5(:) * fc(i-1:i+2))
    end do

! interpolate ghost zones using the explicit method
!
    r(  1) = sum(ce2(:) * fc(  1:  2))
    r(  2) = sum(ce3(:) * fc(  1:  3))
    do i = 3, ng
      r(i) = sum(ce5(:) * fc(i-2:i+2))
    end do
    do i = n - ng, n - 2
      r(i) = sum(ce5(:) * fc(i-2:i+2))
    end do
    r(n-1) = sum(ce3(:) * fc(n-2:  n))
    r(n  ) =              fc(n      )

! solve the tridiagonal system of equations
!
    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))

! apply the monotonicity preserving limiting
!
    call mp_limiting(fc(:), u(:))

! copy the interpolation to the respective vector
!
    fl(1:n) = u(1:n)

!! === right-side interpolation ===
!!
! invert the cell-centered value vector
!
    fi(1:n) = fc(n:1:-1)

! prepare the tridiagonal system coefficients for the interior part
!
    do i = ng, n - ng + 1
      r(i) = sum(ci5(:) * fi(i-1:i+2))
    end do ! i = ng, n - ng + 1

! interpolate ghost zones using the explicit method
!
    r(  1) = sum(ce2(:) * fi(  1:  2))
    r(  2) = sum(ce3(:) * fi(  1:  3))
    do i = 3, ng
      r(i) = sum(ce5(:) * fi(i-2:i+2))
    end do
    do i = n - ng, n - 2
      r(i) = sum(ce5(:) * fi(i-2:i+2))
    end do
    r(n-1) = sum(ce3(:) * fi(n-2:  n))
    r(n  ) =              fi(n      )

! solve the tridiagonal system of equations
!
    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))

! apply the monotonicity preserving limiting
!
    call mp_limiting(fi(:), u(:))

! copy the interpolation to the respective vector
!
    fr(1:n-1) = u(n-1:1:-1)

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_crmp5ld
!
!===============================================================================
!
! subroutine RECONSTRUCT_CRMP7:
! ----------------------------
!
!   Subroutine reconstructs the interface states using the seventh order
!   Compact Reconstruction Monotonicity Preserving (CRMP) method.
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Suresh, A. & Huynh, H. T.,
!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
!          Time Stepping"
!         Journal on Computational Physics,
!         1997, vol. 136, pp. 83-99,
!         http://dx.doi.org/10.1006/jcph.1997.5745
!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
!         "A 5th order monotonicity-preserving upwind compact difference
!          scheme",
!         Science China Physics, Mechanics and Astronomy,
!         Volume 54, Issue 3, pp. 511-522,
!         http://dx.doi.org/10.1007/s11433-010-4220-x
!
!===============================================================================
!
  subroutine reconstruct_crmp7(h, fc, fl, fr)

! include external procedures
!
    use algebra   , only : tridiag

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer :: n, i

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc)) :: fi
    real(kind=8), dimension(size(fc)) :: a, b, c
    real(kind=8), dimension(size(fc)) :: r
    real(kind=8), dimension(size(fc)) :: u

! local parameters
!
    real(kind=8), dimension(3), parameter ::                                   &
                                di  = (/ 2.0d+00, 4.0d+00, 1.0d+00/) / 7.0d+00
    real(kind=8), dimension(5), parameter ::                                   &
                                ci  = (/-1.0d+00, 1.9d+01, 2.39d+02            &
                                              , 1.59d+02, 4.0d+00 /) / 4.2d+02
    real(kind=8), dimension(7), parameter ::                                   &
                                ce7 = (/-3.0d+00, 2.5d+01,-1.01d+02, 3.19d+02  &
                                      , 2.14d+02,-3.8d+01, 4.0d+00/) / 4.2d+02
    real(kind=8), dimension(5), parameter ::                                   &
                                ce5 = (/ 2.0d+00,-1.3d+01, 4.70d+01            &
                                               , 2.70d+01,-3.0d+00/) / 6.0d+01
    real(kind=8), dimension(3), parameter :: &
                                ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00/) / 6.0d+00
    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

! prepare the diagonals of the tridiagonal matrix
!
    do i = 1, ng
      a(i) = 0.0d+00
      b(i) = 1.0d+00
      c(i) = 0.0d+00
    end do
    do i = ng + 1, n - ng - 1
      a(i) = di(1)
      b(i) = di(2)
      c(i) = di(3)
    end do
    do i = n - ng, n
      a(i) = 0.0d+00
      b(i) = 1.0d+00
      c(i) = 0.0d+00
    end do

!! === left-side interpolation ===
!!
! prepare the tridiagonal system coefficients for the interior part
!
    do i = ng, n - ng + 1
      r(i) = sum( ci(:) * fc(i-2:i+2))
    end do

! interpolate ghost zones using the explicit method
!
    r(  1) = sum(ce2(:) * fc(  1:  2))
    r(  2) = sum(ce3(:) * fc(  1:  3))
    r(  3) = sum(ce5(:) * fc(  1:  5))
    do i = 4, ng
      r(i) = sum(ce7(:) * fc(i-3:i+3))
    end do
    do i = n - ng, n - 3
      r(i) = sum(ce7(:) * fc(i-3:i+3))
    end do
    r(n-2) = sum(ce5(:) * fc(n-4:  n))
    r(n-1) = sum(ce3(:) * fc(n-2:  n))
    r(n  ) =              fc(n      )

! solve the tridiagonal system of equations
!
    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))

! apply the monotonicity preserving limiting
!
    call mp_limiting(fc(:), u(:))

! copy the interpolation to the respective vector
!
    fl(1:n) = u(1:n)

!! === right-side interpolation ===
!!
! invert the cell-centered value vector
!
    fi(1:n) = fc(n:1:-1)

! prepare the tridiagonal system coefficients for the interior part
!
    do i = ng, n - ng + 1
      r(i) = sum( ci(:) * fi(i-2:i+2))
    end do ! i = ng, n - ng + 1

! interpolate ghost zones using the explicit method
!
    r(  1) = sum(ce2(:) * fi(  1:  2))
    r(  2) = sum(ce3(:) * fi(  1:  3))
    r(  3) = sum(ce5(:) * fi(  1:  5))
    do i = 4, ng
      r(i) = sum(ce7(:) * fi(i-3:i+3))
    end do
    do i = n - ng, n - 3
      r(i) = sum(ce7(:) * fi(i-3:i+3))
    end do
    r(n-2) = sum(ce5(:) * fi(n-4:  n))
    r(n-1) = sum(ce3(:) * fi(n-2:  n))
    r(n  ) =              fi(n      )

! solve the tridiagonal system of equations
!
    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))

! apply the monotonicity preserving limiting
!
    call mp_limiting(fi(:), u(:))

! copy the interpolation to the respective vector
!
    fr(1:n-1) = u(n-1:1:-1)

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_crmp7
!
!===============================================================================
!
! subroutine RECONSTRUCT_CRMP7LD:
! ------------------------------
!
!   Subroutine reconstructs the interface states using the seventh order
!   low dissipation Compact Reconstruction Monotonicity Preserving (CRMP)
!   method.
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Suresh, A. & Huynh, H. T.,
!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
!          Time Stepping"
!         Journal on Computational Physics,
!         1997, vol. 136, pp. 83-99,
!         http://dx.doi.org/10.1006/jcph.1997.5745
!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
!         "A 5th order monotonicity-preserving upwind compact difference
!          scheme",
!         Science China Physics, Mechanics and Astronomy,
!         Volume 54, Issue 3, pp. 511-522,
!         http://dx.doi.org/10.1007/s11433-010-4220-x
!
!===============================================================================
!
  subroutine reconstruct_crmp7ld(h, fc, fl, fr)

! include external procedures
!
    use algebra   , only : tridiag

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer :: n, i

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc)) :: fi
    real(kind=8), dimension(size(fc)) :: a, b, c
    real(kind=8), dimension(size(fc)) :: r
    real(kind=8), dimension(size(fc)) :: u

! local parameters
!
    real(kind=8), dimension(3), parameter ::                                   &
                              di  = [ 5.0d+00, 1.2d+01, 4.0d+00 ] / 2.1d+01
    real(kind=8), dimension(6), parameter ::                                   &
                              ci  = [-2.00d+00, 4.20d+01, 6.37d+02,            &
                                      5.57d+02, 2.70d+01,-1.00d+00 ] / 1.26d+03
    real(kind=8), dimension(7), parameter ::                                   &
                              ce7 = [-3.0d+00, 2.5d+01,-1.01d+02, 3.19d+02,    &
                                      2.14d+02,-3.8d+01, 4.0d+00 ] / 4.2d+02
    real(kind=8), dimension(5), parameter ::                                   &
                              ce5 = [ 2.0d+00,-1.3d+01, 4.70d+01,              &
                                               2.70d+01,-3.0d+00 ] / 6.0d+01
    real(kind=8), dimension(3), parameter ::                                   &
                              ce3 = [-1.0d+00, 5.0d+00, 2.0d+00 ] / 6.0d+00
    real(kind=8), dimension(2), parameter ::                                   &
                              ce2 = [ 5.0d-01, 5.0d-01 ]
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

! prepare the diagonals of the tridiagonal matrix
!
    do i = 1, ng
      a(i) = 0.0d+00
      b(i) = 1.0d+00
      c(i) = 0.0d+00
    end do
    do i = ng + 1, n - ng - 1
      a(i) = di(1)
      b(i) = di(2)
      c(i) = di(3)
    end do
    do i = n - ng, n
      a(i) = 0.0d+00
      b(i) = 1.0d+00
      c(i) = 0.0d+00
    end do

!! === left-side interpolation ===
!!
! prepare the tridiagonal system coefficients for the interior part
!
    do i = ng, n - ng + 1
      r(i) = sum( ci(:) * fc(i-2:i+3))
    end do

! interpolate ghost zones using the explicit method
!
    r(  1) = sum(ce2(:) * fc(  1:  2))
    r(  2) = sum(ce3(:) * fc(  1:  3))
    r(  3) = sum(ce5(:) * fc(  1:  5))
    do i = 4, ng
      r(i) = sum(ce7(:) * fc(i-3:i+3))
    end do
    do i = n - ng, n - 3
      r(i) = sum(ce7(:) * fc(i-3:i+3))
    end do
    r(n-2) = sum(ce5(:) * fc(n-4:  n))
    r(n-1) = sum(ce3(:) * fc(n-2:  n))
    r(n  ) =              fc(n      )

! solve the tridiagonal system of equations
!
    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))

! apply the monotonicity preserving limiting
!
    call mp_limiting(fc(:), u(:))

! copy the interpolation to the respective vector
!
    fl(1:n) = u(1:n)

!! === right-side interpolation ===
!!
! invert the cell-centered value vector
!
    fi(1:n) = fc(n:1:-1)

! prepare the tridiagonal system coefficients for the interior part
!
    do i = ng, n - ng + 1
      r(i) = sum( ci(:) * fi(i-2:i+3))
    end do ! i = ng, n - ng + 1

! interpolate ghost zones using the explicit method
!
    r(  1) = sum(ce2(:) * fi(  1:  2))
    r(  2) = sum(ce3(:) * fi(  1:  3))
    r(  3) = sum(ce5(:) * fi(  1:  5))
    do i = 4, ng
      r(i) = sum(ce7(:) * fi(i-3:i+3))
    end do
    do i = n - ng, n - 3
      r(i) = sum(ce7(:) * fi(i-3:i+3))
    end do
    r(n-2) = sum(ce5(:) * fi(n-4:  n))
    r(n-1) = sum(ce3(:) * fi(n-2:  n))
    r(n  ) =              fi(n      )

! solve the tridiagonal system of equations
!
    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))

! apply the monotonicity preserving limiting
!
    call mp_limiting(fi(:), u(:))

! copy the interpolation to the respective vector
!
    fr(1:n-1) = u(n-1:1:-1)

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_crmp7ld
!
!===============================================================================
!
! subroutine RECONSTRUCT_OCMP5:
! -----------------------------
!
!   Subroutine reconstructs the interface states using the 5th order Optimized
!   Compact Reconstruction Monotonicity Preserving (CRMP) method.
!
!   Arguments are described in subroutine reconstruct().
!
!   References:
!
!     [1] Myeong-Hwan Ahn, Duck-Joo Lee,
!         "Modified Monotonicity Preserving Constraints for High-Resolution
!          Optimized Compact Scheme",
!         Journal of Scientific Computing,
!         2020, vol. 83, p. 34
!         https://doi.org/10.1007/s10915-020-01221-0
!
!===============================================================================
!
  subroutine reconstruct_ocmp5(h, fc, fl, fr)

    use algebra, only : tridiag

    implicit none

    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

    integer :: n, i

    real(kind=8), dimension(size(fc)) :: fi
    real(kind=8), dimension(size(fc)) :: a, b, c
    real(kind=8), dimension(size(fc)) :: r
    real(kind=8), dimension(size(fc)) :: u

    real(kind=8), dimension(3), parameter ::                                   &
            di  = [ 5.0163016d-01, 1.0d+00, 2.5394716d-01 ]
    real(kind=8), dimension(5), parameter ::                                   &
            ci5 = [- 3.0d+00 * di(1) - 3.0d+00 * di(3) + 2.0d+00,              &
                     2.7d+01 * di(1) + 1.7d+01 * di(3) - 1.3d+01,              &
                     4.7d+01 * di(1) - 4.3d+01 * di(3) + 4.7d+01,              &
                   - 1.3d+01 * di(1) + 7.7d+01 * di(3) + 2.7d+01,              &
                     2.0d+00 * di(1) + 1.2d+01 * di(3) - 3.0d+00 ] / 6.0d+01
    real(kind=8), dimension(5), parameter ::                                   &
            ce5 = [ 2.0d+00,-1.3d+01, 4.7d+01, 2.7d+01,-3.0d+00 ] / 6.0d+01
    real(kind=8), dimension(3), parameter ::                                   &
            ce3 = [-1.0d+00, 5.0d+00, 2.0d+00 ] / 6.0d+00
    real(kind=8), dimension(2), parameter ::                                   &
            ce2 = [ 5.0d-01, 5.0d-01 ]
!
!-------------------------------------------------------------------------------
!
    n = size(fc)

! prepare the diagonals of the tridiagonal matrix
!
    do i = 1, ng
      a(i) = 0.0d+00
      b(i) = 1.0d+00
      c(i) = 0.0d+00
    end do
    do i = ng + 1, n - ng - 1
      a(i) = di(1)
      b(i) = di(2)
      c(i) = di(3)
    end do
    do i = n - ng, n
      a(i) = 0.0d+00
      b(i) = 1.0d+00
      c(i) = 0.0d+00
    end do

!! === left-side interpolation ===
!!
! prepare the right-hand side of the linear system
!
    do i = ng, n - ng + 1
      r(i) = sum(ci5(:) * fc(i-2:i+2))
    end do

! use explicit methods for ghost zones
!
    r(  1) = sum(ce2(:) * fc(  1:  2))
    r(  2) = sum(ce3(:) * fc(  1:  3))
    do i = 3, ng
      r(i) = sum(ce5(:) * fc(i-2:i+2))
    end do
    do i = n - ng, n - 2
      r(i) = sum(ce5(:) * fc(i-2:i+2))
    end do
    r(n-1) = sum(ce3(:) * fc(n-2:  n))
    r(n  ) =              fc(n      )

! solve the tridiagonal system of equations
!
    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))

! apply the monotonicity preserving limiter
!
    call mp_limiting(fc(:), u(:))

! return the interpolated values of the left state
!
    fl(1:n) = u(1:n)

!! === right-side interpolation ===
!!
! invert the cell-centered integrals
!
    fi(1:n) = fc(n:1:-1)

! prepare the right-hand side of the linear system
!
    do i = ng, n - ng + 1
      r(i) = sum(ci5(:) * fi(i-2:i+2))
    end do ! i = ng, n - ng + 1

! use explicit methods for ghost zones
!
    r(  1) = sum(ce2(:) * fi(  1:  2))
    r(  2) = sum(ce3(:) * fi(  1:  3))
    do i = 3, ng
      r(i) = sum(ce5(:) * fi(i-2:i+2))
    end do
    do i = n - ng, n - 2
      r(i) = sum(ce5(:) * fi(i-2:i+2))
    end do
    r(n-1) = sum(ce3(:) * fi(n-2:  n))
    r(n  ) =              fi(n      )

! solve the tridiagonal system of equations
!
    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))

! apply the monotonicity preserving limiter
!
    call mp_limiting(fi(:), u(:))

! return the interpolated values of the right state
!
    fr(1:n-1) = u(n-1:1:-1)

! update the extremum points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_ocmp5
!
!===============================================================================
!
! subroutine PREPARE_GP:
! ---------------------
!
!   Subroutine prepares matrixes for the Gaussian Process (GP) method.
!
!===============================================================================
!
  subroutine prepare_gp(iret)

! include external procedures
!
    use algebra        , only : max_prec, invert
    use constants      , only : pi
    use iso_fortran_env, only : error_unit

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    integer, intent(inout) :: iret

! local variables
!
    logical             :: flag
    integer             :: i, j, ip, jp
    real(kind=max_prec) :: sig, zl, zr, fc

! local arrays for derivatives
!
    real(kind=max_prec), dimension(:,:), allocatable :: cov, agp
    real(kind=max_prec), dimension(:)  , allocatable :: xgp

! local parameters
!
    character(len=*), parameter :: loc = 'INTERPOLATIONS::prepare_gp()'
!
!-------------------------------------------------------------------------------
!
! calculate normal distribution sigma
!
    sig = sqrt(2.0d+00) * sgp

! allocate the convariance matrix and interpolation position vector
!
    allocate(cov(ngp,ngp))
    allocate(agp(ngp,ngp))
    allocate(xgp(ngp))

! prepare the covariance matrix
!
    fc = 0.5d+00 * sqrt(pi) * sig
    i  = 0
    do ip = - mgp, mgp
      i = i + 1
      j = 0
      do jp = - mgp, mgp
        j = j + 1
        zl = (1.0d+00 * (ip - jp) - 0.5d+00) / sig
        zr = (1.0d+00 * (ip - jp) + 0.5d+00) / sig
        cov(i,j) = fc * (erf(zr) - erf(zl))
      end do
    end do

! invert the matrix
!
    call invert(ngp, cov(:,:), agp(:,:), flag)

! continue if the matrix inversion was successful
!
    if (flag) then

! make the inversed matrix symmetric
!
      do j = 1, ngp
        do i = 1, ngp
          if (i /= j) then
            fc       = 0.5d+00 * (agp(i,j) + agp(j,i))
            agp(i,j) = fc
            agp(j,i) = fc
          end if
        end do
      end do

! prepare the interpolation position vector
!
      j = 0
      do jp = - mgp, mgp
        j      = j + 1
        zr     = (0.5d+00 - 1.0d+00 * jp) / sig
        xgp(j) = exp(- zr**2)
      end do

! prepare the interpolation coefficients vector
!
      cgp(1:ngp) = real(matmul(xgp(1:ngp), agp(1:ngp,1:ngp)), kind=8)

    end if

! deallocate the convariance matrix and interpolation position vector
!
    deallocate(cov)
    deallocate(agp)
    deallocate(xgp)

! check if the matrix was inverted successfully
!
    if (.not. flag) then
      write(error_unit,"('[',a,']: ',a)") trim(loc)                            &
                      , "Could not invert the covariance matrix!"
      iret = 202
    end if

!-------------------------------------------------------------------------------
!
  end subroutine prepare_gp
!
!===============================================================================
!
! subroutine RECONSTRUCT_GP:
! -------------------------
!
!   Subroutine reconstructs the interface states using one dimensional
!   high order Gaussian Process (GP) method.
!
!   Arguments are described in subroutine reconstruct().
!
!===============================================================================
!
  subroutine reconstruct_gp(h, fc, fl, fr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8)              , intent(in)  :: h
    real(kind=8), dimension(:), intent(in)  :: fc
    real(kind=8), dimension(:), intent(out) :: fl, fr

! local variables
!
    integer :: n, i

! local arrays for derivatives
!
    real(kind=8), dimension(size(fc)) :: fi
    real(kind=8), dimension(size(fc)) :: u

! local parameters
!
    real(kind=8), dimension(5), parameter ::                                   &
                                ce5 = (/ 2.0d+00,-1.3d+01, 4.70d+01            &
                                               , 2.70d+01,-3.0d+00/) / 6.0d+01
    real(kind=8), dimension(3), parameter :: &
                                ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00/) / 6.0d+00
    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
    n = size(fc)

!! === left-side interpolation ===
!!
! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
    u(  1) = sum(ce2(:) * fc(  1:  2))
    u(  2) = sum(ce3(:) * fc(  1:  3))
    do i = 3, mgp
      u(i) = sum(ce5(:) * fc(i-2:i+2))
    end do
    do i = n + 1 - mgp, n - 2
      u(i) = sum(ce5(:) * fc(i-2:i+2))
    end do
    u(n-1) = sum(ce3(:) * fc(n-2:  n))
    u(n  ) =              fc(n      )

! reconstruct the interface state using the Gauss process interpolation
!
    do i = 1 + mgp, n - mgp
      u(i) = sum(cgp(1:ngp) * (fc(i-mgp:i+mgp) - fc(i))) + fc(i)
    end do

! apply the monotonicity preserving limiting
!
    call mp_limiting(fc(:), u(:))

! copy the interpolation to the respective vector
!
    fl(1:n) = u(1:n)

!! === right-side interpolation ===
!!
! invert the cell-centered value vector
!
    fi(1:n) = fc(n:1:-1)

! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
    u(  1) = sum(ce2(:) * fi(  1:  2))
    u(  2) = sum(ce3(:) * fi(  1:  3))
    do i = 3, mgp
      u(i) = sum(ce5(:) * fi(i-2:i+2))
    end do
    do i = n + 1 - mgp, n - 2
      u(i) = sum(ce5(:) * fi(i-2:i+2))
    end do
    u(n-1) = sum(ce3(:) * fi(n-2:  n))
    u(n  ) =              fi(n      )

! reconstruct the interface state using the Gauss process interpolation
!
    do i = 1 + mgp, n - mgp
      u(i) = sum(cgp(1:ngp) * (fi(i-mgp:i+mgp) - fi(i))) + fi(i)
    end do

! apply the monotonicity preserving limiting
!
    call mp_limiting(fi(:), u(:))

! copy the interpolation to the respective vector
!
    fr(1:n-1) = u(n-1:1:-1)

! update the interpolation of the first and last points
!
    i     = n - 1
    fl(1) = 0.5d+00 * (fc(1) + fc(2))
    fr(i) = 0.5d+00 * (fc(i) + fc(n))
    fl(n) = fc(n)
    fr(n) = fc(n)

!-------------------------------------------------------------------------------
!
  end subroutine reconstruct_gp
!
!===============================================================================
!
! function LIMITER_ZERO:
! ---------------------
!
!   Function returns zero.
!
!   Arguments:
!
!     x    - scaling factor;
!     a, b - the input values;
!
!===============================================================================
!
  function limiter_zero(x, a, b) result(c)

! local variables are not implicit by default
!
    implicit none

! input arguments
!
    real(kind=8), intent(in) :: x, a, b
    real(kind=8)             :: c
!
!-------------------------------------------------------------------------------
!
    c = 0.0d+00

!-------------------------------------------------------------------------------
!
  end function limiter_zero
!
!===============================================================================
!
! function LIMITER_MINMOD:
! -----------------------
!
!   Function returns the minimum module value among two arguments using
!   minmod limiter.
!
!   Arguments:
!
!     x    - scaling factor;
!     a, b - the input values;
!
!===============================================================================
!
  function limiter_minmod(x, a, b) result(c)

! local variables are not implicit by default
!
    implicit none

! input arguments
!
    real(kind=8), intent(in) :: x, a, b
    real(kind=8)             :: c
!
!-------------------------------------------------------------------------------
!
    c = 0.5d+00 * (sign(x, a) + sign(x, b)) * min(abs(a), abs(b))

!-------------------------------------------------------------------------------
!
  end function limiter_minmod
!
!===============================================================================
!
! function LIMITER_MONOTONIZED_CENTRAL:
! ------------------------------------
!
!   Function returns the minimum module value among two arguments using
!   the monotonized central TVD limiter.
!
!   Arguments:
!
!     x    - scaling factor;
!     a, b - the input values;
!
!===============================================================================
!
  function limiter_monotonized_central(x, a, b) result(c)

! local variables are not implicit by default
!
    implicit none

! input arguments
!
    real(kind=8), intent(in) :: x, a, b
    real(kind=8)             :: c
!
!-------------------------------------------------------------------------------
!
    c = (sign(x, a) + sign(x, b)) * min(abs(a), abs(b), 2.5d-01 * abs(a + b))

!-------------------------------------------------------------------------------
!
  end function limiter_monotonized_central
!
!===============================================================================
!
! function LIMITER_SUPERBEE:
! -------------------------
!
!   Function returns the minimum module value among two arguments using
!   superbee limiter.
!
!   Arguments:
!
!     x    - scaling factor;
!     a, b - the input values;
!
!===============================================================================
!
  function limiter_superbee(x, a, b) result(c)

! local variables are not implicit by default
!
    implicit none

! input arguments
!
    real(kind=8), intent(in) :: x, a, b
    real(kind=8)             :: c
!
!-------------------------------------------------------------------------------
!
    c = 0.5d+00 * (sign(x, a) + sign(x, b))                                    &
           * max(min(2.0d+00 * abs(a), abs(b)), min(abs(a), 2.0d+00 * abs(b)))

!-------------------------------------------------------------------------------
!
  end function limiter_superbee
!
!===============================================================================
!
! function LIMITER_VANLEER:
! ------------------------
!
!   Function returns the minimum module value among two arguments using
!   van Leer's limiter.
!
!   Arguments:
!
!     x    - scaling factor;
!     a, b - the input values;
!
!===============================================================================
!
  function limiter_vanleer(x, a, b) result(c)

! local variables are not implicit by default
!
    implicit none

! input arguments
!
    real(kind=8), intent(in) :: x, a, b
    real(kind=8)             :: c
!
!-------------------------------------------------------------------------------
!
    c = a * b
    if (c > 0.0d+00) then
      c = 2.0d+00 * x * c / (a + b)
    else
      c = 0.0d+00
    end if

!-------------------------------------------------------------------------------
!
  end function limiter_vanleer
!
!===============================================================================
!
! function LIMITER_VANALBADA:
! --------------------------
!
!   Function returns the minimum module value among two arguments using
!   van Albada's limiter.
!
!   Arguments:
!
!     x    - scaling factor;
!     a, b - the input values;
!
!===============================================================================
!
  function limiter_vanalbada(x, a, b) result(c)

! local variables are not implicit by default
!
    implicit none

! input arguments
!
    real(kind=8), intent(in) :: x, a, b
    real(kind=8)             :: c
!
!-------------------------------------------------------------------------------
!
    c = x * a * b * (a + b) / max(eps, a * a + b * b)

!-------------------------------------------------------------------------------
!
  end function limiter_vanalbada
!
!===============================================================================
!
! function MINMOD:
! ===============
!
!   Function returns the minimum module value among two arguments.
!
!   Arguments:
!
!     a, b - the input values;
!
!===============================================================================
!
  real(kind=8) function minmod(a, b)

! local variables are not implicit by default
!
    implicit none

! input arguments
!
    real(kind=8), intent(in) :: a, b
!
!-------------------------------------------------------------------------------
!
    minmod = (sign(0.5d+00, a) + sign(0.5d+00, b)) * min(abs(a), abs(b))

    return

!-------------------------------------------------------------------------------
!
  end function minmod
!
!===============================================================================
!
! function MINMOD4:
! ================
!
!   Function returns the minimum module value among four arguments.
!
!   Arguments:
!
!     a, b, c, d - the input values;
!
!===============================================================================
!
  real(kind=8) function minmod4(a, b, c, d)

! local variables are not implicit by default
!
    implicit none

! input arguments
!
    real(kind=8), intent(in) :: a, b, c, d
!
!-------------------------------------------------------------------------------
!
    minmod4 = minmod(minmod(a, b), minmod(c, d))

    return

!-------------------------------------------------------------------------------
!
  end function minmod4
!
!===============================================================================
!
! function MEDIAN:
! ===============
!
!   Function returns the median of three argument values.
!
!   Arguments:
!
!     a, b, c - the input values;
!
!===============================================================================
!
  real(kind=8) function median(a, b, c)

! local variables are not implicit by default
!
    implicit none

! input arguments
!
    real(kind=8), intent(in) :: a, b, c
!
!-------------------------------------------------------------------------------
!
    median = a + minmod(b - a, c - a)

    return

  end function median
!
!===============================================================================
!
! subroutine FIX_POSITIVITY:
! -------------------------
!
!   Subroutine scans the input arrays of the left and right states for
!   negative values. If it finds one, it repeates the state reconstruction
!   the integrated cell-centered values using the zeroth order interpolation.
!
!   Arguments:
!
!     qi     - the cell centered integrals of variable;
!     ql, qr - the left and right states of variable;
!
!===============================================================================
!
  subroutine fix_positivity(qi, ql, qr)

! local variables are not implicit by default
!
    implicit none

! input/output arguments
!
    real(kind=8), dimension(:), intent(in)    :: qi
    real(kind=8), dimension(:), intent(inout) :: ql, qr

! local variables
!
    integer :: n, i, im1, ip1
!
!------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for positivity fix
!
    call start_timer(imf)
#endif /* PROFILE */

! check positivity only if desired
!
    if (positivity) then

! get the input vector size
!
      n = size(qi)

! look for negative values in the states along the vector
!
      do i = 1, n

! check if the left state has a negative value
!
        if (ql(i) <= 0.0d+00) then

! calculate the left neighbour index
!
          im1 = max(1, i - 1)

! limit the states using the zeroth-order reconstruction
!
          ql(i  ) = qi(i)
          qr(im1) = qi(i)

        end if ! fl ≤ 0

! check if the right state has a negative value
!
        if (qr(i) <= 0.0d+00) then

! calculate the right neighbour index
!
          ip1 = min(n, i + 1)

! limit the states using the zeroth-order reconstruction
!
          ql(ip1) = qi(ip1)
          qr(i  ) = qi(ip1)

        end if ! fr ≤ 0

      end do ! i = 1, n

    end if ! positivity == .true.

#ifdef PROFILE
! stop accounting time for positivity fix
!
    call stop_timer(imf)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine fix_positivity
!
!===============================================================================
!
! subroutine CLIP_EXTREMA:
! -----------------------
!
!   Subroutine scans the reconstructed states and check if they didn't leave
!   the allowed limits.  In the case where the limits where exceeded,
!   the states are limited using constant reconstruction.
!
!   Arguments:
!
!     qi     - the cell centered integrals of variable;
!     ql, qr - the left and right states of variable;
!
!===============================================================================
!
  subroutine clip_extrema(qi, ql, qr)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:), intent(in)    :: qi
    real(kind=8), dimension(:), intent(inout) :: ql, qr

! local variables
!
    integer      :: n, i, im1, ip1, ip2
    real(kind=8) :: qmn, qmx
    real(kind=8) :: dql, dqr, dq
!
!------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for extrema clipping
!
    call start_timer(imc)
#endif /* PROFILE */

! get the input vector size
!
    n = size(qi)

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

! calculate indices
!
      im1 = max(1, i - 1)
      ip1 = min(n, i + 1)

! estimate the bounds of the allowed interval for reconstructed states
!
      qmn = min(qi(i), qi(ip1))
      qmx = max(qi(i), qi(ip1))

! check if the left state lays in the allowed range
!
      if (ql(i) < qmn .or. ql(i) > qmx) then

! calculate the left and right derivatives
!
        dql = qi(i  ) - qi(im1)
        dqr = qi(ip1) - qi(i  )

! get the limited slope
!
        dq  = limiter_clip(0.5d+00, dql, dqr)

! calculate new states
!
        ql(i  ) = qi(i) + dq
        qr(im1) = qi(i) - dq

      end if

! check if the right state lays in the allowed range
!
      if (qr(i) < qmn .or. qr(i) > qmx) then

! calculate the missing index
!
        ip2 = min(n, i + 2)

! calculate the left and right derivatives
!
        dql = qi(ip1) - qi(i  )
        dqr = qi(ip2) - qi(ip1)

! get the limited slope
!
        dq  = limiter_clip(0.5d+00, dql, dqr)

! calculate new states
!
        ql(ip1) = qi(ip1) + dq
        qr(i  ) = qi(ip1) - dq

      end if

    end do ! i = 1, n

#ifdef PROFILE
! stop accounting time for extrema clipping
!
    call stop_timer(imc)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine clip_extrema
!
!===============================================================================
!
! subroutine MP_LIMITING:
! ----------------------
!
!   Subroutine applies the monotonicity preserving (MP) limiter to a vector of
!   high order reconstructed interface values.
!
!   Arguments:
!
!     qc - the vector of cell-centered values;
!     qi - the vector of interface values obtained from the high order
!          interpolation as input and its monotonicity limited values as output;
!
!   References:
!
!     [1] Suresh, A. & Huynh, H. T.,
!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
!          Time Stepping"
!         Journal on Computational Physics,
!         1997, vol. 136, pp. 83-99,
!         http://dx.doi.org/10.1006/jcph.1997.5745
!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
!         "A 5th order monotonicity-preserving upwind compact difference
!          scheme",
!         Science China Physics, Mechanics and Astronomy,
!         Volume 54, Issue 3, pp. 511-522,
!         http://dx.doi.org/10.1007/s11433-010-4220-x
!     [3] Myeong-Hwan Ahn, Duck-Joo Lee,
!         "Modified Monotonicity Preserving Constraints for High-Resolution
!          Optimized Compact Scheme",
!         Journal of Scientific Computing,
!         2020, vol. 83, p. 34
!         https://doi.org/10.1007/s10915-020-01221-0
!
!===============================================================================
!
  subroutine mp_limiting(qc, qi)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), dimension(:), intent(in)    :: qc
    real(kind=8), dimension(:), intent(inout) :: qi

! local variables
!
    logical      :: test
    integer      :: n, i, im2, im1, ip1, ip2
    real(kind=8) :: dq, ds, dc0, dc4, dm1, dp1, dml, dmr, bt
    real(kind=8) :: qlc, qmd, qmp, qmn, qmx, qul

! local vectors
!
    real(kind=8), dimension(0:size(qc)+2) :: dm
!
!-------------------------------------------------------------------------------
!
    n = size(qc)

! 1st order derivatives
!
    dm(0  ) = 0.0d+00
    dm(1  ) = 0.0d+00
    dm(2:n) = qc(2:n) - qc(1:n-1)
    dm(n+1) = 0.0d+00
    dm(n+2) = 0.0d+00

! check the monotonicity condition and apply limiting if necessary
!
    do i = 1, n - 1

      ip1 = i + 1

      if (dm(i) * dm(ip1) >= 0.0d+00) then
        dq = kappa * dm(i)
      else
        dq = kbeta * dm(i)
      end if

      qmp  = qc(i) + minmod(dm(ip1), dq)
      ds   = (qi(i) - qc(i)) * (qi(i) - qmp)

      if (ds > eps) then

        im2 = i - 2
        im1 = i - 1
        ip2 = i + 2

        dm1   = dm(i  ) - dm(im1)
        dc0   = dm(ip1) - dm(i  )
        dp1   = dm(ip2) - dm(ip1)
        dc4   = 4.0d+00 * dc0

        dml   = 0.5d+00 * minmod4(dc4 - dm1, 4.0d+00 * dm1 - dc0, dc0, dm1)
        dmr   = 0.5d+00 * minmod4(dc4 - dp1, 4.0d+00 * dp1 - dc0, dc0, dp1)

        qmd   = qc(i) + 0.5d+00 * dm(ip1) - dmr
        qul   = qc(i) +           dq
        qlc   = qc(i) + 0.5d+00 * dq      + dml

        test = qc(i) > max(qc(im1), qc(ip1)) .and.                             &
                       min(qc(im1), qc(ip1)) > max(qc(im2), qc(ip2))
        test = test .or. qc(i) < min(qc(im1), qc(ip1)) .and.                   &
                       max(qc(im1), qc(ip1)) < min(qc(im2), qc(ip2))
        if (test) then
          if (qc(im2) <= qc(ip1) .and. qc(ip2) <= qc(im1)) then
            bt = dc0 / (qc(ip2) + qc(im2) - 2.0d+00 * qc(i))
            if (bt >= 3.0d-01) qlc = qc(im2)
          end if
        end if

        qmx   = max(min(qc(i), qc(ip1), qmd), min(qc(i), qul, qlc))
        qmn   = min(max(qc(i), qc(ip1), qmd), max(qc(i), qul, qlc))

        qi(i) = median(qi(i), qmn, qmx)

      end if

    end do ! i = 1, n - 1

!-------------------------------------------------------------------------------
!
  end subroutine mp_limiting

!===============================================================================
!
end module interpolations