amun-code/sources/evolution.F90

!!******************************************************************************
!!
!!  This file is part of the AMUN source code, a program to perform
!!  Newtonian or relativistic magnetohydrodynamical simulations on uniform or
!!  adaptive mesh.
!!
!!  Copyright (C) 2008-2021 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License as published by
!!  the Free Software Foundation, either version 3 of the License, or
!!  (at your option) any later version.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!
!!******************************************************************************
!!
!! module: EVOLUTION
!!
!!  This module provides an interface for temporal integration with
!!  the stability handling.  New integration methods can be added by
!!  implementing more evolve_* subroutines.
!!
!!******************************************************************************
!
module evolution

#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, ima, imt, imu, imf, iui, imv
#endif /* PROFILE */

! interfaces for procedure pointers
!
  abstract interface
    subroutine evolve_iface()
    end subroutine
    subroutine update_errors_iface(n, m)
      integer, intent(in) :: n, m
    end subroutine
  end interface

! pointer to the temporal integration subroutine
!
  procedure(evolve_iface), pointer, save :: evolve => null()

! pointer to the error update subroutine
!
  procedure(update_errors_iface), pointer, save :: update_errors => null()

! evolution parameters
!
  logical           , save :: error_control = .false.
  character(len=255), save :: name_int  = ""
  character(len=255), save :: name_norm = ""
  integer           , save :: stages    = 2
  integer           , save :: registers = 1
  real(kind=8)      , save :: cfl       = 5.0d-01

! coefficients controlling the decay of scalar potential ѱ
!
  real(kind=8)      , save :: glm_alpha = 5.0d-01

! time variables
!
  integer           , save :: step     = 0
  real(kind=8)      , save :: time     = 0.0d+00
  real(kind=8)      , save :: dt       = 1.0d+00
  real(kind=8)      , save :: dtn      = 1.0d+00
  real(kind=8)      , save :: dth      = 1.0d+00
  real(kind=8)      , save :: dte      = 0.0d+00

! the absolute and relative tolerances, limiting factors, the maximum error,
! the maximum number of passes for the adaptive step,
! the total number of integration iterations, the number of rejected iterations
!
  real(kind=8)      , save :: atol        = 1.0d-04
  real(kind=8)      , save :: rtol        = 1.0d-04
  real(kind=8)      , save :: fac         = 9.0d-01
  real(kind=8)      , save :: facmin      = 1.0d-01
  real(kind=8)      , save :: facmax      = 5.0d+00
  real(kind=8)      , save :: errtol      = 1.0d+00
  real(kind=8)      , save :: chi         = 1.0d+00
  integer           , save :: mrej        = 5
  integer           , save :: niterations = 0
  integer           , save :: nrejections = 0

! errors in three recent steps
!
  real(kind=8), dimension(3), save :: betas, errs

! GLM variables
!
  real(kind=8), dimension(:), allocatable, save :: adecay, aglm

! by default everything is private
!
  private

! declare public subroutines
!
  public :: initialize_evolution, finalize_evolution, print_evolution
  public :: advance, new_time_step

! declare public variables
!
  public :: step, time, dt, dtn, dth, dte, cfl, glm_alpha, registers
  public :: atol, rtol, mrej, niterations, nrejections, errs, errtol

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
  contains
!
!===============================================================================
!!
!!***  PUBLIC SUBROUTINES  *****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine INITIALIZE_EVOLUTION:
! -------------------------------
!
!   Subroutine initializes module EVOLUTION by setting its parameters.
!
!   Arguments:
!
!     verbose - a logical flag turning the information printing;
!     status  - an integer flag for error return value;
!
!===============================================================================
!
  subroutine initialize_evolution(verbose, status)

! include external procedures
!
    use coordinates, only : toplev, adx, ady, adz
    use parameters , only : get_parameter

! local variables are not implicit by default
!
    implicit none

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

! local variables
!
    character(len=255) :: integration = "rk2", error_norm = "l2"
    integer            :: n

! local parameters
!
    character(len=*), parameter :: loc = 'EVOLUTION::initialize_evolution()'
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! set timer descriptions
!
    call set_timer('evolution:: initialization'  , imi)
    call set_timer('evolution:: solution advance', ima)
    call set_timer('evolution:: new time step'   , imt)
    call set_timer('evolution:: solution update' , imu)
    call set_timer('evolution:: flux update'     , imf)
    call set_timer('evolution:: increment update', iui)
    call set_timer('evolution:: variable update' , imv)

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

! reset the status flag
!
    status = 0

! get the integration method and the value of the CFL coefficient
!
    call get_parameter("time_advance", integration)
    call get_parameter("stages"      , stages     )
    call get_parameter("cfl"         , cfl        )
    call get_parameter("glm_alpha"   , glm_alpha  )

! get the absolute and relative tolerances, the maximum number of passes for
! the adaptive step, and limiting factors
!
    call get_parameter("absolute_tolerance", atol)
    call get_parameter("relative_tolerance", rtol)
    call get_parameter("maximum_rejections", mrej)
    call get_parameter("chi"               , chi)
    call get_parameter("factor"            , fac )
    call get_parameter("minimum_factor"    , facmin)
    call get_parameter("maximum_factor"    , facmax)

! select the integration method, check the correctness of the integration
! parameters and adjust the CFL coefficient if necessary
!
    select case(trim(integration))

    case ("euler", "EULER", "rk1", "RK1")

      evolve    => evolve_euler
      registers = 1
      name_int  = "1st order Euler"

    case ("rk2", "RK2")

      evolve    => evolve_rk2
      registers = 2
      name_int  = "2nd order Runge-Kutta"

    case ("ssprk(m,2)", "SSPRK(m,2)")

      error_control = .true.
      evolve    => evolve_ssprk2_m
      registers = 3
      stages    = max(2, min(9, stages))
      cfl       = (stages - 1) * cfl
      write(name_int, "('Optimal 2nd order SSPRK(',i0,',2) with error control')") stages

    case ("rk3", "RK3")

      evolve    => evolve_rk3
      registers = 2
      name_int  = "3rd order Runge-Kutta"

    case ("rk3.4", "RK3.4", "ssprk(4,3)", "SSPRK(4,3)")

      evolve    => evolve_ssprk34
      registers = 2
      cfl       = 2.0d+00 * cfl
      name_int  = "3rd order SSPRK(4,3)"

    case ("rk3.5", "RK3.5", "ssprk(5,3)", "SSPRK(5,3)")

      evolve    => evolve_ssprk35
      registers = 2
      cfl       = 2.65062919143939d+00 * cfl
      name_int  = "3rd order SSPRK(5,3)"

    case ("SSPRK3(2)4", "SSP3(2)4", "ssprk3(2)4", "ssp3(2)4")

      error_control = .true.
      evolve    => evolve_ssp324
      registers = 3
      stages    = 4
      cfl       = 2.0d+00 * cfl
      name_int  = "3ʳᵈ-order 4-step embedded SSP3(2)4 method"

      betas(:) = [ 5.5d-01, -2.7d-01, 5.0d-02 ]

      call get_parameter("beta(1)", betas(1))
      call get_parameter("beta(2)", betas(2))
      call get_parameter("beta(3)", betas(3))

      betas(:) = - betas(:) / 3.0d+00

    case ("rk3.m", "ssprk(m,3)", "SSPRK(m,3)")

      error_control = .true.
      evolve   => evolve_ssprk3_m
      registers = 3
      n = 2
      do while(n**2 <= stages)
        n = n + 1
      end do
      n = n - 1
      stages    = max(4, n**2)
      cfl       = (n - 1) * n * cfl
      write(name_int, "('Optimal 3rd order SSPRK(',i0,',3) with error control')") stages

    case ("rk4.10", "RK4.10", "ssprk(10,4)", "SSPRK(10,4)")

      evolve    => evolve_ssprk4_10
      registers = 2
      cfl       = 6.0d+00 * cfl
      name_int  = "Optimal 4th order SSPRK(10,4)"

    case default

      if (verbose) then
        write(*,*)
        write(*,"(1x,a)") "ERROR!"
        write(*,"(1x,a)") "The selected time advance method is not " //        &
                          "implemented: " // trim(integration)
        write(*,"(1x,a)") "Available methods: 'euler' 'rk2', 'rk3'," //        &
                          " 'ssprk(m,2)', 'ssprk(4,3)', 'ssprk(5,3)'," //      &
                          " 'ssprk(m,3)', 'ssprk(10,4)'."
      end if
      status = 1

    end select

! select error norm
!
    call get_parameter("error_norm", error_norm)

    select case(trim(error_norm))
    case("max", "maximum", "infinity")
      update_errors => update_errors_max
      name_norm = "maximum norm"
    case default
      update_errors => update_errors_l2
      name_norm = "L2-norm"
    end select

! reset recent error history
!
    errs = 1.0d+00

! GLM coefficients for all refinement levels
!
    allocate(adecay(toplev), aglm(toplev))
#if NDIMS == 3
    aglm(:) = - glm_alpha / min(adx(:), ady(:), adz(:))
#else /* NDIMS == 3 */
    aglm(:) = - glm_alpha / min(adx(:), ady(:))
#endif /* NDIMS == 3 */

#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
    call stop_timer(imi)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
  end subroutine initialize_evolution
!
!===============================================================================
!
! subroutine FINALIZE_EVOLUTION:
! -----------------------------
!
!   Subroutine releases memory used by the module.
!
!   Arguments:
!
!     status  - an integer flag for error return value;
!
!===============================================================================
!
  subroutine finalize_evolution(verbose, status)

! include external procedures
!
    use helpers        , only : print_section, print_parameter
    use iso_fortran_env, only : error_unit

! local variables are not implicit by default
!
    implicit none

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

! local parameters
!
    character(len=*), parameter :: loc = 'EVOLUTION::finalize_evolution()'
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for module initialization/finalization
!
    call start_timer(imi)
#endif /* PROFILE */

    if (allocated(adecay)) deallocate(adecay)
    if (allocated(aglm))   deallocate(aglm)

! reset the status flag
!
    status = 0

! nullify pointers
!
    nullify(evolve)
    nullify(update_errors)

! print integration summary
!
    if (niterations > 0) then
      call print_section(verbose, "Integration summary")
      call print_parameter(verbose, "Number of iterations", niterations)
      call print_parameter(verbose, "Number of rejections", nrejections)
    end if

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

!-------------------------------------------------------------------------------
!
  end subroutine finalize_evolution
!
!===============================================================================
!
! subroutine PRINT_EVOLUTION:
! --------------------------
!
!   Subroutine prints module parameters.
!
!   Arguments:
!
!     verbose - a logical flag turning the information printing;
!
!===============================================================================
!
  subroutine print_evolution(verbose)

! import external procedures and variables
!
    use equations, only : magnetized
    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, "Evolution")
      call print_parameter(verbose, "time advance", name_int)
      call print_parameter(verbose, "no. of registers", registers)
      call print_parameter(verbose, "CFL coefficient", cfl)
      if (magnetized) then
        call print_parameter(verbose, "GLM alpha coefficient", glm_alpha)
      end if
      if (error_control) then
        call print_parameter(verbose, "absolute tolerance", atol)
        call print_parameter(verbose, "relative tolerance", rtol)
        call print_parameter(verbose, "error norm"        , name_norm)
        call print_parameter(verbose, "maximum rejections", mrej)
        call print_parameter(verbose, "chi"               , chi)
        call print_parameter(verbose, "factor"            , fac)
        call print_parameter(verbose, "minimum factor"    , facmin)
        call print_parameter(verbose, "maximum factor"    , facmax)
      end if

    end if

!-------------------------------------------------------------------------------
!
  end subroutine print_evolution
!
!===============================================================================
!
! subroutine ADVANCE:
! ------------------
!
!   Subroutine advances the solution by one time step using the selected time
!   integration method.
!
!   Arguments:
!
!     dtnext - the next time step;
!     status - the subroutine call status: 0 for success, otherwise failure;
!
!===============================================================================
!
  subroutine advance(dtnext, status)

! references
!
    use blocks     , only : set_blocks_update
    use coordinates, only : toplev
    use forcing    , only : update_forcing, forcing_enabled
    use mesh       , only : update_mesh

! local variables are not implicit by default
!
    implicit none

! input variables
!
    real(kind=8), intent(in)  :: dtnext
    integer     , intent(out) :: status
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for solution advance
!
    call start_timer(ima)
#endif /* PROFILE */

! find new time step
!
    call new_time_step(dtnext)

! advance the solution using the selected method
!
    call evolve()

! add forcing contribution
!
    if (forcing_enabled) call update_forcing(dt)

! check if we need to perform the refinement step
!
    if (toplev > 1) then

! set all meta blocks to not be updated
!
      call set_blocks_update(.false.)

! check refinement and refine
!
      call update_mesh(status)
      if (status /= 0) go to 100

! update primitive variables
!
      call update_variables(time + dt, 0.0d+00)

! set all meta blocks to be updated
!
      call set_blocks_update(.true.)

    end if ! toplev > 1

#ifdef DEBUG
! check variables for NaNs
!
    call check_variables()
#endif /* DEBUG */

! error entry point
!
    100 continue

#ifdef PROFILE
! stop accounting time for solution advance
!
    call stop_timer(ima)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine advance
!
!===============================================================================
!
! subroutine NEW_TIME_STEP:
! ------------------------
!
!   Subroutine estimates the new time step from the maximum speed in the system
!   and source term constraints.
!
!   Arguments:
!
!     dtnext - next time step;
!
!===============================================================================
!
  subroutine new_time_step(dtnext)

! include external procedures
!
    use equations     , only : maxspeed, cmax, cmax2
#ifdef MPI
    use mpitools      , only : reduce_maximum
#endif /* MPI */

! include external variables
!
    use blocks        , only : block_data, list_data
    use coordinates   , only : adx, ady
#if NDIMS == 3
    use coordinates   , only : adz
#endif /* NDIMS == 3 */
    use sources       , only : viscosity, resistivity

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), intent(in)  :: dtnext

! local pointers
!
    type(block_data), pointer :: pdata

! local variables
!
    integer                   :: mlev
    real(kind=8)              :: cm, dx_min

! local parameters
!
    real(kind=8), parameter   :: eps = tiny(cmax)
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for new time step estimation
!
    call start_timer(imt)
#endif /* PROFILE */

! reset the maximum speed, and the highest level
!
    cmax   = eps
    mlev   = 1

! assign pdata with the first block on the data block list
!
    pdata => list_data

! iterate over all data blocks in order to find the maximum speed among them
! and the highest level which is required to obtain the minimum spacial step
!
    do while (associated(pdata))

! update the maximum level
!
      mlev = max(mlev, pdata%meta%level)

! get the maximum speed for the current block
!
      cm = maxspeed(pdata%q(:,:,:,:))

! update the maximum speed
!
      cmax = max(cmax, cm)

! assign pdata to the next block
!
      pdata => pdata%next

    end do ! over data blocks

#ifdef MPI
! reduce maximum speed and level over all processors
!
    call reduce_maximum(cmax)
    call reduce_maximum(mlev)
#endif /* MPI */

! calculate the squared cmax
!
    cmax2 = cmax * cmax

! find the smallest spacial step
!
#if NDIMS == 2
    dx_min = min(adx(mlev), ady(mlev))
#endif /* NDIMS == 2 */
#if NDIMS == 3
    dx_min = min(adx(mlev), ady(mlev), adz(mlev))
#endif /* NDIMS == 3 */

! calculate the new time step
!
    dth = cfl * dx_min / max(cmax                                              &
                        + 2.0d+00 * max(viscosity, resistivity) / dx_min, eps)

! round the time
!
    dtn = dth
    if (error_control .and. dte > 0.0d+00) dtn = min(dtn, dte)
    if (dtnext > 0.0d+00) dtn = min(dtn, dtnext)
    dt = dtn

#ifdef PROFILE
! stop accounting time for new time step estimation
!
    call stop_timer(imt)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine new_time_step
!
!===============================================================================
!!
!!***  PRIVATE SUBROUTINES  ****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine EVOLVE_EULER:
! -----------------------
!
!   Subroutine advances the solution by one time step using the 1st order
!   Euler integration method.
!
!   References:
!
!     [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P.,
!         "Numerical Recipes in Fortran",
!         Cambridge University Press, Cambridge, 1992
!
!===============================================================================
!
  subroutine evolve_euler()

    use blocks    , only : block_data, list_data
    use boundaries, only : boundary_fluxes
    use equations , only : ibp, cmax
    use sources   , only : update_sources

    implicit none

    type(block_data), pointer :: pdata

    real(kind=8) :: tm, dtm
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
    call start_timer(imu)
#endif /* PROFILE */

    tm  = time + dt
    dtm = dt

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dt * pdata%du(:,:,:,:)

      pdata%u => pdata%uu(:,:,:,:,1)

      pdata => pdata%next
    end do

    if (ibp > 0) then
      adecay(:) = exp(aglm(:) * cmax * dt)

      pdata => list_data
      do while (associated(pdata))

        pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)

        pdata => pdata%next
      end do
    end if

    call update_variables(tm, dtm)

#ifdef PROFILE
    call stop_timer(imu)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine evolve_euler
!
!===============================================================================
!
! subroutine EVOLVE_RK2:
! ---------------------
!
!   Subroutine advances the solution by one time step using the 2nd order
!   Runge-Kutta time integration method.
!
!   References:
!
!     [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P.,
!         "Numerical Recipes in Fortran",
!         Cambridge University Press, Cambridge, 1992
!
!===============================================================================
!
  subroutine evolve_rk2()

    use blocks    , only : block_data, list_data
    use boundaries, only : boundary_fluxes
    use equations , only : ibp, cmax
    use sources   , only : update_sources

    implicit none

    type(block_data), pointer :: pdata

    real(kind=8) :: tm, dtm
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
    call start_timer(imu)
#endif /* PROFILE */

!= 1st step: U(1) = U(n) + dt * F[U(n)]
!
    tm  = time + dt
    dtm = dt

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + dt * pdata%du(:,:,:,:)

      pdata%u => pdata%uu(:,:,:,:,2)

      pdata => pdata%next
    end do

    call update_variables(tm, dtm)

!= 2nd step: U(n+1) = 1/2 U(n) + 1/2 U(1) + 1/2 dt * F[U(1)]
!
    tm  = time + dt
    dtm = 0.5d+00 * dt

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,1) = 0.5d+00 * (pdata%uu(:,:,:,:,1)                     &
                                     + pdata%uu(:,:,:,:,2)                     &
                               +  dt * pdata%du(:,:,:,:))

      pdata%u => pdata%uu(:,:,:,:,1)

      pdata => pdata%next
    end do

    if (ibp > 0) then
      adecay(:) = exp(aglm(:) * cmax * dt)

      pdata => list_data
      do while (associated(pdata))

        pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)

        pdata => pdata%next
      end do
    end if

    call update_variables(tm, dtm)

#ifdef PROFILE
    call stop_timer(imu)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine evolve_rk2
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK2_M:
! --------------------------
!
!   Subroutine advances the solution by one time step using the 2nd order
!   m-stage Strong Stability Preserving Runge-Kutta time integration method.
!   Up to 9 stages are allowed, due to stability problems with more stages.
!
!   References:
!
!     [1] Gottlieb, S. and Gottlieb, L.-A., J.
!         "Strong Stability Preserving Properties of Runge-Kutta Time
!          Discretization Methods for Linear Constant Coefficient Operators",
!         Journal of Scientific Computing,
!         2003, vol. 18, no. 1, pp. 83-109
!
!===============================================================================
!
  subroutine evolve_ssprk2_m()

    use blocks     , only : block_data, list_data, get_dblocks
    use boundaries , only : boundary_fluxes
    use coordinates, only : nc => ncells, nb, ne
    use equations  , only : errors, nf, ibp, cmax
#ifdef MPI
    use mpitools   , only : reduce_maximum
#endif /* MPI */
    use sources    , only : update_sources

    implicit none

    type(block_data), pointer :: pdata

    logical      :: test
    integer      :: n, l, nrej
    real(kind=8) :: tm, dtm, ds, umax, emax
    real(kind=8) :: fc, fcmn, fcmx

    logical     , save :: first = .true.
    real(kind=8), save :: ft, fl, fr, gt, gl

    real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr

    real(kind=8), parameter :: k1 = -2.9d-01, k2 = 1.05d-01, k3 = -5.0d-02
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
    call start_timer(imu)
#endif /* PROFILE */

    if (first) then

      ft = 1.0d+00 / (stages - 1)
      fl = 1.0d+00 / stages
      fr = 1.0d+00 - fl

      gt = fl - ft
      gl = fl

      first = .false.

    end if

    fc   = fac
    fcmn = facmin
    fcmx = facmax

    l = get_dblocks()
#if NDIMS == 3
    allocate(lerr(l,nf,nc,nc,nc))
#else /* NDIMS == 3 */
    allocate(lerr(l,nf,nc,nc, 1))
#endif /* NDIMS == 3 */

    test = .true.
    nrej = 0

    do while(test)

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

      ds = ft * dt

!= 1st step: U(1) = U(n), U(2) = U(n)
!
      pdata => list_data
      do while (associated(pdata))

        pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)
        pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,1)

        pdata%u => pdata%uu(:,:,:,:,2)

        pdata => pdata%next
      end do

!= 2nd step: U(1) = U(1) + dt/(m-1) F[U(1)], for i = 1, ..., m-1
!
      do n = 1, stages - 1

        tm  = time + n * ds
        dtm = ds

        pdata => list_data
        do while (associated(pdata))

          call update_increment(pdata)
          call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

          pdata => pdata%next
        end do

        call boundary_fluxes()

        l = 1
        pdata => list_data
        do while (associated(pdata))

          pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)

#if NDIMS == 3
          lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * pdata%du(:,nb:ne,nb:ne,nb:ne)
#else /* NDIMS == 3 */
          lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * pdata%du(:,nb:ne,nb:ne,  :  )
#endif /* NDIMS == 3 */
          l = l + 1

          pdata => pdata%next
        end do

        call update_variables(tm, dtm)

      end do ! n = 1, stages - 1

!= the final step: U(1) = (m-1)/m U(1) + 1/m U(2) + dt/m F[U(1)]
!
      tm  = time + dt
      dtm = fl * dt

      pdata => list_data
      do while (associated(pdata))

        call update_increment(pdata)
        call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

        pdata => pdata%next
      end do

      call boundary_fluxes()

      l = 1
      pdata => list_data
      do while (associated(pdata))

        pdata%uu(:,:,:,:,2) = fr *  pdata%uu(:,:,:,:,2)                        &
                            + fl * (pdata%uu(:,:,:,:,3)                        &
                            + dt *  pdata%du(:,:,:,:))

#if NDIMS == 3
        lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * pdata%du(:,nb:ne,nb:ne,nb:ne)
#else /* NDIMS == 3 */
        lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * pdata%du(:,nb:ne,nb:ne,  :  )
#endif /* NDIMS == 3 */
        l = l + 1

        pdata => pdata%next
      end do

      call update_variables(tm, dtm)

! find umax
!
      umax = 0.0d+00
      pdata => list_data
      do while (associated(pdata))
#if NDIMS == 3
        umax = max(umax, maxval(abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1))),         &
                         maxval(abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,2))))
#else /* NDIMS == 3 */
        umax = max(umax, maxval(abs(pdata%uu(:,nb:ne,nb:ne,  :  ,1))),         &
                         maxval(abs(pdata%uu(:,nb:ne,nb:ne,  :  ,2))))
#endif /* NDIMS == 3 */
        pdata => pdata%next
      end do

! get the maximum error of each variable
!
      do l = 1, nf
        errors(l) = dt * maxval(abs(lerr(:,l,:,:,:)))
      end do

#ifdef MPI
      call reduce_maximum(umax)
      call reduce_maximum(errors)
#endif /* MPI */

! calculate tolerance and time step
!
      emax   = maxval(errors) / (atol + rtol * umax)

      if (emax <= 1.0d+00 .or. nrej >= mrej) then
        test = .false.

        errs(3) = errs(2)
        errs(2) = errs(1)
        errs(1) = emax

        errtol  = emax

        dte = dt * min(fcmx, max(fcmn,                                         &
                                 fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))

        fcmx = facmax
      else
        errs(1) = emax

        dte = dt * min(fcmx, max(fcmn,                                         &
                                 fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))
        dt  = dte

        fcmx = fac

        nrej        = nrej        + 1 ! rejection count in the current step
        nrejections = nrejections + 1
      end if

      niterations = niterations + 1

    end do

    if (allocated(lerr)) deallocate(lerr)

!= final step: U(n+1) = U(1)
!
    tm  = time + dt
    dtm = ft * dt

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)

      pdata%u => pdata%uu(:,:,:,:,1)

      pdata => pdata%next
    end do

    if (ibp > 0) then
      adecay(:) = exp(aglm(:) * cmax * dt)

      pdata => list_data
      do while (associated(pdata))

        pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)

        pdata => pdata%next
      end do
    end if

    call update_variables(tm, dtm)

#ifdef PROFILE
    call stop_timer(imu)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine evolve_ssprk2_m
!
!===============================================================================
!
! subroutine EVOLVE_RK3:
! ---------------------
!
!   Subroutine advances the solution by one time step using the 3rd order
!   Runge-Kutta time integration method.
!
!   References:
!
!     [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P.,
!         "Numerical Recipes in Fortran",
!         Cambridge University Press, Cambridge, 1992
!
!
!===============================================================================
!
  subroutine evolve_rk3()

    use blocks    , only : block_data, list_data
    use boundaries, only : boundary_fluxes
    use equations , only : ibp, cmax
    use sources   , only : update_sources

    implicit none

    type(block_data), pointer :: pdata

    real(kind=8) :: ds
    real(kind=8) :: tm, dtm

    real(kind=8), parameter :: f21 = 3.0d+00 / 4.0d+00, f22 = 1.0d+00 / 4.0d+00
    real(kind=8), parameter :: f31 = 1.0d+00 / 3.0d+00, f32 = 2.0d+00 / 3.0d+00
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
    call start_timer(imu)
#endif /* PROFILE */

!= 1st substep: U(1) = U(n) + dt F[U(n)]
!
    ds  = dt
    tm  = time + ds
    dtm = ds

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:)

      pdata%u => pdata%uu(:,:,:,:,2)

      pdata => pdata%next
    end do

    call update_variables(tm, dtm)

!= 2nd step: U(2) = 3/4 U(n) + 1/4 U(1) + 1/4 dt F[U(1)]
!
    ds  = f22 * dt
    tm  = time + 0.5d+00 * dt
    dtm = ds

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,2) = f21 * pdata%uu(:,:,:,:,1)                          &
                          + f22 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)

      pdata => pdata%next
    end do

    call update_variables(tm, dtm)

!= 3rd step: U(n+1) = 1/3 U(n) + 2/3 U(2) + 2/3 dt F[U(2)]
!
    ds  = f32 * dt
    tm  = time + dt
    dtm = ds

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,1) = f31 * pdata%uu(:,:,:,:,1)                          &
                          + f32 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)

      pdata%u => pdata%uu(:,:,:,:,1)

      pdata => pdata%next
    end do

    if (ibp > 0) then
      adecay(:) = exp(aglm(:) * cmax * dt)

      pdata => list_data
      do while (associated(pdata))

        pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)

        pdata => pdata%next
      end do
    end if

    call update_variables(tm, dtm)

#ifdef PROFILE
    call stop_timer(imu)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine evolve_rk3
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK34:
! -------------------------
!
!   Subroutine advances the solution by one time step using the 3rd order
!   4-stage Strong Stability Preserving Runge-Kutta time integration method.
!
!   References:
!
!     [1] Ruuth, S. J.,
!         "Global Optimization of Explicit Strong-Stability-Preserving
!          Runge-Kutta methods",
!         Mathematics of Computation,
!         2006, vol. 75, no. 253, pp. 183-207
!
!===============================================================================
!
  subroutine evolve_ssprk34()

    use blocks    , only : block_data, list_data
    use boundaries, only : boundary_fluxes
    use equations , only : ibp, cmax
    use sources   , only : update_sources

    implicit none

    type(block_data), pointer :: pdata

    real(kind=8) :: ds
    real(kind=8) :: tm, dtm
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
    call start_timer(imu)
#endif /* PROFILE */

!= 1st step: U(1) = U(n) + 1/2 dt F[U(n)]
!
    ds  = dt / 2.0d+00
    tm  = time + ds
    dtm = ds

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:)

      pdata%u => pdata%uu(:,:,:,:,2)

      pdata => pdata%next
    end do

    call update_variables(tm, dtm)

!= 2nd step: U(2) = U(2) + 1/2 dt F[U(1)]
!
    tm  = time + dt

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dtm * pdata%du(:,:,:,:)

      pdata => pdata%next
    end do

    call update_variables(tm, dtm)

!= 3rd step: U(3) = 2/3 U(n) + 1/3 (U(2) + 1/2 dt F[U(2)])
!
    tm  = time + ds

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,2) = (2.0d+00 * pdata%uu(:,:,:,:,1)                     &
                                     + pdata%uu(:,:,:,:,2)                     &
                                     + dtm * pdata%du(:,:,:,:)) / 3.0d+00

      pdata => pdata%next
    end do

    call update_variables(tm, dtm)

!= the final step: U(n+1) = U(3) + 1/2 dt F[U(3)]
!
    tm  = time + dt

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)

      pdata%u => pdata%uu(:,:,:,:,1)

      pdata => pdata%next
    end do

    if (ibp > 0) then
      adecay(:) = exp(aglm(:) * cmax * dt)

      pdata => list_data
      do while (associated(pdata))

        pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)

        pdata => pdata%next
      end do
    end if

    call update_variables(tm, dtm)

#ifdef PROFILE
    call stop_timer(imu)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine evolve_ssprk34
!
!===============================================================================
!
! subroutine EVOLVE_SSP324:
! --------------------------
!
!   Subroutine advances the solution by one time step using the 3ʳᵈ-order
!   4-stage embedded Strong Stability Preserving Runge-Kutta time integration
!   method SSP3(2)4 with the error control.
!
!   References:
!
!     [1] Ranocha, H., Dalcin, L., Parsani, M., Ketcheson, D. I.,
!         "Optimized Runge-Kutta Methods with Automatic Step Size Control
!          for Compressible Computational Fluid Dynamics",
!         2021, arXiv:2104.06836v2
!     [2] Conde, S., Fekete, I., Shadid, J. N.,
!         "Embedded error estimation and adaptive step-size control for
!          optimal explicit strong stability preserving Runge–Kutta methods"
!         2018, arXiv:1806.08693v1
!     [3] Gottlieb, S., Ketcheson, D., Shu, C.-W.,
!         "Strong stability preserving Runge-Kutta and multistep time
!          discretizations",
!         2011, World Scientific Publishing
!
!===============================================================================
!
  subroutine evolve_ssp324()

    use blocks     , only : block_data, list_data
    use boundaries , only : boundary_fluxes
    use coordinates, only : nb, ne
    use equations  , only : errors, ibp, cmax
    use sources    , only : update_sources

    implicit none

    type(block_data), pointer :: pdata

    logical         :: test
    integer         :: nrej, i
    real(kind=8)    :: tm, dtm, dh, fc

    real(kind=8), parameter :: onethird = 1.0d+00 / 3.0d+00,                   &
                               twothird = 2.0d+00 / 3.0d+00
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
    call start_timer(imu)
#endif /* PROFILE */

    test = .true.
    nrej = 0

! at the entry point we assume the previous solution of conserved variables U(n)
! is stored in pdata%u(:,:,:,:,1) and the primitive variables are already
! updated from this array;
!
! pdata%uu(:,:,:,:,1) - the previous solution, U(n)
! pdata%uu(:,:,:,:,2) - the new 3rd order solution, U(1)
! pdata%uu(:,:,:,:,3) - the new 2nd order solution, U(2)
!
    do while(test)

! initiate the fractional time step
!
      dh = 5.0d-01 * dt

!= preparation step: U(1) = U(n)
!
      pdata => list_data
      do while (associated(pdata))

        pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)

        pdata%u => pdata%uu(:,:,:,:,2)

        pdata => pdata%next
      end do

!= 1st to 3rd steps: U(1) = U(1) + ½ dt F[U(1)], for i = 1...3
!
      do i = 1, 3

        tm  = time + (i - 1) * dh
        dtm = dh

        pdata => list_data
        do while (associated(pdata))

          call update_increment(pdata)
          call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

          pdata => pdata%next
        end do

        call boundary_fluxes()

        pdata => list_data
        do while (associated(pdata))

          pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:)

          pdata => pdata%next
        end do

        call update_variables(tm, dtm)

      end do ! i = 1, 3

!= 4th step: U(1) = ⅔ U(n) + ⅓ U(1), U(2) = ⅓ U(n) + ⅔ U(1)
!
      pdata => list_data
      do while (associated(pdata))

        pdata%uu(:,:,:,:,3) = onethird * pdata%uu(:,:,:,:,1)                   &
                            + twothird * pdata%uu(:,:,:,:,2)
        pdata%uu(:,:,:,:,2) = twothird * pdata%uu(:,:,:,:,1)                   &
                            + onethird * pdata%uu(:,:,:,:,2)

        pdata => pdata%next
      end do

      call update_variables(tm, dtm)

!= 5th step: U(1) = U(1) + ½ dt F[U(1)] <- 3ʳᵈ-order candidate
!
      tm  = time + dh
      dtm = dh

      pdata => list_data
      do while (associated(pdata))

        call update_increment(pdata)
        call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

        pdata => pdata%next
      end do

      call boundary_fluxes()

      pdata => list_data
      do while (associated(pdata))

        pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:)

        pdata => pdata%next
      end do

      call update_variables(tm, dtm)

!= 6th step: U(2) = ½ (U(1) + U(2)) <- 2ⁿᵈ-order approximation
!
      pdata => list_data
      do while (associated(pdata))

        pdata%uu(:,:,:,:,3) = 5.0d-01 * (pdata%uu(:,:,:,:,3)                   &
                                       + pdata%uu(:,:,:,:,2))

        pdata => pdata%next
      end do

!= 7th step: decay the magnetic divergence potential of both solutions
!
      if (ibp > 0) then
        adecay(:) = exp(aglm(:) * cmax * dt)

        pdata => list_data
        do while (associated(pdata))

          pdata%uu(ibp,:,:,:,2) = adecay(pdata%meta%level)                     &
                                                       * pdata%uu(ibp,:,:,:,2)
          pdata%uu(ibp,:,:,:,3) = adecay(pdata%meta%level)                     &
                                                       * pdata%uu(ibp,:,:,:,3)

          pdata => pdata%next
        end do
      end if

! update the vector of errors
!
      call update_errors(2, 3)

! calculate the tolerance and estimate the next time step due to the error
!
      errtol  = maxval(errors)
      errs(1) = errtol
      fc      = product(errs(:)**betas(:))
      fc      = 1.0d+00 + chi * atan((fc - 1.0d+00) / chi)
      dte     = dt * fc

      if (fc > fac .or. nrej >= mrej) then
        test = .false.

        errs(3) = errs(2)
        errs(2) = errs(1)

      else
        dt  = dte

        nrej        = nrej        + 1 ! rejection count in the current step
        nrejections = nrejections + 1

! since the solution was rejected, we have to revert the primitive variables
! to the previous state
!
        pdata => list_data
        do while (associated(pdata))

          pdata%u => pdata%uu(:,:,:,:,1)

          pdata => pdata%next
        end do

        call update_variables(tm, dtm)

      end if

      niterations = niterations + 1

    end do

!= final step: U(n+1) = U(1) - update the accepted solution
!
    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)

      pdata%u => pdata%uu(:,:,:,:,1)

      pdata => pdata%next
    end do

#ifdef PROFILE
    call stop_timer(imu)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine evolve_ssp324
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK35:
! -------------------------
!
!   Subroutine advances the solution by one time step using the 3rd order
!   5-stage Strong Stability Preserving Runge-Kutta time integration method.
!
!   References:
!
!     [1] Ruuth, S. J.,
!         "Global Optimization of Explicit Strong-Stability-Preserving
!          Runge-Kutta methods",
!         Mathematics of Computation,
!         2006, vol. 75, no. 253, pp. 183-207
!
!===============================================================================
!
  subroutine evolve_ssprk35()

    use blocks    , only : block_data, list_data
    use boundaries, only : boundary_fluxes
    use equations , only : ibp, cmax
    use sources   , only : update_sources

    implicit none

    type(block_data), pointer :: pdata

    real(kind=8) :: ds
    real(kind=8) :: tm, dtm

    real(kind=8), parameter :: b1  = 3.77268915331368d-01
    real(kind=8), parameter :: b3  = 2.42995220537396d-01
    real(kind=8), parameter :: b4  = 2.38458932846290d-01
    real(kind=8), parameter :: b5  = 2.87632146308408d-01
    real(kind=8), parameter :: a31 = 3.55909775063327d-01
    real(kind=8), parameter :: a33 = 6.44090224936673d-01
    real(kind=8), parameter :: a41 = 3.67933791638137d-01
    real(kind=8), parameter :: a44 = 6.32066208361863d-01
    real(kind=8), parameter :: a53 = 2.37593836598569d-01
    real(kind=8), parameter :: a55 = 7.62406163401431d-01
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
    call start_timer(imu)
#endif /* PROFILE */

!= 1st step: U(1) = U(n) + b1 dt F[U(n)]
!
    ds  = b1 * dt
    tm  = time + ds
    dtm = ds

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:)

      pdata%u => pdata%uu(:,:,:,:,2)

      pdata => pdata%next
    end do

    call update_variables(tm, dtm)

!= 2nd step: U(2) = U(1) + b1 dt F[U(1)]
!
    tm  = time + 2.0d+00 * ds
    dtm = ds

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)

      pdata => pdata%next
    end do

    call update_variables(tm, dtm)

!= 3rd step: U(3) = a31 U(n) + a33 U(2) + b3 dt F[U(2)]
!
    ds  = b3 * dt
    tm  = time + (2.0d+00 * a33 * b1 + b3) * dt
    dtm = ds

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,2) = a31 * pdata%uu(:,:,:,:,1)                          &
                          + a33 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)

      pdata => pdata%next
    end do

    call update_variables(tm, dtm)

!= 4th step: U(4) = a41 U(n) + a44 U(3) + b4 dt F[U(3)]
!
    ds  = b4 * dt
    tm  = time + ((2.0d+00 * b1 * a33 + b3) * a44 + b4) * dt
    dtm = ds

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,1) = a41 * pdata%uu(:,:,:,:,1)                          &
                          + a44 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)

      pdata%u => pdata%uu(:,:,:,:,1)

      pdata => pdata%next
    end do

    call update_variables(tm, dtm)

!= the final step: U(n+1) = a53 U(2) + a55 U(4) + b5 dt F[U(4)]
!
    ds  = b5 * dt
    tm  = time + dt
    dtm = ds

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,1) = a53 * pdata%uu(:,:,:,:,2)                          &
                          + a55 * pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:)

      pdata => pdata%next
    end do

    if (ibp > 0) then
      adecay(:) = exp(aglm(:) * cmax * dt)

      pdata => list_data
      do while (associated(pdata))

        pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)

        pdata => pdata%next
      end do
    end if

    call update_variables(tm, dtm)

#ifdef PROFILE
    call stop_timer(imu)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine evolve_ssprk35
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK3_M:
! --------------------------
!
!   Subroutine advances the solution by one time step using the 3rd order
!   m-stage Strong Stability Preserving Runge-Kutta time integration method.
!
!   References:
!
!     [1] Gottlieb, S., Ketcheson, D., Shu C. - W.,
!         "Strong stability preserving Runge-Kutta and multistep
!          time discretizations",
!         World Scientific Publishing, 2011
!
!===============================================================================
!
  subroutine evolve_ssprk3_m()

    use blocks     , only : block_data, list_data, get_dblocks
    use boundaries , only : boundary_fluxes
    use coordinates, only : nc => ncells, nb, ne
    use equations  , only : errors, nf, ibp, cmax
#ifdef MPI
    use mpitools   , only : reduce_maximum
#endif /* MPI */
    use sources    , only : update_sources

    implicit none

    type(block_data), pointer :: pdata

    logical     , save :: first = .true.
    integer     , save :: n, n1, n2, n3, n4
    real(kind=8), save :: r, f1, f2, g1, g2

    logical      :: test
    integer      :: i, l, nrej
    real(kind=8) :: tm, dtm, ds, umax, emax
    real(kind=8) :: fc, fcmn, fcmx

    real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr

    real(kind=8), parameter :: k1 = -5.8d-01 / 3.0d+00, k2 = 7.0d-02,          &
                               k3 = -1.0d-01 / 3.0d+00
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
    call start_timer(imu)
#endif /* PROFILE */

    if (first) then

      n  = int(sqrt(1.0d+00 * stages))
      n1 = (n - 1) * (n - 2) / 2
      n2 = n1 + 1
      n3 = n * (n + 1) / 2 - 1
      n4 = n * (n + 1) / 2 + 1

      r  = (1.0d+00 * (2 * n - 1))
      f1 = (1.0d+00 *  n     ) / r
      f2 = (1.0d+00 * (n - 1)) / r
      r  = 1.0d+00 * (stages - n)

      g1 = 1.0d+00 / (    stages - n) - 1.0d+00 / stages
      g2 = 1.0d+00 / (2 * stages - n) - 1.0d+00 / stages

      first = .false.

    end if

    fc   = fac
    fcmn = facmin
    fcmx = facmax

    l = get_dblocks()
#if NDIMS == 3
    allocate(lerr(l,nf,nc,nc,nc))
#else /* NDIMS == 3 */
    allocate(lerr(l,nf,nc,nc, 1))
#endif /* NDIMS == 3 */

    test = .true.
    nrej = 0

    do while(test)

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

      ds = dt / r

!= 1st step: U(1) = U(n)
!
      pdata => list_data
      do while (associated(pdata))

        pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)

        pdata%u => pdata%uu(:,:,:,:,2)

        pdata => pdata%next
      end do

!= 2ns step: U(1) = U(1) + dt/r F[U(1)], for i = 1, ..., (n-1)*(n-2)/2
!
      do i = 1, n1

        tm  = time + i * ds
        dtm = ds

        pdata => list_data
        do while (associated(pdata))

          call update_increment(pdata)
          call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

          pdata => pdata%next
        end do

        call boundary_fluxes()

        l = 1
        pdata => list_data
        do while (associated(pdata))

          pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)

#if NDIMS == 3
          lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne,nb:ne)
#else /* NDIMS == 3 */
          lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne,  :  )
#endif /* NDIMS == 3 */
          l = l + 1

          pdata => pdata%next
        end do

        call update_variables(tm, dtm)

      end do ! n = 1, n1

!= 3rd step: U(2) = U(1)
!
! iterate over all data blocks
!
      pdata => list_data
      do while (associated(pdata))

        pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,2)

        pdata => pdata%next

      end do

!= 4th step: U(1) = U(1) + dt/r F[U(1)], for i = (n-1)*(n-2)/2+1, ..., n*(n+1)/2-1
!
      do i = n2, n3

        tm  = time + i * ds
        dtm = ds

        pdata => list_data
        do while (associated(pdata))

          call update_increment(pdata)
          call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

          pdata => pdata%next
        end do

        call boundary_fluxes()

        l = 1
        pdata => list_data
        do while (associated(pdata))

          pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)

#if NDIMS == 3
          lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne,nb:ne)
#else /* NDIMS == 3 */
          lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne,  :  )
#endif /* NDIMS == 3 */
          l = l + 1

          pdata => pdata%next
        end do

        call update_variables(tm, dtm)

      end do ! i = n2, n3

!= 5th step: U(1) = n * U(2) + (n - 1) * ( U(1) + dt/r F[U(1)] ) / (2 * n - 1)
!
      tm  = time + dt
      dtm = ds

      pdata => list_data
      do while (associated(pdata))

        call update_increment(pdata)
        call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

        pdata => pdata%next
      end do

      call boundary_fluxes()

      l = 1
      pdata => list_data
      do while (associated(pdata))

        pdata%uu(:,:,:,:,2) = f1 *  pdata%uu(:,:,:,:,3)                        &
                            + f2 * (pdata%uu(:,:,:,:,2)                        &
                            + ds *  pdata%du(:,:,:,:))

#if NDIMS == 3
        lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne,nb:ne)
#else /* NDIMS == 3 */
        lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne,  :  )
#endif /* NDIMS == 3 */
        l = l + 1

        pdata => pdata%next
      end do

      call update_variables(tm, dtm)

!= 6th step: U(1) = U(1) + dt/r F[U(1)], for i = n*(n+1)/2, ..., m
!
      do i = n4, stages

        tm  = time + (i - n) * ds
        dtm = ds

        pdata => list_data
        do while (associated(pdata))

          call update_increment(pdata)
          call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

          pdata => pdata%next
        end do

        call boundary_fluxes()

        l = 1
        pdata => list_data
        do while (associated(pdata))

          pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)

#if NDIMS == 3
          lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne,nb:ne)
#else /* NDIMS == 3 */
          lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne,  :  )
#endif /* NDIMS == 3 */
          l = l + 1

          pdata => pdata%next
        end do

        call update_variables(tm, dtm)

      end do ! i = n4, stages

! find umax
!
      umax = 0.0d+00
      pdata => list_data
      do while (associated(pdata))
#if NDIMS == 3
        umax = max(umax, maxval(abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1))),         &
                         maxval(abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,2))))
#else /* NDIMS == 3 */
        umax = max(umax, maxval(abs(pdata%uu(:,nb:ne,nb:ne,  :  ,1))),         &
                         maxval(abs(pdata%uu(:,nb:ne,nb:ne,  :  ,2))))
#endif /* NDIMS == 3 */
        pdata => pdata%next
      end do

! get the maximum error of each variable
!
      do l = 1, nf
        errors(l) = dt * maxval(abs(lerr(:,l,:,:,:)))
      end do

#ifdef MPI
      call reduce_maximum(umax)
      call reduce_maximum(errors)
#endif /* MPI */

! calculate tolerance and time step
!
      emax = maxval(errors) / (atol + rtol * umax)

      if (emax <= 1.0d+00 .or. nrej >= mrej) then
        test = .false.

        errs(3) = errs(2)
        errs(2) = errs(1)
        errs(1) = emax

        errtol  = emax

        dte = dt * min(fcmx, max(fcmn,                                         &
                                 fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))

        fcmx = facmax
      else
        errs(1) = emax

        dte = dt * min(fcmx, max(fcmn,                                         &
                                 fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))
        dt  = dte

        fcmx = fac

        nrej        = nrej        + 1 ! rejection count in the current step
        nrejections = nrejections + 1
      end if

      niterations = niterations + 1

    end do

    if (allocated(lerr)) deallocate(lerr)

!= final step: U(n+1) = U(1)
!
    tm  = time + dt
    dtm = dt / r

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)

      pdata%u => pdata%uu(:,:,:,:,1)

      pdata => pdata%next
    end do

    call update_variables(tm, dtm)

    if (ibp > 0) then
      adecay(:) = exp(aglm(:) * cmax * dt)

      pdata => list_data
      do while (associated(pdata))

        pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)

        pdata => pdata%next
      end do
    end if

#ifdef PROFILE
    call stop_timer(imu)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine evolve_ssprk3_m
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK4_10:
! ---------------------------
!
!   Subroutine advances the solution by one time step using the 4rd order
!   10-stage Strong Stability Preserving Runge-Kutta time integration method.
!
!   References:
!
!     [1] Gottlieb, S., Ketcheson, D., Shu C. - W.,
!         "Strong stability preserving Runge-Kutta and multistep
!          time discretizations",
!         World Scientific Publishing, 2011
!
!===============================================================================
!
  subroutine evolve_ssprk4_10()

    use blocks    , only : block_data, list_data
    use boundaries, only : boundary_fluxes
    use equations , only : ibp, cmax
    use sources   , only : update_sources

    implicit none

    type(block_data), pointer :: pdata

    integer      :: n
    real(kind=8) :: tm, dtm

    real(kind=8), dimension(9) :: ds

    real(kind=8), dimension(9), parameter :: c =                               &
                   (/ 1.0d+00, 2.0d+00, 3.0d+00, 4.0d+00, 2.0d+00,             &
                               3.0d+00, 4.0d+00, 5.0d+00, 6.0d+00 /) / 6.0d+00
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
    call start_timer(imu)
#endif /* PROFILE */

    ds(:) = c(:) * dt

!= 1st step: U(2) = U(1)
!
    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)

      pdata => pdata%next
    end do

!= 2nd step: U(1) = [1 + dt/6 L] U(1), for i = 1, ..., 5
!
    do n = 1, 5

      tm  = time + ds(n)
      dtm = ds(1)

      pdata => list_data
      do while (associated(pdata))

        call update_increment(pdata)
        call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

        pdata => pdata%next
      end do

      call boundary_fluxes()

      pdata => list_data
      do while (associated(pdata))

        pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:)

        pdata => pdata%next
      end do

      call update_variables(tm, dtm)

    end do ! n = 1, 5

!= 3rd step: U(2) = U(2)/25 + 9/25 U(1), U(1) = 15 U(2) - 5 U(1)
!
    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,2) = (pdata%uu(:,:,:,:,2)                               &
                                      + 9.0d+00 * pdata%uu(:,:,:,:,1)) / 2.5d+01
      pdata%uu(:,:,:,:,1) = 1.5d+01 * pdata%uu(:,:,:,:,2)                      &
                                      - 5.0d+00 * pdata%uu(:,:,:,:,1)

      pdata => pdata%next
    end do

!= 4th step: U(1) = [1 + dt/6 L] U(1), for i = 6, ..., 9
!
! integrate the intermediate steps
!
    do n = 6, 9

      tm  = time + ds(n)
      dtm = ds(1)

      pdata => list_data
      do while (associated(pdata))

        call update_increment(pdata)
        call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

        pdata => pdata%next
      end do

      call boundary_fluxes()

      pdata => list_data
      do while (associated(pdata))

        pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:)

        pdata => pdata%next
      end do

      call update_variables(tm, dtm)

    end do ! n = 6, 9

!= the final step: U(n+1) = U(2) + 3/5 U(1) + 1/10 dt F[U(1)]
!
    tm  = time + dt
    dtm = dt / 1.0d+01

    pdata => list_data
    do while (associated(pdata))

      call update_increment(pdata)
      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))

      pdata => pdata%next
    end do

    call boundary_fluxes()

    pdata => list_data
    do while (associated(pdata))

      pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)                                &
                          + 6.0d-01 * pdata%uu(:,:,:,:,1)                      &
                          +     dtm * pdata%du(:,:,:,:)

      pdata => pdata%next
    end do

    if (ibp > 0) then
      adecay(:) = exp(aglm(:) * cmax * dt)

      pdata => list_data
      do while (associated(pdata))

        pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)

        pdata => pdata%next
      end do
    end if

    call update_variables(tm, dtm)

#ifdef PROFILE
    call stop_timer(imu)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine evolve_ssprk4_10
!
!===============================================================================
!
! subroutine UPDATE_INCREMENT:
! ---------------------------
!
!   Subroutine calculates the conservative variable increment from
!   directional fluxes.
!
!   Arguments:
!
!     pdata - the point to data block storing the directional fluxes;
!
!===============================================================================
!
  subroutine update_increment(pdata)

! include external variables
!
    use blocks     , only : block_data
    use coordinates, only : nbl, ne, neu, nn => bcells
    use coordinates, only : adx, ady, adxi, adyi
#if NDIMS == 3
    use coordinates, only : adz, adzi
#endif /* NDIMS == 3 */
    use equations  , only : nf, ns
    use equations  , only : idn, isl, isu
    use schemes    , only : update_flux

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    type(block_data), pointer, intent(inout) :: pdata

! local variables
!
    integer      :: i  , j  , k = 1, p
    integer      :: im1, ip1
    integer      :: jm1, jp1
#if NDIMS == 3
    integer      :: km1, kp1
#endif /* NDIMS == 3 */
    real(kind=8) :: df

    real(kind=8), dimension(NDIMS) :: dh, dhi
#if NDIMS == 3
    real(kind=8), dimension(nf,nn,nn,nn,3) :: f
#else /* NDIMS == 3 */
    real(kind=8), dimension(nf,nn,nn, 1,2) :: f
#endif /* NDIMS == 3 */
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for the increment update
!
    call start_timer(iui)
#endif /* PROFILE */

! prepare the coordinate intervals
!
    dh(1)  = adx(pdata%meta%level)
    dh(2)  = ady(pdata%meta%level)
#if NDIMS == 3
    dh(3)  = adz(pdata%meta%level)
#endif /* NDIMS == 3 */
    dhi(1) = adxi(pdata%meta%level)
    dhi(2) = adyi(pdata%meta%level)
#if NDIMS == 3
    dhi(3) = adzi(pdata%meta%level)
#endif /* NDIMS == 3 */

! calculate the interface fluxes
!
    call update_flux(dh(1:NDIMS), pdata%q(:,:,:,:), f(:,:,:,:,1:NDIMS))

! store block interface fluxes
!
    pdata%fx(:,1,:,:) = f(:,nbl,:,:,1)
    pdata%fx(:,2,:,:) = f(:,ne ,:,:,1)
    pdata%fy(:,:,1,:) = f(:,:,nbl,:,2)
    pdata%fy(:,:,2,:) = f(:,:,ne ,:,2)
#if NDIMS == 3
    pdata%fz(:,:,:,1) = f(:,:,:,nbl,3)
    pdata%fz(:,:,:,2) = f(:,:,:,ne ,3)
#endif /* NDIMS == 3 */

! calculate the variable update from the directional fluxes
!
#if NDIMS == 3
    do k = nbl, neu
      km1 = k - 1
#endif /* NDIMS == 3 */
      do j = nbl, neu
        jm1 = j - 1
        do i = nbl, neu
          im1 = i - 1

#if NDIMS == 3
          pdata%du(1:nf,i,j,k) = - dhi(1) * (f(:,i,j,k,1) - f(:,im1,j,k,1))    &
                                 - dhi(2) * (f(:,i,j,k,2) - f(:,i,jm1,k,2))    &
                                 - dhi(3) * (f(:,i,j,k,3) - f(:,i,j,km1,3))
#else /* NDIMS == 3 */
          pdata%du(1:nf,i,j,k) = - dhi(1) * (f(:,i,j,k,1) - f(:,im1,j,k,1))    &
                                 - dhi(2) * (f(:,i,j,k,2) - f(:,i,jm1,k,2))
#endif /* NDIMS == 3 */
        end do ! i = nbl, neu
      end do ! j = nbl, neu
#if NDIMS == 3
    end do ! k = nbl, neu
#endif /* NDIMS == 3 */

! update passive scalars
!
    if (ns > 0) then

! reset passive scalar increments
!
      pdata%du(isl:isu,:,:,:) = 0.0d+00

#if NDIMS == 3
      do k = nbl - 1, neu
        kp1 = k + 1
#endif /* NDIMS == 3 */
        do j = nbl - 1, neu
          jp1 = j + 1
          do i = nbl - 1, neu
            ip1 = i + 1

! X-face
!
            if (f(idn,i,j,k,1) >= 0.0d+00) then
              do p = isl, isu
                df = dhi(1) * f(idn,i,j,k,1) * pdata%q(p,i  ,j,k)
                pdata%du(p,i  ,j,k) = pdata%du(p,i  ,j,k) - df
                pdata%du(p,ip1,j,k) = pdata%du(p,ip1,j,k) + df
              end do
            else
              do p = isl, isu
                df = dhi(1) * f(idn,i,j,k,1) * pdata%q(p,ip1,j,k)
                pdata%du(p,i  ,j,k) = pdata%du(p,i  ,j,k) - df
                pdata%du(p,ip1,j,k) = pdata%du(p,ip1,j,k) + df
              end do
            end if

! Y-face
!
            if (f(idn,i,j,k,2) >= 0.0d+00) then
              do p = isl, isu
                df = dhi(2) * f(idn,i,j,k,2) * pdata%q(p,i,j  ,k)
                pdata%du(p,i,j  ,k) = pdata%du(p,i,j  ,k) - df
                pdata%du(p,i,jp1,k) = pdata%du(p,i,jp1,k) + df
              end do
            else
              do p = isl, isu
                df = dhi(2) * f(idn,i,j,k,2) * pdata%q(p,i,jp1,k)
                pdata%du(p,i,j  ,k) = pdata%du(p,i,j  ,k) - df
                pdata%du(p,i,jp1,k) = pdata%du(p,i,jp1,k) + df
              end do
            end if

#if NDIMS == 3
! Z-face
!
            if (f(idn,i,j,k,3) >= 0.0d+00) then
              do p = isl, isu
                df = dhi(3) * f(idn,i,j,k,3) * pdata%q(p,i,j,k  )
                pdata%du(p,i,j,k  ) = pdata%du(p,i,j,k  ) - df
                pdata%du(p,i,j,kp1) = pdata%du(p,i,j,kp1) + df
              end do
            else
              do p = isl, isu
                df = dhi(3) * f(idn,i,j,k,3) * pdata%q(p,i,j,kp1)
                pdata%du(p,i,j,k  ) = pdata%du(p,i,j,k  ) - df
                pdata%du(p,i,j,kp1) = pdata%du(p,i,j,kp1) + df
              end do
            end if
#endif /* NDIMS == 3 */

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

#ifdef PROFILE
! stop accounting time for the increment update
!
    call stop_timer(iui)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine update_increment
!
!===============================================================================
!
! subroutine UPDATE_VARIABLES:
! ---------------------------
!
!   Subroutine iterates over all data blocks and converts the conservative
!   variables to their primitive representation.
!
!   Arguments:
!
!     tm  - time at the moment of update;
!     dtm - time step since the last update;
!
!===============================================================================
!
  subroutine update_variables(tm, dtm)

! include external procedures
!
    use blocks        , only : set_neighbors_update
    use boundaries    , only : boundary_variables
    use equations     , only : update_primitive_variables
    use equations     , only : fix_unphysical_cells, correct_unphysical_states
    use shapes        , only : update_shapes

! include external variables
!
    use blocks        , only : block_meta, list_meta
    use blocks        , only : block_data, list_data

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    real(kind=8), intent(in) :: tm, dtm

! local pointers
!
    type(block_meta), pointer :: pmeta
    type(block_data), pointer :: pdata
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for variable update
!
    call start_timer(imv)
#endif /* PROFILE */

! update primitive variables in the changed blocks
!
    pdata => list_data
    do while (associated(pdata))
      pmeta => pdata%meta

      if (pmeta%update) call update_primitive_variables(pdata%u, pdata%q)

      pdata => pdata%next
    end do

! update boundaries
!
    call boundary_variables(tm, dtm)

! apply shapes in blocks which need it
!
    pdata => list_data
    do while (associated(pdata))
      pmeta => pdata%meta

      if (pmeta%update) call update_shapes(pdata, tm, dtm)

      pdata => pdata%next
    end do

! correct unphysical states if detected
!
    if (fix_unphysical_cells) then

! if an unphysical cell appeared in a block while updating its primitive
! variables it could be propagated to its neighbors through boundary update;
! mark all neighbors of such a block to be verified and corrected for
! unphysical cells too
!
      pmeta => list_meta
      do while (associated(pmeta))

        if (pmeta%update) call set_neighbors_update(pmeta)

        pmeta => pmeta%next
      end do

! verify and correct, if necessary, unphysical cells in recently updated blocks
!
      pdata => list_data
      do while (associated(pdata))
        pmeta => pdata%meta

        if (pmeta%update)                                                      &
              call correct_unphysical_states(step, pmeta%id, pdata%q, pdata%u)

        pdata => pdata%next
      end do
    end if

#ifdef PROFILE
! stop accounting time for variable update
!
    call stop_timer(imv)
#endif /* PROFILE */

!-------------------------------------------------------------------------------
!
  end subroutine update_variables
!
!===============================================================================
!
! subroutine UPDATE_ERRORS_L2:
! ---------------------------
!
!   Subroutine updates the errors with respect to the absolute and relative
!   tolerances using the L2 norm. The errors are calculated per variable.
!
!   Arguments:
!
!     n - the index of the higher order solution
!     m - the index of the lower order solution
!
!===============================================================================
!
  subroutine update_errors_l2(n, m)

    use blocks     , only : block_data, list_data, get_nleafs
    use coordinates, only : ncells, nb, ne
    use equations  , only : nf, errors
#ifdef MPI
    use mpitools   , only : reduce_sum
#endif /* MPI */

    implicit none

    integer, intent(in) :: n, m

    integer :: l
    real(kind=8) :: fnorm

    type(block_data), pointer :: pdata

!-------------------------------------------------------------------------------
!
    errors(:) = 0.0d+00
    fnorm     = 1.0d+00 / (get_nleafs() * ncells**NDIMS)

    pdata => list_data
    do while (associated(pdata))
      do l = 1, nf
#if NDIMS == 3
        errors(l) = errors(l) +                                                &
                      sum((abs(pdata%uu(l,nb:ne,nb:ne,nb:ne,n)                 &
                             - pdata%uu(l,nb:ne,nb:ne,nb:ne,m)) /              &
                      (atol + rtol *                                           &
                           max(abs(pdata%uu(l,nb:ne,nb:ne,nb:ne,n)),           &
                               abs(pdata%uu(l,nb:ne,nb:ne,nb:ne,m)))))**2)
#else /* NDIMS == 3 */
        errors(l) = errors(l) +                                                &
                      sum((abs(pdata%uu(l,nb:ne,nb:ne,  :  ,n)                 &
                             - pdata%uu(l,nb:ne,nb:ne,  :  ,m)) /              &
                      (atol + rtol *                                           &
                           max(abs(pdata%uu(l,nb:ne,nb:ne,  :  ,n)),           &
                               abs(pdata%uu(l,nb:ne,nb:ne,  :  ,m)))))**2)
#endif /* NDIMS == 3 */
      end do

      pdata => pdata%next
    end do

#ifdef MPI
    call reduce_sum(errors)
#endif /* MPI */

    errors(:) = sqrt(fnorm * errors(:))

!-------------------------------------------------------------------------------
!
  end subroutine update_errors_l2
!
!===============================================================================
!
! subroutine UPDATE_ERRORS_MAX:
! ----------------------------
!
!   Subroutine updates the errors with respect to the absolute and relative
!   tolerances using the maximum norm. The errors are calculated per variable.
!
!   Arguments:
!
!     n - the index of the higher order solution
!     m - the index of the lower order solution
!
!===============================================================================
!
  subroutine update_errors_max(n, m)

    use blocks     , only : block_data, list_data
    use coordinates, only : nb, ne
    use equations  , only : nf, errors
#ifdef MPI
    use mpitools   , only : reduce_maximum
#endif /* MPI */

    implicit none

    integer, intent(in) :: n, m

    integer :: l

    type(block_data), pointer :: pdata

!-------------------------------------------------------------------------------
!
    errors(:) = 0.0d+00

    pdata => list_data
    do while (associated(pdata))
      do l = 1, nf
#if NDIMS == 3
        errors(l) = max(errors(l),                                             &
                        maxval(abs(pdata%uu(l,nb:ne,nb:ne,nb:ne,n)             &
                                 - pdata%uu(l,nb:ne,nb:ne,nb:ne,m)) /          &
                      (atol + rtol *                                           &
                           max(abs(pdata%uu(l,nb:ne,nb:ne,nb:ne,n)),           &
                               abs(pdata%uu(l,nb:ne,nb:ne,nb:ne,m))))))
#else /* NDIMS == 3 */
        errors(l) = max(errors(l),                                             &
                        maxval(abs(pdata%uu(l,nb:ne,nb:ne,  :  ,n)             &
                                 - pdata%uu(l,nb:ne,nb:ne,  :  ,m)) /          &
                      (atol + rtol *                                           &
                           max(abs(pdata%uu(l,nb:ne,nb:ne,  :  ,n)),           &
                               abs(pdata%uu(l,nb:ne,nb:ne,  :  ,m))))))
#endif /* NDIMS == 3 */
      end do

      pdata => pdata%next
    end do

#ifdef MPI
    call reduce_maximum(errors)
#endif /* MPI */

!-------------------------------------------------------------------------------
!
  end subroutine update_errors_max
#ifdef DEBUG
!
!===============================================================================
!
! subroutine CHECK_VARIABLES:
! --------------------------
!
!   Subroutine iterates over all data blocks and converts the conservative
!   variables to their primitive representation.
!
!
!===============================================================================
!
  subroutine check_variables()

! include external procedures
!
    use coordinates    , only : nn => bcells
    use equations      , only : nv, pvars, cvars
    use ieee_arithmetic, only : ieee_is_nan

! include external variables
!
    use blocks        , only : block_meta, list_meta
    use blocks        , only : block_data, list_data

! local variables are not implicit by default
!
    implicit none

! local variables
!
    integer :: i, j, k = 1, p

! local pointers
!
    type(block_meta), pointer :: pmeta
    type(block_data), pointer :: pdata
!
!-------------------------------------------------------------------------------
!
! associate the pointer with the first block on the data block list
!
    pdata => list_data

! iterate over all data blocks
!
    do while (associated(pdata))

! associate pmeta with the corresponding meta block
!
      pmeta => pdata%meta

! check if there are NaNs in primitive variables
!
#if NDIMS == 3
      do k = 1, nn
#endif /* NDIMS == 3 */
        do j = 1, nn
          do i = 1, nn
            do p = 1, nv
              if (ieee_is_nan(pdata%u(p,i,j,k))) then
                print *, 'U NaN:', cvars(p), pdata%meta%id, i, j, k
              end if
              if (ieee_is_nan(pdata%q(p,i,j,k))) then
                print *, 'Q NaN:', pvars(p), pdata%meta%id, i, j, k
              end if
            end do ! p = 1, nv
          end do ! i = 1, nn
        end do ! j = 1, nn
#if NDIMS == 3
      end do ! k = 1, nn
#endif /* NDIMS == 3 */

! assign pointer to the next block
!
      pdata => pdata%next

    end do

!-------------------------------------------------------------------------------
!
  end subroutine check_variables
#endif /* DEBUG */

!===============================================================================
!
end module evolution