!!******************************************************************************
!!
!!  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-2019 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/>.
!!
!!******************************************************************************
!!
!! program: AMUN
!!
!!  AMUN is a code to perform numerical simulations in a fluid approximation on
!!  adaptive mesh for Newtonian and relativistic environments with or without
!!  magnetic field.
!!
!!
!!******************************************************************************
!
program amun

! include external subroutines used in this module
!
  use iso_fortran_env, only : error_unit
  use blocks         , only : initialize_blocks, finalize_blocks, get_nleafs
  use blocks         , only : build_leaf_list
  use boundaries     , only : initialize_boundaries, finalize_boundaries
  use boundaries     , only : print_boundaries, boundary_variables
  use coordinates    , only : initialize_coordinates, finalize_coordinates
  use coordinates    , only : print_coordinates
  use coordinates    , only : nn => bcells
  use domains        , only : initialize_domains, finalize_domains
  use equations      , only : initialize_equations, finalize_equations
  use equations      , only : print_equations
  use equations      , only : nv
  use evolution      , only : initialize_evolution, finalize_evolution
  use evolution      , only : print_evolution
  use evolution      , only : advance, new_time_step
  use evolution      , only : step, time, dt
  use gravity        , only : initialize_gravity, finalize_gravity
  use helpers        , only : print_welcome, print_section, print_parameter
  use integrals      , only : initialize_integrals, finalize_integrals
  use integrals      , only : store_integrals
  use interpolations , only : initialize_interpolations, finalize_interpolations
  use interpolations , only : print_interpolations
  use io             , only : initialize_io, finalize_io, print_io
  use io             , only : restart_snapshot_number, restart_from_snapshot
  use io             , only : read_snapshot_parameter
  use io             , only : read_restart_snapshot, write_restart_snapshot
  use io             , only : write_snapshot, next_tout, precise_snapshots
  use mesh           , only : initialize_mesh, finalize_mesh
  use mesh           , only : generate_mesh, store_mesh_stats
  use mpitools       , only : initialize_mpitools, finalize_mpitools
#ifdef MPI
  use mpitools       , only : bcast_integer_variable, reduce_sum_real_array
#endif /* MPI */
  use mpitools       , only : master, nprocs, nproc, check_status
  use operators      , only : initialize_operators, finalize_operators
  use parameters     , only : read_parameters, finalize_parameters
#ifdef MPI
  use parameters     , only : redistribute_parameters
#endif /* MPI */
  use parameters     , only : get_parameter
  use problems       , only : initialize_problems, finalize_problems
  use random         , only : initialize_random, finalize_random
  use refinement     , only : initialize_refinement, finalize_refinement
  use refinement     , only : print_refinement
  use schemes        , only : initialize_schemes, finalize_schemes
  use schemes        , only : print_schemes
  use shapes         , only : initialize_shapes, finalize_shapes, print_shapes
  use sources        , only : initialize_sources, finalize_sources
  use sources        , only : print_sources
  use timers         , only : initialize_timers, finalize_timers
  use timers         , only : start_timer, stop_timer, set_timer, get_timer
  use timers         , only : get_timer_total, timer_enabled, timer_description
  use timers         , only : get_count, ntimers
  use user_problem   , only : initialize_user_problem, finalize_user_problem

! module variables are not implicit by default
!
  implicit none

! the initialization success and status flags
!
  logical               :: initialization_succeeded
  logical               :: proceed = .true.
  integer               :: status  = 0
  integer               :: quit

! the run number (for restarted runs)
!
  integer               :: nrun    = 0

! default parameters
!
  character(len=64)     :: problem = "none"
  character(len=32)     :: eqsys   = "hydrodynamic"
  character(len=32)     :: eos     = "adiabatic"
  integer               :: ncells  = 8
  integer               :: nghosts = 2
  integer               :: toplev  = 1
  integer, dimension(3) :: bdims   = 1
  integer               :: nmax    = huge(1), ndat = 1
  real(kind=8)          :: xmin    = 0.0d+00, xmax    = 1.0d+00
  real(kind=8)          :: ymin    = 0.0d+00, ymax    = 1.0d+00
  real(kind=8)          :: zmin    = 0.0d+00, zmax    = 1.0d+00
  real(kind=8)          :: tmax    = 0.0d+00, trun = 9.999d+03, tsav = 3.0d+01
  real(kind=8)          :: dtnext  = 0.0d+00

! timer indices
!
  integer               :: iin, iev, itm
#ifdef PROFILE
  integer               :: ipr, ipi
#endif /* PROFILE */

! iteration and time variables
!
  integer               :: i, ed, eh, em, es, ec
  integer               :: nsteps = 1
  character(len=80)     :: fmt, tmp

  real(kind=8)          :: tbeg, thrs
  real(kind=8)          :: tm_curr, tm_exec, tm_conv, tm_last = 0.0d+00

#ifdef INTEL
! the type of the function SIGNAL should be defined for Intel compiler
!
  integer(kind=4)       :: signal
#endif /* INTEL */

#ifdef SIGNALS
! references to functions handling signals
!
#ifdef GNU
  intrinsic signal
#endif /* GNU */
  external  terminate

! signal definitions
!
  integer, parameter    :: SIGERR = -1
  integer, parameter    :: SIGINT =  2, SIGABRT =  6, SIGTERM = 15
#endif /* SIGNALS */

! an array to store execution times
!
  real(kind=8), dimension(ntimers) :: tm

! common block
!
  common /termination/ quit
!
!-------------------------------------------------------------------------------
!
! initialize local flags
!
  quit = 0

! initialize module TIMERS
!
  call initialize_timers()

! set timer descriptions
!
  call set_timer('INITIALIZATION', iin)
  call set_timer('EVOLUTION'     , iev)
  call set_timer('TERMINATION'   , itm)

#ifdef SIGNALS
! assign function terminate() to handle signals
!
#ifdef GNU
  if (signal(SIGINT , terminate)     == SIGERR) status = 1
  if (signal(SIGABRT, terminate)     == SIGERR) status = 2
  if (signal(SIGTERM, terminate)     == SIGERR) status = 3
#endif /* GNU */
#ifdef INTEL
  if (signal(SIGINT , terminate, -1) == SIGERR) status = 1
  if (signal(SIGABRT, terminate, -1) == SIGERR) status = 2
  if (signal(SIGTERM, terminate, -1) == SIGERR) status = 3
#endif /* INTEL */

! in the case of problems with signal handler assignment, quit the program
!
  if (status /= 0) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem initializing the signal handler!"
    call finalize_timers()
    call exit(1)
  end if
#endif /* SIGNALS */

! initialize module MPITOOLS
!
  call initialize_mpitools(status)

  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing MPITOOLS module!"
    end if
    go to 4000
  end if

! set the initiallization success flag
!
  initialization_succeeded = .false.

! start time accounting for the initialization
!
  call start_timer(iin)

! print welcome messages
!
  call print_welcome(master)
  call print_section(master, "Parallelization")
#ifdef MPI
  call print_parameter(master, "MPI processes", nprocs)
#endif /* MPI */

! initialize and read parameters from the parameter file
!
  if (master) call read_parameters(status)

  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem reading parameters!"
    end if
    go to 3900
  end if

#ifdef MPI
! redistribute parameters among all processors
!
  call redistribute_parameters(status)

  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem broadcasting parameters!"
    end if
    go to 3900
  end if
#endif /* MPI */

! initialize IO early to handle restart snapshots if necessary
!
  call initialize_io(master, status)

  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing IO module!"
    end if
    go to 3800
  end if

! get the run number
!
  nrun = max(1, restart_snapshot_number() + 1)

! if the run is from a restarted job, read the fixed parameters from
! the restart snapshot, otherwise, read them from the parameter file
!
  if (restart_from_snapshot()) then
    call read_snapshot_parameter("problem", problem, status)
    if (check_status(status /= 0)) go to 3800
    call read_snapshot_parameter("eqsys"  , eqsys  , status)
    if (check_status(status /= 0)) go to 3800
    call read_snapshot_parameter("eos"    , eos    , status)
    if (check_status(status /= 0)) go to 3800
    call read_snapshot_parameter("ncells" , ncells , status)
    if (check_status(status /= 0)) go to 3800
    call read_snapshot_parameter("maxlev" , toplev , status)
    if (check_status(status /= 0)) go to 3800
    call read_snapshot_parameter("rdims"  , bdims  , status)
    if (check_status(status /= 0)) go to 3800
    call read_snapshot_parameter("xmin"   , xmin   , status)
    if (check_status(status /= 0)) go to 3800
    call read_snapshot_parameter("xmax"   , xmax   , status)
    if (check_status(status /= 0)) go to 3800
    call read_snapshot_parameter("ymin"   , ymin   , status)
    if (check_status(status /= 0)) go to 3800
    call read_snapshot_parameter("ymax"   , ymax   , status)
    if (check_status(status /= 0)) go to 3800
#if NDIMS == 3
    call read_snapshot_parameter("zmin"   , zmin   , status)
    if (check_status(status /= 0)) go to 3800
    call read_snapshot_parameter("zmax"   , zmax   , status)
    if (check_status(status /= 0)) go to 3800
#endif /* NDIMS == 3 */
  else
    call get_parameter("problem"          , problem )
    call get_parameter("equation_system"  , eqsys   )
    call get_parameter("equation_of_state", eos     )
    call get_parameter("ncells"           , ncells  )
    call get_parameter("maxlev"           , toplev  )
    call get_parameter("xblocks"          , bdims(1))
    call get_parameter("yblocks"          , bdims(2))
#if NDIMS == 3
    call get_parameter("zblocks"          , bdims(3))
#endif /* NDIMS == 3 */
    call get_parameter("xmin"             , xmin    )
    call get_parameter("xmax"             , xmax    )
    call get_parameter("ymin"             , ymin    )
    call get_parameter("ymax"             , ymax    )
#if NDIMS == 3
    call get_parameter("zmin"             , zmin    )
    call get_parameter("zmax"             , zmax    )
#endif /* NDIMS == 3 */
  end if

! get the number of ghost zones
!
  call get_parameter("nghosts", nghosts)

! get the execution termination parameters
!
  call get_parameter("nmax", nmax)
  call get_parameter("tmax", tmax)
  call get_parameter("trun", trun)
  call get_parameter("tsav", tsav)

! correct the run time by the save time
!
  trun   = trun - tsav / 6.0d+01

! initialize dtnext
!
  dtnext = 2.0d+00 * tmax

! get integral calculation interval
!
  call get_parameter("ndat", ndat)

! print the problem name here, so in initialize_user_problem() we can print
! parameters
!
  call print_section(master, "Problem")
  call print_parameter(master, "problem name", trim(problem))

! initialize the remaining modules
!
  call initialize_random(1, 0)
  if (check_status(status /= 0)) go to 3700
  call initialize_equations(eqsys, eos, master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing EQUATIONS module!"
    end if
    go to 3600
  end if
  call initialize_interpolations(nghosts, master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing INTERPOLATIONS module!"
    end if
    go to 3500
  end if
  call initialize_coordinates(ncells, nghosts, toplev, bdims, xmin, xmax,      &
                                       ymin, ymax, zmin, zmax, master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing COORDINATES module!"
    end if
    go to 3400
  end if
  call initialize_blocks(nv, nv, nn, master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing BLOCKS module!"
    end if
    go to 3300
  end if
  call initialize_operators(master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing OPERATORS module!"
    end if
    go to 3200
  end if
  call initialize_sources(master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing SOURCES module!"
    end if
    go to 3100
  end if
  call initialize_user_problem(problem, master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing USER_PROBLEM module!"
    end if
    go to 3000
  end if
  call initialize_problems(problem, master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing PROBLEMS module!"
    end if
    go to 2900
  end if
  call initialize_domains(problem, master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing DOMAINS module!"
    end if
    go to 2800
  end if
  call initialize_boundaries(master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing BOUNDARIES module!"
    end if
    go to 2700
  end if
  call initialize_refinement(master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing REFINEMENT module!"
    end if
    go to 2600
  end if
  call initialize_mesh(nrun, master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing MESH module!"
    end if
    go to 2500
  end if
  call initialize_shapes(master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing SHAPES module!"
    end if
    go to 2400
  end if
  call initialize_gravity(master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing GRAVITY module!"
    end if
    go to 2300
  end if
  call initialize_schemes(master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing SCHEMES module!"
    end if
    go to 2200
  end if
  call initialize_evolution(master, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing EVOLUTION module!"
    end if
    go to 2100
  end if
  call initialize_integrals(master, nrun, status)
  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Problem initializing INTEGRALS module!"
    end if
    go to 2000
  end if

! print module information
!
  call print_equations(master)
  call print_sources(master)
  call print_coordinates(master)
  call print_boundaries(master)
  call print_shapes(master)
  call print_refinement(master)
  call print_evolution(master)
  call print_schemes(master)
  call print_interpolations(master)
  call print_io(master)

! check if we initiate new problem or restart previous job
!
  if (restart_from_snapshot()) then

! reconstruct the meta and data block structures from a given restart file
!
    call read_restart_snapshot(status)

    if (check_status(status /= 0)) then
      if (master) then
        write(error_unit,"('[AMUN::program]: ', a)")                           &
                         "Problem reading restart snapshot!"
      end if
      go to 1000
    end if

! update the list of leafs
!
    call build_leaf_list(status)

    if (check_status(status /= 0)) then
      if (master) then
        write(error_unit,"('[AMUN::program]: ', a)")                           &
                         "Problem building the list of leafs!"
      end if
      go to 1000
    end if

! update boundaries
!
    call boundary_variables(time, dtnext)

  else

! generate the initial mesh, refine that mesh to the desired level according to
! the initialized problem
!
    call generate_mesh(status)

    if (check_status(status /= 0)) then
      if (master) then
        write(error_unit,"('[AMUN::program]: ', a)")                           &
                         "Problem generating the initial mesh!"
      end if
      go to 1000
    end if

! update boundaries
!
    call boundary_variables(0.0d+00, dtnext)

! calculate new timestep
!
    call new_time_step(dtnext)

  end if

! store mesh statistics
!
  call store_mesh_stats(step, time)

! store integrals and data to a file
!
  if (.not. restart_from_snapshot()) then

    call store_integrals()
    call write_snapshot()

  end if

! stop time accounting for the initialization
!
  call stop_timer(iin)

! set the initiallization success flag
!
  initialization_succeeded = .true.

! start time accounting for the evolution
!
  call start_timer(iev)

! get current time in seconds
!
  if (master) &
    tbeg = time

! print progress info on master processor
!
  if (master) then

! initialize estimated remaining time of calculations
!
    ed = 9999
    eh =   23
    em =   59
    es =   59

! print progress info
!
    write(*,*)
    write(*,"(1x,a)"   ) "Evolving the system:"
    write(*,'(4x,a4,5x,a4,11x,a2,12x,a6,7x,a3)') 'step', 'time', 'dt'          &
                                                 , 'blocks', 'ETA'
#ifdef INTEL
    write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' //      &
            ',1i2.2,"s",15x,a1,$)')                                            &
                        step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
#else /* INTEL */
    write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' //      &
            ',1i2.2,"s",15x,a1)',advance="no")                                 &
                        step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
#endif /* INTEL */

  end if

! main loop
!
  do while((nsteps <= nmax) .and. (time < tmax) .and. proceed)

! get the next snapshot time
!
    if (precise_snapshots) dtnext = next_tout() - time

! performe one step evolution
!
    call advance(dtnext, status)

    if (check_status(status /= 0)) then
      if (master) then
        write(error_unit,"('[AMUN::program]: ', a)")                           &
                         "Advancing to next step failed!"
      end if
      go to 1000
    end if

! advance the iteration number and time
!
    time   = time   + dt
    step   = step   +  1
    nsteps = nsteps +  1

! get current time in seconds
!
    tm_curr = get_timer_total()

! compute elapsed time
!
    thrs = tm_curr / 3.6d+03

! store mesh statistics
!
    call store_mesh_stats(step, time)

! store integrals
!
    call store_integrals()

! write down the restart snapshot
!
    call write_restart_snapshot(thrs, nrun, status)

    if (check_status(status /= 0)) then
      if (master) then
        write(error_unit,"('[AMUN::program]: ', a)")                           &
                         "Writing restart snapshot failed!"
      end if
      go to 1000
    end if

! store data
!
    call write_snapshot()

! check if the time exceeds execution time limit
!
    if (thrs > trun) quit = 100

! print progress info to console, but not too often
!
    if (master) then
      if (time >= tmax .or. (tm_curr - tm_last) >= 1.0d+00) then

! calculate days, hours, seconds
!
        ec = int(tm_curr * (tmax - time) / max(1.0d-03, time - tbeg), kind = 4)
        es = max(0, min(863999999, ec))
        ed = es / 86400
        es = es - 86400 * ed
        eh = es /  3600
        es = es -  3600 * eh
        em = es /    60
        es = es -    60 * em

#ifdef INTEL
        write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' //  &
                ',1i2.2,"s",15x,a1,$)')                                        &
                        step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
#else /* INTEL */
        write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' //  &
                ',1i2.2,"s",15x,a1)',advance="no")                             &
                        step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
#endif /* INTEL */

! update the timestamp
!
        tm_last = tm_curr

      end if
    end if

! check if the loop should be kept going
!
    proceed = .not. check_status(quit /= 0)

  end do ! main loop

! add one empty line
!
  if (master) write(*,*)

! stop time accounting for the evolution
!
  call stop_timer(iev)

! write down the restart snapshot
!
  call write_restart_snapshot(1.0d+16, nrun, status)

  if (check_status(status /= 0)) then
    if (master) then
      write(error_unit,"('[AMUN::program]: ', a)")                             &
                       "Writing restart snapshot failed!"
    end if
  end if

! a label to go to if there are any problems, but since all modules have been
! initialized, we have to finalize them first
!
  1000 continue

! start time accounting for the termination
!
  call start_timer(itm)

! finalize modules
!
  2000 continue
  call finalize_integrals(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing INTEGRALS module!"
  end if
  2100 continue
  call finalize_evolution(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing EVOLUTION module!"
  end if
  2200 continue
  call finalize_schemes(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing SCHEMES module!"
  end if
  2300 continue
  call finalize_gravity(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing GRAVITY module!"
  end if
  2400 continue
  call finalize_shapes(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing SHAPES module!"
  end if
  2500 continue
  call finalize_mesh(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing MESH module!"
  end if
  2600 continue
  call finalize_refinement(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing REFINEMENT module!"
  end if
  2700 continue
  call finalize_boundaries(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing BOUNDARIES module!"
  end if
  2800 continue
  call finalize_domains(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing DOMAINS module!"
  end if
  2900 continue
  call finalize_problems(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing PROBLEMS module!"
  end if
  3000 continue
  call finalize_user_problem(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing USER_PROBLEM module!"
  end if
  3100 continue
  call finalize_sources(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing SOURCES module!"
  end if
  3200 continue
  call finalize_operators(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing OPERATORS module!"
  end if
  3300 continue
  call finalize_blocks(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing BLOCKS module!"
  end if
  3400 continue
  call finalize_coordinates(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing COORDINATES module!"
  end if
  3500 continue
  call finalize_interpolations(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing INTERPOLATIONS module!"
  end if
  3600 continue
  call finalize_equations(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing EQUATIONS module!"
  end if
  3700 continue
  call finalize_random()
  3800 continue
  call finalize_io(status)
  if (check_status(status /= 0) .and. master) then
    write(error_unit,"('[AMUN::program]: ', a)")                               &
                     "Problem finalizing IO module!"
  end if
  3900 continue
  call finalize_parameters()

! print execution times only if the initialization was fine
!
  if (initialization_succeeded) then

! stop time accounting for the termination
!
    call stop_timer(itm)

! get total time
!
    tm(1) = get_timer_total()

! get subtasks timers
!
    do i = 2, ntimers
      tm(i) = get_timer(i)
    end do

#ifdef MPI
! sum up timers from all processes
!
    call reduce_sum_real_array(ntimers, tm(:), status)
#endif /* MPI */

! print timings only on master processor
!
    if (master) then

! calculate the conversion factor
!
      tm_conv = 1.0d+02 / tm(1)

! print one empty line
!
      write (*,'(a)') ''

! print the execution times
!
      write (tmp,"(a)") "(2x,a32,1x,':',3x,1f16.3,' secs = ',f6.2,' %')"

      write (*,'(1x,a)') 'EXECUTION TIMINGS'
      do i = 2, ntimers
        if (timer_enabled(i)) then
          if (get_count(i) > 0) then
            write (*,tmp) timer_description(i), tm(i), tm_conv * tm(i)
          end if ! timer counter > 0
        end if ! enabled
      end do

! print the execution times
!
      write (tmp,"(a)") "(1x,a14,20x,':',3x,1f16.3,' secs = ',f6.2,' %')"
      write (*,tmp) 'TOTAL CPU TIME', tm(1)         , 1.0d+02
      write (*,tmp) 'TIME PER STEP ', tm(1) / nsteps, 1.0d+02 / nsteps
#ifdef MPI
      write (*,tmp) 'TIME PER CPU  ', tm(1) / nprocs, 1.0d+02 / nprocs
#endif /* MPI */

! get the execution time
!
      tm_exec = get_timer_total()

! convert the execution time to days, hours, minutes, and seconds and print it
!
      es = int(tm_exec)
      ec = int(1000 * (tm_exec - es))
      ed = es / 86400
      es = es - 86400 * ed
      eh = es /  3600
      es = es -  3600 * eh
      em = es /    60
      es = es -    60 * em
      write (tmp,"(a)") "(1x,a14,20x,':',3x,i14,'d'" //                        &
                                       ",i3.2,'h',i3.2,'m',i3.2,'.',i3.3,'s')"
      write (*,tmp) 'EXECUTION TIME', ed, eh, em, es, ec

    end if
  end if

! finalize module MPITOOLS
!
  4000 continue
  call finalize_mpitools(status)

! finalize module TIMERS
!
  call finalize_timers()

!-------------------------------------------------------------------------------
!
end program

#ifdef SIGNALS
!
!===============================================================================
!
! function TERMINATE:
! ------------------
!
!   Function sets variable iterm after receiving a signal.
!
!
!===============================================================================
!
integer(kind=4) function terminate(sig_num)

  implicit none

! input arguments
!
  integer(kind=4), intent(in) :: sig_num
  integer                     :: iterm

! common block
!
  common /termination/ iterm

!-------------------------------------------------------------------------------
!
#ifdef INTEL
  iterm     = sig_num
#else /* INTEL */
  iterm     = 15
#endif /* INTEL */
  terminate = 1

!-------------------------------------------------------------------------------
!
end
#endif /* SIGNALS */

!===============================================================================
!