! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2011-2017 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2011-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2011-2012 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2011-2014 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2011 Jan Hueckelheim <j.hueckelheim@grs-sim.de>
! Copyright (c) 2012-2014 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016-2017 Raphael Haupt <raphael.haupt@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ***************************************************************************** !
!> In this module, the control structure for computing each time step is
!! established.
!!
!! Various different control structures can be defined and set by function
!! pointers.
!! The current main routine is recursive_multilevel which recursively computes
!! the new time step for all the schemes defined on all levels, starting from
!! the coarsest.
!!
?? include 'header/lbm_macros.inc'
module mus_control_module

  ! include treelm modules
  use env_module,         only: rk, labelLen
  use tem_timer_module,   only: tem_startTimer, tem_stopTimer
  use tem_general_module, only: tem_general_type
  use tem_debug_module,   only: main_debug, dbgUnit
  use tem_time_module,    only: tem_time_advance, tem_time_dump
  use tem_logging_module, only: logUnit
  use tem_varSys_module,  only: tem_varSys_type
  use tem_stencil_module, only: tem_stencilHeader_type

  ! include musubi modules
  use mus_pdf_module,                only: pdf_data_type, mus_swap_now_next
  use mus_param_module,              only: mus_param_type
  use mus_geom_module,               only: mus_geom_type
  use mus_bc_general_module,         only: set_boundary
  use mus_aux_module,                only: check_flow_status, &
    &                                      mus_update_relaxParams
  use mus_scheme_type_module,        only: mus_scheme_type
  use mus_scheme_header_module,      only: mus_scheme_header_type
  use mus_field_module,              only: mus_field_type
  use mus_source_var_module,         only: mus_source_type
  use mus_source_module,             only: mus_apply_sourceTerms
  use mus_debug_module,              only: dump_debug_Info
  use mus_IBM_module,                only: mus_inamuro_IBM, mus_buildBuffIBM
  use mus_timer_module,              only: mus_timerHandles, nStages
  use mus_physics_module,            only: mus_convertFac_type
  use mus_auxField_module,           only: mus_auxFieldVar_type, &
    &                                      mus_proc_calcAuxField
  use mus_derVarPos_module,          only: mus_derVarPos_type

  implicit none

  private

  public :: mus_control_type
  public :: mus_init_control


  !> Datatype containing mapping of control routines to function pointers
  type mus_control_type
    procedure( computation ), pointer :: do_computation => null()
  end type mus_control_type

  abstract interface
    !> Interface describes the main control routine which does computation
    !! set boundary and check flow status
    subroutine computation( me, scheme, geometry, params, iLevel)
      import :: mus_scheme_type, mus_geom_type, mus_param_type, &
        &       mus_control_type 
      !> self control type
      class( mus_control_type ) :: me
      !> container for the scheme
      type( mus_scheme_type ), intent(inout) :: scheme
      !> geometry infomation
      type( mus_geom_type ), intent(inout) :: geometry
      !> global parameters
      type( mus_param_type ), intent(inout) :: params
      !> Level counter variable
      integer, intent(in) :: iLevel
    end subroutine computation
  end interface

  integer, save :: iStage = 0
  logical, save :: running = .false.

contains

! ***************************************************************************** !
  !> This routines sets the function pointer to main control routine
  !!
  !! control routine is chosen based on the type of simulation.
  !! like single level, multi-level because multilevel requires
  !! recursive routine
  !!
  !! in the lua file you can now select
  !! `control_routine = '...'`
  !! where possible values are currently
  !!
  !! - `benchmark`: strongly reduced control routine with only single level,
  !!                no sources, etc.
  !! mainly for benchmarking
  !!
  !! - `multiLevel`: full multilevel, multiLevel routine
  !! - if nothing is given, the full multilevel, multiLevel routine is chosen
  subroutine mus_init_control( controlRoutine, me, minLevel, maxLevel )
    ! --------------------------------------------------------------------------
    character(len=labelLen), intent(in) :: controlRoutine
    !> contains function pointer to point control routine
    type( mus_control_type ), intent(out) :: me
    integer, intent(in) :: minLevel, maxLevel
    ! --------------------------------------------------------------------------

    ! Then select according to special need
    if ( trim( controlRoutine ) == 'benchmark' ) then
      me%do_computation => do_benchmark
      write(logUnit(5),"(A)") "Select benchmark control routine."
    elseif ( trim( controlRoutine ) == 'debug' ) then
      ! me%do_computation => do_recursive_debug
      write(logUnit(5),"(A)") "Select debug recursive control routine."
    else ! No special requirement from user, then select by single or multi levels
      if ( minLevel == maxLevel ) then
        me%do_computation => do_fast_singleLevel
        write(logUnit(5),"(A)") "Select fast single level control routine."
      else
        me%do_computation => do_recursive_multiLevel
        write(logUnit(5),"(A)") "Select recursive multi level control routine."
      end if
    end if

  end subroutine mus_init_control
! ***************************************************************************** !


! ***************************************************************************** !
  !> Main control routine: Update the time step for all levels.
  !! Main steps:
  !!   * do BC at iLevel
  !!   * do compute kernel at iLevel
  !!   * do do_computeFinerIntpAndExchange at iLevel if iLevel < maxLevel
  !!     * do recursive at iLevel+1
  !!     * intp My Coarser ghost (iLevel) from Finer (iLevel+1)
  !!     * do exchange bufferFromFiner at iLevel
  !!   * exchange buffer at iLevel
  !!   * exchange bufferFromCoarser at iLevel if iLevel > minLevel
  !!   * do do_intpCoarserAndExchange at iLevel if iLevel < maxLevel
  !!     * intp Finer Ghost (iLevel+1) from my coarser (iLevel)
  !!     * exchange bufferFromCoarser at iLevel+1
  !!
  recursive subroutine do_recursive_multiLevel( me, scheme, geometry, params, &
    &                                           iLevel)
    ! --------------------------------------------------------------------------
    !> self control type
    class( mus_control_type ) :: me
    !> container for the scheme
    type( mus_scheme_type ), intent(inout) :: scheme
    !> geometry infomation
    type( mus_geom_type ), intent(inout)      :: geometry
    !> global parameters
    type( mus_param_type ), intent(inout)  :: params
    !> the current level
    integer, intent(in) :: iLevel
    ! --------------------------------------------------------------------------
    integer :: now, next
    ! --------------------------------------------------------------------------

    ! swap double buffer index for current level
    call mus_swap_now_next( scheme%pdf( iLevel ) )
    now  = scheme%pdf(iLevel)%nNow
    next = scheme%pdf(iLevel)%nNext

    ! -------------------------------------------------------------------------
    ! set boundary for each field in current scheme
    call set_boundary( field       = scheme%field,                  &
      &                pdf         = scheme%pdf(iLevel),            &
      &                state       = scheme%state(iLevel)%val,      &
      &                levelDesc   = scheme%levelDesc(iLevel),      &
      &                tree        = geometry%tree,                 &
      &                iLevel      = iLevel,                        &
      &                nBCs        = geometry%boundary%nBCtypes,    &
      &                params      = params,                        &
      &                layout      = scheme%layout,                 &
      &                physics     = params%physics,                &
      &                varSys      = scheme%varSys,                 &
      &                mixture     = scheme%mixture,                &
      &                derVarPos   = scheme%derVarPos,              &
      &                globBC      = scheme%globBC                  )
    ! -------------------------------------------------------------------------

    ! --------------------------------------------------------------------------
    ! Compute auxField from pre-collision state for fluid and ghostFromCoarser
    ! and exchange them if turbulence is active
    call mus_calcAuxFieldAndExchange(                             &
      & auxField          = scheme%auxField(iLevel),              &
      & calcAuxField      = scheme%calcAuxField,                  &
      & state             = scheme%state(iLevel)%val(:, now),     &
      & pdfData           = scheme%pdf(iLevel),                   &
      & nFields           = scheme%nFields,                       &
      & field             = scheme%field(:),                      &
      & globSrc           = scheme%globSrc,                       &
      & stencil           = scheme%layout%fStencil,               &
      & varSys            = scheme%varSys,                        &
      & derVarPos         = scheme%derVarPos,                     &
      & general           = params%general,                       &
      & phyConvFac        = params%physics%fac(iLevel),           &
      & iLevel            = iLevel,                               &   
      & minLevel          = geometry%tree%global%minLevel,        &
      & schemeHeader      = scheme%header                         )   
    ! --------------------------------------------------------------------------

    ! --------------------------------------------------------------------------
    ! Update parameters, relaxation time .etc
    call mus_update_relaxParams( scheme, iLevel, params%general%simControl%now, &
      &                          params%physics, params%lattice )
    ! --------------------------------------------------------------------------

    ! -------------------------------------------------------------------------
    call start_stageTimer()
    ! write(logUnit(5), "(A,I0)") 'Compute on level: ', iLevel
    ! Compute current scheme of current level
    call tem_startTimer( timerHandle =  mus_timerHandles%compute(iLevel) )

!$omp parallel
    call scheme%compute(                                   &
      &  fieldProp = scheme%field(:)%fieldProp,            &
      &  inState   = scheme%state(iLevel)%val(:,Now),      &
      &  outState  = scheme%state(iLevel)%val(:,Next),     &
      &  auxField  = scheme%auxField(ilevel)%val(:),       &
      &  neigh     = scheme%pdf(iLevel)%neigh(:),          &
      &  nElems    = scheme%pdf(iLevel)%nSize,             &
      &  nSolve    = scheme%pdf(iLevel)%nElems_solve,      &
      &  level     = iLevel,                               &
      &  layout    = scheme%layout,                        &
      &  params    = params,                               &
      &  derVarPos = scheme%derVarPos,                     &
      &  varSys    = scheme%varSys                         )
!$omp end parallel

    call tem_stopTimer( timerHandle =  mus_timerHandles%compute(iLevel) )
    ! -------------------------------------------------------------------------

    ! --------------------------------------------------------------------------
    call mus_apply_sourceTerms( field      = scheme%field(:),                &
      &                         nFields    = scheme%nFields,                 &
      &                         globSrc    = scheme%globSrc,                 &
      &                         pdf        = scheme%pdf(iLevel),             &
      &                         varSys     = scheme%varSys,                  &
      &                         iLevel     = iLevel,                         &
      &                         time       = params%general%simControl%now,  &
      &                         state      = scheme%state(iLevel)%val,       &
      &                         auxField   = scheme%auxField(iLevel)%val(:), &
      &                         derVarPos  = scheme%derVarPos(:),            &
      &                         phyConvFac = params%physics%fac(iLevel),     &
      &                         tree       = geometry%tree                   )
    ! --------------------------------------------------------------------------

    ! when not on finest level, go to next level
    if( iLevel < geometry%tree%global%maxLevel ) then
      ! Do computation on finer level (L+1)
      ! Fill my coarser element (L) from finer (L+1)
      call do_computeFinerIntpAndExchange( me, scheme, geometry, params,   &
        &                                  iLevel                          )
    end if

    call stop_stageTimer()
! -------------------------------------------------------------------------
    ! Communicate the halo elements of each scheme on current level
    call tem_startTimer( timerHandle =  mus_timerHandles%comm(iLevel) )
    ! communicate halo elements
    call params%general%commPattern%exchange_real(            &
      &  send         = scheme%levelDesc(iLevel)%sendbuffer,  &
      &  recv         = scheme%levelDesc(iLevel)%recvbuffer,  &
      &  state        = scheme%state(iLevel)%val(:,Next),     &
      &  message_flag = iLevel,                               &
      &  comm         = params%general%proc%comm              )

    ! communicate turbulent viscosity, required for interpolation
    if (trim(scheme%header%kind) == 'fluid' .or. &
      & trim(scheme%header%kind) == 'fluid_incompressible') then
      if (scheme%field(1)%fieldProp%fluid%turbulence%active) then
        call params%general%commPattern%exchange_real(                &
          & recv         = scheme%field(1)%fieldProp%fluid%turbulence &
          &                      %dataOnLvl(iLevel)%recvbuffer,       &
          & send         = scheme%field(1)%fieldProp%fluid%turbulence &
          &                      %dataOnLvl(iLevel)%sendbuffer,       &
          & state        = scheme%field(1)%fieldProp%fluid%turbulence &
          &                      %dataOnLvl(iLevel)%visc(:),          &
          & message_flag = iLevel+100,                               &
          & comm         = params%general%proc%comm                  )
      end if
    end if

    call tem_stopTimer( timerHandle =  mus_timerHandles%comm(iLevel) )
! -------------------------------------------------------------------------

    ! communicate my finer ghost element from coarse level
    if ( iLevel > geometry%tree%global%minLevel ) then
      call tem_startTimer( timerHandle                                &
        &                  = mus_timerHandles%commFromCoarser(iLevel) )
      call params%general%commPattern%exchange_real(                  &
        & send    = scheme%levelDesc(iLevel)%sendbufferFromCoarser,   &
        & recv    = scheme%levelDesc(iLevel)%recvbufferFromCoarser,   &
        & state   = scheme%state(iLevel)%val(:, Next),                &
        & message_flag = iLevel,                                      &
        & comm    = params%general%proc%comm                          )

      ! communicate turbulent viscosity, required for interpolation
      if (trim(scheme%header%kind) == 'fluid' .or. &
        & trim(scheme%header%kind) == 'fluid_incompressible') then
        if (scheme%field(1)%fieldProp%fluid%turbulence%active) then
          call params%general%commPattern%exchange_real(                     &
            & recv         = scheme%field(1)%fieldProp%fluid%turbulence      &
            &                      %dataOnLvl(iLevel)%recvBufferFromCoarser, &
            & send         = scheme%field(1)%fieldProp%fluid%turbulence      &
            &                      %dataOnLvl(iLevel)%sendBufferFromCoarser, &
            & state        = scheme%field(1)%fieldProp%fluid%turbulence      &
            &                      %dataOnLvl(iLevel)%visc(:),               &
            & message_flag = iLevel+200,                                     &
            & comm         = params%general%proc%comm                        )
        end if 
      end if 

      call tem_stopTimer( timerHandle                                 &
        &                 = mus_timerHandles%commFromCoarser(iLevel)  )
    end if

    ! --------------------------------------------------------------------------
    ! Interpolate the ghost elements on the finer level(L+1) with data provided
    ! from current level(L).
    ! This only needs to be done once in two finer time steps
    if ( iLevel < geometry%tree%global%maxLevel ) then
      call do_intpCoarserAndExchange( scheme, params, iLevel )
    end if ! if not on finest level
    ! --------------------------------------------------------------------------

    ! update the time counters. MH: checked. Please dont move
    ! Increasing with the smallest time step (maxLevel)
    ! KM: time is advanced here since new time is required to update sources.
    if( iLevel == geometry%tree%global%maxLevel ) then
      call tem_time_advance( me     = params%general%simControl%now, &
        &                    sim_dt = params%physics%dtLvl( geometry &
        &                                   %tree%global%maxLevel )  )
    endif

    ! check flow status including perform checks, tracking, restart and geomInc
    ! only at minLevel to complete one multilevel cycle
    ! this check can also be done using mus_time_modulo(now, reqInterval)
    if( iLevel == geometry%tree%global%minLevel ) then
      ! check if some action has to be taken based on timeControl:
      ! tracking, global reduction operations, restart
      call check_flow_status( scheme   = scheme,                    &
        &                     general  = params%general,            &
        &           restart_triggered  = params%restart_triggered,  &
        &                    geometry  = geometry )
      iStage = 0
      running = .false.
    end if
  end subroutine do_recursive_multiLevel
! ***************************************************************************** !


! ***************************************************************************** !
  !> This routine does:
  !! 1. computation twice at finer level (iLevel+1),
  !! 2. interpolate my coarse ghost element (iLevel) from finer level (iLevel+1)
  !! 3. exchange the data of my coarse ghost elements between process
  recursive subroutine do_computeFinerIntpAndExchange( me, scheme, geometry,   &
    &                                        params, iLevel                    )
    ! --------------------------------------------------------------------------
    class( mus_control_type ) :: me !< contains control routine function pointer
    !> containers for the different schemes
    type( mus_scheme_type ), intent(inout), target  :: scheme
    !> geometry infomation
    type( mus_geom_type ), intent(inout)    :: geometry
    !> global parameters
    type( mus_param_type ), intent(inout)  :: params
    integer, intent(in) :: iLevel     !< Level counter variable
    ! --------------------------------------------------------------------------
    integer :: thisLevelNext, nextLevelNext
    !> diffusive or acoustic nesting loop
    integer :: iNestingLoop
    ! --------------------------------------------------------------------------

    ! Perform the number of nested time steps on the finer level L+1
    ! according to scaling type :
    !   diffusive: nNesting = 4
    !   acoustic:  nNesting = 2
    do iNestingLoop = 1, params%nNesting
      call me%do_computation( scheme, geometry, params, iLevel+1 )
    end do

    ! Interpolate my coarse ghost elements from the finer ones
    thisLevelNext = scheme%pdf(iLevel)%nNext
    nextLevelNext = scheme%pdf(iLevel+1)%nNext

    call start_stageTimer()
    ! write(logUnit(5), "(A,I0)") 'intpFromFiner on level: ', iLevel
    ! -------------------------------------------------------------------------
    call tem_startTimer( timerHandle =  mus_timerHandles%intpFromFiner(iLevel) )
    ! -------------------------------------------------------------------------
!$omp parallel
    call scheme%intp%fillMineFromFiner%do_intp(                             &
      &     fieldProp   = scheme%field(:)%fieldProp,                        &
      &     sState      = scheme%state(iLevel+1)%val(:,nextLevelNext),      &
      &     snSize      = scheme%pdf(iLevel+1)%nSize,                       &
      &     tnSize      = scheme%pdf(iLevel)%nSize,                         &
      &     tState      = scheme%state(iLevel)%val(:,thisLevelNext),        &
      &     tAuxField   = scheme%auxField(iLevel)%val(:),                   &
      &     tLevelDesc  = scheme%levelDesc(iLevel),                         &
      &     level       = iLevel,                                           &
      &     nTargets    = scheme%levelDesc(iLevel)%intpFromFiner%nVals,     &
      &     targetList  = scheme%levelDesc(iLevel)%intpFromFiner%val,       &
      &     layout      = scheme%layout,                                    &
      &     physics     = params%physics,                                   &
      &     varSys      = scheme%varSys,                                    &
      &     derVarPos   = scheme%derVarPos(:),                              &
      &     time        = params%general%simControl%now                     )
!$omp end parallel
    call tem_stopTimer( timerHandle = mus_timerHandles%intpFromFiner(iLevel) )
    ! -------------------------------------------------------------------------

    ! if ( main_debug%checkEachAlgorithmicStep ) then
    !   buffer=' after fill Mine FromFiner       '
    !   call dump_debug_info( buffer, scheme, params, iLevel, 1 ,pdf=scheme%pdf(iLevel))
    ! end if
    call stop_stageTimer()
    ! -------------------------------------------------------------------------
    ! Exchange the coarse ghost elements
    call tem_startTimer( timerHandle = mus_timerHandles%commFromFiner(iLevel) )
    call params%general%commPattern%exchange_real(                      &
      &    send    = scheme%levelDesc(iLevel)%sendbufferFromFiner,      &
      &    recv    = scheme%levelDesc(iLevel)%recvbufferFromFiner,      &
      &    state   = scheme%state(iLevel)%val(:,thisLevelNext),         &
      &    message_flag   = iLevel,                                     &
      &    comm    = params%general%proc%comm                  )

    ! exchange velocity halo fromFiner, required to compute velocity
    ! gradient
    call params%general%commPattern%exchange_real(                    &
      &  send         = scheme%auxField(iLevel)%sendBufferFromFiner,  &
      &  recv         = scheme%auxField(iLevel)%recvBufferFromFiner , &
      &  state        = scheme%auxField(iLevel)%val(:),               &
      &  message_flag = iLevel+100,                                   &
      &  comm         = params%general%proc%comm                      )

    if (trim(scheme%header%kind) == 'fluid' .or. &
      & trim(scheme%header%kind) == 'fluid_incompressible') then
      if (scheme%field(1)%fieldProp%fluid%turbulence%active) then
        ! communicate turbulent viscosity, required for interpolation
        call params%general%commPattern%exchange_real(                    &
          &  recv         = scheme%field(1)%fieldProp%fluid%turbulence    &
          &                       %dataOnLvl(iLevel)%recvBufferFromFiner, &
          &  send         = scheme%field(1)%fieldProp%fluid%turbulence    &
          &                       %dataOnLvl(iLevel)%sendBufferFromFiner, &
          &  state        = scheme%field(1)%fieldProp%fluid%turbulence    &
          &                       %dataOnLvl(iLevel)%visc(:),             &
          &  message_flag = iLevel+200,                                   &
          &  comm         = params%general%proc%comm                      )
      end if
    end if

    call tem_stopTimer( timerHandle = mus_timerHandles%commFromFiner(iLevel) )
    ! -------------------------------------------------------------------------

    ! if ( main_debug%checkEachAlgorithmicStep ) then
    !   buffer = ' after exch Mine FromFiner'
    !   call dump_debug_info( buffer, scheme, params, iLevel, 1 ,pdf=scheme%pdf(iLevel))
    ! end if

  end subroutine do_computeFinerIntpandExchange
! ***************************************************************************** !


! ***************************************************************************** !
  !> This routine utilizes fluid elements on my level (L) to fill finer
  !! ghost elements on next level (L+1).
  !! Then it exchanges the datas of finer ghost elements (L+1) between process.
  subroutine do_intpCoarserAndExchange( scheme, params, iLevel )
    ! --------------------------------------------------------------------------
    !> containers for the different schemes
    type( mus_scheme_type ), target, intent(inout)  :: scheme
    !> global parameters
    type( mus_param_type ), intent(in) :: params
    !> Level counter variable
    integer, intent(in) :: iLevel
    ! --------------------------------------------------------------------------
    integer :: thisLevelNext, nextLevelNext
    integer :: iOrder
    ! --------------------------------------------------------------------------

    thisLevelNext = scheme%pdf(iLevel)%nNext
    nextLevelNext = scheme%pdf(iLevel+1)%nNext

    call start_stageTimer()
    ! write(logUnit(5), "(A,I0)") 'intpFromCoarser on level: ', iLevel+1
    ! Fill finer ghost elements (L+1) by low order interpolation
    call tem_startTimer( timerHandle = mus_timerHandles%intpFromCoarser(iLevel+1) )

    ! Fill finer ghost elements (L+1) by starting from lower order 
    ! to higher order interpolation
!$omp parallel
    do iOrder = 0, scheme%intp%config%order
      call scheme%intp%fillFinerFromMe(iOrder)%do_intp(              &
        & fieldProp   = scheme%field(:)%fieldProp,                   &
        & sState      = scheme%state(iLevel)%val(:,thisLevelNext),   &
        & snSize      = scheme%pdf(iLevel)%nSize,                    &
        & tnSize      = scheme%pdf(iLevel+1)%nSize,                  &
        & tState      = scheme%state(iLevel+1)%val(:,nextLevelNext), &
        & tAuxField   = scheme%auxField(iLevel+1)%val(:),            &
        & tLevelDesc  = scheme%levelDesc(iLevel+1),                  &
        & level       = iLevel,                                      &
        & nTargets    = scheme%levelDesc(iLevel+1)                   &
        &                     %intpFromCoarser(iOrder)%nVals,        &
        & targetList  = scheme%levelDesc(iLevel+1)                   &
        &                     %intpFromCoarser(iOrder)%Val,          &
        & layout      = scheme%layout,                               &
        & physics     = params%physics,                              &
        & varSys      = scheme%varSys,                               &
        & derVarPos   = scheme%derVarPos(:),                         &
        & time        = params%general%simControl%now                )
     end do   
!$omp end parallel
      call tem_stopTimer( timerHandle                                  &
        &                 = mus_timerHandles%intpFromCoarser(iLevel+1) )

    ! Debug output
    ! if ( main_debug%checkEachAlgorithmicStep ) then
    !   buffer = ' after fillFinerFromMe'
    !   call dump_debug_info( buffer, scheme, params, iLevel, 1 ,pdf=scheme%pdf(iLevel))
    ! end if

    call stop_stageTimer()
    ! Exchange the fine ghost elements
    call tem_startTimer( timerHandle = mus_timerHandles%commFromCoarser(iLevel+1) )
    call params%general%commPattern%exchange_real(                   &
      &    send = scheme%levelDesc(iLevel+1)%SendBufferFromCoarser,  &
      &    recv = scheme%levelDesc(iLevel+1)%RecvBufferFromCoarser,  &
      &    state   = scheme%state(iLevel+1)%val(:,nextLevelNext),    &
      &    message_flag   = iLevel,                                  &
      &    comm    = params%general%proc%comm                        )

    ! exchange velocity halo fromCoarser, required to compute velocity
    ! gradient
    call params%general%commPattern%exchange_real(                        &
      &  send         = scheme%auxField(iLevel+1)%sendBufferFromCoarser,  &
      &  recv         = scheme%auxField(iLevel+1)%recvBufferFromCoarser , &
      &  state        = scheme%auxField(iLevel+1)%val(:),                 &
      &  message_flag = iLevel+100,                                       &
      &  comm         = params%general%proc%comm                          )

    if (trim(scheme%header%kind) == 'fluid' .or. &
      & trim(scheme%header%kind) == 'fluid_incompressible') then
      if (scheme%field(1)%fieldProp%fluid%turbulence%active) then
        ! communicate turbulent viscosity, required for interpolation
        call params%general%commPattern%exchange_real(                       &
          & recv         = scheme%field(1)%fieldProp%fluid%turbulence        &
          &                      %dataOnLvl(iLevel+1)%recvBufferFromCoarser, &
          & send         = scheme%field(1)%fieldProp%fluid%turbulence        &
          &                      %dataOnLvl(iLevel+1)%sendBufferFromCoarser, &
          & state        = scheme%field(1)%fieldProp%fluid%turbulence        &
          &                      %dataOnLvl(iLevel+1)%visc(:),               &
          & message_flag = iLevel+200,                                       &
          & comm         = params%general%proc%comm                          )
      end if
    end if
    call tem_stopTimer( timerHandle = mus_timerHandles%commFromCoarser(iLevel+1) )

    ! Debug output
    ! if ( main_debug%checkEachAlgorithmicStep ) then
    !   buffer = '  after exchFinerFromMe'
    !   call dump_debug_info( buffer, scheme, params, iLevel, 1 ,pdf=scheme%pdf(iLevel))
    ! end if

  end subroutine do_intpCoarserandExchange
! ***************************************************************************** !


! **************************************************************************** !
  !> Control routine for an optimized workflow with reduced functionality.
  !!
  !! No sources, no multilevel, no multiLevel.
  !! Use for benchmarking
  !!
  subroutine do_fast_singleLevel( me, scheme, geometry, params, iLevel )
    ! --------------------------------------------------------------------------
    !> self control type
    !! dummy variable in this routine, required by interface
    class( mus_control_type ) :: me
    !> container for the scheme
    type( mus_scheme_type ), intent(inout)  :: scheme
    !> geometry infomation
    type( mus_geom_type ), intent(inout)    :: geometry
    !> global parameters
    type( mus_param_type ), intent(inout)   :: params
    !> Level counter variable
    integer, intent(in) :: iLevel
    ! --------------------------------------------------------------------------
    integer :: loc_level
    integer :: now, next
    ! --------------------------------------------------------------------------

    loc_level = iLevel

!$omp single
    ! swap double buffer index for current level
    call mus_swap_now_next( scheme%pdf( iLevel ) )
    now  = scheme%pdf(iLevel)%nNow
    next = scheme%pdf(iLevel)%nNext

    ! --------------------------------------------------------------------------
    !set boundary for each field in current scheme
    call set_boundary( field       = scheme%field,                  &
      &                pdf         = scheme%pdf(iLevel),            &
      &                state       = scheme%state(iLevel)%val,      &
      &                levelDesc   = scheme%levelDesc(iLevel),      &
      &                tree        = geometry%tree,                 &
      &                iLevel      = iLevel,                        &
      &                nBCs        = geometry%boundary%nBCtypes,    &
      &                params      = params,                        &
      &                layout      = scheme%layout,                 &
      &                physics     = params%physics,                &
      &                varSys      = scheme%varSys,                 &
      &                mixture     = scheme%mixture,                &
      &                derVarPos   = scheme%derVarPos,              &
      &                globBC      = scheme%globBC                  )
    ! --------------------------------------------------------------------------

    ! --------------------------------------------------------------------------
    ! Compute auxField from pre-collision state for fluid and ghostFromCoarser
    ! and exchange them if turbulence is active
    call mus_calcAuxFieldAndExchange(                             &
      & auxField          = scheme%auxField(iLevel),              &
      & calcAuxField      = scheme%calcAuxField,                  &
      & state             = scheme%state(iLevel)%val(:, now),     &
      & pdfData           = scheme%pdf(iLevel),                   &
      & nFields           = scheme%nFields,                       &
      & field             = scheme%field(:),                      &
      & globSrc           = scheme%globSrc,                       &
      & stencil           = scheme%layout%fStencil,               &
      & varSys            = scheme%varSys,                        &
      & derVarPos         = scheme%derVarPos,                     &
      & general           = params%general,                       &
      & phyConvFac        = params%physics%fac(iLevel),           &
      & iLevel            = iLevel,                               &
      & minLevel          = geometry%tree%global%minLevel,        &
      & schemeHeader      = scheme%header                         )   
    ! --------------------------------------------------------------------------

    ! --------------------------------------------------------------------------
    ! Update parameters, relaxation time .etc
    call mus_update_relaxParams( scheme, iLevel, params%general%simControl%now, &
      &                          params%physics, params%lattice )
    ! --------------------------------------------------------------------------

    ! -------------------------------------------------------------------------
    ! Compute current scheme of current level
    call tem_startTimer( timerHandle =  mus_timerHandles%compute(iLevel) )
!$omp end single

!$omp parallel
    call scheme%compute(                                            &
      &           fieldProp = scheme%field(:)%fieldProp,            &
      &           inState   = scheme%state(iLevel)%val(:,Now),      &
      &           outState  = scheme%state(iLevel)%val(:,Next),     &
      &           auxField  = scheme%auxField(ilevel)%val(:),       &
      &           neigh     = scheme%pdf(iLevel)%neigh(:),          & 
      &           nElems    = scheme%pdf(iLevel)%nSize,             &
      &           nSolve    = scheme%pdf(iLevel)%nElems_solve  ,    &
      &           level     = iLevel,                               &
      &           layout    = scheme%layout,                        &
      &           params    = params,                               &
      &           derVarPos = scheme%derVarPos,                     &
      &           varSys    = scheme%varSys                         )
!$omp end parallel

!$omp single
    call tem_stopTimer( timerHandle =  mus_timerHandles%compute(iLevel) )
    ! --------------------------------------------------------------------------

    ! --------------------------------------------------------------------------
    call mus_apply_sourceTerms( field      = scheme%field(:),                &
      &                         nFields    = scheme%nFields,                 &
      &                         globSrc    = scheme%globSrc,                 &
      &                         pdf        = scheme%pdf(iLevel),             &
      &                         varSys     = scheme%varSys,                  &
      &                         iLevel     = iLevel,                         &
      &                         time       = params%general%simControl%now,  &
      &                         state      = scheme%state(iLevel)%val,       &
      &                         auxField   = scheme%auxField(iLevel)%val(:), &
      &                         derVarPos  = scheme%derVarPos(:),            &      
      &                         phyConvFac = params%physics%fac(iLevel),     &
      &                         tree       = geometry%tree                   )
    ! --------------------------------------------------------------------------

    ! -------------------------------------------------------------------------
    ! ... check if at least one of the IBMs is active
    if ( geometry%globIBM%nIBMs > 0 ) then
      call mus_buildBuffIBM(                              &
        &       me          = geometry%globIBM%IBM,       &
        &       commPattern = params%general%commPattern, &
        &       globTree    = geometry%tree,              &
        &       params      = params,                     &
        &       layout      = scheme%layout,              &
        &       levelDesc   = scheme%levelDesc(iLevel),   &
        &       iLevel      = loc_Level                   )
    end if
    ! -------------------------------------------------------------------------

    ! Communicate the halo elements of each scheme on current level
    call tem_startTimer( timerHandle =  mus_timerHandles%comm(iLevel) )
    call params%general%commPattern%exchange_real(             &
      &    send    = scheme%levelDesc(iLevel)%sendbuffer,      &
      &    recv    = scheme%levelDesc(iLevel)%recvbuffer,      &
      &    state   = scheme%state(iLevel)%val(:,Next),         &
      &    message_flag   = iLevel,                            &
      &    comm    = params%general%proc%comm                  )

    call tem_stopTimer( timerHandle =  mus_timerHandles%comm(iLevel) )

    ! Increasing with the smallest time step (maxLevel)
    ! KM: time is advanced here since new time is required to update sources.
    call tem_time_advance( me = params%general%simControl%now,   &
      &                    sim_dt = params%physics%dtLvl(iLevel ))

    ! update the immersed boundaries if available
    ! ... and over the schemes
    ! ... check if at least one of the IBMs is active
    if( geometry%globIBM%nIBMs > 0 )then
      call mus_inamuro_IBM(                                     &
        &      me          = geometry%globIBM%IBM,              &
        &      commPattern = params%general%commPattern,        &
        &      globTree    = geometry%tree,                     &
        &      general     = params%general,                    &
        &      pdf         = scheme%pdf(iLevel),                &
        &      layout      = scheme%layout,                     &
        &      levelDesc   = scheme%levelDesc(iLevel),          &
        &      globSys     = scheme%varSys,                     &
        &      stateVarMap = scheme%stateVarMap%varPos%val(:),  &
        &      convFac     = params%physics%fac(iLevel),        &
        &      iField      = 1,                                 &
        &      state       = scheme%state(iLevel)%val,          &
        &      iLevel      = loc_Level                          )
    end if

    ! check if some action has to be taken based on timeControl:
    ! tracking, global reduction operations, restart
    call check_flow_status( scheme   = scheme,                    &
      &                     general  = params%general,            &
      &           restart_triggered  = params%restart_triggered,  &
      &                    geometry  = geometry )
!$omp end single

  end subroutine do_fast_singleLevel
! ***************************************************************************** !

! ***************************************************************************** !
  subroutine do_benchmark( me, scheme, geometry, params, iLevel )
    ! --------------------------------------------------------------------------
    !> self control type
    !! dummy variable in this routine, required by interface
    class( mus_control_type ) :: me
    !> containers for the different schemes
    type( mus_scheme_type ), intent(inout)  :: scheme
    !> geometry infomation
    type( mus_geom_type ), intent(inout)    :: geometry
    !> global parameters
    type( mus_param_type ), intent(inout)   :: params
    !> Level counter variable
    integer, intent(in) :: iLevel
    ! --------------------------------------------------------------------------
    integer :: now, next
    ! --------------------------------------------------------------------------

!$omp single
    ! swap double buffer index for current level
    call mus_swap_now_next( scheme%pdf( iLevel ) )
    now  = scheme%pdf(iLevel)%nNow
    next = scheme%pdf(iLevel)%nNext

    ! --------------------------------------------------------------------------
    !set boundary for each field in current scheme
    call set_boundary( field       = scheme%field,                  &
      &                pdf         = scheme%pdf(iLevel),            &
      &                state       = scheme%state(iLevel)%val,      &
      &                levelDesc   = scheme%levelDesc(iLevel),      &
      &                tree        = geometry%tree,                 &
      &                iLevel      = iLevel,                        &
      &                nBCs        = geometry%boundary%nBCtypes,    &
      &                params      = params,                        &
      &                layout      = scheme%layout,                 &
      &                physics     = params%physics,                &
      &                varSys      = scheme%varSys,                 &
      &                mixture     = scheme%mixture,                &
      &                derVarPos   = scheme%derVarPos,              &
      &                globBC      = scheme%globBC                  )
    ! --------------------------------------------------------------------------

    ! --------------------------------------------------------------------------
    ! Compute auxField from pre-collision state for fluid and ghostFromCoarser
    ! and exchange them if turbulence is active
    call mus_calcAuxFieldAndExchange(                             &
      & auxField          = scheme%auxField(iLevel),              &
      & calcAuxField      = scheme%calcAuxField,                  &
      & state             = scheme%state(iLevel)%val(:, now),     &
      & pdfData           = scheme%pdf(iLevel),                   &
      & nFields           = scheme%nFields,                       &
      & field             = scheme%field(:),                      &
      & globSrc           = scheme%globSrc,                       &
      & stencil           = scheme%layout%fStencil,               &
      & varSys            = scheme%varSys,                        &
      & derVarPos         = scheme%derVarPos,                     &
      & general           = params%general,                       &
      & phyConvFac        = params%physics%fac(iLevel),           &
      & iLevel            = iLevel,                               &   
      & minLevel          = geometry%tree%global%minLevel,        &
      & schemeHeader      = scheme%header                         )   
    ! --------------------------------------------------------------------------

    ! --------------------------------------------------------------------------
    ! Update parameters, relaxation time .etc
    call mus_update_relaxParams( scheme, iLevel, params%general%simControl%now, &
      &                          params%physics, params%lattice )
    ! --------------------------------------------------------------------------

    ! --------------------------------------------------------------------------
    ! Compute current scheme of current level
    call tem_startTimer( timerHandle = mus_timerHandles%compute(iLevel) )
!$omp end single

!$omp parallel
    call scheme%compute(                                        &
      &           fieldProp = scheme%field(:)%fieldProp,        &
      &           inState   = scheme%state(iLevel)%val(:,Now),  &
      &           outState  = scheme%state(iLevel)%val(:,Next), &
      &           auxField  = scheme%auxField(ilevel)%val(:),   &
      &           neigh     = scheme%pdf(iLevel)%neigh(:),      &
      &           nElems    = scheme%pdf(iLevel)%nSize,         &
      &           nSolve    = scheme%pdf(iLevel)%nElems_solve,  &
      &           level     = iLevel,                           &
      &           layout    = scheme%layout,                    &
      &           params    = params,                           &
      &           derVarPos = scheme%derVarPos,                 &
      &           varSys    = scheme%varSys                     )
!$omp end parallel

!$omp single
    call tem_stopTimer( timerHandle = mus_timerHandles%compute(iLevel) )
    ! --------------------------------------------------------------------------

    ! --------------------------------------------------------------------------
    call mus_apply_sourceTerms( field      = scheme%field(:),                &
      &                         nFields    = scheme%nFields,                 &
      &                         globSrc    = scheme%globSrc,                 &
      &                         pdf        = scheme%pdf(iLevel),             &
      &                         varSys     = scheme%varSys,                  &
      &                         iLevel     = iLevel,                         &
      &                         time       = params%general%simControl%now,  &
      &                         state      = scheme%state(iLevel)%val,       &
      &                         auxField   = scheme%auxField(iLevel)%val(:), &
      &                         derVarPos  = scheme%derVarPos(:),            &      
      &                         phyConvFac = params%physics%fac(iLevel),     &
      &                         tree       = geometry%tree                   )
    ! --------------------------------------------------------------------------

    ! Communicate the halo elements of each scheme on current level
    call tem_startTimer( timerHandle = mus_timerHandles%comm(iLevel) )
    call params%general%commPattern%exchange_real(                  &
      &         send    = scheme%levelDesc(iLevel)%sendbuffer,      &
      &         recv    = scheme%levelDesc(iLevel)%recvbuffer,      &
      &         state   = scheme%state(iLevel)%val(:,Next),         &
      &         message_flag   = iLevel,                            &
      &         comm    = params%general%proc%comm                  )

    ! -------------------------------------------------------------------------

    call tem_stopTimer( timerHandle = mus_timerHandles%comm(iLevel) )

    ! Increasing with the smallest time step (maxLevel)
    call tem_time_advance( me = params%general%simControl%now,   &
      &                    sim_dt = params%physics%dtLvl(iLevel ))

    ! check if some action has to be taken based on timeControl:
    ! tracking, global reduction operations, restart
    call check_flow_status( scheme   = scheme,                    &
      &                     general  = params%general,            &
      &           restart_triggered  = params%restart_triggered,  &
      &                    geometry  = geometry )
!$omp end single
  end subroutine do_benchmark
! ***************************************************************************** !

! ***************************************************************************** !
  subroutine start_stageTimer( )

    if ( .not. running ) then
      iStage = mod( iStage, nStages ) + 1
      call tem_startTimer( timerHandle =  mus_timerHandles%stage(iStage) )
      running = .true.
    end if

  end subroutine start_stageTimer
! ***************************************************************************** !

! ***************************************************************************** !
  subroutine stop_stageTimer( )

    if ( running ) then
      call tem_stopTimer( timerHandle =  mus_timerHandles%stage(iStage) )
      running = .false.
    end if

  end subroutine stop_stageTimer
! ***************************************************************************** !


  ! ************************************************************************** !
  !> This routine compute auxField variable from pre-collision pdf and exchange
  !! halos
  subroutine mus_calcAuxFieldAndExchange(auxField, calcAuxField, state, &
    &  pdfData, nFields, field, globSrc, stencil, varSys, derVarPos,    &
    &  phyConvFac, general, iLevel, minLevel, schemeHeader)
    !---------------------------------------------------------------------------
    !> auxilary field array
    type(mus_auxFieldVar_type), intent(inout) :: auxField
    !> function pointer to calculate auxField
    procedure(mus_proc_calcAuxField), pointer, intent(in) :: calcAuxField
    !> state array
    real(kind=rk), intent(in) :: state(:)
    !> contains neigh array and nElems on current level
    type(pdf_data_type), intent(in) :: pdfData
    !> Number of fields
    integer, intent(in) :: nFields
    !> contains sources of all fields
    type(mus_field_type), intent(in)  :: field(nFields)
    !> global source
    type(mus_source_type), intent(in)  :: globSrc
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> contains auxField position of all fields in varSys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    !> physics conversion factors for this level
    type(mus_convertFac_type), intent(in) :: phyConvFac
    !> contains commPattern, MPI communicator and simControl
    type(tem_general_type), intent(in) :: general
    !> current level
    integer, intent(in) :: iLevel
    !> minlevel
    integer, intent(in) :: minLevel
    !> scheme header
    type(mus_scheme_header_type), intent(in) :: schemeHeader
    !---------------------------------------------------------------------------
    integer :: nSolve, iField, iSrc
    !---------------------------------------------------------------------------
    ! calculate auxField only for fluids and ghostfromcoarser (buffer ghost)
    ! ghostFromFiner are interpolated
    nSolve = pdfData%nElems_fluid + pdfData%nElems_ghostFromCoarser
    ! calculate velocity for fluid and ghost elements
    call calcAuxField(auxField  = auxField%val(:), &
      &               state     = state,           & 
      &               neigh     = pdfData%neigh,   &
      &               nSize     = pdfData%nSize,   &
      &               nSolve    = nSolve,          &
      &               iLevel    = iLevel,          &
      &               stencil   = stencil,         &
      &               varSys    = varSys,          &
      &               derVarPos = derVarPos        )

    ! update auxField with source term  
    ! Field add field source and then add globSrc
    do iField = 1, nFields
      do iSrc = 1, field(iField)%source%varDict%nVals
        call field(iField)%source%method(iSrc)%addSrcToAuxField( &
          & auxField   = auxField%val(:),                        &
          & iLevel     = iLevel,                                 &
          & time       = general%simControl%now,                 &
          & varSys     = varSys,                                 &
          & phyConvFac = phyConvFac,                             &
          & derVarPos  = derVarPos                               )
      end do
    end do

    ! apply global source
    do iSrc = 1, globSrc%varDict%nVals
      call globSrc%method(iSrc)%addSrcToAuxField( &
        & auxField   = auxField%val(:),           &
        & iLevel     = iLevel,                    &
        & time       = general%simControl%now,    &
        & varSys     = varSys,                    &
        & phyConvFac = phyConvFac,                &
        & derVarPos  = derVarPos                  )
    end do

    ! communicate velocity field. Requires for tubulence to compute ShearRate
    ! from velocity gradient.
    ! exchange velocity halo on current level
    call general%commpattern%exchange_real(   &
      &  send         = auxField%sendBuffer,  &
      &  recv         = auxField%recvBuffer , &
      &  state        = auxField%val(:),      &
      &  message_flag = iLevel+100,           &
      &  comm         = general%proc%comm     )

  !    ! communicate ghost halos fromFiner
  !    if (iLevel < maxLevel) then
  !       call pattern%exchange_real(                       &
  !        &  send         = auxField%sendBufferFromFiner,  &
  !        &  recv         = auxField%recvBufferFromFiner , &
  !        &  state        = auxField%val(:),               &
  !        &  message_flag = iLevel+101,                    &
  !        &  comm         = comm                           )
  !    end if

    ! communicate ghost halos from coarser
    if (iLevel > minLevel) then
       call general%commpattern%exchange_real(             &
        &  send         = auxField%sendBufferFromCoarser,  &
        &  recv         = auxField%recvBufferFromCoarser , &
        &  state        = auxField%val(:),                 &
        &  message_flag = iLevel+101,                      &
        &  comm         = general%proc%comm                )
    end if

  end subroutine mus_calcAuxFieldAndExchange
  ! ************************************************************************** !


end module mus_control_module
! ***************************************************************************** !
