!!******************************************************************************
!!
!!  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) 2017-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: PROBLEM
!!
!!  This module handles all operations related to the simulated problem, i.e.
!!  the initialization of problem parameters and related modules (equations,
!!  mesh, evolution, schemes, time statistics), the problem execution and
!!  termination.
!!
!!*****************************************************************************
!
module problem

  implicit none

! this flags determines if the problem is being resumed
!
  logical, save :: resumed = .false.

! the number of resumed run
!
  integer, save :: nrun    = 0

! the problem's permament parameters
!
  character(len=64)           , save :: name    = "none"
  character(len=32)           , save :: eqsys   = "hydrodynamic"
  character(len=32)           , save :: eos     = "adiabatic"
  character(len=32)           , save :: rngtype = "same"
  integer                     , save :: ncells  = 8
  integer                     , save :: maxlev  = 1
  integer     , dimension(3)  , save :: bdims   = 1
  real(kind=8), dimension(3,2), save :: dbnds

! the remaining problem's parameters
!
  integer                     , save :: nghosts = 2
  integer                     , save :: nmax    = huge(1)
  integer                     , save :: ndat    = 1
  real(kind=8)                , save :: tmax    = 0.0d+00
  real(kind=8)                , save :: trun    = 9.9d+99
  real(kind=8)                , save :: tsav    = 3.0d+01

! the minimum workspace size
!
  integer                     , save :: nwork   = 4194304 ! 32MiB

  private
  public :: initialize_problem, finalize_problem, execute_problem
  public :: print_problem_info
  public :: resumed, name, nrun, eqsys, eos, ncells, nghosts, maxlev
  public :: bdims, dbnds, rngtype, nmax, ndat, tmax, trun, tsav, nwork

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
  contains
!
!===============================================================================
!!
!!***  PUBLIC SUBROUTINES  *****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine INITIALIZE_PROBLEM:
! -----------------------------
!
!   Subroutine initializes this module and all related modules.
!
!   Arguments:
!
!     verbose - the verbose flag;
!     status  - the subroutine call status;
!
!===============================================================================
!
  subroutine initialize_problem(verbose, status)

    use blocks         , only : initialize_blocks
    use boundaries     , only : initialize_boundaries
    use coordinates    , only : initialize_coordinates, nn => bcells
    use equations      , only : initialize_equations, nf, nv
    use evolution      , only : initialize_evolution, registers
    use interpolations , only : initialize_interpolations
    use iso_fortran_env, only : error_unit
    use problems       , only : initialize_problems

    implicit none

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

    character(len=*), parameter :: loc = 'PROBLEM::initialize_problem()'

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

    call restore_parameters(status)
    if (status /=0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not restore the problem parameters!"
      return
    end if

    call initialize_equations(eqsys, eos, verbose, status)
    if (status /= 0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not initialize equations!"
      return
    end if

    call initialize_interpolations(nghosts, verbose, status)
    if (status /= 0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not initialize interpolations!"
      return
    end if
    call initialize_coordinates(ncells, nghosts, maxlev, bdims,                &
                                dbnds(1,1), dbnds(1,2),                        &
                                dbnds(2,1), dbnds(2,2),                        &
                                dbnds(3,1), dbnds(3,2), verbose, status)
    if (status /= 0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not initialize coordinates!"
      return
    end if

    call initialize_evolution(verbose, status)
    if (status /= 0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not initialize evolution!"
      return
    end if

    call initialize_blocks(nv, nf, nn, registers, status)
    if (status /= 0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not initialize blocks!"
      return
    end if

    call initialize_problems(name, nrun, verbose, status)
    if (status /= 0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not initialize problems!"
      return
    end if

    call initialize_boundaries(status)
    if (status /= 0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not initialize boundaries!"
      return
    end if

    nwork = 3 * NDIMS * nf * (ncells + 2 * nghosts)**NDIMS

!-------------------------------------------------------------------------------
!
  end subroutine initialize_problem
!
!===============================================================================
!
! subroutine FINALIZE_PROBLEM:
! ---------------------------
!
!   Subroutine frees all the memory used by this module and finalizes related
!   modules.
!
!   Arguments:
!
!     verbose - the verbose flag;
!     status  - the subroutine call status;
!
!===============================================================================
!
  subroutine finalize_problem(verbose, status)

    use blocks         , only : finalize_blocks
    use boundaries     , only : finalize_boundaries
    use coordinates    , only : finalize_coordinates
    use equations      , only : finalize_equations
    use evolution      , only : finalize_evolution
    use interpolations , only : finalize_interpolations
    use iso_fortran_env, only : error_unit
    use problems       , only : finalize_problems

    implicit none

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

    character(len=*), parameter :: loc = 'PROBLEM::finalize_problem()'

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

    call finalize_boundaries(status)
    if (status /= 0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not finalize boundaries!"
      return
    end if

    call finalize_problems(status)
    if (status /= 0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not finalize problems!"
      return
    end if

    call finalize_blocks(status)
    if (status /= 0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not finalize blocks!"
      return
    end if

    call finalize_evolution(verbose, status)
    if (status /= 0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not finalize evolution!"
      return
    end if

    call finalize_coordinates(status)
    if (status /= 0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not finalize coordinates!"
      return
    end if

    call finalize_interpolations(status)
    if (status /= 0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not finalize interpolations!"
      return
    end if

    call finalize_equations(status)
    if (status /= 0) then
      write(error_unit,"('[',a,']: ',a)") trim(loc),                           &
                       "Could not finalize equations!"
      return
    end if

!-------------------------------------------------------------------------------
!
  end subroutine finalize_problem
!
!===============================================================================
!
! subroutine PRINT_PROBLEM_INFO:
! -----------------------------
!
!   Subroutine prints all problem related information.
!
!   Arguments:
!
!     verbose - the verbose flag;
!
!===============================================================================
!
  subroutine print_problem_info(verbose)

    use boundaries     , only : print_boundaries
    use coordinates    , only : print_coordinates
    use equations      , only : print_equations
    use evolution      , only : print_evolution
    use helpers        , only : print_section, print_parameter
    use interpolations , only : print_interpolations

    implicit none

    logical, intent(in) :: verbose

!-------------------------------------------------------------------------------
!
    call print_section(verbose, "Problem")
    call print_parameter(verbose, "simulated problem", name)
    call print_parameter(verbose, "resumed", resumed, "yes")

    call print_equations(verbose)
    call print_interpolations(verbose)
    call print_coordinates(verbose)
    call print_evolution(verbose)
    call print_boundaries(verbose)

!-------------------------------------------------------------------------------
!
  end subroutine print_problem_info
!
!===============================================================================
!
! subroutine EXECUTE_PROBLEM:
! ---------------------------
!
!   Subroutine executes the problem by performing its time evolution.
!
!   Arguments:
!
!     status - the subroutine call status;
!
!===============================================================================
!
  subroutine execute_problem(status)

    implicit none

    integer, intent(out) :: status

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

!-------------------------------------------------------------------------------
!
  end subroutine execute_problem
!
!===============================================================================
!!
!!***  PRIVATE SUBROUTINES  ****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine RESTORE_PARAMETERS:
! -----------------------------
!
!   Subroutine reads the parameters from the parameters' file, then
!   restores them if the job is resumed.
!
!   Arguments:
!
!     status - the subroutine call status;
!
!===============================================================================
!
  subroutine restore_parameters(status)

    use io        , only : restart_from_snapshot, restart_snapshot_number
    use io        , only : read_snapshot_parameter
    use parameters, only : get_parameter

    implicit none

    integer, intent(out) :: status

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

    resumed = restart_from_snapshot()

    nrun = max(1, restart_snapshot_number() + 1)

    dbnds(:,1) = 0.0d+00
    dbnds(:,2) = 1.0d+00

    if (resumed) then
      call read_snapshot_parameter("problem", name, status)
      if (status /=0) return
      call read_snapshot_parameter("eqsys"  , eqsys, status)
      if (status /=0) return
      call read_snapshot_parameter("eos"    , eos  , status)
      if (status /=0) return
      call read_snapshot_parameter("ncells" , ncells, status)
      if (status /=0) return
      call read_snapshot_parameter("maxlev" , maxlev, status)
      if (status /=0) return
      call read_snapshot_parameter("xblocks", bdims(1), status)
      if (status /=0) return
      call read_snapshot_parameter("yblocks", bdims(2), status)
      if (status /=0) return
      call read_snapshot_parameter("zblocks", bdims(3), status)
      if (status /=0) return
      call read_snapshot_parameter("xmin", dbnds(1,1), status)
      if (status /=0) return
      call read_snapshot_parameter("xmax", dbnds(1,2), status)
      if (status /=0) return
      call read_snapshot_parameter("ymin", dbnds(2,1), status)
      if (status /=0) return
      call read_snapshot_parameter("ymax", dbnds(2,2), status)
      if (status /=0) return
#if NDIMS == 3
      call read_snapshot_parameter("zmin", dbnds(3,1), status)
      if (status /=0) return
      call read_snapshot_parameter("zmax", dbnds(3,2), status)
      if (status /=0) return
#endif /* NDIMS == 3 */
    else
      call get_parameter("problem", name)
      call get_parameter("equation_system"  , eqsys)
      call get_parameter("equation_of_state", eos)
      call get_parameter("ncells", ncells)
      call get_parameter("maxlev", maxlev)
      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", dbnds(1,1))
      call get_parameter("xmax", dbnds(1,2))
      call get_parameter("ymin", dbnds(2,1))
      call get_parameter("ymax", dbnds(2,2))
#if NDIMS == 3
      call get_parameter("zmin", dbnds(3,1))
      call get_parameter("zmax", dbnds(3,2))
#endif /* NDIMS == 3 */
    end if

    call get_parameter("rngtype", rngtype)
    call get_parameter("tmax"   , tmax)
    call get_parameter("trun"   , trun)
    call get_parameter("tsav"   , tsav)

    trun = trun - tsav / 6.0d+01

!-------------------------------------------------------------------------------
!
  end subroutine restore_parameters

!===============================================================================
!
end module problem