! Copyright (c) 2011-2014 Jens Zudrop <j.zudrop@grs-sim.de>
! Copyright (c) 2011 Metin Cakircali <m.cakircali@grs-sim.de>
! Copyright (c) 2011-2018 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2011 Gaurang Phadke <g.phadke@grs-sim.de>
! Copyright (c) 2011-2012 Laura Didinger <l.didinger@grs-sim.de>
! Copyright (c) 2012 Melven Zoellner <yameta@freenet.de>
! Copyright (c) 2012 Vyacheslav Korchagin <v.korchagin@grs-sim.de>
! Copyright (c) 2013-2017 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2013-2016 Nikhil Anand <nikhil.anand@uni-siegen.de>
! Copyright (c) 2013-2018, 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016-2017 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016-2017 Tobias Girresser <tobias.girresser@student.uni-siegen.de>
! Copyright (c) 2016-2017 Daniel Petró <daniel.petro@student.uni-siegen.de>
! Copyright (c) 2016-2019 Neda Ebrahimi Pour <neda.epour@uni-siegen.de>
! Copyright (c) 2017 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
!
! Permission to use, copy, modify, and distribute this software for any
! purpose with or without fee is hereby granted, provided that the above
! copyright notice and this permission notice appear in all copies.
!
! THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHORS DISCLAIM ALL WARRANTIES
! WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
! MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
! ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
! WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
! ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
! OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
! **************************************************************************** !

!> author: Jens Zudrop
!! Initialization of Ateles
!!
!! This module collects the various initialization tasks, that need to be done,
!! to set up the simulation.
module atl_initialize_module
  use, intrinsic :: iso_c_binding,   only: c_loc, C_F_POINTER
  use env_module,                    only: rk, pathLen
  use aotus_module,                  only: flu_state, close_config

  use treelmesh_module,              only: treelmesh_type
  use tem_aux_module,                only: tem_open_distconf
  use tem_logging_module,            only: logUnit
  use tem_variable_module,           only: tem_variable_type,                  &
    &                                      tem_variable_load
  use tem_derived_module,            only: tem_varSys_append_luaVar
  use tem_varSys_module,             only: tem_varSys_solverData_evalElem_type
  use tem_meshInfo_module,           only: tem_varSys_append_meshInfoVar
  use tem_stencil_module,            only: tem_stencilHeader_type,   &
    &                                      init,                     &
    &                                      d3q6_cxDir,               &
    &                                      tem_stencil_map_toTreelmDef
  use tem_precice_module,            only: precice_available, &
    &                                      tem_precice_create
  use tem_operation_var_module,      only: tem_opVar_reduction_transient_init

  use tem_spacetime_fun_module,      only: tem_create_subTree_of_st_funList
  use tem_convergence_module,        only: tem_init_convergence

  use ply_dynArray_project_module,   only: dyn_projectionArray_type, init
  use ply_poly_project_module,       only: ply_poly_project_type, &
    &                                      ply_fill_project_list
  use ply_sampled_tracking_module,   only: ply_sampled_track_init

  use atl_bc_header_module,          only: atl_store_bcVarPos
  use atl_bc_state_module,           only: atl_bc_state_set_fromPoint
  use atl_equation_module,           only: atl_equations_type
  use atl_equation_init_module,      only: atl_init_equation
  use atl_godunovFlux_module,        only: atl_flux_initGodunov
  use atl_container_module,          only: atl_init_elem_container, &
    &                                      atl_element_container_type
  use atl_space_basis,               only: atl_init_spacebasis
  use atl_initial_condition_module,  only: atl_load_initial_condition
  use atl_solver_param_module,       only: atl_solver_param_type
  use atl_restart_module,            only: atl_initRestart, &
    &                                      atl_readRestart
  use atl_source_module,             only: atl_initialize_sources, &
    &                                      atl_allocate_sourcedata
  use atl_source_types_module,       only: atl_init_source_type
  use atl_scheme_module,             only: atl_modg_scheme_prp,    &
    &                                      atl_modg_2d_scheme_prp, &
    &                                      atl_modg_1d_scheme_prp, &
    &                                      atl_scheme_type
  use atl_materialIni_module,        only: atl_init_materialParams
  use atl_materialPrp_module,        only: atl_init_material_type
  use atl_load_project_module,       only: atl_load_projection
  use atl_varSys_module,             only: atl_varSys_solverData_type, &
    &                                      atl_init_varSys_solverData, &
    &                                      atl_set_stFun_getElement,   &
    &                                      atl_create_fortranVar,      &
    &                                      atl_varSys_load_user
  use atl_operator_module,           only: atl_set_opVar_getElement
  use atl_boundary_module,           only: atl_fill_BCIndex
  use atl_penalization_module,       only: atl_init_penalization
  use atl_weights_module,            only: atl_initWeights
  use atl_precice_module,            only: atl_init_precice

  implicit none

  private

  public :: atl_initialize


contains


  ! ****************************************************************************
  !> Routine to initialize the complete Ateles solver.
  !!
  !! It reads the configuration from the given Lua script, and
  !! creates the necessary data structures for the computation.
  !!
  !! For a detailed description see [Initialization](../page/init.html)
  subroutine atl_initialize( params, container, equation, tree, varSys_data, &
    &                        timestepLimit, poly_proj_list                   )
    ! --------------------------------------------------------------------------
    type(atl_solver_param_type), intent(inout) :: params
    !> Container of mesh data.
    !! This includes information about the mesh parts, the scheme, but also the
    !! initial kernel states (i.e. initial conditions).
    type(atl_element_container_type), intent(inout), target :: container

    !> Description on the equation system to solve.
    type(atl_Equations_type), intent(inout), target :: equation

    !> Mesh data in treelmesh format.
    type(treelmesh_type), intent(inout) :: tree

    !> Data infomation of the variable System
    type(atl_varSys_solverData_type), intent(inout), target :: varSys_data

    !> timesteplimitation due to the precice timestep
    real(kind=rk), intent(inout), optional :: timestepLimit

    !> unique list for projection methods
    type(ply_poly_project_type), intent(inout), allocatable, target &
      & :: poly_proj_list(:)
    ! -----local variables------------------------------------------------------
    !> unique list for projection methods
    type(dyn_projectionArray_type) :: dynprojectArray
    !> position pointing if projection method used for ic
    integer, allocatable :: ic_prj_pos_list(:)
    ! The number of degrees of freedoms to write to/read from
    ! a tracking or restart file
    integer :: nDofsIO
    integer :: scheme_dim
    integer :: stateShape(3)
    type(tem_stencilHeader_type) :: stencil
    !> This structure takes the information about source terms from the equation
    !! system and stores them to evaluate the present source terms in the
    !! configuration file. After adding the active source terms to the variable
    !! system, this structure is not needed anymore.
    type(atl_init_source_type) :: initSource
    type(atl_init_material_type) :: initMaterial
    ! Name of the file, the configuration is read from.
    character(len=PathLen) :: configFile
    ! --------------------------------------------------------------------------
    type(tem_variable_type), allocatable :: userVars(:)
    type(tem_variable_type), allocatable :: atlVars(:)
    ! --------------------------------------------------------------------------
    integer, allocatable :: vError(:)
    type(tem_varSys_solverData_evalElem_type) :: solverData_evalElem
    type(flu_State) :: eqn_conf
    ! --------------------------------------------------------------------------

    solverData_evalElem%solver_bundle = c_loc(varSys_data)
    solverData_evalElem%stFun_setter => atl_set_stfun_getElement
    solverData_evalElem%opVar_setter => atl_set_opVar_getElement

    call atl_bc_state_set_fromPoint(solverData_evalElem)

    configFile = trim(params%general%solver%configFile)
    ! Create the appropriate D3Q6 stencil
    call init( me     = stencil,   &
      &        QQN    = 6,         &
      &        QQ     = 6,         &
      &        useAll = .false.,   &
      &        nDims  = 3,         &
      &        label  = 'd3q6',    &
      &        cxDir  = d3q6_cxDir )
    call tem_stencil_map_toTreelmDef(stencil)

    ! **************************************************************************
    ! EQUATION
    !
    ! Initialize the equation description
    call atl_init_equation( equation     = equation,                      &
      &                     conf         = params%general%solver%conf(1), &
      &                     varSys_Data  = varSys_data,                   &
      &                     initSource   = initSource,                    &
      &                     initMaterial = initMaterial                   )

    if ( (trim(equation%eq_kind) == 'unknown') ) then

      ! If there was no equation system provided in the configuration,
      ! try to get it from the restart file, when available.
      if ( params%general%restart%controller%readRestart ) then

        call tem_open_distconf(                                              &
          & L        = eqn_conf,                                             &
          & filename = trim(params%general%restart%controller%readFileName), &
          & proc     = params%general%proc                                   )

        call atl_init_equation( equation     = equation,            &
          &                     conf         = eqn_conf,            &
          &                     varSys_Data  = varSys_data,         &
          &                     initSource   = initSource,          &
          &                     initMaterial = initMaterial         )

        call close_config(eqn_conf)

        if ( (trim(equation%eq_kind) /= 'unknown') ) then
          write(logunit(1),*) 'Using the equation table from the restart!'
        end if

      end if

    else
      write(logunit(1),*) 'Using the equation table from the configuration.'
    end if

    if ( (trim(equation%eq_kind) == 'unknown') ) &
      &  write(logunit(1),*) 'WARNING: no equation table found !!'


    select case(trim(equation%eq_kind))
    case( 'euler', 'euler_2d', 'euler_1d',                      &
      &   'navier_stokes', 'navier_stokes_2d',                  &
      &   'filtered_navier_stokes', 'filtered_navier_stokes_2d' )
      call atl_flux_initGodunov(equation%euler%isen_coef)
    end select


    ! Init precice if pecice is used
    if (precice_available) then
      write(logUnit(1),*) 'Create precice environment'
      call tem_precice_create( rank      = params%general%proc%rank,     &
        &                      comm_size = params%general%proc%comm_size )
    end if

    ! **************************************************************************
    ! ELEMENT CONTAINER
    !
    ! Initialize the simulation container (including the scheme)
    call atl_init_elem_container(                                        &
      & me              = container,                                     &
      & equation        = equation,                                      &
      & conf            = params%general%solver%conf(1),                 &
      & tree            = tree,                                          &
      & time            = params%general%simControl%now,                 &
      & readRestart     = params%general%restart%controller%readRestart, &
      & proc            = params%general%proc,                           &
      & commPattern     = params%general%commPattern,                    &
      & boundary        = params%boundary                                )


    ! **************************************************************************
    ! PROJECTION
    !
    ! Now the general and all individuall projection tables are loaded and
    ! added to the DA for projection
    call atl_load_projection(                                      &
      & minLevel           = tree%global%minlevel,                 &
      & maxLevel           = tree%global%maxlevel,                 &
      & equation           = equation,                             &
      & poly_proj_pos      = container%cubes%poly_proj_pos,        &
      & dynprojectArray    = dynProjectArray,                      &
      & conf               = params%general%solver%conf(1),        &
      & scheme_list        = container%cubes%scheme_list,          &
      & source_projPos     = container%cubes%source%poly_proj_pos, &
      & boundary_list      = container%cubes%boundary_list,        &
      & boundary_list_stab = container%cubes%boundary_stab_list,   &
      & ic_pos_list        = ic_prj_pos_list,                      &
      & material_list      = container%cubes%material_list         )

    ! After loading all information needed for individuall projection method
    ! we can finally fill up the projection body
    !Figure out the dimension of the scheme to reduce overhead
    select case(container%cubes%scheme_list(tree%global%minLevel)%scheme)
      case (atl_modg_scheme_prp)
         scheme_dim = 3
      case (atl_modg_2d_scheme_prp)
         scheme_dim = 2
      case (atl_modg_1d_scheme_prp)
         scheme_dim = 1
    end select

    call ply_fill_project_list( proj_list        = poly_proj_list,  &
      &                         dyn_projectArray = dynProjectArray, &
      &                         scheme_dim       = scheme_dim       )


    call atl_alloc_temp( minLevel    = tree%global%minlevel,       &
      &                  maxLevel    = tree%global%maxlevel,       &
      &                  equation    = equation,                   &
      &                  proj_list   = poly_proj_list,             &
      &                  scheme_list = container%cubes%scheme_list )

    ! **************************************************************************
    ! VARIABLE SYSTEM
    !
    ! Load the user defined variables from the configuration file and store them
    ! in a temporary list.
    call tem_variable_load( me             = userVars,                      &
      &                     conf           = params%general%solver%conf(1), &
      &                     load_solvervar = atl_varSys_load_user,          &
      &                     vError         = vError                         )

    ! take the temporary list with the lua variables and try to add them to the
    ! variable system.
    call tem_varSys_append_luaVar(                      &
      & luaVar                   = userVars,            &
      & varSys                   = equation%varSys,     &
      & st_funList               = equation%stFunList,  &
      & solverData_evalElem      = solverData_evalElem  )

    ! Create some default ateles variables.
    ! For example: Lot of boundary condition variables requires constant zero as
    ! spacetime-function so create a variable "zero_const" and refer to this
    ! name in boundary condition variables
    call atl_create_fortranVar( me = atlVars )

    ! Append hard coded ateles spacetime function variable to the
    ! variable system.
    ! @TODO KM: Do no add atlVars into stFunList
    call tem_varSys_append_luaVar(                      &
      & luaVar                   = atlVars,             &
      & varSys                   = equation%varSys,     &
      & st_funList               = equation%stFunList,  &
      & solverData_evalElem      = solverData_evalElem  )

    ! Append variables with information on the mesh.
    call tem_varSys_append_meshInfoVar(equation%varSys)

    ! Create a subtree for each space time function based on the defined shape.
    ! This subtree can be used to create a list of elements that are affected
    ! by the space time function. This simplifies the evaluation of the
    ! particular space time function as there is no need to loop over all
    ! elements and check whether they are affected by the space time function.
    ! Instead, the space time function has a list of affected elements to loop
    ! over.
    call tem_create_subTree_of_st_funList( &
      & me      = equation%stFunList,      &
      & tree    = tree,                    &
      & bc_prop = params%boundary          )

    ! Initialize the global instance of atl_varSys_data
    call atl_init_varSys_solverdata(                      &
      & this           = varSys_data,                     &
      & equation       = equation,                        &
      & tree           = tree,                            &
      & schemeList     = container%cubes%scheme_list,     &
      & statedataList  = container%cubes%stateData_list,  &
      & kerneldataList = container%cubes%kernelData_list, &
      & meshList       = container%cubes%mesh_list,       &
      & polyProjList   = poly_proj_list,                  &
      & levelPointer   = container%cubes%levelPointer,    &
      & polyProjPos    = container%cubes%poly_proj_pos    )

    ! **************************************************************************
    ! BOUNDARY CODITIONS
    !
    ! Store position of a variable name defined for each boundary
    ! variable in tem_bc_state_type
    call atl_store_bcVarPos( bc     = container%cubes%bc, &
      &                      varSys = equation%varSys     )

    ! For BC we use the seuptIndices and getValOfIndex routines to access
    ! variables in the varSys. This need to be done after bc are initialized
    ! in init_cube_container and after the projection is
    ! initialized
    call atl_fill_BCIndex( tree           = tree,                          &
      &                    bc             = container%cubes%bc,            &
      &                    boundary_list  = container%cubes%boundary_list, &
      &                    nDim           = equation%nDimensions,          &
      &                    varSys         = equation%varSys,               &
      &                    poly_proj_list = poly_proj_list,                &
      &                    mesh_list      = container%cubes%mesh_list      )

    ! **************************************************************************
    ! SOURCES
    !
    call atl_initialize_sources(                        &
      & source         = container%cubes%source,        &
      & initSource     = initSource,                    &
      & conf           = params%general%solver%conf(1), &
      & equation       = equation,                      &
      & poly_proj_list = poly_proj_list,                &
      & mesh_list      = container%cubes%mesh_list,     &
      & tree           = tree,                          &
      & varSys_data    = varSys_data                    )

    stateShape = shape( container%cubes                                    &
      &                          %statedata_list(container%cubes%minLevel) &
      &                          %state                                    )

    call atl_allocate_sourceData( source      = container%cubes%source, &
      &                           nDofs       = stateShape(2),          &
      &                           nComponents = stateShape(3)           )

    ! **************************************************************************
    ! LB weights measurement setup
    !
    call atl_initWeights( tree )

    ! **************************************************************************
    ! MATERIAL
    !
    call atl_init_materialParams(                             &
      & equation             = equation,                      &
      & tree                 = tree,                          &
      & varSys_data          = varSys_data,                   &
      & material_list        = container%cubes%material_list, &
      & mesh_list            = container%cubes%mesh_list,     &
      & scheme_list          = container%cubes%scheme_list,   &
      & boundary_list        = container%cubes%boundary_list, &
      & time                 = params%general%simControl%now, &
      & conf                 = params%general%solver%conf(1), &
      & proc                 = params%general%proc,           &
      & commPattern          = params%general%commPattern,    &
      & poly_proj_list       = poly_proj_list,                &
      & levelpointer         = container%cubes%levelpointer,  &
      & initMaterial         = initMaterial                   )

    ! **************************************************************************
    ! PENALIZATION
    !
    call atl_init_penalization(                                        &
      & tree                  = tree,                                  &
      & penalizationdata_list = container%cubes%penalizationdata_list, &
      & scheme_list           = container%cubes%scheme_list,           &
      & equation              = equation,                              &
      & mesh_list             = container%cubes%mesh_list              )
!PV!      ,             &
!PV!      & penalization_list     = container%cubes%penalization_list,     &
!PV!      & conf                  = params%general%solver%conf(1),         &
!PV!      & levelpointer          = container%cubes%levelpointer           )


    ! The number of degrees of freedoms to write out in tracking or
    ! restart
    nDofsIO = maxval(container%cubes%scheme_list(:)%nDofs)

    ! Fill the spatial polynomial basis, that is needed for integrations
    ! and point evaluations.
    call atl_init_spacebasis( scheme_list = container%cubes%scheme_list, &
      &                       minlevel    = tree%global%minlevel,        &
      &                       maxlevel    = tree%global%maxlevel         )

    ! **************************************************************************
    ! RESTART
    !
    call atl_initRestart( restart    = params%general%restart, &
      &                   equation   = equation,               &
      &                   solver     = params%general%solver,  &
      &                   mesh       = container%cubes,        &
      &                   tree       = tree,                   &
      &                   nDofsIO    = nDofsIO,                &
      &                   scheme_dim = scheme_dim              )

    ! **************************************************************************
    ! INITIAL CONDITION
    !
    if (params%general%restart%controller%readRestart) then
      ! Actually read the restart data as initial condition
      call atl_readRestart( mesh           = container%cubes,        &
        &                   restart        = params%general%restart, &
        &                   tree           = tree,                   &
        &                   proc           = params%general%proc,    &
        &                   equation       = equation                )
    else
      ! No restart, get the initial condition from the configuration
      call atl_load_initial_condition(                    &
        & cube_container = container%cubes,               &
        & equation       = equation,                      &
        & conf           = params%general%solver%conf(1), &
        & tree           = tree,                          &
        & luaFile        = trim(configFile),              &
        & prj_pos        = ic_prj_pos_list,               &
        & poly_proj_list = poly_proj_list                 )
    end if


    ! **************************************************************************
    ! PRECICE
    !
    !> COUPLING USING EQUIDISTANT AND NON-EQUIDISTANT POINTS
    !! Reminder: For the equidistant and non-equidistant points the call for
    !! the setup_indices is included in the
    !! atl_write_equiPoints/atl_writenonequiPoints routine.
    if (precice_available) then
      call atl_init_precice(                         &
        & params        = params,                    &
        & equation      = equation,                  &
        & tree          = tree,                      &
        & mesh_list     = container%cubes%mesh_list, &
        & timestepLimit = timestepLimit              )
    end if

    ! **************************************************************************
    ! TRACKING
    !
    call ply_sampled_track_init( me      = params%plySampleTrack, &
      &                          mesh    = tree,                  &
      &                          solver  = params%general%solver, &
      &                          varSys  = equation%varSys,       &
      &                          bc      = params%boundary,       &
      &                          stencil = stencil,               &
      &                          proc    = params%general%proc,   &
      &                          ndofs   = ndofsIO,               &
      &                          ndims   = equation%nDimensions   )

    ! **************************************************************************
    ! CONVERGENCE
    !
    if ( params%general%simControl%abortCriteria%steady_state ) then
      call tem_init_convergence( me       = params%general%simControl    &
        &                                                 %abortCriteria &
        &                                                 %convergence,  &
        &                        tree     = tree,                        &
        &                        bc_prop  = params%boundary,             &
        &                        globProc = params%general%proc,         &
        &                        varSys   = equation%varSys,             &
        &                        ndofs    = ndofsIO                      )
    end if

    ! Initialize reduction_transient operation variables
    call tem_opVar_reduction_transient_init(      &
      & varSys         = equation%varSys,         &
      & nElems         = tree%nElems,             &
      & redTransVarMap = equation%redTransVarmap, &
      & nDofs          = nDofsIO                  )

  end subroutine atl_initialize
  ! ****************************************************************************


  ! ****************************************************************************
  !> Allocate temporary arrays.
  subroutine atl_alloc_temp( minLevel, maxLevel, proj_list, &
    &                        scheme_list, equation   )
  ! ----------------------------------------------------------------------------
    integer, intent(in) :: minLevel
    integer, intent(in) :: maxLevel
    type(ply_poly_project_type), intent(inout), allocatable :: proj_list(:)
    type(atl_equations_type), intent(in) :: equation
    type(atl_scheme_type), intent(inout) :: scheme_list(minlevel:maxLevel)
    ! --------------------------------------------------------------------------

    ! allocate the temporary arrays required for the calculation of fluxes
    ! Note : They are preallocated so that the allocation is not needed at the
    ! time they are used ( in the thread parallel region)
    select case(scheme_list(minLevel)%scheme)
      case (atl_modg_scheme_prp)

        allocate( scheme_list(minlevel)%temp_over(      &
          & maxVal(proj_list(:)%body_3d%oversamp_dofs), &
          & equation%temp%nScal,                        &
          & equation%temp%oversamp )                    )
        allocate( scheme_list(minlevel)%temp_nodal(   &
          & maxval(proj_list(:)%body_3d%nquadpoints), &
          & equation%temp%nScal,                      &
          & equation%temp%nodal )                     )
        allocate( scheme_list(minlevel)%temp_modal( &
          & maxVal(proj_list(:)%body_3d%ndofs),     &
          & equation%temp%nScal,                    &
          & equation%temp%modal )                   )

        scheme_list(minlevel)%temp_over = 0.0_rk
        scheme_list(minlevel)%temp_nodal = 0.0_rk
        scheme_list(minlevel)%temp_modal = 0.0_rk

      case (atl_modg_2d_scheme_prp)

        allocate( scheme_list(minlevel)%temp_over(      &
          & maxVal(proj_list(:)%body_2d%oversamp_dofs), &
          & equation%temp%nScal,                        &
          & equation%temp%oversamp )                    )
        allocate( scheme_list(minlevel)%temp_nodal(    &
          & maxval(proj_list(:)%body_2d%nquadpoints),  &
          & equation%temp%nScal,                       &
          & equation%temp%nodal )                      )
        allocate( scheme_list(minlevel)%temp_modal( &
          & maxVal(proj_list(:)%body_2d%ndofs),     &
          & equation%temp%nScal,                    &
          & equation%temp%modal )                   )

        scheme_list(minlevel)%temp_over = 0.0_rk
        scheme_list(minlevel)%temp_nodal = 0.0_rk
        scheme_list(minlevel)%temp_modal = 0.0_rk

      case (atl_modg_1d_scheme_prp)

        allocate( scheme_list(minlevel)%temp_over(       &
          & maxVal(proj_list(:)%body_1d%oversamp_dofs ), &
          & equation%temp%nScal,                         &
          & equation%temp%oversamp )                     )
        allocate( scheme_list(minlevel)%temp_nodal(   &
          & maxval(proj_list(:)%body_1d%nquadpoints), &
          & equation%temp%nScal,                      &
          & equation%temp%nodal )                     )
        allocate( scheme_list(minlevel)%temp_modal( &
          & maxVal(proj_list(:)%body_1d%ndofs),     &
          & equation%temp%nScal,                    &
          & equation%temp%modal )                   )

        scheme_list(minlevel)%temp_over = 0.0_rk
        scheme_list(minlevel)%temp_nodal = 0.0_rk
        scheme_list(minlevel)%temp_modal = 0.0_rk

    end select

  end subroutine atl_alloc_temp
  ! ****************************************************************************

end module atl_initialize_module
