!!******************************************************************************
!!
!!  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

  implicit none

! 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.
  logical           , save :: fsal          = .false.
  logical           , save :: fsal_update   = .true.
  character(len=255), save :: name_int  = ""
  character(len=255), save :: name_norm = ""
  integer           , save :: order     = 2
  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 :: dth      = 1.0d+00
  real(kind=8)      , save :: dte      = 0.0d+00
  real(kind=8)      , save :: dtp      = huge(dt) ! dt for precise snapshots

! 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

! coefficients of Runge-Kutta methods
!
  real(kind=8), dimension(:)  , allocatable, save :: c, bt, dl, bh
  real(kind=8), dimension(:,:), allocatable, save :: gm

! 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 :: initialize_time_step, new_time_step, advance

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

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

    use coordinates, only : toplev, adx, ady, adz
    use forcing    , only : initialize_forcing
    use gravity    , only : initialize_gravity
    use helpers    , only : print_message
    use mpitools   , only : check_status
    use parameters , only : get_parameter
    use schemes    , only : initialize_schemes
    use shapes     , only : initialize_shapes
    use sources    , only : initialize_sources

    implicit none

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

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

    character(len=*), parameter :: loc = 'EVOLUTION::initialize_evolution()'

!-------------------------------------------------------------------------------
!
    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
      order     = 1
      registers = 2
      name_int  = "1st order Euler"

    case ("rk2", "RK2")

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

    case ("rk3", "RK3")

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

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

      evolve    => evolve_ssprk34
      order     = 3
      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
      order     = 3
      registers = 3
      cfl       = 2.65062919143939d+00 * cfl
      name_int  = "3rd order SSPRK(5,3)"

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

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

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

      error_control = .true.
      evolve    => evolve_ssprk2_m
      order     = 2
      registers = 3
      stages    = max(2, min(9, stages))
      cfl       = (stages - 1) * cfl
      write(name_int, "('2ⁿᵈ-order ',i0,'-step embedded " //                   &
                      "SSPRK(',i0,',2) method')") stages, stages

      allocate(c(stages), dl(stages), bt(stages), bh(stages))

      do n = 1, stages
        c(n) = n - 1.0d+00
      end do
      c (:)      = c(:)    / (stages - 1)
      bt(:)      = 1.0d+00 / (stages - 1)
      bt(stages) = 1.0d+00 /  stages
      bh(:)      = 1.0d+00 /  stages
      bh(1)      = (stages + 1.0d+00) / stages**2
      bh(stages) = (stages - 1.0d+00) / stages**2
      dl(1)      = (stages - 1.0d+00) / stages
      dl(2)      = 1.0d+00 - dl(1)

      betas(:) = [ 5.8d-01, -2.1d-01, 1.0d-01 ]

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

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

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

      error_control = .true.
      evolve    => evolve_ssp324
      order     = 3
      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 ("SSPRK3(2)m", "SSP3(2)m", "ssprk3(2)m", "ssp3(2)m", "SSPRK(m,3)", "ssprk(m,3)")

      error_control = .true.
      evolve   => evolve_ssprk32m
      order     = 3
      registers = 4
      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, "('3ʳᵈ-order order ',i0,'-step embedded SSPRK3(2)n² method')") stages

      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

      if (stages == 4 .and. verbose) then
        write(*,*)
        write(*,*) "ADVICE:"
        write(*,*) "The method 'SSPRK3(2)m' is primarily designed for n²,"  // &
                   " where n > 2. Since you intend to use 4 stages with"    // &
                   " this method, it is highly advised to use method"       // &
                   " 'SSPRK3(2)4', since it slightly better optimized and"  // &
                   " significantly used less memory."
      end if

    case ("RK3(2)5", "rk3(2)5")

!   References:
!
!     [1] Ranocha, H., Dalcin, L., Parsani, M., Ketcheson, D. I.
!         "Coefficients of Optimized Runge-Kutta Methods with Automatic
!          Step Size Control for Compressible Computational Fluid Dynamics",
!         2021, DOI:10.5281/zenodo.4671927,
!         https://github.com/ranocha/Optimized-RK-CFD
!
      error_control = .true.
      evolve    => evolve_3sstarplus
      registers = 4
      order     = 3
      stages    = 5
      cfl       = 2.5d+00 * cfl
      name_int  = "Optimized 3ʳᵈ-order 5-step embedded RK3(2)5"

      allocate(c(stages+1), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
      c (:)   = [ 0.000000000000000000000000000000000000d+00,                  &
                  2.300285062878154351930669430512780706d-01,                  &
                  4.050049049262914975700372321130661410d-01,                  &
                  8.947823877926760224705450466361360720d-01,                  &
                  7.235108137218888081489570284485201518d-01,                  &
                  1.000000000000000000000000000000000000d+00 ]
      dl(:)   = [ 1.000000000000000000000000000000000000d+00,                  &
                  3.407687209321455242558804921815861422d-01,                  &
                  3.414399280584625023244387687873774697d-01,                  &
                  7.229302732875589702087936723400941329d-01,                  &
                  0.000000000000000000000000000000000000d+00 ]
      bt(:)   = [ 2.300285062878154351930669430512780706d-01,                  &
                  3.021457892454169700189445968126242994d-01,                  &
                  8.025601039472704213300183888573974531d-01,                  &
                  4.362158997637629844305216319994356355d-01,                  &
                  1.129268494470295369172265188216779157d-01 ]
      bh(:)   = [ 1.046363371354093758897668305991705199d-01,                  &
                  9.520431574956758809511173383346476348d-02,                  &
                  4.482446645568668405072421350300379357d-01,                  &
                  2.449030295461310135957132640369862245d-01,                  &
                  1.070116530120251819121660365003405564d-01,                  &
                  0.000000000000000000000000000000000000d+00 ]
      gm(1,:) = [ 0.000000000000000000000000000000000000d+00,                  &
                  2.587669070352079020144955303389306026d-01,                  &
                 -1.324366873994502973977035353758550057d-01,                  &
                  5.055601231460399101814291350373559483d-02,                  &
                  5.670552807902877312521811889846000976d-01 ]
      gm(2,:) = [ 1.000000000000000000000000000000000000d+00,                  &
                  5.528418745102160639901976698795928733d-01,                  &
                  6.731844400389673824374042790213570079d-01,                  &
                  2.803103804507635075215805236096803381d-01,                  &
                  5.521508873507393276457754945308880998d-01 ]
      gm(3,:) = [ 0.000000000000000000000000000000000000d+00,                  &
                  0.000000000000000000000000000000000000d+00,                  &
                  0.000000000000000000000000000000000000d+00,                  &
                  2.752585813446636957256614568573008811d-01,                  &
                 -8.950548709279785077579454232514633376d-01 ]

      betas(:) = [ 6.4d-01, -3.1d-01, 4.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(2)5F", "rk3(2)5f")

!   References:
!
!     [1] Ranocha, H., Dalcin, L., Parsani, M., Ketcheson, D. I.
!         "Coefficients of Optimized Runge-Kutta Methods with Automatic
!          Step Size Control for Compressible Computational Fluid Dynamics",
!         2021, DOI:10.5281/zenodo.4671927,
!         https://github.com/ranocha/Optimized-RK-CFD
!
      error_control = .true.
      fsal          = .true.
      evolve    => evolve_3sstarplus
      registers = 4
      order     = 3
      stages    = 5
      cfl       = 2.5d+00 * cfl
      name_int  = "Optimized 3ʳᵈ-order 5-step embedded RK3(2)5 FSAL"

      allocate(c(stages+1), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
      c (:)   = [ 0.000000000000000000000000000000000000d+00,                  &
                  2.300298624518076223899418286314123354d-01,                  &
                  4.050046072094990912268498160116125481d-01,                  &
                  8.947822893693433545220710894560512805d-01,                  &
                  7.235136928826589010272834603680114769d-01,                  &
                  1.000000000000000000000000000000000000d+00 ]
      dl(:)   = [ 1.000000000000000000000000000000000000d+00,                  &
                  3.407655879334525365094815965895763636d-01,                  &
                  3.414382655003386206551709871126405331d-01,                  &
                  7.229275366787987419692007421895451953d-01,                  &
                  0.000000000000000000000000000000000000d+00 ]
      bt(:)   = [ 2.300298624518076223899418286314123354d-01,                  &
                  3.021434166948288809034402119555380003d-01,                  &
                  8.025606185416310937583009085873554681d-01,                  &
                  4.362158943603440930655148245148766471d-01,                  &
                  1.129272530455059129782111662594436580d-01 ]
      bh(:)   = [ 9.484166705035703392326247283838082847d-02,                  &
                  1.726371339430353766966762629176676070d-01,                  &
                  3.998243189084371024483169698618455770d-01,                  &
                  1.718016807580178450618829007973835152d-01,                  &
                  5.881914422155740300718268359027168467d-02,                  &
                  1.020760551185952388626787099944507877d-01 ]
      gm(1,:) = [ 0.000000000000000000000000000000000000d+00,                  &
                  2.587771979725733308135192812685323706d-01,                  &
                 -1.324380360140723382965420909764953437d-01,                  &
                  5.056033948190826045833606441415585735d-02,                  &
                  5.670532000739313812633197158607642990d-01 ]
      gm(2,:) = [ 1.000000000000000000000000000000000000d+00,                  &
                  5.528354909301389892439698870483746541d-01,                  &
                  6.731871608203061824849561782794643600d-01,                  &
                  2.803103963297672407841316576323901761d-01,                  &
                  5.521525447020610386070346724931300367d-01 ]
      gm(3,:) = [ 0.000000000000000000000000000000000000d+00,                  &
                  0.000000000000000000000000000000000000d+00,                  &
                  0.000000000000000000000000000000000000d+00,                  &
                  2.752563273304676380891217287572780582d-01,                  &
                 -8.950526174674033822276061734289327568d-01 ]

      betas(:) = [ 7.0d-01, -2.3d-01, 0.0d+00 ]

      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 ("RK4(3)9", "rk4(3)9")

!   References:
!
!     [1] Ranocha, H., Dalcin, L., Parsani, M., Ketcheson, D. I.
!         "Coefficients of Optimized Runge-Kutta Methods with Automatic
!          Step Size Control for Compressible Computational Fluid Dynamics",
!         2021, DOI:10.5281/zenodo.4671927,
!         https://github.com/ranocha/Optimized-RK-CFD
!
      error_control = .true.
      evolve    => evolve_3sstarplus
      registers = 4
      order     = 4
      stages    = 9
      cfl       = 3.5d+00 * cfl
      name_int  = "Optimized 4ᵗʰ-order 9-step embedded RK4(3)9"

      allocate(c(stages+1), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
      c (:)   = [  0.000000000000000000000000000000000000d+00,                 &
                   2.836343531977826022543660465926414772d-01,                 &
                   5.484073767552486705240014599676811834d-01,                 &
                   3.687229456675706936558667052479014150d-01,                 &
                  -6.806119916032093175251948474173648331d-01,                 &
                   3.518526451892056368706593492732753284d-01,                 &
                   1.665941920204672094647868254892387293d+00,                 &
                   9.715276989307335935187466054546761665d-01,                 &
                   9.051569554420045339601721625247585643d-01,                 &
                   1.000000000000000000000000000000000000d+00 ]
      dl(:)   = [  1.000000000000000000000000000000000000d+00,                 &
                   1.262923854387806460989545005598562667d+00,                 &
                   7.574967177560872438940839460448329992d-01,                 &
                   5.163591158111222863455531895152351544d-01,                 &
                  -2.746333792042827389548936599648122146d-02,                 &
                  -4.382674653941770848797864513655752318d-01,                 &
                   1.273587103668392811985704533534301656d+00,                 &
                  -6.294740045442794829622796613103492913d-01,                 &
                   0.000000000000000000000000000000000000d+00 ]
      bt(:)   = [  2.836343531977826022543660465926414772d-01,                 &
                   9.736497978646965372894268287659773644d-01,                 &
                   3.382358566377620380505126936670933370d-01,                 &
                  -3.584937820217850715182820651063453804d-01,                 &
                  -4.113955814725134294322006403954822487d-03,                 &
                   1.427968962196019024010757034274849198d+00,                 &
                   1.808467712038743032991177525728915926d-02,                 &
                   1.605771316794521018947553625079465692d-01,                 &
                   2.952226811394310028003810072027839487d-01 ]
      bh(:)   = [  4.550655927970944948340364817140593012d-02,                 &
                   1.175968310492638562142460384341959193d-01,                 &
                   3.658257330515213200375475084421083608d-02,                 &
                  -5.311555834355629559010061596928357525d-03,                 &
                   5.178250012713127329531367677410650996d-03,                 &
                   4.954639022118682638697706200022961443d-01,                 &
                  -5.999303132737865921441409466809521699d-03,                 &
                   9.405093434568315929035250835218733824d-02,                 &
                   2.169318087627035072893925375820310602d-01,                 &
                   0.000000000000000000000000000000000000d+00 ]
      gm(1,:) = [  0.000000000000000000000000000000000000d+00,                 &
                  -4.655641301259180308677051498071354582d+00,                 &
                  -7.720264924836063859141482018013692338d-01,                 &
                  -4.024423213419724605695005429153112050d+00,                 &
                  -2.129685246739018613087466942802498152d-02,                 &
                  -2.435022519234470128602335652131234586d+00,                 &
                   1.985627480986167686791439120784668251d-02,                 &
                  -2.810790112885283952929218377438668784d-01,                 &
                   1.689434895835535695524003319503844110d-01 ]
      gm(2,:) = [  1.000000000000000000000000000000000000d+00,                 &
                   2.499262752607825957145627300817258023d+00,                 &
                   5.866820365436136799319929406678132638d-01,                 &
                   1.205141365412670762568835277881144391d+00,                 &
                   3.474793796700868848597960521248007941d-01,                 &
                   1.321346140128723105871355808477092220d+00,                 &
                   3.119636324379370564023292317172847140d-01,                 &
                   4.351419055894087609560896967082486864d-01,                 &
                   2.359698299440788299161958168555704234d-01 ]
      gm(3,:) = [  0.000000000000000000000000000000000000d+00,                 &
                   0.000000000000000000000000000000000000d+00,                 &
                   0.000000000000000000000000000000000000d+00,                 &
                   7.621037111138170045618771082985664430d-01,                 &
                  -1.981182159087218433914909510116664154d-01,                 &
                  -6.228960706317566993192689455719570179d-01,                 &
                  -3.752246993432626328289874575355102038d-01,                 &
                  -3.355436539000946543242869676125143358d-01,                 &
                  -4.560963110717484359015342341157302403d-02 ]

      betas(:) = [ 2.5d-01, -1.2d-01, 0.0d+00 ]

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

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

    case ("RK4(3)9F", "rk4(3)9f")

!   References:
!
!     [1] Ranocha, H., Dalcin, L., Parsani, M., Ketcheson, D. I.
!         "Coefficients of Optimized Runge-Kutta Methods with Automatic
!          Step Size Control for Compressible Computational Fluid Dynamics",
!         2021, DOI:10.5281/zenodo.4671927,
!         https://github.com/ranocha/Optimized-RK-CFD
!
      error_control = .true.
      fsal          = .true.
      evolve    => evolve_3sstarplus
      registers = 4
      order     = 4
      stages    = 9
      cfl       = 3.5d+00 * cfl
      name_int  = "Optimized 4ᵗʰ-order 9-step embedded RK4(3)9 FSAL"

      allocate(c(stages+1), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
      c (:)   = [  0.000000000000000000000000000000000000d+00,                 &
                   2.836343005184365275160654678626695428d-01,                 &
                   5.484076570002894365286665352032296535d-01,                 &
                   3.687228761669438493478872632332010073d-01,                 &
                  -6.806126440140844191258463830024463902d-01,                 &
                   3.518526124230705801739919476290327750d-01,                 &
                   1.665941994879593315477304663913129942d+00,                 &
                   9.715279295934715835299192116436237065d-01,                 &
                   9.051569840159589594903399929316959062d-01,                 &
                   1.000000000000000000000000000000000000d+00 ]
      dl(:)   = [  1.000000000000000000000000000000000000d+00,                 &
                   1.262923876648114432874834923838556100d+00,                 &
                   7.574967189685911558308119415539596711d-01,                 &
                   5.163589453140728104667573195005629833d-01,                 &
                  -2.746327421802609557034437892013640319d-02,                 &
                  -4.382673178127944142238606608356542890d-01,                 &
                   1.273587294602656522645691372699677063d+00,                 &
                  -6.294740283927400326554066998751383342d-01,                 &
                   0.000000000000000000000000000000000000d+00 ]
      bt(:)   = [  2.836343005184365275160654678626695428d-01,                 &
                   9.736500104654741223716056170419660217d-01,                 &
                   3.382359225242515288768487569778320563d-01,                 &
                  -3.584943611106183357043212309791897386d-01,                 &
                  -4.113944068471528211627210454497620358d-03,                 &
                   1.427968894048586363415504654313371031d+00,                 &
                   1.808470948394314017665968411915568633d-02,                 &
                   1.605770645946802213926893453819236685d-01,                 &
                   2.952227015964591648775833803635147962d-01 ]
      bh(:)   = [  2.483675912451591196775756814283216443d-02,                 &
                   1.866327774562103796990092260942180726d-01,                 &
                   5.671080795936984495604436622517631183d-02,                 &
                  -3.447695439149287702616943808570747099d-03,                 &
                   3.602245056516636472203469198006404016d-03,                 &
                   4.545570622145088936800484247980581766d-01,                 &
                  -2.434665289427612407531544765622888855d-04,                 &
                   6.642755361103549971517945063138312147d-02,                 &
                   1.613697079523505006226025497715177578d-01,                 &
                   4.955424859358438183052504342394102722d-02 ]
      gm(1,:) = [  0.000000000000000000000000000000000000d+00,                 &
                  -4.655641447335068552684422206224169103d+00,                 &
                  -7.720265099645871829248487209517314217d-01,                 &
                  -4.024436690519806086742256154738379161d+00,                 &
                  -2.129676284018530966221583708648634733d-02,                 &
                  -2.435022509790109546199372365866450709d+00,                 &
                   1.985627297131987000579523283542615256d-02,                 &
                  -2.810791146791038566946663374735713961d-01,                 &
                   1.689434168754859644351230590422137972d-01 ]
      gm(2,:) = [  1.000000000000000000000000000000000000d+00,                 &
                   2.499262792574495009336242992898153462d+00,                 &
                   5.866820377718875577451517985847920081d-01,                 &
                   1.205146086523094569925592464380295241d+00,                 &
                   3.474793722186732780030762737753849272d-01,                 &
                   1.321346060965113109321230804210670518d+00,                 &
                   3.119636464694193615946633676950358444d-01,                 &
                   4.351419539684379261368971206040518552d-01,                 &
                   2.359698130028753572503744518147537768d-01 ]
      gm(3,:) = [  0.000000000000000000000000000000000000d+00,                 &
                   0.000000000000000000000000000000000000d+00,                 &
                   0.000000000000000000000000000000000000d+00,                 &
                   7.621006678721315291614677352949377871d-01,                 &
                  -1.981182504339400567765766904309673119d-01,                 &
                  -6.228959218699007450469629366684127462d-01,                 &
                  -3.752248380775956442989480369774937099d-01,                 &
                  -3.355438309135169811915662336248989661d-01,                 &
                  -4.560955005031121479972862973705108039d-02 ]

      betas(:) = [ 3.8d-01, -1.8d-01, 1.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(:) / 4.0d+00

    case ("RK5(4)10", "rk5(4)10")

!   References:
!
!     [1] Ranocha, H., Dalcin, L., Parsani, M., Ketcheson, D. I.
!         "Coefficients of Optimized Runge-Kutta Methods with Automatic
!          Step Size Control for Compressible Computational Fluid Dynamics",
!         2021, DOI:10.5281/zenodo.4671927,
!         https://github.com/ranocha/Optimized-RK-CFD
!
      error_control = .true.
      evolve    => evolve_3sstarplus
      registers = 4
      order     = 5
      stages    = 10
      cfl       = 4.2d+00 * cfl
      name_int  = "Optimized 5ᵗʰ-order 10-step embedded RK5(4)10"

      allocate(c(stages+1), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
      c (:)   = [  0.000000000000000000000000000000000000d+00,                 &
                   2.597883575710995826783320802193635406d-01,                 &
                   9.904573115730917688557891428202061598d-02,                 &
                   2.155511882303785204133426661931565216d-01,                 &
                   5.007950078421880417512789524851012021d-01,                 &
                   5.592251914858131230054392022144328176d-01,                 &
                   5.449986973408778242805929551952000165d-01,                 &
                   7.615224662599497796472095353126697300d-01,                 &
                   8.427062083059167761623893618875787414d-01,                 &
                   9.152209807185253394871325258038753352d-01,                 &
                   1.000000000000000000000000000000000000d+00 ]
      dl(:)   = [  1.000000000000000000000000000000000000d+00,                 &
                  -1.331778409133849616712007380176762548d-01,                 &
                   8.260422785246030254485064732649153253d-01,                 &
                   1.513700430513332405798616943654007796d+00,                 &
                  -1.305810063177048110528482211982726539d+00,                 &
                   3.036678789342507704281817524408221954d+00,                 &
                  -1.449458267074592489788800461540171106d+00,                 &
                   3.834313873320957483471400258279635203d+00,                 &
                   4.122293971923324492772059928094971199d+00,                 &
                   0.000000000000000000000000000000000000d+00 ]
      bt(:)   = [  2.597883575710995826783320802193635406d-01,                 &
                   1.777008800169541694837687556103565007d-02,                 &
                   2.481636637328140606807905234325691851d-01,                 &
                   7.941736827560429420202759490815682546d-01,                 &
                   3.885391296871822541486945325814526190d-01,                 &
                   1.455051664264339366757555740296587660d-01,                 &
                   1.587517379462528932413419955691782412d-01,                 &
                   1.650605631567659573994022720500446501d-01,                 &
                   2.118093299943235065178000892467421832d-01,                 &
                   1.559392340339606299335442956580114440d-01 ]
      bh(:)   = [  5.734588484676193812418453938089759359d-02,                 &
                   1.971447518039733870541652912891291496d-02,                 &
                   7.215296605683716720707226840456658773d-02,                 &
                   1.739659489807939956977075317768151880d-01,                 &
                   3.703693600445487815015171515640585668d-01,                 &
                  -1.215599039055065009827765147821222534d-01,                 &
                   1.180372945491121604465067725859678821d-01,                 &
                   4.155688823364870056536983972605056553d-02,                 &
                   1.227886627910379901351569893551486490d-01,                 &
                   1.456284232223684285998448928597043056d-01,                 &
                   0.000000000000000000000000000000000000d+00 ]
      gm(1,:) = [  0.000000000000000000000000000000000000d+00,                 &
                   4.043660078504695837542588769963326988d-01,                 &
                  -8.503427464263185087039788184485627962d-01,                 &
                  -6.950894167072419998080989313353063399d+00,                 &
                   9.238765225328278557805080247596562995d-01,                 &
                  -2.563178039957404359875124580586147888d+00,                 &
                   2.545744869966347362604059848503340890d-01,                 &
                   3.125831733863168874151935287174374515d-01,                 &
                  -7.007114800567584871263283872289072079d-01,                 &
                   4.839620970980726631935174740648996010d-01 ]
      gm(2,:) = [  1.000000000000000000000000000000000000d+00,                 &
                   6.871467069752345566001768382316915820d-01,                 &
                   1.093024760468898686510433898645775908d+00,                 &
                   3.225975382330161123625348062949430509d+00,                 &
                   1.041153700841396427100436517666787823d+00,                 &
                   1.292821488864702752767390075072674807d+00,                 &
                   7.391462769297006312785029455392854586d-01,                 &
                   1.239129257039300081860496157739352186d-01,                 &
                   1.842753479366766790220633908793933781d-01,                 &
                   5.712788942697077644959290025755003720d-02 ]
      gm(3,:) = [  0.000000000000000000000000000000000000d+00,                 &
                   0.000000000000000000000000000000000000d+00,                 &
                   0.000000000000000000000000000000000000d+00,                 &
                  -2.393405159342139386425044844626597490d+00,                 &
                  -1.902854422095986544338294743445530533d+00,                 &
                  -2.820042210583207174321941694153843259d+00,                 &
                  -1.832698464130564949123807896975136336d+00,                 &
                  -2.199094510750697865007677774395365522d-01,                 &
                  -4.082430660384876496971887725512427800d-01,                 &
                  -1.377669791121207993339861855818881150d-01 ]

      betas(:) = [ 4.7d-01, -2.0d-01, 6.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(:) / 5.0d+00

    case ("RK5(4)10F", "rk5(4)10f")

!   References:
!
!     [1] Ranocha, H., Dalcin, L., Parsani, M., Ketcheson, D. I.
!         "Coefficients of Optimized Runge-Kutta Methods with Automatic
!          Step Size Control for Compressible Computational Fluid Dynamics",
!         2021, DOI:10.5281/zenodo.4671927,
!         https://github.com/ranocha/Optimized-RK-CFD
!
      error_control = .true.
      fsal          = .true.
      evolve    => evolve_3sstarplus
      registers = 4
      order     = 5
      stages    = 10
      cfl       = 4.2d+00 * cfl
      name_int  = "Optimized 5ᵗʰ-order 10-step embedded RK5(4)10 FSAL"

      allocate(c(stages+1), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
      c (:)   = [  0.000000000000000000000000000000000000d+00,                 &
                   2.597883554788674084039539165398464630d-01,                 &
                   9.904573247592460887087003212056568980d-02,                 &
                   2.155511890524058691860390281856497503d-01,                 &
                   5.007950088969676776844289399972611534d-01,                 &
                   5.592251911688643533787800688765883636d-01,                 &
                   5.449986978853637084972622392134732553d-01,                 &
                   7.615224694532590139829150720490417596d-01,                 &
                   8.427062083267360939805493320684741215d-01,                 &
                   9.152209805057669959657927210873423883d-01,                 &
                   1.000000000000000000000000000000000000d+00 ]
      dl(:)   = [  1.000000000000000000000000000000000000d+00,                 &
                  -1.331778419508803397033287009506932673d-01,                 &
                   8.260422814750207498262063505871077303d-01,                 &
                   1.513700425755728332485300719652378197d+00,                 &
                  -1.305810059935023735972298885749903694d+00,                 &
                   3.036678802924163246003321318996156380d+00,                 &
                  -1.449458274398895177922690618003584514d+00,                 &
                   3.834313899176362315089976408899373409d+00,                 &
                   4.122293760012985409330881631526514714d+00,                 &
                   0.000000000000000000000000000000000000d+00 ]
      bt(:)   = [  2.597883554788674084039539165398464630d-01,                 &
                   1.777008889438867858759149597539211023d-02,                 &
                   2.481636629715501931294746189266601496d-01,                 &
                   7.941736871152005775821844297293296135d-01,                 &
                   3.885391285642019129575902994397298066d-01,                 &
                   1.455051657916305055730603387469193768d-01,                 &
                   1.587517385964749337690916959584348979d-01,                 &
                   1.650605617880053419242434594242509601d-01,                 &
                   2.118093284937153836908655490906875007d-01,                 &
                   1.559392342362059886106995325687547506d-01 ]
      bh(:)   = [ -2.019255440012066080909442770590267512d-02,                 &
                   2.737903480959184339932730854141598275d-02,                 &
                   3.028818636145965534365173822296811090d-01,                 &
                  -3.656843880622222190071445247906780540d-02,                 &
                   3.982664774676767729863101188528827405d-01,                 &
                  -5.715959421140685436681459970502471634d-02,                 &
                   9.849855103848558320961101178888983150d-02,                 &
                   6.654601552456084978615342374581437947d-02,                 &
                   9.073479542748112726465375642050504556d-02,                 &
                   8.432289325330803924891866923939606351d-02,                 &
                   4.529095628204896774513180907141004447d-02 ]
      gm(1,:) = [  0.000000000000000000000000000000000000d+00,                 &
                   4.043660121685749695640462197806189975d-01,                 &
                  -8.503427289575839690883191973980814832d-01,                 &
                  -6.950894175262117526410215315179482885d+00,                 &
                   9.238765192731084931855438934978371889d-01,                 &
                  -2.563178056509891340215942413817786020d+00,                 &
                   2.545744879365226143946122067064118430d-01,                 &
                   3.125831707411998258746812355492206137d-01,                 &
                  -7.007114414440507927791249989236719346d-01,                 &
                   4.839621016023833375810172323297465039d-01 ]
      gm(2,:) = [  1.000000000000000000000000000000000000d+00,                 &
                   6.871467028161416909922221357014564412d-01,                 &
                   1.093024748914750833700799552463885117d+00,                 &
                   3.225975379607193001678365742708874597d+00,                 &
                   1.041153702510101386914019859778740444d+00,                 &
                   1.292821487912164945157744726076279306d+00,                 &
                   7.391462755788122847651304143259254381d-01,                 &
                   1.239129251371800313941948224441873274d-01,                 &
                   1.842753472370123193132193302369345580d-01,                 &
                   5.712788998796583446479387686662738843d-02 ]
      gm(3,:) = [  0.000000000000000000000000000000000000d+00,                 &
                   0.000000000000000000000000000000000000d+00,                 &
                   0.000000000000000000000000000000000000d+00,                 &
                  -2.393405133244194727221124311276648940d+00,                 &
                  -1.902854422421760920850597670305403139d+00,                 &
                  -2.820042207399977261483046412236557428d+00,                 &
                  -1.832698465277380999601896111079977378d+00,                 &
                  -2.199094483084671192328083958346519535d-01,                 &
                  -4.082430635847870963724591602173546218d-01,                 &
                  -1.377669797880289713535665985132703979d-01 ]

      betas(:) = [ 4.5d-01, -1.3d-01, 0.0d+00 ]

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

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

    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)'."
        write(*,"(1x,a)") "Available embedded methods: 'ssprk3(2)4',"       // &
                          " 'ssprk3(2)m', 'rk3(2)5', 'rk4(3)9', 'rk5(4)10',"// &
                          "'rk3(2)5f', 'rk4(3)9f', 'rk5(4)10f'."
      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 */

! initialize other modules
!
    call initialize_schemes(verbose, status)
    if (status /=0) &
      call print_message(loc, "Could not initialize module SCHEMES!")
    if (check_status(status /= 0)) return

    call initialize_sources(verbose, status)
    if (status /=0) &
      call print_message(loc, "Could not initialize module SOURCES!")
    if (check_status(status /= 0)) return

    call initialize_forcing(verbose, status)
    if (status /=0) &
      call print_message(loc, "Could not initialize module FORCING!")
    if (check_status(status /= 0)) return

    call initialize_gravity(verbose, status)
    if (status /=0) &
      call print_message(loc, "Could not initialize module GRAVITY!")
    if (check_status(status /= 0)) return

    call initialize_shapes(status)
    if (status /=0) &
      call print_message(loc, "Could not initialize module SHAPES!")
    if (check_status(status /= 0)) return

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

    use forcing, only : finalize_forcing
    use gravity, only : finalize_gravity
    use helpers, only : print_section, print_parameter, print_message
    use schemes, only : finalize_schemes
    use shapes , only : finalize_shapes
    use sources, only : finalize_sources

    implicit none

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

    character(len=*), parameter :: loc = 'EVOLUTION::finalize_evolution()'

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

    call finalize_shapes(status)
    if (status /=0) &
      call print_message(loc, "Could not finalize module SHAPES!")

    call finalize_gravity(status)
    if (status /=0) &
      call print_message(loc, "Could not finalize module GRAVITY!")

    call finalize_forcing(status)
    if (status /=0) &
      call print_message(loc, "Could not finalize module FORCING!")

    call finalize_sources(status)
    if (status /=0) &
      call print_message(loc, "Could not finalize module SOURCES!")

    call finalize_schemes(status)
    if (status /=0) &
      call print_message(loc, "Could not finalize module SCHEMES!")

    if (allocated(adecay)) deallocate(adecay)
    if (allocated(aglm))   deallocate(aglm)
    if (allocated(c))      deallocate(c)
    if (allocated(dl))     deallocate(dl)
    if (allocated(bt))     deallocate(bt)
    if (allocated(bh))     deallocate(bh)
    if (allocated(gm))     deallocate(gm)

! 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

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

    use equations, only : magnetized
    use forcing  , only : print_forcing
    use helpers  , only : print_section, print_parameter
    use schemes  , only : print_schemes
    use shapes   , only : print_shapes
    use sources  , only : print_sources

    implicit none

    logical, intent(in) :: verbose

!-------------------------------------------------------------------------------
!
    call print_sources(verbose)
    call print_forcing(verbose)
    call print_shapes(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

    call print_schemes(verbose)

!-------------------------------------------------------------------------------
!
  end subroutine print_evolution
!
!===============================================================================
!
! subroutine ADVANCE:
! ------------------
!
!   Subroutine advances the solution by one time step using the selected time
!   integration method.
!
!   Arguments:
!
!     status - the subroutine call status: 0 for success, otherwise failure;
!
!===============================================================================
!
  subroutine advance(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
!
    integer, intent(out) :: status
!
!-------------------------------------------------------------------------------
!
! 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, status)
      if (status /= 0) go to 100

! 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

!-------------------------------------------------------------------------------
!
  end subroutine advance
!
!===============================================================================
!
! subroutine INITIALIZE_TIME_STEP:
! -------------------------------
!
!   Subroutine estimates the initial time step.
!
!   References:
!
!     [1] Courant, R., Friedrichs, K., Lewy, H.,
!         "Über die partiellen Differenzengleichungen der mathematischen Physik",
!         Mathematische Annalen, 1928, volume 100, pages 32–74,
!         DOI: 10.1007/BF01448839
!     [2] Hairer, E., Nørsett, S. P., Wanner, G.,
!         "Solving Ordinary Differential Equations I: Nonstiff Problems",
!         Springer-Verlag Berlin Heidelberg, 1993,
!         ISBN: 978-3-540-56670-0, DOI: 10.1007/978-3-540-78862-1
!
!===============================================================================
!
  subroutine initialize_time_step()

    use blocks     , only : block_data, list_data, get_nleafs
    use coordinates, only : nc => ncells, nb, ne, adx, ady
#if NDIMS == 3
    use coordinates, only : adz
#endif /* NDIMS == 3 */
    use equations  , only : nf, maxspeed, cmax, cmax2
    use helpers    , only : print_message
#ifdef MPI
    use mpitools   , only : reduce_maximum, reduce_sum
#endif /* MPI */
    use sources    , only : viscosity, resistivity, update_sources

    implicit none

    type(block_data), pointer  :: pdata

    logical                    :: flag
    integer                    :: mlev, n, status
    real(kind=8)               :: dx_min, fnorm, h0, h1
    real(kind=8), dimension(3) :: d

    real(kind=8), dimension(:,:,:,:), allocatable :: sc, df

    real(kind=8), parameter :: eps = tiny(cmax)

    character(len=*), parameter :: loc = 'EVOLUTION:initialize_time_step()'
!
!-------------------------------------------------------------------------------
!
    cmax   = eps
    mlev   = 1

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

      mlev = max(mlev, pdata%meta%level)
      cmax = max(cmax, maxspeed(pdata%q(:,:,:,:)))

      pdata => pdata%next
    end do

#ifdef MPI
    call reduce_maximum(cmax)
    call reduce_maximum(mlev)
#endif /* MPI */

    cmax2 = cmax * cmax

#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 */

! determine the time step due to the stability
!
    dth = dx_min / max(cmax                                                    &
                        + 2.0d+00 * max(viscosity, resistivity) / dx_min, eps)
    dt  = cfl * dth

! determine the time step due to the error
!
    if (error_control) then

#if NDIMS == 3
      allocate(sc(nf,nc,nc,nc),df(nf,nc,nc,nc), stat=status)
#else /* NDIMS == 3 */
      allocate(sc(nf,nc,nc, 1),df(nf,nc,nc, 1), stat=status)
#endif /* NDIMS == 3 */
#ifdef MPI
      call reduce_maximum(status)
#endif /* MPI */
      if (status > 0) then
        call print_message(loc,                                                &
                       "Cannot allocate memory for the time step estimation!")
        go to 100
      end if

      fnorm = 1.0d+00 / (get_nleafs() * nc**NDIMS)

      d(:) = 0.0d+00

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

        call update_increment(pdata)
        call update_sources(pdata, 0.0d+00, 0.0d+00, pdata%du(:,:,:,:))

#if NDIMS == 3
        sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1))

        d(1) = d(1) + sum((pdata%uu(:,nb:ne,nb:ne,nb:ne,1) / sc(:,:,:,:))**2)
        d(2) = d(2) + sum((pdata%du(:,nb:ne,nb:ne,nb:ne)   / sc(:,:,:,:))**2)
#else /* NDIMS == 3 */
        sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne,  :  ,1))

        d(1) = d(1) + sum((pdata%uu(:,nb:ne,nb:ne,  :  ,1) / sc(:,:,:,:))**2)
        d(2) = d(2) + sum((pdata%du(:,nb:ne,nb:ne,  :  )   / sc(:,:,:,:))**2)
#endif /* NDIMS == 3 */

        pdata => pdata%next
      end do

#ifdef MPI
      call reduce_sum(d(1:2))
#endif /* MPI */

      d(1:2) = sqrt(fnorm * d(1:2))

      if (minval(d(1:2)) >= 1.0d-05) then
        h0 = 1.0d-02 * d(1) / d(2)
      else
        h0 = 1.0d-06
      end if

      n    = 10
      flag = .true.

      do while(flag)

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

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

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

          pdata => pdata%next
        end do

        call update_variables(time + h0, h0, status)

        if (status /= 0) then
          h0   = 2.5d-01 * h0
          flag = n > 0
          n    = n - 1
        else
          flag = .false.
        end if

      end do

      if (status /= 0) then
        call print_message(loc, "Could not estimate the initial time step " // &
                         "due to the error. Setting it to the CFL time step.")
        dte = dt
      else
        pdata => list_data
        do while (associated(pdata))

#if NDIMS == 3
          df(:,:,:,:) = pdata%du(:,nb:ne,nb:ne,nb:ne)
          sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1))
#else /* NDIMS == 3 */
          df(:,:,:,:) = pdata%du(:,nb:ne,nb:ne,  :  )
          sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne,  :  ,1))
#endif /* NDIMS == 3 */

          call update_increment(pdata)
          call update_sources(pdata, time + h0, 0.0d+00, pdata%du(:,:,:,:))

#if NDIMS == 3
          df(:,:,:,:) = df(:,:,:,:) - pdata%du(:,nb:ne,nb:ne,nb:ne)
#else /* NDIMS == 3 */
          df(:,:,:,:) = df(:,:,:,:) - pdata%du(:,nb:ne,nb:ne,  :  )
#endif /* NDIMS == 3 */
          d(3) = d(3) + sum((df(:,:,:,:) / sc(:,:,:,:))**2)

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

          pdata => pdata%next
        end do

#ifdef MPI
        call reduce_sum(d(3:3))
#endif /* MPI */

        d(3) = sqrt(fnorm * d(3)) / h0

        if (max(d(2), d(3)) > 1.0d-15) then
          h1 = (1.0d-02 / max(d(2), d(3)))**(1.0d+00 / (order + 1))
        else
          h1 = max(1.0d-06, 1.0d-03 * h0)
        end if

        dte = min(1.0d+02 * h0, h1)
        dt  = min(dt, dte)

      end if

      if (allocated(sc)) deallocate(sc)
      if (allocated(df)) deallocate(df)

      100 continue

    end if

!-------------------------------------------------------------------------------
!
  end subroutine initialize_time_step
!
!===============================================================================
!
! subroutine NEW_TIME_STEP:
! ------------------------
!
!   Subroutine estimates the new time step from the maximum speed in the system
!   and source term constraints.
!
!   References:
!
!     [1] Courant, R., Friedrichs, K., Lewy, H.,
!         "Über die partiellen Differenzengleichungen der mathematischen Physik",
!         Mathematische Annalen, 1928, volume 100, pages 32–74,
!         DOI: 10.1007/BF01448839
!
!===============================================================================
!
  subroutine new_time_step()

    use blocks        , only : block_data, list_data
    use coordinates   , only : adx, ady
#if NDIMS == 3
    use coordinates   , only : adz
#endif /* NDIMS == 3 */
    use equations     , only : maxspeed, cmax, cmax2
#ifdef MPI
    use mpitools      , only : reduce_maximum
#endif /* MPI */
    use sources       , only : viscosity, resistivity

    implicit none

    type(block_data), pointer :: pdata

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

    real(kind=8), parameter   :: eps = tiny(cmax)
!
!-------------------------------------------------------------------------------
!
    cmax   = eps
    mlev   = 1

    pdata => list_data

    do while (associated(pdata))

      mlev = max(mlev, pdata%meta%level)
      cm   = maxspeed(pdata%q(:,:,:,:))
      cmax = max(cmax, cm)

      pdata => pdata%next

    end do

#ifdef MPI
    call reduce_maximum(cmax)
    call reduce_maximum(mlev)
#endif /* MPI */

    cmax2 = cmax * cmax

#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 */

    dth = dx_min / max(cmax                                                    &
                        + 2.0d+00 * max(viscosity, resistivity) / dx_min, eps)
    dt  = cfl * dth
    dt  = min(dt, dtp)
    if (error_control) dt = min(dt, dte)

!-------------------------------------------------------------------------------
!
  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

    integer      :: status
    real(kind=8) :: t
!
!-------------------------------------------------------------------------------
!
    status = 1

    do while(status /= 0)

      t = time + dt

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

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

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

        pdata => pdata%next
      end do

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

        call update_increment(pdata)
        call update_sources(pdata, t, dt, pdata%du(:,:,:,:))

        pdata => pdata%next
      end do

      call boundary_fluxes()

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

        pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dt * 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(t, dt, status)

      if (status /= 0) dt = 2.5d-01 * dt

    end do

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

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

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

      pdata => pdata%next
    end do

!-------------------------------------------------------------------------------
!
  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

    integer      :: status
    real(kind=8) :: t, ds
!
!-------------------------------------------------------------------------------
!
    status = 1

    do while(status /= 0)

      ds = 5.0d-01 * dt

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

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

        call update_increment(pdata)
        call update_sources(pdata, t, 0.0d+00, 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(t, dt, status)

      if (status /= 0) go to 100

!= 2nd step: U(n+1) = 1/2 U(n) + 1/2 U(1) + 1/2 dt * F[U(1)]
!
      pdata => list_data
      do while (associated(pdata))

        call update_increment(pdata)
        call update_sources(pdata, t, ds, pdata%du(:,:,:,:))

        pdata => pdata%next
      end do

      call boundary_fluxes()

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

        pdata%uu(:,:,:,:,2) = 5.0d-01 * (pdata%uu(:,:,:,:,1)                   &
                                       + pdata%uu(:,:,:,:,2))                  &
                                  + 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(t, ds, status)

      100 continue

      if (status /= 0) dt = 2.5d-01 * dt

    end do

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

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

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

      pdata => pdata%next
    end do

!-------------------------------------------------------------------------------
!
  end subroutine evolve_rk2
!
!===============================================================================
!
! 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

    integer      :: status
    real(kind=8) :: t, ds

    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
!
!-------------------------------------------------------------------------------
!
    status = 1

    do while(status /= 0)

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

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

        call update_increment(pdata)
        call update_sources(pdata, t, 0.0d+00, 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(t, dt, status)

      if (status /= 0) go to 100

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

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

        call update_increment(pdata)
        call update_sources(pdata, t, ds, 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(t, ds, status)

      if (status /= 0) go to 100

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

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

        call update_increment(pdata)
        call update_sources(pdata, t, ds, pdata%du(:,:,:,:))

        pdata => pdata%next
      end do

      call boundary_fluxes()

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

        pdata%uu(:,:,:,:,2) = 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(t, dt, status)

      100 continue

      if (status /= 0) dt = 2.5d-01 * dt

    end do

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

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

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

      pdata => pdata%next
    end do

!-------------------------------------------------------------------------------
!
  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

    integer      :: status
    real(kind=8) :: t, ds
!
!-------------------------------------------------------------------------------
!
    status = 1

    do while(status /= 0)

!= 1st step: U(1) = U(n) + 1/2 dt F[U(n)]
!
      ds = 5.0d-01 * dt
      t  = time + ds

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

        call update_increment(pdata)
        call update_sources(pdata, t, 0.0d+00, 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(t, ds, status)

      if (status /= 0) go to 100

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

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

        call update_increment(pdata)
        call update_sources(pdata, t, ds, 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(t, ds, status)

      if (status /= 0) go to 100

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

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

        call update_increment(pdata)
        call update_sources(pdata, t, ds, 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)                   &
                                       + ds * pdata%du(:,:,:,:)) / 3.0d+00

        pdata => pdata%next
      end do

      call update_variables(t, ds, status)

      if (status /= 0) go to 100

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

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

        call update_increment(pdata)
        call update_sources(pdata, t, ds, 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

      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(t, dt, status)

      100 continue

      if (status /= 0) dt = 2.5d-01 * dt

    end do

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

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

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

      pdata => pdata%next
    end do

!-------------------------------------------------------------------------------
!
  end subroutine evolve_ssprk34
!
!===============================================================================
!
! 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

    integer      :: status
    real(kind=8) :: t, ds

    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
!
!-------------------------------------------------------------------------------
!
    status = 1

    do while(status /= 0)

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

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

        call update_increment(pdata)
        call update_sources(pdata, t, 0.0d+00, 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(t, ds, status)

      if (status /= 0) go to 100

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

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

        call update_increment(pdata)
        call update_sources(pdata, t, ds, 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(t, ds, status)

      if (status /= 0) go to 100

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

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

        call update_increment(pdata)
        call update_sources(pdata, t, ds, 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(t, ds, status)

      if (status /= 0) go to 100

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

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

        call update_increment(pdata)
        call update_sources(pdata, t, ds, pdata%du(:,:,:,:))

        pdata => pdata%next
      end do

      call boundary_fluxes()

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

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

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

        pdata => pdata%next
      end do

      call update_variables(t, ds, status)

      if (status /= 0) go to 100

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

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

        call update_increment(pdata)
        call update_sources(pdata, t, ds, pdata%du(:,:,:,:))

        pdata => pdata%next
      end do

      call boundary_fluxes()

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

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

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

        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(t, dt, status)

      100 continue

      if (status /= 0) dt = 2.5d-01 * dt

    end do

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

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

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

      pdata => pdata%next
    end do

!-------------------------------------------------------------------------------
!
  end subroutine evolve_ssprk35
!
!===============================================================================
!
! 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      :: i, status
    real(kind=8) :: t

    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
!
!-------------------------------------------------------------------------------
!
    ds(:) = c(:) * dt

    status = 1

    do while(status /= 0)

!= 1st step: U(2) = U(1)
!
      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) = [1 + dt/6 L] U(1), for i = 1, ..., 5
!
      do i = 1, 5

        t = time + ds(i)

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

          call update_increment(pdata)
          call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:))

          pdata => pdata%next
        end do

        call boundary_fluxes()

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

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

          pdata => pdata%next
        end do

        call update_variables(t, ds(1), status)

        if (status /= 0) go to 100

      end do ! i = 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(:,:,:,:,3) = (pdata%uu(:,:,:,:,3)                             &
                            + 9.0d+00 * pdata%uu(:,:,:,:,2)) / 2.5d+01
        pdata%uu(:,:,:,:,2) = 1.5d+01 * pdata%uu(:,:,:,:,3)                    &
                            - 5.0d+00 * pdata%uu(:,:,:,:,2)

        pdata => pdata%next
      end do

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

        t = time + ds(i)

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

          call update_increment(pdata)
          call update_sources(pdata, t, ds(1), pdata%du(:,:,:,:))

          pdata => pdata%next
        end do

        call boundary_fluxes()

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

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

          pdata => pdata%next
        end do

        call update_variables(t, ds(1), status)

        if (status /= 0) go to 100

      end do ! i = 6, 9

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

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

        call update_increment(pdata)
        call update_sources(pdata, t, 1.0d-01 * dt, pdata%du(:,:,:,:))

        pdata => pdata%next
      end do

      call boundary_fluxes()

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

        pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,3)                              &
                            + 6.0d-01 * pdata%uu(:,:,:,:,2)                    &
                            + 1.0d-01 * dt * 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(t, dt, status)

      100 continue

      if (status /= 0) dt = 2.5d-01 * dt

    end do

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

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

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

      pdata => pdata%next
    end do

!-------------------------------------------------------------------------------
!
  end subroutine evolve_ssprk4_10
!
!===============================================================================
!
! EMBEDDED METHODS
!
!===============================================================================
!
! 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] 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
!     [2] 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
    use boundaries , only : boundary_fluxes
    use coordinates, only : nb, ne
    use equations  , only : errors, ibp, cmax
    use helpers    , only : print_message
    use sources    , only : update_sources

    implicit none

    type(block_data), pointer :: pdata

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

    character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssprk2_m()'

!-------------------------------------------------------------------------------
!
    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 2nd order solution, U(1)
! pdata%uu(:,:,:,:,3) - the new 1st order solution, U(2)
!
    do while(test)

!= preparatory step: U(1) = U(n), U(2) = 0, U(3) = 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

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

        dtm = c(i) * dt
        tm  = time + dtm

        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)                            &
                              + bt(i) * dt * pdata%du(:,:,:,:)
          pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,3)                            &
                              + bh(i) * dt * pdata%du(:,:,:,:)

          pdata => pdata%next
        end do

        call update_variables(tm, dtm, status)
        if (status /= 0) go to 100

      end do ! i = 1, stages - 1

!= 2nd step: U(1) = (m-1)/m U(1) + 1/m U(2) + dt/m F[U(1)]
!
      i   = stages
      dtm = c(i) * dt
      tm  = time + dtm

      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) = dl(1) * pdata%uu(:,:,:,:,2)                      &
                            + dl(2) * pdata%uu(:,:,:,:,1)                      &
                            + bt(i) * dt * pdata%du(:,:,:,:)
        pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,3)                              &
                            + bh(i) * dt * pdata%du(:,:,:,:)

        pdata => pdata%next
      end do

      call update_variables(tm, dtm, status)
      if (status /= 0) go to 100

!= 3rd 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

      100 continue

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

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

      else
        if (status == 0) then
          dt   = dte
          nrej = nrej + 1 ! rejection count in the current step
        else
          dt     = 2.5d-01 * dt
        end if
        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, status)
        if (status /= 0) &
          call print_message(loc, "Possible bug: the revert to " //            &
                                  "the previous state found it unphysical.")

      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

!-------------------------------------------------------------------------------
!
  end subroutine evolve_ssprk2_m
!
!===============================================================================
!
! 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 helpers    , only : print_message
    use sources    , only : update_sources

    implicit none

    type(block_data), pointer :: pdata

    logical         :: test
    integer         :: nrej, i, status
    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

    character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssp324()'

!-------------------------------------------------------------------------------
!
    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, status)
        if (status /= 0) go to 100

      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, status)
      if (status /= 0) go to 100

!= 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, status)
      if (status /= 0) go to 100

!= 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

      100 continue

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

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

      else
        if (status == 0) then
          dt   = dte
          nrej = nrej + 1 ! rejection count in the current step
        else
          dt     = 2.5d-01 * dt
        end if
        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, status)
        if (status /= 0) &
          call print_message(loc, "Possible bug: the revert to " //            &
                                  "the previous state found it unphysical.")
      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

!-------------------------------------------------------------------------------
!
  end subroutine evolve_ssp324
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK32M:
! --------------------------
!
!   Subroutine advances the solution by one time step using the 3ʳᵈ-order
!   m-stage (m=n²) embedded Strong Stability Preserving Runge-Kutta time
!   integration method SSP3(2)m with the error control.
!
!   References:
!
!     [1] 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
!     [2] Gottlieb, S., Ketcheson, D., Shu, C.-W.,
!         "Strong stability preserving Runge-Kutta and multistep time
!          discretizations",
!         2011, World Scientific Publishing
!
!===============================================================================
!
  subroutine evolve_ssprk32m()

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

    implicit none

    type(block_data), pointer :: pdata

    logical     , save :: first = .true.
    integer     , save :: i1, i2, i3, i4, n, m
    real(kind=8), save :: f1, f2, f3, f4

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

    character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssprk32m()'

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

      n = 2
      do while(n**2 < stages)
        n = n + 1
      end do
      m  = stages - n
      i1 = (n - 1) * (n - 2) / 2
      i2 = i1 + 1
      i3 = n * (n + 1) / 2
      i4 = i3 + 1

      fc = real(2 * n - 1, kind=8)
      f1 = real(n - 1, kind=8) / fc
      f2 = real(n    , kind=8) / fc
      f4 = real(m    , kind=8) / stages
      f3 = 1.0d+00 - f4

      first = .false.

    end if

    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 intermediate solution, U(2)
! pdata%uu(:,:,:,:,4) - the new 2nd order solution, U(3)
!
    do while(test)

! initiate the fractional time step
!
      dh = dt / m

!= 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 step: U(1) = U(1) + dt/r F[U(1)], for i = 1...(n-1)*(n-2)/2
!
      do i = 1, i1

        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, status)
        if (status /= 0) go to 100

      end do ! i = 1, i1

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

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

        pdata => pdata%next
      end do

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

        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, status)
        if (status /= 0) go to 100

      end do ! i = i2, i3

!= 4th step: U(3) = (1 - r/s) * U(n) + r/s * U(1),
!            U(1) = [(n - 1) * U(1) + n * U(2)] / (2*n - 1),
!            U(3) = U(3) - r/s * U(1)
!
      pdata => list_data
      do while (associated(pdata))

        pdata%uu(:,:,:,:,4) = f3 * pdata%uu(:,:,:,:,1)                         &
                            + f4 * pdata%uu(:,:,:,:,2)
        pdata%uu(:,:,:,:,2) = f1 * pdata%uu(:,:,:,:,2)                         &
                            + f2 * pdata%uu(:,:,:,:,3)
        pdata%uu(:,:,:,:,4) =      pdata%uu(:,:,:,:,4)                         &
                            - f4 * pdata%uu(:,:,:,:,2)

        pdata => pdata%next
      end do

      call update_variables(tm, dtm, status)
      if (status /= 0) go to 100

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

        tm  = time + (i - 1 - n) * 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, status)
        if (status /= 0) go to 100

      end do ! i = i4, stages

!= 7th step: U(3) = U(3) + r/s * U(1)
!
      pdata => list_data
      do while (associated(pdata))

        pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4) + f4 * pdata%uu(:,:,:,:,2)

        pdata => pdata%next
      end do

!= 8th 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,:,:,:,4) = adecay(pdata%meta%level)                     &
                                                       * pdata%uu(ibp,:,:,:,4)

          pdata => pdata%next
        end do
      end if

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

! 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

      100 continue

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

        errs(3) = errs(2)
        errs(2) = errs(1)
      else
        if (status == 0) then
          dt   = dte
          nrej = nrej + 1 ! rejection count in the current step
        else
          dt     = 2.5d-01 * dt
        end if
        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, status)
        if (status /= 0) &
          call print_message(loc, "Possible bug: the revert to " //            &
                                  "the previous state found it unphysical.")
      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

!-------------------------------------------------------------------------------
!
  end subroutine evolve_ssprk32m
!
!===============================================================================
!
! subroutine EVOLVE_3SSTARPLUS:
! ----------------------------
!
!   Subroutine advances the solution by one time step using the minimum
!   storage implementation of 3S*+ methods.
!
!   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.06836
!     [2] Gottlieb, S., Ketcheson, D., Shu C. - W.,
!         "Strong stability preserving Runge-Kutta and multistep
!          time discretizations",
!         World Scientific Publishing, 2011
!
!===============================================================================
!
  subroutine evolve_3sstarplus()

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

    implicit none

    type(block_data), pointer :: pdata

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

    character(len=*), parameter :: loc = 'EVOLUTION:evolve_3sstarplus()'

!-------------------------------------------------------------------------------
!
    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 intermediate solution, U(2)
! pdata%uu(:,:,:,:,4) - the new 2nd order solution, U(3)
!
    do while(test)

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

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

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

        pdata => pdata%next
      end do

!= 1st stage
!
      if (fsal_update) then
        pdata => list_data
        do while (associated(pdata))

          call update_increment(pdata)
          call update_sources(pdata, time, 0.0d+00, pdata%du(:,:,:,:))

          pdata => pdata%next
        end do

        call boundary_fluxes()
      end if

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

        pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)                              &
                            + bt(1) * dt * pdata%du(:,:,:,:)
        pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4)                              &
                            + bh(1) * dt * pdata%du(:,:,:,:)

        pdata => pdata%next
      end do

      dtm = c(2) * dt
      call update_variables(time, dtm, status)
      if (status /= 0) go to 100

!= remaining stages
!
      do i = 2, stages

          dtm = c(i) * dt
          tm  = time + dtm

          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(:,:,:,:,3) =           pdata%uu(:,:,:,:,3)                  &
                              +   dl(i) * pdata%uu(:,:,:,:,2)
          pdata%uu(:,:,:,:,2) = gm(1,i) * pdata%uu(:,:,:,:,2)                  &
                              + gm(2,i) * pdata%uu(:,:,:,:,3)                  &
                              + gm(3,i) * pdata%uu(:,:,:,:,1)                  &
                              + bt(i) * dt * pdata%du(:,:,:,:)
          pdata%uu(:,:,:,:,4) =           pdata%uu(:,:,:,:,4)                  &
                              + bh(i) * dt * pdata%du(:,:,:,:)

          pdata => pdata%next
        end do

        dtm = c(i+1) * dt
        call update_variables(tm, dtm, status)
        if (status /= 0) go to 100

      end do ! i = 2, stages

!= extra 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,:,:,:,4) = adecay(pdata%meta%level)                     &
                                                       * pdata%uu(ibp,:,:,:,4)

          pdata => pdata%next
        end do
      end if

!= extra step: FSAL
!
      if (fsal) then

        tm  = time + dt

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

          call update_increment(pdata)
          call update_sources(pdata, tm, 0.0d+00, pdata%du(:,:,:,:))

          pdata => pdata%next
        end do

        call boundary_fluxes()

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

          pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4)                            &
                              + bh(stages+1) * dt * pdata%du(:,:,:,:)

          pdata => pdata%next
        end do

      end if

!= error estimation
!
      call update_errors(2, 4)

! 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

      100 continue

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

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

        if (fsal) fsal_update = .true.
      else
        if (status == 0) then
          dt   = dte
          nrej = nrej + 1 ! rejection count in the current step
        else
          dt     = 2.5d-01 * dt
        end if
        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, status)
        if (status /= 0) &
          call print_message(loc, "Possible bug: the revert to " //            &
                                  "the previous state found it unphysical.")

        if (fsal) fsal_update = .false.
      end if

      niterations = niterations + 1

    end do

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

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

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

      pdata => pdata%next
    end do

    tm = time + dt
    call update_variables(tm, dt, status)
    if (status /= 0) &
      call print_message(loc, "Possible bug: the new state has been " //       &
                              "accepted, nevertheless the variable "  //       &
                              "update found it unphysical.")

!-------------------------------------------------------------------------------
!
  end subroutine evolve_3sstarplus
!
!===============================================================================
!
! 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)

    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 helpers    , only : print_message
    use schemes    , only : update_flux
    use workspace  , only : resize_workspace, work, work_in_use

    implicit none

    type(block_data), pointer, intent(inout) :: pdata

    logical, save :: first = .true.

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

    real(kind=8), dimension(NDIMS) :: dh, dhi

    real(kind=8), dimension(:,:,:,:,:)  , pointer, save :: f
    real(kind=8), dimension(:,:,:,:,:,:), pointer, save :: s

    character(len=*), parameter :: loc = 'EVOLUTION:update_increment()'

!-------------------------------------------------------------------------------
!
    if (first) then
      i = NDIMS * nf * nn**NDIMS
      j = 3 * i

      call resize_workspace(j, status)
      if (status /= 0) then
        call print_message(loc, "Could not resize the workspace!")
        go to 100
      end if

#if NDIMS == 3
      f(1:nf,1:nn,1:nn,1:nn,1:3)     => work(  1:i)
      s(1:nf,1:nn,1:nn,1:nn,1:2,1:3) => work(i+1:j)
#else /* NDIMS == 3 */
      f(1:nf,1:nn,1:nn,1: 1,1:2)     => work(  1:i)
      s(1:nf,1:nn,1:nn,1: 1,1:2,1:2) => work(i+1:j)
#endif /* NDIMS == 3 */

      first = .false.
    end if

    if (work_in_use) &
      call print_message(loc, "Workspace is being used right now! " //         &
                              "Corruptions can occur!")

    work_in_use = .true.

    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(:,:,:,:),                            &
                                    s(:,:,:,:,:,:), f(:,:,:,:,:))

! store the 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

      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

    work_in_use = .false.

    100 continue

!-------------------------------------------------------------------------------
!
  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;
!     status - status flag indicating if the update was successful;
!
!===============================================================================
!
  subroutine update_variables(tm, dtm, status)

    use blocks    , only : block_meta, list_meta
    use blocks    , only : block_data, list_data
    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
#ifdef MPI
    use mpitools  , only : reduce_maximum
#endif /* MPI */
    use shapes    , only : update_shapes

    implicit none

    real(kind=8), intent(in)  :: tm, dtm
    integer     , intent(out) :: status

    type(block_meta), pointer :: pmeta
    type(block_data), pointer :: pdata

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

    pdata => list_data
    do while (associated(pdata))
      pmeta => pdata%meta

      if (pmeta%update) then
        call update_primitive_variables(pdata%u, pdata%q, status)
        if (status /= 0) go to 100
      end if

      pdata => pdata%next
    end do

    100 continue
#ifdef MPI
    call reduce_maximum(status)
#endif /* MPI */
    if (status /= 0) go to 200

    call boundary_variables(tm, dtm)

    pdata => list_data
    do while (associated(pdata))
      pmeta => pdata%meta

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

      pdata => pdata%next
    end do

    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

    200 continue

!-------------------------------------------------------------------------------
!
  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