! Copyright (c) 2022 Kannan Masilamani <kannan.masilamani@dlr.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.
! ***************************************************************************** !
!> author: Kannan Masilamani
!! Module containing subroutines for building MUSUBI specific source 
!! variables
!!
module mus_source_type_module
  ! include treelm modules
  use env_module,               only: rk, long_k, labelLen
  use tem_tools_module,         only: upper_to_lower
  use tem_aux_module,           only: tem_abort
  use tem_varMap_module,        only: tem_possible_variable_type, &
    &                                 init, append, truncate,     &
    &                                 tem_variable_loadMapping
  use tem_stencil_module,       only: tem_stencilHeader_type
  use tem_varSys_module,        only: tem_varSys_type
  use tem_stringKeyValuePair_module, only: grw_stringKeyValuePairArray_type, &
    &                                      init, append, truncate
  use tem_stencil_module,       only: tem_stencilHeader_type
  use tem_time_module,          only: tem_time_type
  use treelmesh_module,         only: treelmesh_type
  use tem_logging_module,       only: logUnit
  use tem_shape_module,         only: tem_shape_type, tem_load_shape
  use tem_subTree_type_module,  only: tem_subTree_type

  ! include musubi modules
  use mus_scheme_header_module, only: mus_scheme_header_type
  use mus_physics_module,       only: mus_convertFac_type
  use mus_derVarPos_module,     only: mus_derVarPos_type
  use mus_absorbLayer_module,   only: mus_absorbLayer_type,        &
    &                                 mus_absorbLayer_dynAvg_type, &
    &                                 mus_load_absorbLayer

  ! include aotus modules
  use aotus_module,     only: flu_State, aot_get_val, aoterr_Fatal, &
    &                         aoterr_NonExistent, aoterr_WrongType
  use aot_table_module, only: aot_table_open, aot_table_close, aot_get_val

  implicit none
  private

  public :: mus_source_type
  public :: mus_source_op_type
  public :: mus_turbChannelForce_type
  public :: mus_HRRCorrectionTerm_type

  public :: mus_create_poss_srcVar
  public :: mus_load_source_var
  public :: mus_source_cleanup

  public :: mus_applySrc_dummy
  public :: mus_addSrcToAuxField_dummy
  public :: mus_updateSrcVar_dummy
  
  ! ************************************************************************** !
  !> Stores correction term for HRR_bgk
  type mus_HRRCorrectionTerm_type
    !> phi
    real(kind=rk), allocatable :: phi(:)
  end type mus_HRRCorrectionTerm_type
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Contains information to compute average bulk velocity for dynamic_force.
  !! In turb_channel_force, the force term is adapted according to difference
  !! between reference bulk velocity and simulated plane average bulk velocity
  !! to avoid linear increase in simulated bulk velocits.
  !! For more information:
  !! https://www.wias-berlin.de/people/john/ELECTRONIC_PAPERS/JR07.IJNMF.pdf
  type mus_turbChannelForce_type
    !> tracking shapes
    type(tem_shape_type)  :: geom_utau(1) !1
    type(tem_shape_type)  :: geom_umean(1) !1

    !> sub-tree resulting from the elements within the tracking shape
    !! The sub-tree also holds the sub-communicator
    !! This data needs to be UPDATED after balance
    type(tem_subTree_type) :: subTree_utau !one for utau and umean
    type(tem_subTree_type) :: subTree_umean !one for utau and umean

    !> Reference bulk velocity in physical unit
    real(kind=rk) :: refVelBulk

    !> Characteristic height in physical unit
    real(kind=rk) :: refHeight

    !> Stream-wise direction to compute average velocity
    !! x=1, y=2, z=3
    integer :: flow_direction

    !> Dynamic Force term for turbulent channel in physical unit [m/s^2]
    !! F_dyn = (refVelBulk-avgVelXBulk) * refVelBulk / refHeight
    real(kind=rk) :: forceDyn(3)

    !> Global number of elements in defined shape 
    integer :: nElemsGlobal_utau
    integer :: nElemsGlobal_umean
  end type mus_turbChannelForce_type
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Contains source elements position in state array and idx to access 
  !! data variable refered in config file.
  !! This type is defined for each level
  type mus_source_elems_type
    !> Number of source elements on this level.
    !! nFluids + nGhosts
    integer :: nElems

    !> Position of elements in state array to apply source terms.
    !! Position in state array is same as position in total list
    !! Size: nElems
    integer, allocatable :: posInTotal(:)

    !> Index to access point data type to retrieve values from variable 
    !! refered for source variable
    integer, allocatable :: idx(:)

    !> source field value obtained from ST_fun data variable.
    !! Filled only for elements where source is active i.e. elements in
    !! posInTotal. 
    !! size: nElems*nComponents
    !! \todo KM: might be not neccessary
!KM!    real(kind=rk), allocatable :: val(:)

    !> Contains time average values of density and velocity for dynamic 
    !! absorblayer.
    !! \todo KM: 02042021 Introduce method_data c_ptr and point to 
    !! dynAvg for absorbLayer and change intent(inout) to intent(in) in
    !! proc_addSrcToAuxField.
    type(mus_absorbLayer_dynAvg_type) :: dynAvg
    
    ! source term for HRR_bgk
    type(mus_HRRCorrectionTerm_type) :: HRR_Corr
  end type mus_source_elems_type
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Description contains list of elements on which source is active and
  !! function pointer to update source
  type mus_source_op_type
    !> Position of this source term variable in the varSys
    integer :: srcTerm_varPos

    !> Position of data variable provided in config file in the varSys
    integer :: data_varPos

    !> Contains source elements position in state array/total list for
    !! each level
    type(mus_source_elems_type), allocatable :: elemLvl(:)

    !> Function to update state with source term
    procedure(proc_apply_source), pointer :: applySrc => null()

    !> to use source field array
    !KM!logical :: useSrcFieldVal = .false.

    !> name of the source variable
    character(len=labelLen) :: varname

    !> Function pointer to append source field to auxilary variable
    procedure(proc_addSrcToAuxField), pointer :: addSrcToAuxField => null()

    !> Function pointer to update source variable which are dependent on 
    !! auxField.
    procedure(proc_updateSourceVar), pointer :: updateSourceVar => null()

    !> Order of approximation for source like force, electric_field,
    !! charge_density.
    !! Order = 1, uses force term in BE approximated by forward Euler method
    !! Order = 2, uses force term in BE approximated by Trapezoidal method.
    !! For order 2, macroscopic source is also added to auxField.
    !! For fluid, fluid_incompressible, multispecies_liquid: source is added
    !! to momentum and for poisson: source is added to potential.
    !! Default: order = 2.
    integer :: order

    !> Additional config information for absorbLayer
    type(mus_absorbLayer_type) :: absLayer

    !> Contains information to compute average bulk velocity for dynamic_force
    type(mus_turbChannelForce_type) :: turbChanForce
  end type mus_source_op_type
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> Description of musubi source type
  type mus_source_type
    !> Contains source elements position in tree%treeID and
    !! function pointer to update source
    !! Size: varDict%nVals
    type(mus_source_op_type), allocatable :: method(:)

    !> Dictionary of source variable with 
    !! varDict%val()%key is the name of source variable and
    !! varDict%val()%value is the name of variable provided for the key
    type(grw_stringKeyValuePairArray_type) :: varDict
  end type mus_source_type
  ! *************************************************************************** !


  ! *************************************************************************** !
  abstract interface
    !> Abstract interface to update state with source terms
    subroutine proc_apply_source( fun, inState, outState, neigh, auxField, &
      & nPdfSize, iLevel, varSys, time, phyConvFac, derVarPos )
      import :: rk, mus_source_op_type, tem_varSys_type, tem_time_type, &
        &       mus_convertFac_type, mus_derVarPos_type

      !> Description of method to update source
      class(mus_source_op_type), intent(in) :: fun

      !> input  pdf vector
      !! \todo KM: instate is passed to compute auxField.
      !! Since auxField is precomputed from instate and passed to this routine.
      !! instate can be removed
      real(kind=rk), intent(in)  ::  inState(:)

      !> output pdf vector
      real(kind=rk), intent(inout) :: outState(:)

      !> connectivity Array corresponding to state vector
      integer,intent(in) :: neigh(:)

      !> auxField array
      real(kind=rk), intent(in) :: auxField(:)

      !> number of elements in state Array
      integer, intent(in) :: nPdfSize

      !> current level
      integer, intent(in) :: iLevel

      !> variable system
      type(tem_varSys_type), intent(in) :: varSys

      !> Point in time at which to evaluate the variable.
      type(tem_time_type), intent(in)  :: time

      !> Physics conversion factor for current level
      type(mus_convertFac_type), intent(in) :: phyConvFac

      !> position of derived quantities in varsys
      type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    end subroutine proc_apply_source

    !> Interface to add source to auxField vars in source_op_type for 
    !! all nSolve elements (nFluids+nGhostFromCoarser+nGhostFromFiner). 
    !! Halo elements are exchanged
    subroutine proc_addSrcToAuxField(fun, auxField, iLevel, time, varSys, &
      & phyConvFac, derVarPos)
      import :: rk, tem_varSys_type, mus_derVarPos_type, mus_convertFac_type, &
        &       tem_time_type, mus_source_op_type

      !> Description of method to update source
      class(mus_source_op_type), intent(inout) :: fun
      !> output auxField array
      real(kind=rk), intent(inout)          :: auxField(:)
      !> current level
      integer, intent(in)                   :: iLevel
      !> current timing information
      type(tem_time_type), intent(in)       :: time
      !> variable system definition
      type(tem_varSys_type), intent(in)     :: varSys
      !> Physics conversion factor for current level
      type(mus_convertFac_type), intent(in) :: phyConvFac
      !> position of derived quantities in varsys
      type(mus_derVarPos_type), intent(in)  :: derVarPos(:)
    end subroutine proc_addSrcToAuxField 

    !> Interface to update source variable which has dependency on auxField.
    !! Applied on all nSolve elements (nFluids+nGhostFromCoarser+nGhostFromFiner).
    !! Halo elements are exchanged
    !! This should be called after adding sorce term to state so that both
    !! auxField and apply_source uses same source value in one multilevel cycle.
    subroutine proc_updateSourceVar(fun, auxField, iLevel, varSys, phyConvFac, &
      & derVarPos)
      import :: rk, tem_varSys_type, mus_derVarPos_type, mus_convertFac_type, &
        &       treelmesh_type, mus_source_op_type

      !> Description of method to update source
      class(mus_source_op_type), intent(inout) :: fun
      !> input auxField array on current level 
      real(kind=rk), intent(in)          :: auxField(:)
      !> current level
      integer, intent(in) :: iLevel
      !> variable system definition
      type(tem_varSys_type), intent(in)     :: varSys
      !> Physics conversion factor on current level
      type(mus_convertFac_type), intent(in) :: phyConvFac
      !> position of derived quantities in varsys
      type(mus_derVarPos_type), intent(in)  :: derVarPos(:)
    end subroutine proc_updateSourceVar


  end interface
  ! *************************************************************************** !

contains

  ! *************************************************************************** !
  !> Routine initialize possible source variable depends on scheme kind
  subroutine mus_create_poss_srcVar(poss_srcVar, schemeHeader)
    ! --------------------------------------------------------------------------!
    !> possible source variables
    type(tem_possible_variable_type), intent(out) :: poss_srcVar

    !> Identifier of the scheme
    type(mus_scheme_header_type), intent(in) :: schemeHeader
    ! --------------------------------------------------------------------------!
    integer :: QQ
    ! --------------------------------------------------------------------------!
    ! used only for HRR Correction term
    select case (trim(schemeHeader%layout))
    case ('d2q9')
      QQ = 9
    case ('d3q19')
      QQ = 19
    case ('d3q27')
      QQ = 27
    case default
      ! Musubi will be aborted later, other layouts not supported!!!!
      QQ = -1
    end select

    write(logUnit(10),*) 'Creating possible source terms '
    call init(me = poss_srcVar, length = 2 )

    select case(trim(schemeHeader%kind))
    case ('fluid', 'fluid_incompressible')
      ! body force
      call append(me          = poss_srcVar, &
        &         varName     = 'force',     &
        &         nComponents = 3            )
      ! dynamic body force for turbulent channel
      call append(me          = poss_srcVar,                &
        &         varName     = 'turb_channel_force_accel', &
        &         nComponents = 3                           )
      ! absorb layer, STfun should return sponge_strength
      ! Target pressure and velocity are defined under absorb_layer_target
      call append(me          = poss_srcVar,    &
        &         varName     = 'absorb_layer', &
        &         nComponents = 1               )

      ! Absorn layer for inlet boundary with dynamic pressure and
      ! constant velocity.
      call append(me          = poss_srcVar,          &
        &         varName     = 'absorb_layer_inlet', &
        &         nComponents = 1                     )
      ! Absorn layer for outlet boundary with constant pressure and
      ! dynamic velocity.
      call append(me          = poss_srcVar,           &
        &         varName     = 'absorb_layer_outlet', &
        &         nComponents = 1                     )
      ! HRR correction term, might be usefull also for other schemes??
      call append(me          = poss_srcVar,      &
        &         varName     = 'hrr_correction', &
        &         nComponents = QQ                )

    case ('passive_scalar')
      call append(me          = poss_srcVar,       &
        &         varName     = 'equal_injection', &
        &         nComponents = 1                  )

      call append(me          = poss_srcVar, &
        &         varName     = 'injection', &
        &         nComponents = 1            )

    case ('nernst_planck')
      call append(me          = poss_srcVar,      &
        &         varName     = 'electric_field', &
        &         nComponents = 3                 )

    case ('multispecies_liquid')
      call append(me          = poss_srcVar,      &
        &         varName     = 'electric_field', &
        &         nComponents = 3                 )

      call append(me          = poss_srcVar, &
        &         varName     = 'force',     &
        &         nComponents = 3            )

    case ('poisson')
      call append(me          = poss_srcVar,      &
        &         varName     = 'charge_density', &
        &         nComponents = 1                 )

    case default
      write(logUnit(1),*) 'No possible source term defined for scheme kind:' &
        &                 //trim(schemeHeader%kind)
    end select
    call truncate(poss_srcVar)

  end subroutine mus_create_poss_srcVar
  ! *************************************************************************** !



  ! ***************************************************************************!
  !> Routine load musubi source terms for given key.
  !! key is glob_source or source
  subroutine mus_load_source_var(me, possVars, conf, parent, key, varSys)
    ! --------------------------------------------------------------------------!
    !> Source variable type to initialize
    type(mus_source_type), intent(out) :: me
    !> possible source variables
    type(tem_possible_variable_type), intent(in) :: possVars
    !> flu state
    type( flu_State ) :: conf
    !> parent handle if scheme table is defined
    integer, intent(in), optional :: parent
    !> key to load source term
    character(len=*), intent(in) :: key
    !> Global variable system
    type(tem_varSys_type), intent(inout) :: varSys
    ! --------------------------------------------------------------------------!
    integer :: iSrc, iError, srchandle
    character(len=labelLen) :: varname_order
    ! --------------------------------------------------------------------------!
    ! initialize growing array stringKeyValuePair
    call init( me = me%varDict )

    ! load the global source variables
    call tem_variable_loadMapping( possVars = possVars,   &
      &                            conf     = conf,       &
      &                            parent   = parent,     &
      &                            key      = key,        &
      &                            varDict  = me%varDict, &
      &                            varSys   = varSys      )

    ! truncate varDict
    call truncate( me = me%varDict )

    ! If source is defined then reopen source table to load order of
    ! approximatation to use for defined source variable
    if (me%varDict%nVals > 0) then
      ! allocate source method
      allocate(me%method(me%varDict%nVals))

      call aot_table_open( L       = conf,      &
        &                  parent  = parent,    &
        &                  thandle = srchandle, &
        &                  key     = key        )

      do iSrc = 1, me%varDict%nVals
         varname_order = trim(me%varDict%val(iSrc)%key)//'_order'
         call aot_get_val( L       = conf,                  &
           &               thandle = srchandle,             &
           &               key     = varname_order,         &
           &               val     = me%method(iSrc)%order, &
           &               default = 2,                     &
           &               ErrCode = iError                 )

         if (btest(iError, aoterr_Fatal)) then
           write(logUnit(1),*)'FATAL Error occured, while retrieving source' &
             &                //'order'
           if (btest(iError, aoterr_WrongType)) then
             write(logUnit(1),*)'Variable has wrong type!'
             write(logUnit(1),*)'STOPPING'
             call tem_abort()
           endif
         end if
         write(logUnit(1),'(A,I0)') ' Order of approximation for '          &
           &                        //trim(me%varDict%val(iSrc)%key)//': ', &
           &                        me%method(iSrc)%order

         ! Load additional information for absorblayer
         select case (trim(me%varDict%val(iSrc)%key)) 
         case ('turb_channel_force_accel')
           call load_turbChanForce( me     = me%method(iSrc)%turbChanForce, &
             &                      conf   = conf,                          &
             &                      key    = 'turb_channel_force_dynamic',  &
             &                      parent = srchandle                      )
         case ('absorb_layer')
           call mus_load_absorbLayer( me       = me%method(iSrc)%absLayer &
             &                                     %config,               &
             &                        conf     = conf,                    &
             &                        key      = 'absorb_layer_target',   &
             &                        parent   = srchandle,               &
             &                        loadPres = .true.,                  &
             &                        loadVel  = .true.                   )
         case ('absorb_layer_inlet')
           call mus_load_absorbLayer( me       = me%method(iSrc)%absLayer     &
             &                                     %config,                   &
             &                        conf     = conf,                        &
             &                        key      = 'absorb_layer_inlet_target', &
             &                        parent   = srchandle,                   &
             &                        loadPres = .false.,                     &
             &                        loadVel  = .true.                       )
         case ('absorb_layer_outlet')
           call mus_load_absorbLayer( me       = me%method(iSrc)%absLayer      &
             &                                     %config,                    &
             &                        conf     = conf,                         &
             &                        key      = 'absorb_layer_outlet_target', &
             &                        parent   = srchandle,                    &
             &                        loadPres = .true.,                       &
             &                        loadVel  = .false.                       )
         end select
      end do
    else
      allocate(me%method(0))
    end if  

  end subroutine mus_load_source_var
  ! ***************************************************************************!

  ! ************************************************************************** !
  !> Load shape, bulk velocity and height for turbulent channel force
  subroutine load_turbChanForce(me, conf, key, parent)
    ! -------------------------------------------------------------------------!
    !> Turbulent channel force 
    type(mus_turbChannelForce_type), intent(out) :: me
    !> flu state
    type( flu_State ) :: conf
    !> Table name to load target states
    character(len=*), intent(in) :: key
    !> parent source handle
    integer, intent(in) :: parent
    ! -------------------------------------------------------------------------!
    integer :: turbForce_handle, iError
    character(len=1) :: flow_direction
    ! -------------------------------------------------------------------------!
    ! -------------------------------------------------------------------------!
 
    call aot_table_open( L       = conf,             &
      &                  parent  = parent,           &
      &                  thandle = turbForce_handle, &
      &                  key     = trim(key)         )

    ! Load Shape, bulk velocity and height for turbulent channel force
    write(logUnit(1),*) ' * Turbulent channel force:'
    call aot_get_val( L       = conf,                &
      &               thandle = turbForce_handle,    &
      &               key     = 'ref_velocity_bulk', &
      &               val     = me%refVelBulk,       &
      &               ErrCode = iError               )
    if (btest(iError, aoterr_Fatal)) then
      write(logUnit(1), *) 'Error loading ref_velocity_bulk!'
      call tem_abort()
    end if
    write(logUnit(1),*) '    ref_velocity_bulk =', me%refVelBulk

    call aot_get_val( L       = conf,             &
      &               thandle = turbForce_handle, &
      &               key     = 'ref_height',     &
      &               val     = me%refHeight,     &
      &               ErrCode = iError            )
    if (btest(iError, aoterr_Fatal)) then
      write(logUnit(1), *) 'Error loading ref_height!'
      call tem_abort()
    end if

    call aot_get_val( L       = conf,              &
      &               thandle = turbForce_handle,  &
      &               key     = 'flow_direction', &
      &               val     = flow_direction,   &
      &               ErrCode = iError            )
    if (btest(iError, aoterr_Fatal)) then
      write(logUnit(1), *) 'Error loading flow_direction!'
      call tem_abort()
    end if

    write(logUnit(1),*) '    flow_direction ='//flow_direction
    select case(flow_direction)
    case ('x', 'X')
      me%flow_direction = 1
    case ('y', 'Y')
      me%flow_direction = 2
    case ('z', 'Z')
      me%flow_direction = 3
    case default
      call tem_abort('Unknown flow_direction')
    end select

    ! load geometry
    call tem_load_shape( conf    = conf,             &
      &                  parent  = turbForce_handle, &
      &                  key     = 'shape_utau',     & 
      &                  me      = me%geom_utau(1)   )

    call tem_load_shape( conf    = conf,             &
      &                  parent  = turbForce_handle, &
      &                  key     = 'shape_umean',    & 
      &                  me      = me%geom_umean(1)  )

    if (size(me%geom_utau(1)%canoND) /= 1) then
      write(logUnit(1),*) 'Error: Requires single shape for turb_channel_force'
      call tem_abort()
    end if

    if (size(me%geom_umean(1)%canoND) /= 1) then
      write(logUnit(1),*) 'Error: Requires single shape for turb_channel_force'
      call tem_abort()
    end if


  end subroutine load_turbChanForce
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> This routine act as a destructor for source type.
  !! The arrays allocated in mus_init_sourceTerms are destroyed here
  subroutine mus_source_cleanup(me)
    ! --------------------------------------------------------------------------
    type(mus_source_type), intent(inout)   :: me
    ! --------------------------------------------------------------------------
    integer :: iSrc
    ! --------------------------------------------------------------------------
    ! KM: DO NOT DESTROY VARDICT IN SOURCE_TYPE AS IT CONTAINS CONFIG INFO
    do iSrc = 1, me%varDict%nVals 
      deallocate(me%method(iSrc)%elemLvl)
    end do  
  end subroutine mus_source_cleanup
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Dummy routine for apply source
  subroutine mus_applySrc_dummy( fun, inState, outState, neigh, auxField, &
    & nPdfSize, iLevel, varSys, time, phyConvFac, derVarPos )
    ! --------------------------------------------------------------------------
    !> Description of method to update source
    class(mus_source_op_type), intent(in) :: fun
    !> input  pdf vector
    !! \todo KM: instate is passed to compute auxField.
    !! Since auxField is precomputed from instate and passed to this routine.
    !! instate can be removed
    real(kind=rk), intent(in)  ::  inState(:)
    !> output pdf vector
    real(kind=rk), intent(inout) :: outState(:)
    !> connectivity Array corresponding to state vector
    integer,intent(in) :: neigh(:)
    !> auxField array
    real(kind=rk), intent(in) :: auxField(:)
    !> number of elements in state Array
    integer, intent(in) :: nPdfSize
    !> current level
    integer, intent(in) :: iLevel
    !> variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> Point in time at which to evaluate the variable.
    type(tem_time_type), intent(in)  :: time
    !> Physics conversion factor for current level
    type(mus_convertFac_type), intent(in) :: phyConvFac
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! --------------------------------------------------------------------------
    ! abort only if source is active and function pointer is not assinged
    write(logUnit(6),*) 'WARNING: Dummy routine for applySrc'
  end subroutine mus_applySrc_dummy
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Dummy routine for add source to auxField
  subroutine mus_addSrcToAuxField_dummy(fun, auxField, iLevel, time, varSys, &
    &                                   phyConvFac, derVarPos)
    ! ------------------------------------------------------------------------ !
    !> Description of method to update source
    class(mus_source_op_type), intent(inout) :: fun
    !> output auxField array
    real(kind=rk), intent(inout)         :: auxField(:)
    !> current level
    integer, intent(in)                :: iLevel
    !> current timing information
    type(tem_time_type), intent(in)    :: time
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> Physics conversion factor for current level
    type(mus_convertFac_type), intent(in) :: phyConvFac
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ------------------------------------------------------------------------ !
    ! abort only if source is active and function pointer is not assinged
    write(logUnit(6),*) 'WARNING: Dummy routine for addSrcToAuxField'
  end subroutine mus_addSrcToAuxField_dummy
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Dummy routine for update source variable
  subroutine mus_updateSrcVar_dummy(fun, auxField, iLevel, varSys, phyConvFac, &
    &                               derVarPos)
    ! ------------------------------------------------------------------------ !
    !> Description of method to update source
    class(mus_source_op_type), intent(inout) :: fun
    !> input auxField array on current level 
    real(kind=rk), intent(in)          :: auxField(:)
    !> current level
    integer, intent(in) :: iLevel
    !> variable system definition
    type(tem_varSys_type), intent(in)     :: varSys
    !> Physics conversion factor on current level
    type(mus_convertFac_type), intent(in) :: phyConvFac
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in)  :: derVarPos(:)
    ! --------------------------------------------------------------------------
    write(logUnit(6),*) 'WARNING: Dummy routine for updateSourceVar'
  end subroutine mus_updateSrcVar_dummy
  ! ************************************************************************** !


end module mus_source_type_module
! **************************************************************************** !
