tem_spongelayer_module.f90 Source File


This file depends on

sourcefile~~tem_spongelayer_module.f90~~EfferentGraph sourcefile~tem_spongelayer_module.f90 tem_spongelayer_module.f90 sourcefile~tem_geometry_module.f90 tem_geometry_module.f90 sourcefile~tem_spongelayer_module.f90->sourcefile~tem_geometry_module.f90 sourcefile~tem_logging_module.f90 tem_logging_module.f90 sourcefile~tem_spongelayer_module.f90->sourcefile~tem_logging_module.f90 sourcefile~tem_aux_module.f90 tem_aux_module.f90 sourcefile~tem_spongelayer_module.f90->sourcefile~tem_aux_module.f90 sourcefile~env_module.f90 env_module.f90 sourcefile~tem_spongelayer_module.f90->sourcefile~env_module.f90 sourcefile~treelmesh_module.f90 treelmesh_module.f90 sourcefile~tem_spongelayer_module.f90->sourcefile~treelmesh_module.f90 sourcefile~tem_param_module.f90 tem_param_module.f90 sourcefile~tem_spongelayer_module.f90->sourcefile~tem_param_module.f90 sourcefile~tem_geometry_module.f90->sourcefile~tem_logging_module.f90 sourcefile~tem_geometry_module.f90->sourcefile~env_module.f90 sourcefile~tem_geometry_module.f90->sourcefile~treelmesh_module.f90 sourcefile~tem_geometry_module.f90->sourcefile~tem_param_module.f90 sourcefile~tem_property_module.f90 tem_property_module.f90 sourcefile~tem_geometry_module.f90->sourcefile~tem_property_module.f90 sourcefile~tem_subtree_type_module.f90 tem_subTree_type_module.f90 sourcefile~tem_geometry_module.f90->sourcefile~tem_subtree_type_module.f90 sourcefile~tem_tools_module.f90 tem_tools_module.f90 sourcefile~tem_geometry_module.f90->sourcefile~tem_tools_module.f90 sourcefile~tem_float_module.f90 tem_float_module.f90 sourcefile~tem_geometry_module.f90->sourcefile~tem_float_module.f90 sourcefile~tem_topology_module.f90 tem_topology_module.f90 sourcefile~tem_geometry_module.f90->sourcefile~tem_topology_module.f90 sourcefile~tem_debug_module.f90 tem_debug_module.f90 sourcefile~tem_geometry_module.f90->sourcefile~tem_debug_module.f90 sourcefile~tem_logging_module.f90->sourcefile~env_module.f90 sourcefile~tem_aux_module.f90->sourcefile~tem_logging_module.f90 sourcefile~tem_aux_module.f90->sourcefile~env_module.f90 sourcefile~tem_lua_requires_module.f90 tem_lua_requires_module.f90 sourcefile~tem_aux_module.f90->sourcefile~tem_lua_requires_module.f90 sourcefile~tem_aux_module.f90->sourcefile~tem_tools_module.f90 sourcefile~tem_comm_env_module.f90 tem_comm_env_module.f90 sourcefile~tem_aux_module.f90->sourcefile~tem_comm_env_module.f90 sourcefile~tem_revision_module.f90 tem_revision_module.f90 sourcefile~tem_aux_module.f90->sourcefile~tem_revision_module.f90 sourcefile~treelmesh_module.f90->sourcefile~tem_logging_module.f90 sourcefile~treelmesh_module.f90->sourcefile~tem_aux_module.f90 sourcefile~treelmesh_module.f90->sourcefile~env_module.f90 sourcefile~treelmesh_module.f90->sourcefile~tem_property_module.f90 sourcefile~treelmesh_module.f90->sourcefile~tem_tools_module.f90 sourcefile~tem_global_module.f90 tem_global_module.f90 sourcefile~treelmesh_module.f90->sourcefile~tem_global_module.f90 sourcefile~tem_sparta_module.f90 tem_sparta_module.f90 sourcefile~treelmesh_module.f90->sourcefile~tem_sparta_module.f90 sourcefile~treelmesh_module.f90->sourcefile~tem_topology_module.f90 sourcefile~tem_param_module.f90->sourcefile~env_module.f90 sourcefile~tem_property_module.f90->sourcefile~env_module.f90 sourcefile~tem_prophead_module.f90 tem_prophead_module.f90 sourcefile~tem_property_module.f90->sourcefile~tem_prophead_module.f90 sourcefile~tem_lua_requires_module.f90->sourcefile~env_module.f90 sourcefile~tem_subtree_type_module.f90->sourcefile~tem_logging_module.f90 sourcefile~tem_subtree_type_module.f90->sourcefile~tem_aux_module.f90 sourcefile~tem_subtree_type_module.f90->sourcefile~env_module.f90 sourcefile~tem_subtree_type_module.f90->sourcefile~treelmesh_module.f90 sourcefile~tem_subtree_type_module.f90->sourcefile~tem_property_module.f90 sourcefile~tem_subtree_type_module.f90->sourcefile~tem_global_module.f90 sourcefile~tem_tools_module.f90->sourcefile~env_module.f90 sourcefile~tem_global_module.f90->sourcefile~tem_logging_module.f90 sourcefile~tem_global_module.f90->sourcefile~tem_aux_module.f90 sourcefile~tem_global_module.f90->sourcefile~env_module.f90 sourcefile~tem_global_module.f90->sourcefile~tem_prophead_module.f90 sourcefile~tem_sparta_module.f90->sourcefile~tem_logging_module.f90 sourcefile~tem_sparta_module.f90->sourcefile~tem_aux_module.f90 sourcefile~tem_sparta_module.f90->sourcefile~env_module.f90 sourcefile~tem_sparta_module.f90->sourcefile~tem_float_module.f90 sourcefile~tem_float_module.f90->sourcefile~env_module.f90 sourcefile~tem_topology_module.f90->sourcefile~env_module.f90 sourcefile~tem_debug_module.f90->sourcefile~tem_logging_module.f90 sourcefile~tem_debug_module.f90->sourcefile~env_module.f90 sourcefile~tem_debug_module.f90->sourcefile~tem_tools_module.f90 sourcefile~tem_prophead_module.f90->sourcefile~env_module.f90

Files dependent on this one

sourcefile~~tem_spongelayer_module.f90~~AfferentGraph sourcefile~tem_spongelayer_module.f90 tem_spongelayer_module.f90 sourcefile~tem_spatial_module.f90 tem_spatial_module.f90 sourcefile~tem_spatial_module.f90->sourcefile~tem_spongelayer_module.f90 sourcefile~tem_ini_condition_module.f90 tem_ini_condition_module.f90 sourcefile~tem_ini_condition_module.f90->sourcefile~tem_spatial_module.f90 sourcefile~tem_spacetime_var_module.f90 tem_spacetime_var_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_spatial_module.f90 sourcefile~tem_spacetime_fun_module.f90 tem_spacetime_fun_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_variable_module.f90 tem_variable_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_variable_module.f90 sourcefile~tem_spacetime_fun_module.f90->sourcefile~tem_spatial_module.f90 sourcefile~tem_varmap_module.f90 tem_varMap_module.f90 sourcefile~tem_varmap_module.f90->sourcefile~tem_spacetime_var_module.f90 sourcefile~tem_varmap_module.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_logical_operator_test.f90 tem_logical_operator_test.f90 sourcefile~tem_logical_operator_test.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_derived_module.f90 tem_derived_module.f90 sourcefile~tem_logical_operator_test.f90->sourcefile~tem_derived_module.f90 sourcefile~tem_logical_operator_test.f90->sourcefile~tem_variable_module.f90 sourcefile~tem_varsys_derivevar_test.f90 tem_varSys_deriveVar_test.f90 sourcefile~tem_varsys_derivevar_test.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_variable_evaltype_test.f90 tem_variable_evaltype_test.f90 sourcefile~tem_variable_evaltype_test.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_variable_evaltype_test.f90->sourcefile~tem_derived_module.f90 sourcefile~tem_variable_evaltype_test.f90->sourcefile~tem_variable_module.f90 sourcefile~tem_variable_combine_test.f90 tem_variable_combine_test.f90 sourcefile~tem_variable_combine_test.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_variable_combine_test.f90->sourcefile~tem_derived_module.f90 sourcefile~tem_variable_combine_test.f90->sourcefile~tem_variable_module.f90 sourcefile~tem_face_module.f90 tem_face_module.f90 sourcefile~tem_face_module.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_derived_module.f90->sourcefile~tem_spacetime_var_module.f90 sourcefile~tem_derived_module.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_derived_module.f90->sourcefile~tem_variable_module.f90 sourcefile~tem_variable_extract_test.f90 tem_variable_extract_test.f90 sourcefile~tem_variable_extract_test.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_variable_extract_test.f90->sourcefile~tem_derived_module.f90 sourcefile~tem_variable_extract_test.f90->sourcefile~tem_variable_module.f90 sourcefile~tem_bc_module.f90 tem_bc_module.f90 sourcefile~tem_bc_module.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_bc_module.f90->sourcefile~tem_varmap_module.f90 sourcefile~tem_varsys_test.f90 tem_varSys_test.f90 sourcefile~tem_varsys_test.f90->sourcefile~tem_spacetime_var_module.f90 sourcefile~tem_varsys_test.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_varsys_test.f90->sourcefile~tem_variable_module.f90 sourcefile~tem_surfacedata_module.f90 tem_surfaceData_module.f90 sourcefile~tem_surfacedata_module.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_varsys_stfunvar_test.f90 tem_varSys_stfunVar_test.f90 sourcefile~tem_varsys_stfunvar_test.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_varsys_stfunvar_test.f90->sourcefile~tem_derived_module.f90 sourcefile~tem_varsys_stfunvar_test.f90->sourcefile~tem_variable_module.f90 sourcefile~tem_varsys_opvar_test.f90 tem_varSys_opVar_test.f90 sourcefile~tem_varsys_opvar_test.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_varsys_opvar_test.f90->sourcefile~tem_derived_module.f90 sourcefile~tem_varsys_opvar_test.f90->sourcefile~tem_variable_module.f90 sourcefile~tem_variable_module.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_spacetime_fun_test.f90 tem_spacetime_fun_test.f90 sourcefile~tem_spacetime_fun_test.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_parallel_singlelevel_test.f90 tem_parallel_singlelevel_test.f90 sourcefile~tem_parallel_singlelevel_test.f90->sourcefile~tem_face_module.f90 sourcefile~tem_serial_multilevel_2_test.f90 tem_serial_multilevel_2_test.f90 sourcefile~tem_serial_multilevel_2_test.f90->sourcefile~tem_face_module.f90 sourcefile~tem_depend_module.f90 tem_depend_module.f90 sourcefile~tem_depend_module.f90->sourcefile~tem_varmap_module.f90 sourcefile~tem_restart_module.f90 tem_restart_module.f90 sourcefile~tem_restart_module.f90->sourcefile~tem_varmap_module.f90 sourcefile~tem_tracking_module.f90 tem_tracking_module.f90 sourcefile~tem_tracking_module.f90->sourcefile~tem_varmap_module.f90 sourcefile~tem_serial_singlelevel_test.f90 tem_serial_singlelevel_test.f90 sourcefile~tem_serial_singlelevel_test.f90->sourcefile~tem_face_module.f90 sourcefile~tem_operation_var_module.f90 tem_operation_var_module.f90 sourcefile~tem_operation_var_module.f90->sourcefile~tem_varmap_module.f90 sourcefile~tem_operation_var_module.f90->sourcefile~tem_variable_module.f90 sourcefile~hvs_output_module.f90 hvs_output_module.f90 sourcefile~hvs_output_module.f90->sourcefile~tem_varmap_module.f90 sourcefile~tem_convergence_module.f90 tem_convergence_module.f90 sourcefile~tem_convergence_module.f90->sourcefile~tem_varmap_module.f90

Contents


Source Code

! Copyright (c) 2013-2015 Nikhil Anand <nikhil.anand@uni-siegen.de>
! Copyright (c) 2013,2015-2016,2021-2022 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2014 Jens Zudrop <j.zudrop@grs-sim.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2018 Neda Ebrahimi Pour <neda.epour@uni-siegen.de>
! Copyright (c) 2019 Harald Klimach <harald.klimach@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 COPYRIGHT HOLDERS AND CONTRIBUTORS "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 THE COPYRIGHT HOLDER 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.
! ****************************************************************************** !
!> This module gathers the various spatial spongeLayer definitions
!!
module tem_spongeLayer_module

  ! include treelm modules
  use env_module,          only: rk, long_k, labelLen
  use tem_param_module,    only: PI
  use tem_aux_module,      only: tem_abort
  use tem_logging_module,  only: logUnit
  use treelmesh_module,    only: treelmesh_type
  use tem_geometry_module, only: tem_BaryOfId

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

  implicit none
  private

  public :: tem_spongeLayer_plane_type, tem_load_spongeLayer_plane
  public :: tem_spongeLayer_radial_type, tem_load_spongeLayer_radial
  public :: tem_spongeLayer_box_type, tem_load_spongeLayer_box
  public :: tem_spongelayer_plane_for
  public :: tem_spongelayer_box_for
  public :: tem_spongelayer_box2d_for
  public :: tem_spongeLayer_radial_for
  public :: tem_viscSpongeLayer_plane_for
  public :: tem_viscSpongeLayer_box_for
  public :: tem_viscSpongeLayer_box2d_for
  public :: tem_viscSpongeLayer_radial_for


  !> This type contains base data defined for all sponge layers
  type spongeLayer_base_type
    !> Thickness of the sponge layer.
    !! For planar sponge thickness is defined implicitly in place_normal
    real(kind=rk) :: thickness
    !> Damp factor or strength for the sponge Layer
    real(kind=rk) :: dampFactor
    !> damping exponent for the sponge layer
    real(kind=rk) :: dampExponent
    !> damping profile
    character(len=labelLen) :: dampProfile
    !> target states.
    !! For viscous sponge, viscosity is stored and multiplied with sponge
    !! strength
    real(kind=rk), allocatable :: targetState(:)
  end type spongeLayer_base_type

  !> This type contains data to define spongeLayer plane
  type, extends(spongeLayer_base_type) :: tem_spongeLayer_plane_type
    !> Sponge Plane origin
    real(kind=rk) :: origin(3)
    !> Sponge Plane normal
    real(kind=rk) :: normal(3)
  end type tem_spongeLayer_plane_type

  !> This type contains data to define spongeLayer box
  type, extends(spongeLayer_base_type) :: tem_spongeLayer_box_type
    !> Box origin, bottom left corner of sponge layer
    real(kind=rk) :: origin(3)
    !> Length of box in each dimension 
    real(kind=rk) :: extent(3)
    !> To create sponge box with rounded corners
    logical :: rounded_corner
    !> Corner radius for rounded box
    real(kind=rk) :: corner_radius
  end type tem_spongeLayer_box_type

  !> This type contains data to define spongeLayer radial
  type, extends(spongeLayer_base_type) :: tem_spongeLayer_radial_type
    !> Sponge radial origin
    real(kind=rk) :: origin(3)
    !> Sponge inner radius i.e. sponge start.
    !! Outer radius is computed by adding thickness to inner radius.
    real(kind=rk) :: radius
  end type tem_spongeLayer_radial_type

  !> Interface for sponge layer plane
  interface tem_spongeLayer_plane_for
    module procedure spongeLayer_plane_scalar_for_coord
    module procedure spongeLayer_plane_scalar_for_treeIDs
    module procedure spongeLayer_plane_vector_for_coord
    module procedure spongeLayer_plane_vector_for_treeIDs
  end interface tem_spongeLayer_plane_for

  !> Interface for sponge layer box
  interface tem_spongeLayer_box_for
    module procedure spongeLayer_box_scalar_for_coord
    module procedure spongeLayer_box_scalar_for_treeIDs
    module procedure spongeLayer_box_vector_for_coord
    module procedure spongeLayer_box_vector_for_treeIDs
  end interface tem_spongeLayer_box_for

  !> Interface for sponge layer box 2d
  interface tem_spongeLayer_box2d_for
    module procedure spongeLayer_box2d_scalar_for_coord
    module procedure spongeLayer_box2d_scalar_for_treeIDs
    module procedure spongeLayer_box2d_vector_for_coord
    module procedure spongeLayer_box2d_vector_for_treeIDs
  end interface tem_spongeLayer_box2d_for

  !> Interface for sponge layer radial
  interface tem_spongeLayer_radial_for
    module procedure spongeLayer_radial_scalar_for_coord
    module procedure spongeLayer_radial_scalar_for_treeIDs
    module procedure spongeLayer_radial_vector_for_coord
    module procedure spongeLayer_radial_vector_for_treeIDs
  end interface tem_spongeLayer_radial_for

  !> Interface for viscous sponge layer plane
  interface tem_viscSpongeLayer_plane_for
    module procedure viscSpongeLayer_plane_for_coord
    module procedure viscSpongeLayer_plane_for_treeIDs
  end interface tem_viscSpongeLayer_plane_for

  !> Interface for viscous sponge layer box
  interface tem_viscSpongeLayer_box_for
    module procedure viscSpongeLayer_box_for_coord
    module procedure viscSpongeLayer_box_for_treeIDs
  end interface tem_viscSpongeLayer_box_for

  !> Interface for viscous sponge layer box
  interface tem_viscSpongeLayer_box2d_for
    module procedure viscSpongeLayer_box2d_for_coord
    module procedure viscSpongeLayer_box2d_for_treeIDs
  end interface tem_viscSpongeLayer_box2d_for

  !> Interface for viscous sponge layer radial
  interface tem_viscSpongeLayer_radial_for
    module procedure viscSpongeLayer_radial_for_coord
    module procedure viscSpongeLayer_radial_for_treeIDs
  end interface tem_viscSpongeLayer_radial_for


contains


  ! -------------------------------------------------------------------------- !
  !> This subroutine load data for standard plane sponge layer
  !! Example:
  !!
  !!```lua
  !! spatial = {
  !!   -- supported options: 'spongelayer_plane', 'spongelayer_plane_1d', 
  !!   --                    'spongelayer_plane_2d', 'viscous_spongelayer_plane'
  !!   predefined = 'spongelayer',
  !!   origin = {0.0,0.0,0.0},
  !!   normal = {1.0, 0.0, 0.0},
  !!   thickness = 0.5,
  !!   damp_profile = 'linear', --'exponential', 'polynomial_n5', 'polynomial_n6'
  !!   damp_factor = 0.5,
  !!   damp_exponent = 1.0,
  !!   target_state = {
  !!     Default: density, velocityX, velocityY, velocityZ and pressure
  !!     density = 1.0,
  !!     pressure = 1.0,
  !!     velocityX = 0.0, velocityY = 0.0, velocityZ = 0.0
  !! }
  !!```
  subroutine tem_load_spongeLayer_plane(me, conf, thandle, ndim, nComp, &
    &                                   stateName)
    ! --------------------------------------------------------------------------
    !> Plane spongeLayer data type
    type(tem_spongeLayer_plane_type), intent(out) :: me
    !> lua state type
    type(flu_State) :: conf
    !> aotus parent handle
    integer, intent(in) :: thandle
    !> number of Dimension for nonViscous sponges
    integer, intent(in) :: nDim
    !> Number of component of St-Fun variable under which this spatial function
    !! is defined
    integer, intent(in) :: nComp
    !> Load stateName from target_state table
    character(len=*), intent(in), optional :: stateName
    ! --------------------------------------------------------------------------
    integer :: vError(3), errfatal(3)
    real(kind=rk) :: thickness
    ! --------------------------------------------------------------------------
    errfatal = aotErr_Fatal
    ! Plane_origin
    call aot_get_val( L       = conf,      &
      &               thandle = thandle,   &
      &               key     = 'origin',  &
      &               val     = me%origin, &
      &               ErrCode = vError     )
    if (any(btest( vError, errFatal )) ) then
      write(logUnit(1),*) 'ERROR reading the plane_origin of sponge layer. ' &
        &              // 'It should have 3 entries for each coordinate.'
      call tem_abort()
    end if

    write(logUnit(1),*) ' * Origin =', me%origin

    ! Plane_normal
    call aot_get_val( L       = conf,      &
      &               thandle = thandle,   &
      &               key     = 'normal',  &
      &               val     = me%normal, &
      &               ErrCode = vError     )
    if (any(btest( vError, errFatal )) ) then
      write(logUnit(1),*) 'ERROR reading the plane_normal of sponge layer. ' &
        &              // 'It should have 3 entries for each coordinate.'
      call tem_abort()
    end if
    write(logUnit(1),*) ' * Normal =', me%normal

    ! Compute thickness from normal and use it only if thickness is not
    ! defined seperately.
    thickness = sqrt(me%normal(1)**2 + me%normal(2)**2 + me%normal(3)**2)

    ! Normalize the normal
    me%normal = me%normal/thickness
    write(logUnit(1),*) ' * Normalized normal =', me%normal

    ! Load base information required for sponge layer definition like
    ! damp_factor, damp_exponent and target_state
    call load_spongeLayer( conf      = conf,                     &
      &                    thandle   = thandle,                  &
      &                    me        = me%spongeLayer_base_type, &
      &                    nDim      = nDim,                     &
      &                    nComp     = nComp,                    &
      &                    stateName = stateName,                &
      &                    thickness = thickness                 )

  end subroutine tem_load_spongeLayer_plane
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This subroutine load data for radial sponge layer
  !! Example:
  !!
  !!```lua
  !! spatial = {
  !!   --supported options: 'spongelayer_radial','sponge_radial_2d',
  !!   --                   'viscous_spongelayer_radial',
  !!   --                   'viscous_spongelayer_radial_2d'
  !!   predefined = 'viscous_spongelayer_radial', 
  !!   origin = {0.0,0.0,0.0},
  !!   radius = 1.0, -- Sponge start
  !!   thickness = 0.3,
  !!   damp_profile = 'linear', --'exponential', 'polynomial_n5', 'polynomial_n6'
  !!   damp_factor = 0.5,
  !!   damp_exponent = 1.0,
  !!   target_state = {
  !!     Default: density, velocityX, velocityY, velocityZ and pressure
  !!     viscosity = 1e-3
  !! }
  !!```
  subroutine tem_load_spongeLayer_radial(me, conf, thandle, nDim, nComp, &
    &                                    stateName)
    ! --------------------------------------------------------------------------
    !> Radial spongeLayer data type
    type(tem_spongeLayer_radial_type), intent(out) :: me
    !> lua state type
    type(flu_State) :: conf
    !> aotus parent handle
    integer, intent(in) :: thandle
    !> number of Dimension for nonViscous sponges
    integer, intent(in) :: nDim
    !> Number of component of St-Fun variable under which this spatial function
    !! is defined
    integer, intent(in) :: nComp
    !> Load stateName from target_state table
    character(len=*), intent(in), optional :: stateName
    ! --------------------------------------------------------------------------
    integer :: iError
    integer :: vError(3), errfatal(3)
    ! --------------------------------------------------------------------------

    errfatal = aotErr_Fatal
    ! Sponge origin
    call aot_get_val( L       = conf,      &
      &               thandle = thandle,   &
      &               key     = 'origin',  &
      &               val     = me%origin, &
      &               ErrCode = vError     )
    if (any(btest( vError, errFatal )) ) then
      write(logUnit(1),*) 'ERROR reading the sponge origin, ' &
        &              // 'origin is not well defined. '   &
        &              // 'It should have 3 entries for each coordinate.'
      call tem_abort()
    end if
    write(logUnit(1),*) ' * Origin =', me%origin

    ! Sponge inner radius
    call aot_get_val( L       = conf,      &
      &               thandle = thandle,   &
      &               key     = 'radius',  &
      &               val     = me%radius, &
      &               ErrCode = iError     )
    if (btest(iError, aotErr_Fatal)) then
      write(logUnit(1),*) 'ERROR reading the sponge inner radius. '
      call tem_abort()
    end if
    write(logUnit(1),*) ' * Inner radius =', me%radius

    ! Load base information required for sponge layer definition like
    ! damp_factor, damp_exponent and target_state
    call load_spongeLayer( conf      = conf,                     &
      &                    thandle   = thandle,                  &
      &                    me        = me%spongeLayer_base_type, &
      &                    nDim      = nDim,                     &
      &                    nComp     = nComp,                    &
      &                    stateName = stateName                 )

  end subroutine tem_load_spongeLayer_radial
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This subroutine load data for sponge layer box
  !! Example:
  !!
  !!```lua
  !! spatial = {
  !!   --supported options: 'spongelayer_box', 'spongelayer_box_2d'
  !!   --                   'viscous_spongelayer_box',
  !!   predefined = 'spongelayer_box',
  !!   origin = {0.0,0.0,0.0},
  !!   extent = {1.0, 2.0, 3.0},
  !!   thickness = 0.3,
  !!   damp_profile = 'linear', --'exponential', 'polynomial_n5', 'polynomial_n6'
  !!   damp_factor = 0.5,
  !!   damp_exponent = 1.0,
  !!   target_state = {
  !!     Default: density, velocityX, velocityY, velocityZ and pressure
  !!     density = 1.0,
  !!     pressure = 1.0,
  !!     velocityX = 0.0, velocityY = 0.0, velocityZ = 0.0
  !! }
  !!```
  subroutine tem_load_spongeLayer_box(me, conf, thandle, ndim, nComp, stateName)
    ! --------------------------------------------------------------------------
    !> Box spongeLayer data type
    type(tem_spongeLayer_box_type), intent(out) :: me
    !> lua state type
    type(flu_State) :: conf
    !> aotus parent handle
    integer, intent(in) :: thandle
    !> number of Dimension for nonViscous sponges
    integer, intent(in) :: nDim
    !> Number of component of St-Fun variable under which this spatial function
    !! is defined
    integer, intent(in) :: nComp
    !> Load stateName from target_state table
    character(len=*), intent(in), optional :: stateName
    ! --------------------------------------------------------------------------
    integer :: vError(3), errfatal(3), iError
    real(kind=rk) :: min_halfextent
    ! --------------------------------------------------------------------------
    errfatal = aotErr_Fatal
    ! Box origin
    call aot_get_val( L       = conf,      &
      &               thandle = thandle,   &
      &               key     = 'origin',  &
      &               val     = me%origin, &
      &               ErrCode = vError     )
    if (any(btest( vError, errFatal )) ) then
      write(logUnit(1),*) 'ERROR reading the origin of box sponge layer. ' &
        &              // 'It should have 3 entries for each coordinate.'
      call tem_abort()
    end if

    write(logUnit(1),*) ' * Origin =', me%origin

    ! Box extent
    call aot_get_val( L       = conf,      &
      &               thandle = thandle,   &
      &               key     = 'extent',  &
      &               val     = me%extent, &
      &               ErrCode = vError     )
    if (any(btest( vError, errFatal )) ) then
      write(logUnit(1),*) 'ERROR reading the extent of box sponge layer. ' &
        &              // 'It should have 3 entries for each dimension.'
      call tem_abort()
    end if

    write(logUnit(1),*) ' * Extent =', me%extent

    ! Additional parameters for box with rounded corner
    call aot_get_val( L       = conf,              &
      &               thandle = thandle,           &
      &               key     = 'rounded_corner',  &
      &               val     = me%rounded_corner, &
      &               default = .false.,           &
      &               ErrCode = iError             )
    if (btest( iError, aotErr_Fatal )) then
      write(logUnit(1),*) 'WARNING: reading the rounded_corner of box '&
        &               //'sponge layer. Setting it to false'
      me%rounded_corner = .false.
    end if
    write(logUnit(1),*) ' * Rounded_corner =', me%rounded_corner

    if (me%rounded_corner) then
      call aot_get_val( L       = conf,              &
        &               thandle = thandle,           &
        &               key     = 'corner_radius',   &
        &               val     = me%corner_radius,  &
        &               ErrCode = iError             )
      if (btest( iError, aotErr_Fatal )) then
        write(logUnit(1),*) 'ERROR reading the corner_radius of box '&
          &               //'sponge layer with rounded corner.'
        call tem_abort()
      end if
      ! if corner radius is greater than minimum half extent of the box
      ! then set corner_radius to minimum of half extent
      min_halfextent = minval(me%extent(1:nDim))*0.5_rk
      if (me%corner_radius > min_halfextent) then
        write(logUnit(1),*) 'WARNING: corner_radius is greater than half of '
        write(logUnit(1),*) '  min of box extent. '
        write(logUnit(1),*) 'Setting it to half of min extent'
        me%corner_radius = min_halfextent
      end if
      write(logUnit(1),*) ' * Corner radius =', me%corner_radius
    end if

    ! Load base information required for sponge layer definition like
    ! damp_factor, damp_exponent and target_state
    call load_spongeLayer( conf      = conf,                     &
      &                    thandle   = thandle,                  &
      &                    me        = me%spongeLayer_base_type, &
      &                    nDim      = nDim,                     &
      &                    nComp     = nComp,                    &
      &                    stateName = stateName                 )

  end subroutine tem_load_spongeLayer_box
  ! -------------------------------------------------------------------------- !


  ! -------------------------------------------------------------------------- !
  !> This routine load base info for sponge layer
  subroutine load_spongeLayer(conf, thandle, me, ndim, nComp, stateName, &
    &                         thickness)
    ! --------------------------------------------------------------------------
    !> lua state type
    type(flu_State) :: conf
    !> aotus parent handle
    integer, intent(in) :: thandle
    !> base spongeLayer data type
    type(spongeLayer_base_type), intent(out) :: me
    !> number of spatial dimensions
    integer, intent(in) :: ndim
    !> Number of component of St-Fun variable under which this spatial function
    !! is defined
    integer, intent(in) :: nComp
    !> Load stateName from target_state table
    character(len=*), intent(in), optional :: stateName
    !> Thickness computed from sponge layer plane normal. Use this thickness
    !! If thickness is not defined.
    real(kind=rk), intent(in), optional :: thickness
    ! --------------------------------------------------------------------------
    integer :: iError, ts_handle
    ! --------------------------------------------------------------------------
    ! Thickness
    call aot_get_val( L       = conf,         &
      &               thandle = thandle,      &
      &               key     = 'thickness',  &
      &               val     = me%thickness, &
      &               ErrCode = iError        )
    if (btest( iError, aotErr_Fatal )) then
      if (present(thickness)) then
        me%thickness = thickness
      else
        write(logUnit(1),*) 'ERROR reading the thickness of sponge layer.'
        write(logUnit(1),*) 'Thickness is required to calculate sponge end.'
        call tem_abort()
      end if
    end if

    write(logUnit(1),*) ' * Thickness =', me%thickness

    !damp_factor
    call aot_get_val( L       = conf,          &
      &               thandle = thandle,       &
      &               key     = 'damp_factor', &
      &               val     = me%dampFactor, &
      &               ErrCode = iError         )
    if (btest( iError, aotErr_Fatal )) then
      write(logUnit(1),*) 'ERROR reading the damp_factor of sponge layer.'
      call tem_abort()
    end if
    write(logUnit(1),*) ' * Damp_factor =', me%dampFactor

    !damp_profile
    call aot_get_val( L       = conf,           &
      &               thandle = thandle,        &
      &               key     = 'damp_profile', &
      &               val     = me%dampProfile, &
      &               ErrCode = iError          )

    ! Viscous sponge works only with exponential profile so no need load 
    ! damp_profile
    if (present(stateName)) then
      if (trim(stateName) == 'viscosity') then
        me%dampProfile = 'exponential'
        iError = 0
      end if
    end if

    if (btest( iError, aotErr_Fatal )) then
      write(logUnit(1),*) 'ERROR reading the damp_profile of sponge layer.'
      call tem_abort()
    end if
    write(logUnit(1),*) ' * Damp_profile =', me%dampProfile

    select case (trim(me%dampProfile))
    case ('linear')
      me%dampExponent = 1.0
    case ('exponential')
      !damp_exponent
      call aot_get_val( L       = conf,            &
        &               thandle = thandle,         &
        &               key     = 'damp_exponent', &
        &               val     = me%dampExponent, &
        &               ErrCode = iError,          &
        &               default = 1.0_rk           )

      if (btest( iError, aotErr_Fatal )) then
        write(logUnit(1),*) 'ERROR reading the damp_exponent of exponential ' &
          &                 //'sponge layer.'
        call tem_abort()
      end if
    case ('polynomial_n5', 'polynomial_n6')
      me%dampExponent = 1.0
    case default
      write(logUnit(1),*) 'ERROR unknown damp_profile for sponge layer.'
      write(logUnit(1),*) 'Supported options: '
      write(logUnit(1),*) '  linear, exponential, polynomial_n5, polynomial_n6'
      call tem_abort()
    end select

    write(logUnit(1),*) ' * Damp_exponent =', me%dampExponent


    ! Load stateName provided by caller function
    if ( present(stateName) .and. nComp == 1) then
      allocate(me%targetState(1))
      call aot_table_open( L       = conf,          &
        &                  parent  = thandle,       &
        &                  thandle = ts_handle,     &
        &                  key     = 'target_state' )

      call aot_get_val( L       = conf,               &
        &               thandle = ts_handle,          &
        &               key     = trim(stateName),    &
        &               val     = me%targetState(1),  &
        &               ErrCode = iError              )

      if (btest(iError, aotErr_Fatal)) then
        write(logUnit(1),*) 'ERROR reading the target state: '//trim(stateName)
        call tem_abort()
      end if
      call aot_table_close(conf, thandle)

      write(logUnit(1),*) ' * Target state:'
      write(logUnit(1),*) '   '//trim(stateName)//' =', me%targetState(1)
    else if (nComp > 1) then
      call load_defaultTargetState( conf        = conf,          &
        &                           parent      = thandle,       &
        &                           nDim        = nDim,          &
        &                           nComp       = nComp,         &
        &                           targetState = me%targetState )
    else
      write(logUnit(1),*) 'WARNING: nComp = 1 so no target states are loaded'
    end if

  end subroutine load_spongeLayer
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This routine loads state names from target_state table
  subroutine load_defaultTargetState(conf, parent, nDim, nComp, targetState)
    ! --------------------------------------------------------------------------
    !> lua state type
    type(flu_State) :: conf
    !> aotus parent handle
    integer, intent(in) :: parent
    !> Number of dimension
    integer, intent(in) :: nDim
    !> Number of component of St-Fun variable under which this spatial function
    !! is defined
    integer, intent(in) :: nComp
    !> Target state value
    real(kind=rk), allocatable, intent(out) :: targetState(:)
    ! --------------------------------------------------------------------------
    integer :: thandle, iState, nState
    integer :: iError
    !> List of stateNames
    character(len=labelLen), allocatable :: stateName(:)
    ! --------------------------------------------------------------------------
    nState = 0
    select case (nDim)
    case (3)
      nState = 5
      allocate(stateName(nState))
      stateName = [ 'density  ', 'velocityX', &
        &           'velocityY', 'velocityZ', &
        &           'pressure '               ]
    case (2)
      nState = 4
      allocate(stateName(nState))
      stateName = [ 'density  ', 'velocityX', &
        &           'velocityY', 'pressure '  ]

    case (1)
      nState = 3
      allocate(stateName(nState))
      stateName = [ 'density  ', 'velocityX', &
        &           'pressure '               ]
    end select

    ! nState must be nComp - 1 to return all target states when evaluating
    ! the space-time function
    if (nComp /= nState+1) then
      write(logUnit(1),*) 'Error: Expected ncomponents = ', nState+1
      write(logUnit(1),*) '  Defined ncomponents in st-fun is ',  ncomp
      call tem_abort()
    end if

    ! target_state
    allocate(targetState(nState))
    call aot_table_open( L       = conf,          &
      &                  parent  = parent,        &
      &                  thandle = thandle,       &
      &                  key     = 'target_state' )

    write(logUnit(1),*) ' * Target state:'
    do iState = 1, nState
      call aot_get_val( L       = conf,                    &
           &            thandle = thandle,                 &
           &            key     = trim(stateName(iState)), &
           &            val     = targetState(iState),     &
           &            ErrCode = iError                   )

      if (btest(iError, aoterr_Fatal)) then
        write(*,*) 'FATAL Error occured, when loading target state: ' &
          &      // trim(stateName(iState))// '! Aborting'
        call tem_abort()
      end if

      write(logUnit(1),*) '   '//trim(stateName(iState))//' =', &
        &                 targetState(iState)
    end do

    call aot_table_close(conf, thandle)

  end subroutine load_defaultTargetState
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !
  !                            SPONGE LAYER PLANE
  !
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function returns the sigma for the planar shape spongelayer
  function spongelayer_plane_scalar_for_coord(me, coord, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_plane_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    select case (trim(me%dampProfile))
    case ('linear', 'exponential')
      res(:) = spongeLayer_plane_expon_for_coord(me, coord, n)
    case ('polynomial_n5')
      res(:) = spongeLayer_plane_polyn5_for_coord(me, coord, n)
    case ('polynomial_n6')
      res(:) = spongeLayer_plane_polyn6_for_coord(me, coord, n)
    end select

  end function spongelayer_plane_scalar_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function returns the sigma for the spongelayer computed from
  !! exponential function
  function spongelayer_plane_expon_for_coord(me, coord, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_plane_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(3), normal(3), vec(3), proj_len
    ! --------------------------------------------------------------------------
    origin(:) = me%origin
    ! Normal was normalized in load routine so multiply thickness to convert
    ! normal to user defined normal
    normal(:) = me%normal * me%thickness

    do i = 1,n
      vec(:) = coord(i,:) - origin(:)
      proj_len = (vec(1)*normal(1)+ vec(2)*normal(2)+vec(3)*normal(3))/   &
                 (normal(1)**2 + normal(2)**2 + normal(3)**2)
      sigma = me%dampFactor*((proj_len)**me%dampExponent)
      if (proj_len > 0 .and. proj_len < 1) then
        res(i) = sigma
      else if (proj_len > 1) then
        res(i) = me%dampFactor
      else
        res(i) = 0.0_rk
      end if

    end do

  end function spongelayer_plane_expon_for_coord
  ! -------------------------------------------------------------------------- !

 
  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the spongeLayer from coord for
  !! the polynomial order n5.
  !! Sponge profile:
  !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$
  !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
  !!
  !! Profile is taken from:
  !! Xu, Hui; Sagaut, Pierre (2013): Analysis of the absorbing layers for the 
  !! weakly-compressible lattice Boltzmann methods. In Journal of Computational 
  !! Physics 245, pp. 14-42. DOI: 10.1016/j.jcp.2013.02.051.
  function spongeLayer_plane_polyn5_for_coord(me, coord, n)  result(res)
    ! --------------------------------------------------------------------------
    !> Spatial sponge layer to evaluate
    type(tem_spongeLayer_plane_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(3), normal(3), vec1(3), vec2(3)
    real(kind=rk) :: proj_len1, proj_len2, const_fac
    ! --------------------------------------------------------------------------
    origin(:) = me%origin
    normal(:) = me%normal
    const_fac = 3125_rk/(256_rk*me%thickness**5)

    do i = 1,n
      vec1(:) = coord(i,:) - origin(:)
      vec2(:) = me%thickness*normal(:) + origin(:) - coord(i,:)
      proj_len1 = vec1(1)*normal(1) + vec1(2)*normal(2) + vec1(3)*normal(3)
      proj_len2 = vec2(1)*normal(1) + vec2(2)*normal(2) + vec2(3)*normal(3)

      if (proj_len1 > 0 .and. proj_len2 > 0) then
        sigma = const_fac * proj_len2 * (proj_len1**4)
        res(i) = sigma*me%dampFactor
      else if (proj_len2 < 0) then ! If coord is beyond thickness
        res(i) = me%dampFactor
      else
        res(i) = 0.0_rk
      end if
    end do

  end function spongeLayer_plane_polyn5_for_coord
  ! -------------------------------------------------------------------------- !


  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the spongeLayer from coord for
  !! the polynomial order n6.
  !! Sponge profile:
  !! $$ \sigma(x) = \sigma_A*729*(L+x0-x)^2*(x-x0)^4)/(16*(L)^6 $$
  !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
  function spongeLayer_plane_polyn6_for_coord(me, coord, n)  result(res)
    ! --------------------------------------------------------------------------
    !> Spatial sponge layer to evaluate
    type(tem_spongeLayer_plane_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(3), normal(3), vec1(3), vec2(3)
    real(kind=rk) :: proj_len1, proj_len2, const_fac
    ! --------------------------------------------------------------------------
    origin(:) = me%origin
    normal(:) = me%normal
    const_fac = 729_rk/(16_rk*me%thickness**6)

    do i = 1,n
      vec1(:) = coord(i,:) - origin(:)
      vec2(:) = me%thickness*normal(:) + origin(:) - coord(i,:)
      proj_len1 = vec1(1)*normal(1) + vec1(2)*normal(2) + vec1(3)*normal(3)
      proj_len2 = vec2(1)*normal(1) + vec2(2)*normal(2) + vec2(3)*normal(3)

      sigma = const_fac * proj_len2**2 * (proj_len1**4)
      if (proj_len1 > 0) then
        res(i) = min(1.0_rk, sigma) * me%dampFactor
      else
        res(i) = 0.0_rk
      end if
    end do

  end function spongeLayer_plane_polyn6_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the spongelayer and fills up
  !!  the res with the target state
  function spongelayer_plane_vector_for_coord(me, nComp, coord, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_plane_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> Number of entrys in each array
    integer, intent(in) :: ncomp
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n,ncomp)
    ! --------------------------------------------------------------------------
    integer :: i
    ! --------------------------------------------------------------------------
    res(:, 1) = spongeLayer_plane_scalar_for_coord(me, coord, n)

    if (ncomp > 1) then
      do i = 1,n
        res(i,2:) = me%targetState(:)
      end do
    end if

  end function spongelayer_plane_vector_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function returns the sigma for the spongelayer from treeids
  function spongelayer_plane_scalar_for_treeIDs(me, treeids, tree, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_plane_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    select case (trim(me%dampProfile))
    case ('linear', 'exponential')
      res(:) = spongeLayer_plane_expon_for_treeIDs(me, treeIDs, tree, n)
    case ('polynomial_n5')
      res(:) = spongeLayer_plane_polyn5_for_treeIDs(me, treeIDs, tree, n)
    case ('polynomial_n6')
      res(:) = spongeLayer_plane_polyn6_for_treeIDs(me, treeIDs, tree, n)
    end select

  end function spongeLayer_plane_scalar_for_treeIDs
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function returns the sigma for the exponential spongelayer from 
  !! treeids
  function spongelayer_plane_expon_for_treeIDs(me, treeids, tree, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_plane_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(3), normal(3), vec(3), coord(3), proj_len
    ! --------------------------------------------------------------------------
    origin(:) = me%origin
    ! Normal was normalized in load routine so multiply thickness to convert
    ! normal to user defined normal
    normal(:) = me%normal * me%thickness

    do i = 1,n
      !barycentric coordinate
      coord = tem_BaryOfId( tree, treeIds(i) )
      vec (:) = coord(:) - origin(:)
      proj_len = (vec(1)*normal(1)+ vec(2)*normal(2)+vec(3)*normal(3))/   &
                 (normal(1)**2 + normal(2)**2 + normal(3)**2)
      sigma = me%dampFactor*((proj_len)**me%dampExponent)
      if (proj_len > 0 .and. proj_len < 1) then
        res(i) = sigma
      else if (proj_len > 1) then
        res(i) = me%dampFactor
      else
        res(i) = 0.0_rk
      end if

    end do

  end function spongeLayer_plane_expon_for_treeIDs
  ! -------------------------------------------------------------------------- !


  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the spongelayer using
  !! polynomial order 5 from treeids.
  !! Sponge profile:
  !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$
  !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
  !!
  !! Profile is taken from:
  !! Xu, Hui; Sagaut, Pierre (2013): Analysis of the absorbing layers for the 
  !! weakly-compressible lattice Boltzmann methods. In Journal of Computational 
  !! Physics 245, pp. 14-42. DOI: 10.1016/j.jcp.2013.02.051.
  function spongeLayer_plane_polyn5_for_treeids(me, treeids, tree, n) &
    & result(res)
    ! --------------------------------------------------------------------------
    !> Spatial absorb layer to evaluate
    type(tem_spongeLayer_plane_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(3), normal(3), vec1(3), vec2(3), coord(3)
    real(kind=rk) :: proj_len1, proj_len2, const_fac
    ! --------------------------------------------------------------------------
    origin(:) = me%origin
    normal(:) = me%normal
    const_fac = 3125_rk/(256_rk*me%thickness**5)

    do i = 1,n
      !barycentric coordinate
      coord = tem_BaryOfId( tree, treeIds(i) )
      vec1(:) = coord(:) - origin(:)
      vec2(:) = me%thickness*normal(:) + origin(:) - coord(:)
      proj_len1 = vec1(1)*normal(1)+ vec1(2)*normal(2)+vec1(3)*normal(3)
      proj_len2 = vec2(1)*normal(1)+ vec2(2)*normal(2)+vec2(3)*normal(3)

      if (proj_len1 > 0 .and. proj_len2 > 0) then
        sigma = const_fac * proj_len2 * (proj_len1**4)
        res(i) = sigma*me%dampFactor
      else if (proj_len2 < 0) then ! If coord is beyond thickness
        res(i) = me%dampFactor
      else
        res(i) = 0.0_rk
      end if

    end do

  end function spongeLayer_plane_polyn5_for_treeids
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the spongeLayer from treeid for
  !! the polynomial order n6.
  !! Sponge profile:
  !! $$ \sigma(x) = \sigma_A*729*(L+x0-x)^2*(x-x0)^4)/(16*(L)^6 $$
  !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
  function spongeLayer_plane_polyn6_for_treeids(me, treeids, tree, n) &
    & result(res)
    ! --------------------------------------------------------------------------
    !> Spatial sponge layer to evaluate
    type(tem_spongeLayer_plane_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(3), normal(3), vec1(3), vec2(3), coord(3)
    real(kind=rk) :: proj_len1, proj_len2, const_fac
    ! --------------------------------------------------------------------------
    origin(:) = me%origin
    normal(:) = me%normal
    const_fac = 729_rk/(16_rk*me%thickness**6)

    do i = 1,n
      !barycentric coordinate
      coord = tem_BaryOfId( tree, treeIds(i) )
      vec1(:) = coord(:) - origin(:)
      vec2(:) = me%thickness*normal(:) + origin(:) - coord(:)
      proj_len1 = vec1(1)*normal(1)+ vec1(2)*normal(2)+vec1(3)*normal(3)
      proj_len2 = vec2(1)*normal(1)+ vec2(2)*normal(2)+vec2(3)*normal(3)

      if (proj_len1 > 0) then
        sigma = const_fac * proj_len2**2 * (proj_len1**4)
        res(i) = min(1.0_rk, sigma) * me%dampFactor
      else
        res(i) = 0.0_rk
      end if

    end do

  end function spongeLayer_plane_polyn6_for_treeids
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the spongelayer and fills up
  !!  the res with the target state
  function spongelayer_plane_vector_for_treeIDs(me, ncomp, treeids, tree, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_plane_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> Number of entrys in each array
    integer, intent(in) :: ncomp
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n,ncomp)
    ! --------------------------------------------------------------------------
    integer :: i
    ! --------------------------------------------------------------------------
    res(:, 1) = spongeLayer_plane_scalar_for_treeIDs(me, treeids, tree, n)

    if (ncomp > 1) then
      do i = 1,n
        res(i,2:) = me%targetState(:)
      end do
    end if

  end function spongelayer_plane_vector_for_treeIDs
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !
  !                            SPONGE LAYER BOX
  !
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function returns the sigma for the box shape spongelayer
  function spongelayer_box_scalar_for_coord(me, coord, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    select case (trim(me%dampProfile))
    case ('linear', 'exponential')
      res(:) = spongeLayer_box_expon_for_coord(me, coord, n)
    case ('polynomial_n5')
      if (me%rounded_corner) then
        res(:) = spongeLayer_box_roundCornerPolyn5_for_coord(me, coord, n)
      else
        res(:) = spongeLayer_box_sharpCornerPolyn5_for_coord(me, coord, n)
      end if
    case ('polynomial_n6')
      if (me%rounded_corner) then
        res(:) = spongeLayer_box_roundCornerPolyn6_for_coord(me, coord, n)
      else
        res(:) = spongeLayer_box_sharpCornerPolyn6_for_coord(me, coord, n)
      end if
    end select

  end function spongelayer_box_scalar_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function returns sigma for the box shape spongelayer for coord for
  !! exponential profile.
  function spongelayer_box_expon_for_coord(me, coord, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(3), extent(3), box_max(3), proj_len
    real(kind=rk) :: coordLoc(3), normal, vec_min(3), vec_max(3)
    real(kind=rk) :: rad, vec_minSqr(3), vec_maxSqr(3)
    ! --------------------------------------------------------------------------
    origin(:) = me%origin
    extent(:) = me%extent
    box_max(:) = origin(:) + extent(:)

    do i = 1,n
      coordLoc = coord(i,:)
      vec_min(:) = coordLoc(:) - origin(:)
      vec_max(:) = coordLoc(:) - box_max(:)
      vec_minSqr(:) = vec_min(:)**2
      vec_maxSqr(:) = vec_max(:)**2
      normal = 1.0_rk
      rad = 0.0_rk

      ! Bottom-South-West: -x,-y,-z
      if (all(coordLoc(:) < origin(:))) then
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(2), vec_minSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2) &
        & .and. coordLoc(3) < origin(3)) then ! Bottom-South-East: x, -y, -z
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2), vec_minSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2) &
        & .and. coordLoc(3) < origin(3)) then ! Bottom-North-West: -x, y, -z
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2), vec_minSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) > box_max(2) &
        & .and. coordLoc(3) < origin(3)) then ! Bottom-North-East: x, y, -z
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2), vec_minSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) < origin(2) &
        & .and. coordLoc(3) > box_max(3)) then ! Top-South-West: -x, -y, z
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2) &
        & .and. coordLoc(3) > box_max(3)) then ! Top-South-East: x, -y, z
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2) &
        & .and. coordLoc(3) > box_max(3)) then ! Top-North-West: -x, y, z
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2), vec_maxSqr(3)) )

      else if (all(coordLoc(:) > box_max(:))) then ! Bottom-North-East: x, y, z
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(2) < origin(2) .and. coordLoc(3) < origin(3)) then
        ! Botton-South: -y,-z
        rad = sqrt( max(vec_minSqr(2), vec_minSqr(3)) )

      else if (coordLoc(2) < origin(2) .and. coordLoc(3) > box_max(3)) then
        ! Top-South: -y, z
        rad = sqrt( max(vec_minSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(2) > box_max(2) .and. coordLoc(3) < origin(3)) then
        ! Botom-North: -y, z
        rad = sqrt( max(vec_maxSqr(2), vec_minSqr(3)) )

      else if (coordLoc(2) > box_max(2) .and. coordLoc(3) > box_max(3)) then
        ! Top-North: y, z
        rad = sqrt( max(vec_maxSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(3) < origin(3)) then
        ! Botton-West: -x,-z
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(3) < origin(3)) then
        ! Bottom-East: x,-z
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(3) > box_max(3)) then
        ! Top-West: -x, z
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(3) > box_max(3)) then
        ! Top-East: x, z
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) < origin(2)) then
        ! South-West: -x,-y
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(2)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2)) then
        ! North-West: -x, y
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2)) then
        ! South-East: x, -y
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) > box_max(2)) then
        ! North-East: x, y
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2)) )

      else if (coordLoc(1) < origin(1)) then ! West: -x
        normal = -1_rk
        rad = vec_min(1)

      else if (coordLoc(2) < origin(2)) then ! South: -y
        normal = -1_rk
        rad = vec_min(2)

      else if (coordLoc(3) < origin(3)) then ! Bottom: -z
        normal = -1_rk
        rad = vec_min(3)

      else if (coordLoc(1) > box_max(1)) then ! East: x
        normal = 1_rk
        rad = vec_max(1)

      else if (coordLoc(2) > box_max(2)) then ! North: y
        normal = 1_rk
        rad = vec_max(2)

      else if (coordLoc(3) > box_max(3)) then ! Top: z
        normal = 1_rk
        rad = vec_max(3)
      end if

      proj_len = rad*normal/me%thickness
      sigma = me%dampFactor*((proj_len)**me%dampExponent)
      if (proj_len > 0 .and. proj_len < 1) then
        res(i) = sigma
      else if (proj_len > 1) then
        res(i) = me%dampFactor
      else
        res(i) = 0.0_rk
      end if

    end do

  end function spongelayer_box_expon_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the box shape sponge layer from 
  !! coord for polynomial n5.
  !! Sponge profile:
  !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$
  !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
  !!
  !! Profile is taken from:
  !! Xu, Hui; Sagaut, Pierre (2013): Analysis of the absorbing layers for the 
  !! weakly-compressible lattice Boltzmann methods. In Journal of Computational 
  !! Physics 245, pp. 14-42. DOI: 10.1016/j.jcp.2013.02.051.
  function spongeLayer_box_sharpCornerPolyn5_for_coord(me, coord, n) &
    & result(res)
    ! --------------------------------------------------------------------------
    !> Spatial sponge layer to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(3), extent(3), coordLoc(3), normal
    real(kind=rk) :: proj_len1, proj_len2
    real(kind=rk) :: vec_min(3), vec_max(3), vec_minSqr(3), vec_maxSqr(3) 
    real(kind=rk) :: rad, const_fac, box_max(3)
    ! --------------------------------------------------------------------------
    origin(:) = me%origin
    extent(:) = me%extent
    box_max(:) = origin(:) + extent(:)

    const_fac = 3125_rk/(256_rk*me%thickness**5)

    do i = 1,n
      coordLoc = coord(i,:)
      vec_min(:) = coordLoc(:) - origin(:)
      vec_max(:) = coordLoc(:) - box_max(:)
      vec_minSqr(:) = vec_min(:)**2
      vec_maxSqr(:) = vec_max(:)**2
      proj_len1 = 0_rk
      proj_len2 = 0_rk
      normal = 1_rk
      rad = 0_rk

      ! Bottom-South-West: -x,-y,-z
      if (all(coordLoc(:) < origin(:))) then
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(2), vec_minSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2) &
        & .and. coordLoc(3) < origin(3)) then ! Bottom-South-East: x, -y, -z
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2), vec_minSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2) &
        & .and. coordLoc(3) < origin(3)) then ! Bottom-North-West: -x, y, -z
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2), vec_minSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) > box_max(2) &
        & .and. coordLoc(3) < origin(3)) then ! Bottom-North-East: x, y, -z
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2), vec_minSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) < origin(2) &
        & .and. coordLoc(3) > box_max(3)) then ! Top-South-West: -x, -y, z
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2) &
        & .and. coordLoc(3) > box_max(3)) then ! Top-South-East: x, -y, z
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2) &
        & .and. coordLoc(3) > box_max(3)) then ! Top-North-West: -x, y, z
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2), vec_maxSqr(3)) )

      else if (all(coordLoc(:) > box_max(:))) then ! Top-North-East: x, y, z
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(2) < origin(2) .and. coordLoc(3) < origin(3)) then
        ! Botton-South: -y,-z
        rad = sqrt( max(vec_minSqr(2), vec_minSqr(3)) )

      else if (coordLoc(2) < origin(2) .and. coordLoc(3) > box_max(3)) then
        ! Top-South: -y, z
        rad = sqrt( max(vec_minSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(2) > box_max(2) .and. coordLoc(3) < origin(3)) then
        ! Botom-North: -y, z
        rad = sqrt( max(vec_maxSqr(2), vec_minSqr(3)) )

      else if (coordLoc(2) > box_max(2) .and. coordLoc(3) > box_max(3)) then
        ! Top-North: y, z
        rad = sqrt( max(vec_maxSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(3) < origin(3)) then
        ! Botton-West: -x,-z
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(3) < origin(3)) then
        ! Bottom-East: x,-z
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(3) > box_max(3)) then
        ! Top-West: -x, z
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(3) > box_max(3)) then
        ! Top-East: x, z
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) < origin(2)) then
        ! South-West: -x,-y
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(2)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2)) then
        ! North-West: -x, y
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2)) then
        ! South-East: x, -y
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) > box_max(2)) then
        ! North-East: x, y
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2)) )

      else if (coordLoc(1) < origin(1)) then ! West: -x
        normal = -1_rk
        rad = vec_min(1)

      else if (coordLoc(2) < origin(2)) then ! South: -y
        normal = -1_rk
        rad = vec_min(2)

      else if (coordLoc(3) < origin(3)) then ! Bottom: -z
        normal = -1_rk
        rad = vec_min(3)

      else if (coordLoc(1) > box_max(1)) then ! East: x
        normal = 1_rk
        rad = vec_max(1)

      else if (coordLoc(2) > box_max(2)) then ! North: y
        normal = 1_rk
        rad = vec_max(2)

      else if (coordLoc(3) > box_max(3)) then ! Top: z
        normal = 1_rk
        rad = vec_max(3)
      end if

      proj_len1 = rad*normal
      proj_len2 = (me%thickness*normal - rad)*normal

      if (proj_len1 > 0 .and. proj_len2 > 0) then
        sigma = const_fac * proj_len2 * (proj_len1**4)
        res(i) = sigma*me%dampFactor
      else if (proj_len2 < 0) then ! If coord is beyond thickness
        res(i) = me%dampFactor
      else
        res(i) = 0.0_rk
      end if
    end do

  end function spongeLayer_box_sharpCornerPolyn5_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the box shape sponge layer from 
  !! coord for polynomial n5.
  !! Sponge profile:
  !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$
  !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
  !!
  !! Profile is taken from:
  !! Xu, Hui; Sagaut, Pierre (2013): Analysis of the absorbing layers for the 
  !! weakly-compressible lattice Boltzmann methods. In Journal of Computational 
  !! Physics 245, pp. 14-42. DOI: 10.1016/j.jcp.2013.02.051.
  function spongeLayer_box_roundCornerPolyn5_for_coord(me, coord, n) &
    & result(res)
    ! --------------------------------------------------------------------------
    !> Spatial sponge layer to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(3), extent(3), coordLoc(3), normal
    real(kind=rk) :: proj_len1, proj_len2
    real(kind=rk) :: vec_corMin(3), vec_corMax(3)
    real(kind=rk) :: vec_corMinSqr(3), vec_corMaxSqr(3) 
    real(kind=rk) :: rad, const_fac, box_max(3)
    real(kind=rk) :: corInRad, corOutRad, corMin(3), corMax(3)
    logical :: isCorner
    ! --------------------------------------------------------------------------
    origin(:) = me%origin
    extent(:) = me%extent
    box_max(:) = origin(:) + extent(:)
    corInRad = me%corner_radius
    corOutRad = corInRad + me%thickness
    corMin(:) = origin(:) + corInRad
    corMax(:) = box_max(:) - corInRad

    const_fac = 3125_rk/(256_rk*me%thickness**5)

    do i = 1,n
      coordLoc = coord(i,:)
      vec_corMin(:) = coordLoc(:) - corMin(:)
      vec_corMax(:) = coordLoc(:) - corMax(:)
      vec_corMinSqr(:) = vec_corMin(:)**2
      vec_corMaxSqr(:) = vec_corMax(:)**2
      proj_len1 = 0_rk
      proj_len2 = 0_rk
      normal = 1_rk
      rad = 0_rk
      isCorner = .true.

      ! Bottom-South-West: -x,-y,-z
      if (all(coordLoc(:) < corMin(:))) then
        rad = sqrt( vec_corMinSqr(1) + vec_corMinSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) < corMin(2) &
        & .and. coordLoc(3) < corMin(3)) then ! Bottom-South-East: x, -y, -z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMinSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) > corMax(2) &
        & .and. coordLoc(3) < corMin(3)) then ! Bottom-North-West: -x, y, -z
        rad = sqrt( vec_corMinSqr(1) + vec_corMaxSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) > corMax(2) &
        & .and. coordLoc(3) < corMin(3)) then ! Bottom-North-East: x, y, -z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMaxSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) < corMin(2) &
        & .and. coordLoc(3) > corMax(3)) then ! Top-South-West: -x, -y, z
        rad = sqrt( vec_corMinSqr(1) + vec_corMinSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) < corMin(2) &
        & .and. coordLoc(3) > corMax(3)) then ! Top-South-East: x, -y, z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMinSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) > corMax(2) &
        & .and. coordLoc(3) > corMax(3)) then ! Top-North-West: -x, y, z
        rad = sqrt( vec_corMinSqr(1) + vec_corMaxSqr(2) + vec_corMaxSqr(3))

      else if (all(coordLoc(:) > corMax(:))) then ! Top-North-East: x, y, z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMaxSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(2) < corMin(2) .and. coordLoc(3) < corMin(3)) then
        ! Botton-South: -y,-z
        rad = sqrt( vec_corMinSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(2) < corMin(2) .and. coordLoc(3) > corMax(3)) then
        ! Top-South: -y, z
        rad = sqrt( vec_corMinSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(2) > corMax(2) .and. coordLoc(3) < corMin(3)) then
        ! Bottom-North: y, -z
        rad = sqrt( vec_corMaxSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(2) > corMax(2) .and. coordLoc(3) > corMax(3)) then
        ! Top-North: y, z
        rad = sqrt( vec_corMaxSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(3) < corMin(3)) then
        ! Botton-West: -x,-z
        rad = sqrt( vec_corMinSqr(1) + vec_corMinSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(3) < corMin(3)) then
        ! Bottom-East: x,-z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMinSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(3) > corMax(3)) then
        ! Top-West: -x, z
        rad = sqrt( vec_corMinSqr(1) + vec_corMaxSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(3) > corMax(3)) then
        ! Top-East: x, z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMaxSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) < corMin(2)) then
        ! South-West: -x,-y
        rad = sqrt(vec_corMinSqr(1) + vec_corMinSqr(2))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) > corMax(2)) then
        ! North-West: -x, y
        rad = sqrt(vec_corMinSqr(1) + vec_corMaxSqr(2))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) < corMin(2)) then
        ! South-East: x, -y
        rad = sqrt(vec_corMaxSqr(1) + vec_corMinSqr(2))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) > corMax(2)) then
        ! North-East: x, y
        rad = sqrt(vec_corMaxSqr(1) + vec_corMaxSqr(2))

      else
        isCorner = .false.
        if (coordLoc(1) < origin(1)) then ! West: -x
          normal = -1_rk
          rad = coordLoc(1) - origin(1)

        else if (coordLoc(2) < origin(2)) then ! South: -y
          normal = -1_rk
          rad = coordLoc(2) - origin(2)

        else if (coordLoc(3) < origin(3)) then ! Bottom: -z
          normal = -1_rk
          rad = coordLoc(3) - origin(3)

        else if (coordLoc(1) > box_max(1)) then ! East: x
          normal = 1_rk
          rad = coordLoc(1) - box_max(1)

        else if (coordLoc(2) > box_max(2)) then ! North: y
          normal = 1_rk
          rad = coordLoc(2) - box_max(2)

        else if (coordLoc(3) > box_max(3)) then ! Top: z
          normal = 1_rk
          rad = coordLoc(3) - box_max(3)
        end if
      end if

      if (isCorner) then
        proj_len1 = rad - corInRad
        proj_len2 = corOutRad - rad
      else
        proj_len1 = rad*normal
        proj_len2 = (me%thickness*normal - rad)*normal
      end if

      if (proj_len1 > 0 .and. proj_len2 > 0) then
        sigma = const_fac * proj_len2 * (proj_len1**4)
        res(i) = sigma*me%dampFactor
      else if (proj_len2 < 0) then ! If coord is beyond thickness
        res(i) = me%dampFactor
      else
        res(i) = 0.0_rk
      end if
    end do

  end function spongeLayer_box_roundCornerPolyn5_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the box shape sponge layer from 
  !! coord for polynomial n6.
  !! Sponge profile:
  !! $$ \sigma(x) = \sigma_A*729*(L+x0-x)^2*(x-x0)^4)/(16*(L)^6 $$
  !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
  function spongeLayer_box_sharpCornerPolyn6_for_coord(me, coord, n) &
    & result(res)
    ! --------------------------------------------------------------------------
    !> Spatial sponge layer to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(3), extent(3), coordLoc(3), normal
    real(kind=rk) :: proj_len1, proj_len2
    real(kind=rk) :: vec_min(3), vec_max(3), vec_minSqr(3), vec_maxSqr(3) 
    real(kind=rk) :: rad, const_fac, box_max(3)
    ! --------------------------------------------------------------------------
    origin(:) = me%origin
    extent(:) = me%extent
    box_max(:) = origin(:) + extent(:)

    const_fac = 729_rk/(16_rk*me%thickness**6)

    do i = 1,n
      coordLoc = coord(i,:)
      vec_min(:) = coordLoc(:) - origin(:)
      vec_max(:) = coordLoc(:) - box_max(:)
      vec_minSqr(:) = vec_min(:)**2
      vec_maxSqr(:) = vec_max(:)**2
      proj_len1 = 0_rk
      proj_len2 = 0_rk
      normal = 1_rk
      rad = 0_rk

      ! Bottom-South-West: -x,-y,-z
      if (all(coordLoc(:) < origin(:))) then
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(2), vec_minSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2) &
        & .and. coordLoc(3) < origin(3)) then ! Bottom-South-East: x, -y, -z
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2), vec_minSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2) &
        & .and. coordLoc(3) < origin(3)) then ! Bottom-North-West: -x, y, -z
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2), vec_minSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) > box_max(2) &
        & .and. coordLoc(3) < origin(3)) then ! Bottom-North-East: x, y, -z
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2), vec_minSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) < origin(2) &
        & .and. coordLoc(3) > box_max(3)) then ! Top-South-West: -x, -y, z
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2) &
        & .and. coordLoc(3) > box_max(3)) then ! Top-South-East: x, -y, z
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2) &
        & .and. coordLoc(3) > box_max(3)) then ! Top-North-West: -x, y, z
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2), vec_maxSqr(3)) )

      else if (all(coordLoc(:) > box_max(:))) then ! Top-North-East: x, y, z
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(2) < origin(2) .and. coordLoc(3) < origin(3)) then
        ! Botton-South: -y,-z
        rad = sqrt( max(vec_minSqr(2), vec_minSqr(3)) )

      else if (coordLoc(2) < origin(2) .and. coordLoc(3) > box_max(3)) then
        ! Top-South: -y, z
        rad = sqrt( max(vec_minSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(2) > box_max(2) .and. coordLoc(3) < origin(3)) then
        ! Botom-North: -y, z
        rad = sqrt( max(vec_maxSqr(2), vec_minSqr(3)) )

      else if (coordLoc(2) > box_max(2) .and. coordLoc(3) > box_max(3)) then
        ! Top-North: y, z
        rad = sqrt( max(vec_maxSqr(2), vec_maxSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(3) < origin(3)) then
        ! Botton-West: -x,-z
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(3) < origin(3)) then
        ! Bottom-East: x,-z
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(3) > box_max(3)) then
        ! Top-West: -x, z
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(3)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(3) > box_max(3)) then
        ! Top-East: x, z
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(3)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) < origin(2)) then
        ! South-West: -x,-y
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(2)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2)) then
        ! North-West: -x, y
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2)) then
        ! South-East: x, -y
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) > box_max(2)) then
        ! North-East: x, y
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2)) )

      else if (coordLoc(1) < origin(1)) then ! West: -x
        normal = -1_rk
        rad = vec_min(1)

      else if (coordLoc(2) < origin(2)) then ! South: -y
        normal = -1_rk
        rad = vec_min(2)

      else if (coordLoc(3) < origin(3)) then ! Bottom: -z
        normal = -1_rk
        rad = vec_min(3)

      else if (coordLoc(1) > box_max(1)) then ! East: x
        normal = 1_rk
        rad = vec_max(1)

      else if (coordLoc(2) > box_max(2)) then ! North: y
        normal = 1_rk
        rad = vec_max(2)

      else if (coordLoc(3) > box_max(3)) then ! Top: z
        normal = 1_rk
        rad = vec_max(3)
      end if

      proj_len1 = rad*normal
      proj_len2 = (me%thickness*normal - rad)*normal

      if (proj_len1 > 0) then
        sigma = const_fac * proj_len2**2 * (proj_len1**4)
        res(i) = min(1.0_rk, sigma) * me%dampFactor
      else
        res(i) = 0.0_rk
      end if
    end do

  end function spongeLayer_box_sharpCornerPolyn6_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the box shape sponge layer from 
  !! coord for polynomial n6.
  !! Sponge profile:
  !! $$ \sigma(x) = \sigma_A*729*(L+x0-x)^2*(x-x0)^4)/(16*(L)^6 $$
  !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
  function spongeLayer_box_roundCornerPolyn6_for_coord(me, coord, n) &
    & result(res)
    ! --------------------------------------------------------------------------
    !> Spatial sponge layer to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(3), extent(3), coordLoc(3), normal
    real(kind=rk) :: proj_len1, proj_len2
    real(kind=rk) :: vec_corMin(3), vec_corMax(3)
    real(kind=rk) :: vec_corMinSqr(3), vec_corMaxSqr(3) 
    real(kind=rk) :: rad, const_fac, box_max(3)
    real(kind=rk) :: corInRad, corOutRad, corMin(3), corMax(3)
    logical :: isCorner
    ! --------------------------------------------------------------------------
    origin(:) = me%origin
    extent(:) = me%extent
    box_max(:) = origin(:) + extent(:)
    corInRad = me%corner_radius
    corOutRad = corInRad + me%thickness
    corMin(:) = origin(:) + corInRad
    corMax(:) = box_max(:) - corInRad

    const_fac = 729_rk/(16_rk*me%thickness**6)

    do i = 1,n
      coordLoc = coord(i,:)
      vec_corMin(:) = coordLoc(:) - corMin(:)
      vec_corMax(:) = coordLoc(:) - corMax(:)
      vec_corMinSqr(:) = vec_corMin(:)**2
      vec_corMaxSqr(:) = vec_corMax(:)**2
      proj_len1 = 0_rk
      proj_len2 = 0_rk
      normal = 1_rk
      rad = 0_rk
      isCorner = .true.

      ! Bottom-South-West: -x,-y,-z
      if (all(coordLoc(:) < corMin(:))) then
        rad = sqrt( vec_corMinSqr(1) + vec_corMinSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) < corMin(2) &
        & .and. coordLoc(3) < corMin(3)) then ! Bottom-South-East: x, -y, -z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMinSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) > corMax(2) &
        & .and. coordLoc(3) < corMin(3)) then ! Bottom-North-West: -x, y, -z
        rad = sqrt( vec_corMinSqr(1) + vec_corMaxSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) > corMax(2) &
        & .and. coordLoc(3) < corMin(3)) then ! Bottom-North-East: x, y, -z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMaxSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) < corMin(2) &
        & .and. coordLoc(3) > corMax(3)) then ! Top-South-West: -x, -y, z
        rad = sqrt( vec_corMinSqr(1) + vec_corMinSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) < corMin(2) &
        & .and. coordLoc(3) > corMax(3)) then ! Top-South-East: x, -y, z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMinSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) > corMax(2) &
        & .and. coordLoc(3) > corMax(3)) then ! Top-North-West: -x, y, z
        rad = sqrt( vec_corMinSqr(1) + vec_corMaxSqr(2) + vec_corMaxSqr(3))

      else if (all(coordLoc(:) > corMax(:))) then ! Top-North-East: x, y, z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMaxSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(2) < corMin(2) .and. coordLoc(3) < corMin(3)) then
        ! Botton-South: -y,-z
        rad = sqrt( vec_corMinSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(2) < corMin(2) .and. coordLoc(3) > corMax(3)) then
        ! Top-South: -y, z
        rad = sqrt( vec_corMinSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(2) > corMax(2) .and. coordLoc(3) < corMin(3)) then
        ! Bottom-North: y, -z
        rad = sqrt( vec_corMaxSqr(2) + vec_corMinSqr(3))

      else if (coordLoc(2) > corMax(2) .and. coordLoc(3) > corMax(3)) then
        ! Top-North: y, z
        rad = sqrt( vec_corMaxSqr(2) + vec_corMaxSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(3) < corMin(3)) then
        ! Botton-West: -x,-z
        rad = sqrt( vec_corMinSqr(1) + vec_corMinSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(3) < corMin(3)) then
        ! Bottom-East: x,-z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMinSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(3) > corMax(3)) then
        ! Top-West: -x, z
        rad = sqrt( vec_corMinSqr(1) + vec_corMaxSqr(3))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(3) > corMax(3)) then
        ! Top-East: x, z
        rad = sqrt( vec_corMaxSqr(1) + vec_corMaxSqr(3))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) < corMin(2)) then
        ! South-West: -x,-y
        rad = sqrt(vec_corMinSqr(1) + vec_corMinSqr(2))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) > corMax(2)) then
        ! North-West: -x, y
        rad = sqrt(vec_corMinSqr(1) + vec_corMaxSqr(2))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) < corMin(2)) then
        ! South-East: x, -y
        rad = sqrt(vec_corMaxSqr(1) + vec_corMinSqr(2))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) > corMax(2)) then
        ! North-East: x, y
        rad = sqrt(vec_corMaxSqr(1) + vec_corMaxSqr(2))

      else
        isCorner = .false.
        if (coordLoc(1) < origin(1)) then ! West: -x
          normal = -1_rk
          rad = coordLoc(1) - origin(1)

        else if (coordLoc(2) < origin(2)) then ! South: -y
          normal = -1_rk
          rad = coordLoc(2) - origin(2)

        else if (coordLoc(3) < origin(3)) then ! Bottom: -z
          normal = -1_rk
          rad = coordLoc(3) - origin(3)

        else if (coordLoc(1) > box_max(1)) then ! East: x
          normal = 1_rk
          rad = coordLoc(1) - box_max(1)

        else if (coordLoc(2) > box_max(2)) then ! North: y
          normal = 1_rk
          rad = coordLoc(2) - box_max(2)

        else if (coordLoc(3) > box_max(3)) then ! Top: z
          normal = 1_rk
          rad = coordLoc(3) - box_max(3)
        end if
      end if

      if (isCorner) then
        proj_len1 = rad - corInRad
        proj_len2 = corOutRad - rad
      else
        proj_len1 = rad*normal
        proj_len2 = (me%thickness*normal - rad)*normal
      end if

      if (proj_len1 > 0) then
        sigma = const_fac * proj_len2**2 * (proj_len1**4)
        res(i) = min(1.0_rk, sigma) * me%dampFactor
      else
        res(i) = 0.0_rk
      end if
    end do

  end function spongeLayer_box_roundCornerPolyn6_for_coord
  ! -------------------------------------------------------------------------- !


  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the spongelayer and fills up
  !!  the res with the target state
  function spongelayer_box_vector_for_coord(me, nComp, coord, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> Number of entrys in each array
    integer, intent(in) :: ncomp
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n,ncomp)
    ! --------------------------------------------------------------------------
    integer :: i
    ! --------------------------------------------------------------------------
    res(:, 1) = spongeLayer_box_scalar_for_coord(me, coord, n)

    if (ncomp > 1) then
      do i = 1,n
        res(i,2:) = me%targetState(:)
      end do
    end if

  end function spongelayer_box_vector_for_coord
  ! -------------------------------------------------------------------------- !


  ! -------------------------------------------------------------------------- !
  !> This function returns the sigma for the spongelayer from treeids
  function spongelayer_box_scalar_for_treeIDs(me, treeids, tree, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    select case (trim(me%dampProfile))
    case ('linear', 'exponential')
      res(:) = spongeLayer_box_expon_for_treeIDs(me, treeIDs, tree, n)
    case ('polynomial_n5')
      if (me%rounded_corner) then
        res(:) = spongeLayer_box_roundCornerPolyn5_for_treeIDs(me, treeIDs, &
          &                                                    tree, n)
      else
        res(:) = spongeLayer_box_sharpCornerPolyn5_for_treeIDs(me, treeIDs, &
          &                                                    tree, n)
      end if
    case ('polynomial_n6')
      if (me%rounded_corner) then
        res(:) = spongeLayer_box_roundCornerPolyn6_for_treeIDs(me, treeIDs, &
          &                                                    tree, n)
      else
        res(:) = spongeLayer_box_sharpCornerPolyn6_for_treeIDs(me, treeIDs, &
          &                                                    tree, n)
      end if
    end select

  end function spongeLayer_box_scalar_for_treeIDs
  ! -------------------------------------------------------------------------- !


  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the box shape spongelayer fom treeid
  !! for exponential profile
  function spongelayer_box_expon_for_treeIDs(me, treeIDs, tree, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: res_i(1), coord(3,1)
    ! --------------------------------------------------------------------------
    do i = 1,n
      !barycentric coordinate
      coord(:,1) = tem_BaryOfId( tree, treeIds(i) )
      res_i = spongeLayer_box_expon_for_coord(me, coord, 1)
      res(i) = res_i(1)
    end do

  end function spongelayer_box_expon_for_treeIDs
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the box shape spongelayer fom treeid
  !! for polynomial n5 profile
  function spongelayer_box_sharpCornerPolyn5_for_treeIDs(me, treeIDs, tree, n) &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: res_i(1), coord(3,1)
    ! --------------------------------------------------------------------------
    do i = 1,n
      !barycentric coordinate
      coord(:,1) = tem_BaryOfId( tree, treeIds(i) )
      res_i = spongeLayer_box_sharpCornerPolyn5_for_coord(me, coord, 1)
      res(i) = res_i(1)
    end do

  end function spongelayer_box_sharpCornerPolyn5_for_treeIDs
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the box shape spongelayer fom treeid
  !! for polynomial n5 profile
  function spongelayer_box_roundCornerPolyn5_for_treeIDs(me, treeIDs, tree, n) &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: res_i(1), coord(3,1)
    ! --------------------------------------------------------------------------
    do i = 1,n
      !barycentric coordinate
      coord(:,1) = tem_BaryOfId( tree, treeIds(i) )
      res_i = spongeLayer_box_roundCornerPolyn5_for_coord(me, coord, 1)
      res(i) = res_i(1)
    end do

  end function spongelayer_box_roundCornerPolyn5_for_treeIDs
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the box shape spongelayer fom treeid
  !! for polynomial n5 profile
  function spongelayer_box_sharpCornerPolyn6_for_treeIDs(me, treeIDs, tree, n) &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: res_i(1), coord(3,1)
    ! --------------------------------------------------------------------------
    do i = 1,n
      !barycentric coordinate
      coord(:,1) = tem_BaryOfId( tree, treeIds(i) )
      res_i = spongeLayer_box_sharpCornerPolyn6_for_coord(me, coord, 1)
      res(i) = res_i(1)
    end do

  end function spongelayer_box_sharpCornerPolyn6_for_treeIDs
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the box shape spongelayer fom treeid
  !! for polynomial n5 profile
  function spongelayer_box_roundCornerPolyn6_for_treeIDs(me, treeIDs, tree, n) &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: res_i(1), coord(3,1)
    ! --------------------------------------------------------------------------
    do i = 1,n
      !barycentric coordinate
      coord(:,1) = tem_BaryOfId( tree, treeIds(i) )
      res_i = spongeLayer_box_roundCornerPolyn6_for_coord(me, coord, 1)
      res(i) = res_i(1)
    end do

  end function spongelayer_box_roundCornerPolyn6_for_treeIDs
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the spongelayer and fills up
  !!  the res with the target state
  function spongelayer_box_vector_for_treeIDs(me, ncomp, treeids, tree, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> Number of entrys in each array
    integer, intent(in) :: ncomp
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n,ncomp)
    ! --------------------------------------------------------------------------
    integer :: i
    ! --------------------------------------------------------------------------
    res(:, 1) = spongeLayer_box_scalar_for_treeIDs(me, treeids, tree, n)

    if (ncomp > 1) then
      do i = 1,n
        res(i,2:) = me%targetState(:)
      end do
    end if

  end function spongelayer_box_vector_for_treeIDs
  ! -------------------------------------------------------------------------- !


  ! -------------------------------------------------------------------------- !
  !
  !                            SPONGE LAYER BOX 2D
  !
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function returns the sigma for the 2d box shape spongelayer
  function spongelayer_box2d_scalar_for_coord(me, coord, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    select case (trim(me%dampProfile))
    case ('linear', 'exponential')
      res(:) = spongeLayer_box2d_expon_for_coord(me, coord, n)
    case ('polynomial_n5')
      if (me%rounded_corner) then
        res(:) = spongeLayer_box2d_roundCornerPolyn5_for_coord(me, coord, n)
      else
        res(:) = spongeLayer_box2d_sharpCornerPolyn5_for_coord(me, coord, n)
      end if
    case ('polynomial_n6')
      if (me%rounded_corner) then
        res(:) = spongeLayer_box2d_roundCornerPolyn6_for_coord(me, coord, n)
      else
        res(:) = spongeLayer_box2d_sharpCornerPolyn6_for_coord(me, coord, n)
      end if
    end select

  end function spongelayer_box2d_scalar_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function returns sigma for the box shape spongelayer for coord for
  !! exponential profile.
  function spongelayer_box2d_expon_for_coord(me, coord, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(2), extent(2), box_max(2), proj_len
    real(kind=rk) :: coordLoc(2), normal, vec_min(2), vec_max(2)
    real(kind=rk) :: rad, vec_minSqr(2), vec_maxSqr(2)
    ! --------------------------------------------------------------------------
    origin(:) = me%origin(1:2)
    extent(:) = me%extent(1:2)
    box_max(:) = origin(:) + extent(:)

    do i = 1,n
      coordLoc = coord(i,1:2)
      vec_min(:) = coordLoc(:) - origin(:)
      vec_max(:) = coordLoc(:) - box_max(:)
      vec_minSqr(:) = vec_min(:)**2
      vec_maxSqr(:) = vec_max(:)**2
      normal = 1.0_rk
      rad = 0.0_rk

      if (coordLoc(1) < origin(1) .and. coordLoc(2) < origin(2)) then
        ! South-West: -x,-y
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(2)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2)) then
        ! North-West: -x, y
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2)) then
        ! South-East: x, -y
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) > box_max(2)) then
        ! North-East: x, y
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2)) )

      else if (coordLoc(1) < origin(1)) then ! West: -x
        normal = -1_rk
        rad = vec_min(1)

      else if (coordLoc(2) < origin(2)) then ! South: -y
        normal = -1_rk
        rad = vec_min(2)

      else if (coordLoc(1) > box_max(1)) then ! East: x
        normal = 1_rk
        rad = vec_max(1)

      else if (coordLoc(2) > box_max(2)) then ! North: y
        normal = 1_rk
        rad = vec_max(2)
      end if

      proj_len = rad*normal/me%thickness
      sigma = me%dampFactor*((proj_len)**me%dampExponent)
      if (proj_len > 0 .and. proj_len < 1) then
        res(i) = sigma
      else if (proj_len > 1) then
        res(i) = me%dampFactor
      else
        res(i) = 0.0_rk
      end if

    end do

  end function spongelayer_box2d_expon_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the 2d box shape sponge layer from 
  !! coord for polynomial n5.
  !! Sponge profile:
  !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$
  !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
  !!
  !! Profile is taken from:
  !! Xu, Hui; Sagaut, Pierre (2013): Analysis of the absorbing layers for the 
  !! weakly-compressible lattice Boltzmann methods. In Journal of Computational 
  !! Physics 245, pp. 14-42. DOI: 10.1016/j.jcp.2013.02.051.
  function spongeLayer_box2d_sharpCornerPolyn5_for_coord(me, coord, n) &
    & result(res)
    ! --------------------------------------------------------------------------
    !> Spatial sponge layer to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(2), extent(2), coordLoc(2), normal
    real(kind=rk) :: proj_len1, proj_len2
    real(kind=rk) :: vec_min(2), vec_max(2), vec_minSqr(2), vec_maxSqr(2) 
    real(kind=rk) :: rad, const_fac, box_max(2)
    ! --------------------------------------------------------------------------
    origin(:) = me%origin(1:2)
    extent(:) = me%extent(1:2)
    box_max(:) = origin(:) + extent(:)

    const_fac = 3125_rk/(256_rk*me%thickness**5)

    do i = 1,n
      coordLoc = coord(i,1:2)
      vec_min(:) = coordLoc(:) - origin(:)
      vec_max(:) = coordLoc(:) - box_max(:)
      vec_minSqr(:) = vec_min(:)**2
      vec_maxSqr(:) = vec_max(:)**2
      proj_len1 = 0_rk
      proj_len2 = 0_rk
      normal = 1_rk
      rad = 0_rk

      if (coordLoc(1) < origin(1) .and. coordLoc(2) < origin(2)) then
        ! South-West: -x,-y
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(2)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2)) then
        ! North-West: -x, y
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2)) then
        ! South-East: x, -y
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) > box_max(2)) then
        ! North-East: x, y
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2)) )

      else if (coordLoc(1) < origin(1)) then ! West: -x
        normal = -1_rk
        rad = vec_min(1)

      else if (coordLoc(2) < origin(2)) then ! South: -y
        normal = -1_rk
        rad = vec_min(2)

      else if (coordLoc(1) > box_max(1)) then ! East: x
        normal = 1_rk
        rad = vec_max(1)

      else if (coordLoc(2) > box_max(2)) then ! North: y
        normal = 1_rk
        rad = vec_max(2)
      end if

      proj_len1 = rad*normal
      proj_len2 = (me%thickness*normal - rad)*normal

      if (proj_len1 > 0 .and. proj_len2 > 0) then
        sigma = const_fac * proj_len2 * (proj_len1**4)
        res(i) = sigma*me%dampFactor
      else if (proj_len2 < 0) then ! If coord is beyond thickness
        res(i) = me%dampFactor
      else
        res(i) = 0.0_rk
      end if
    end do

  end function spongeLayer_box2d_sharpCornerPolyn5_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the 2d box shape sponge layer from 
  !! coord for polynomial n5.
  !! Sponge profile:
  !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$
  !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
  !!
  !! Profile is taken from:
  !! Xu, Hui; Sagaut, Pierre (2013): Analysis of the absorbing layers for the 
  !! weakly-compressible lattice Boltzmann methods. In Journal of Computational 
  !! Physics 245, pp. 14-42. DOI: 10.1016/j.jcp.2013.02.051.
  function spongeLayer_box2d_roundCornerPolyn5_for_coord(me, coord, n) &
    & result(res)
    ! --------------------------------------------------------------------------
    !> Spatial sponge layer to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(2), extent(2), coordLoc(2), normal
    real(kind=rk) :: proj_len1, proj_len2
    real(kind=rk) :: vec_corMin(2), vec_corMax(2)
    real(kind=rk) :: vec_corMinSqr(2), vec_corMaxSqr(2) 
    real(kind=rk) :: rad, const_fac, box_max(2)
    real(kind=rk) :: corInRad, corOutRad, corMin(2), corMax(2)
    logical :: isCorner
    ! --------------------------------------------------------------------------
    origin(:) = me%origin(1:2)
    extent(:) = me%extent(1:2)
    box_max(:) = origin(:) + extent(:)
    corInRad = me%corner_radius
    corOutRad = corInRad + me%thickness
    corMin(:) = origin(:) + corInRad
    corMax(:) = box_max(:) - corInRad

    const_fac = 3125_rk/(256_rk*me%thickness**5)

    do i = 1,n
      coordLoc = coord(i,1:2)
      vec_corMin(:) = coordLoc(:) - corMin(:)
      vec_corMax(:) = coordLoc(:) - corMax(:)
      vec_corMinSqr(:) = vec_corMin(:)**2
      vec_corMaxSqr(:) = vec_corMax(:)**2
      proj_len1 = 0_rk
      proj_len2 = 0_rk
      normal = 1_rk
      rad = 0_rk
      isCorner = .true.

      if (coordLoc(1) < corMin(1) .and. coordLoc(2) < corMin(2)) then
        ! South-West: -x,-y
        rad = sqrt(vec_corMinSqr(1) + vec_corMinSqr(2))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) > corMax(2)) then
        ! North-West: -x, y
        rad = sqrt(vec_corMinSqr(1) + vec_corMaxSqr(2))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) < corMin(2)) then
        ! South-East: x, -y
        rad = sqrt(vec_corMaxSqr(1) + vec_corMinSqr(2))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) > corMax(2)) then
        ! North-East: x, y
        rad = sqrt(vec_corMaxSqr(1) + vec_corMaxSqr(2))

      else
        isCorner = .false.
        if (coordLoc(1) < origin(1)) then ! West: -x
          normal = -1_rk
          rad = coordLoc(1) - origin(1)

        else if (coordLoc(2) < origin(2)) then ! South: -y
          normal = -1_rk
          rad = coordLoc(2) - origin(2)

        else if (coordLoc(1) > box_max(1)) then ! East: x
          normal = 1_rk
          rad = coordLoc(1) - box_max(1)

        else if (coordLoc(2) > box_max(2)) then ! North: y
          normal = 1_rk
          rad = coordLoc(2) - box_max(2)
        end if
      end if

      if (isCorner) then
        proj_len1 = rad - corInRad
        proj_len2 = corOutRad - rad
      else
        proj_len1 = rad*normal
        proj_len2 = (me%thickness*normal - rad)*normal
      end if

      if (proj_len1 > 0 .and. proj_len2 > 0) then
        sigma = const_fac * proj_len2 * (proj_len1**4)
        res(i) = sigma*me%dampFactor
      else if (proj_len2 < 0) then ! If coord is beyond thickness
        res(i) = me%dampFactor
      else
        res(i) = 0.0_rk
      end if
    end do

  end function spongeLayer_box2d_roundCornerPolyn5_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the 2d box shape sponge layer from 
  !! coord for polynomial n6.
  !! Sponge profile:
  !! $$ \sigma(x) = \sigma_A*729*(L+x0-x)^2*(x-x0)^4)/(16*(L)^6 $$
  !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
  function spongeLayer_box2d_sharpCornerPolyn6_for_coord(me, coord, n) &
    & result(res)
    ! --------------------------------------------------------------------------
    !> Spatial sponge layer to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(2), extent(2), coordLoc(2), normal
    real(kind=rk) :: proj_len1, proj_len2
    real(kind=rk) :: vec_min(2), vec_max(2), vec_minSqr(2), vec_maxSqr(2) 
    real(kind=rk) :: rad, const_fac, box_max(2)
    ! --------------------------------------------------------------------------
    origin(:) = me%origin(1:2)
    extent(:) = me%extent(1:2)
    box_max(:) = origin(:) + extent(:)

    const_fac = 729_rk/(16_rk*me%thickness**6)

    do i = 1,n
      coordLoc = coord(i,1:2)
      vec_min(:) = coordLoc(:) - origin(:)
      vec_max(:) = coordLoc(:) - box_max(:)
      vec_minSqr(:) = vec_min(:)**2
      vec_maxSqr(:) = vec_max(:)**2
      proj_len1 = 0_rk
      proj_len2 = 0_rk
      normal = 1_rk
      rad = 0_rk

      if (coordLoc(1) < origin(1) .and. coordLoc(2) < origin(2)) then
        ! South-West: -x,-y
        rad = sqrt( max(vec_minSqr(1), vec_minSqr(2)) )

      else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2)) then
        ! North-West: -x, y
        rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2)) then
        ! South-East: x, -y
        rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2)) )

      else if (coordLoc(1) > box_max(1) .and. coordLoc(2) > box_max(2)) then
        ! North-East: x, y
        rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2)) )

      else if (coordLoc(1) < origin(1)) then ! West: -x
        normal = -1_rk
        rad = vec_min(1)

      else if (coordLoc(2) < origin(2)) then ! South: -y
        normal = -1_rk
        rad = vec_min(2)

      else if (coordLoc(1) > box_max(1)) then ! East: x
        normal = 1_rk
        rad = vec_max(1)

      else if (coordLoc(2) > box_max(2)) then ! North: y
        normal = 1_rk
        rad = vec_max(2)
      end if

      proj_len1 = rad*normal
      proj_len2 = (me%thickness*normal - rad)*normal

      if (proj_len1 > 0) then
        sigma = const_fac * proj_len2**2 * (proj_len1**4)
        res(i) = min(1.0_rk, sigma) * me%dampFactor
      else
        res(i) = 0.0_rk
      end if
    end do

  end function spongeLayer_box2d_sharpCornerPolyn6_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the 2d box shape sponge layer from 
  !! coord for polynomial n6.
  !! Sponge profile:
  !! $$ \sigma(x) = \sigma_A*729*(L+x0-x)^2*(x-x0)^4)/(16*(L)^6 $$
  !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
  function spongeLayer_box2d_roundCornerPolyn6_for_coord(me, coord, n) &
    & result(res)
    ! --------------------------------------------------------------------------
    !> Spatial sponge layer to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(2), extent(2), coordLoc(2), normal
    real(kind=rk) :: proj_len1, proj_len2
    real(kind=rk) :: vec_corMin(2), vec_corMax(2)
    real(kind=rk) :: vec_corMinSqr(2), vec_corMaxSqr(2) 
    real(kind=rk) :: rad, const_fac, box_max(2)
    real(kind=rk) :: corInRad, corOutRad, corMin(2), corMax(2)
    logical :: isCorner
    ! --------------------------------------------------------------------------
    origin(:) = me%origin(1:2)
    extent(:) = me%extent(1:2)
    box_max(:) = origin(:) + extent(:)
    corInRad = me%corner_radius
    corOutRad = corInRad + me%thickness
    corMin(:) = origin(:) + corInRad
    corMax(:) = box_max(:) - corInRad

    const_fac = 729_rk/(16_rk*me%thickness**6)

    do i = 1,n
      coordLoc = coord(i,1:2)
      vec_corMin(:) = coordLoc(:) - corMin(:)
      vec_corMax(:) = coordLoc(:) - corMax(:)
      vec_corMinSqr(:) = vec_corMin(:)**2
      vec_corMaxSqr(:) = vec_corMax(:)**2
      proj_len1 = 0_rk
      proj_len2 = 0_rk
      normal = 1_rk
      rad = 0_rk
      isCorner = .true.

      if (coordLoc(1) < corMin(1) .and. coordLoc(2) < corMin(2)) then
        ! South-West: -x,-y
        rad = sqrt(vec_corMinSqr(1) + vec_corMinSqr(2))

      else if (coordLoc(1) < corMin(1) .and. coordLoc(2) > corMax(2)) then
        ! North-West: -x, y
        rad = sqrt(vec_corMinSqr(1) + vec_corMaxSqr(2))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) < corMin(2)) then
        ! South-East: x, -y
        rad = sqrt(vec_corMaxSqr(1) + vec_corMinSqr(2))

      else if (coordLoc(1) > corMax(1) .and. coordLoc(2) > corMax(2)) then
        ! North-East: x, y
        rad = sqrt(vec_corMaxSqr(1) + vec_corMaxSqr(2))

      else
        isCorner = .false.
        if (coordLoc(1) < origin(1)) then ! West: -x
          normal = -1_rk
          rad = coordLoc(1) - origin(1)

        else if (coordLoc(2) < origin(2)) then ! South: -y
          normal = -1_rk
          rad = coordLoc(2) - origin(2)

        else if (coordLoc(1) > box_max(1)) then ! East: x
          normal = 1_rk
          rad = coordLoc(1) - box_max(1)

        else if (coordLoc(2) > box_max(2)) then ! North: y
          normal = 1_rk
          rad = coordLoc(2) - box_max(2)
        end if
      end if

      if (isCorner) then
        proj_len1 = rad - corInRad
        proj_len2 = corOutRad - rad
      else
        proj_len1 = rad*normal
        proj_len2 = (me%thickness*normal - rad)*normal
      end if

      if (proj_len1 > 0) then
        sigma = const_fac * proj_len2**2 * (proj_len1**4)
        res(i) = min(1.0_rk, sigma) * me%dampFactor
      else
        res(i) = 0.0_rk
      end if
    end do

  end function spongeLayer_box2d_roundCornerPolyn6_for_coord
  ! -------------------------------------------------------------------------- !


  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the spongelayer and fills up
  !!  the res with the target state
  function spongelayer_box2d_vector_for_coord(me, nComp, coord, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> Number of entrys in each array
    integer, intent(in) :: ncomp
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n,ncomp)
    ! --------------------------------------------------------------------------
    integer :: i
    ! --------------------------------------------------------------------------
    res(:, 1) = spongeLayer_box2d_scalar_for_coord(me, coord, n)

    if (ncomp > 1) then
      do i = 1,n
        res(i,2:) = me%targetState(:)
      end do
    end if

  end function spongelayer_box2d_vector_for_coord
  ! -------------------------------------------------------------------------- !


  ! -------------------------------------------------------------------------- !
  !> This function returns the sigma for the spongelayer from treeids
  function spongelayer_box2d_scalar_for_treeIDs(me, treeids, tree, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    select case (trim(me%dampProfile))
    case ('linear', 'exponential')
      res(:) = spongeLayer_box2d_expon_for_treeIDs(me, treeIDs, tree, n)
    case ('polynomial_n5')
      if (me%rounded_corner) then
        res(:) = spongeLayer_box2d_roundCornerPolyn5_for_treeIDs(me, treeIDs, &
          &                                                      tree, n)
      else
        res(:) = spongeLayer_box2d_sharpCornerPolyn5_for_treeIDs(me, treeIDs, &
          &                                                      tree, n)
      end if
    case ('polynomial_n6')
      if (me%rounded_corner) then
        res(:) = spongeLayer_box2d_roundCornerPolyn6_for_treeIDs(me, treeIDs, &
          &                                                      tree, n)
      else
        res(:) = spongeLayer_box2d_sharpCornerPolyn6_for_treeIDs(me, treeIDs, &
          &                                                      tree, n)
      end if
    end select

  end function spongeLayer_box2d_scalar_for_treeIDs
  ! -------------------------------------------------------------------------- !


  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the box shape spongelayer fom treeid
  !! for exponential profile
  function spongelayer_box2d_expon_for_treeIDs(me, treeIDs, tree, n)  &
    &                           result(res)
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: res_i(1), coord(3,1)
    ! --------------------------------------------------------------------------
    do i = 1,n
      !barycentric coordinate
      coord(:,1) = tem_BaryOfId( tree, treeIds(i) )
      res_i = spongeLayer_box2d_expon_for_coord(me, coord, 1)
      res(i) = res_i(1)
    end do

  end function spongelayer_box2d_expon_for_treeIDs
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the box shape spongelayer fom treeid
  !! for polynomial n5 profile
  function spongelayer_box2d_sharpCornerPolyn5_for_treeIDs(me, treeIDs, tree, &
    &                                                      n) result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: res_i(1), coord(3,1)
    ! --------------------------------------------------------------------------
    do i = 1,n
      !barycentric coordinate
      coord(:,1) = tem_BaryOfId( tree, treeIds(i) )
      res_i = spongeLayer_box2d_sharpCornerPolyn5_for_coord(me, coord, 1)
      res(i) = res_i(1)
    end do

  end function spongelayer_box2d_sharpCornerPolyn5_for_treeIDs
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the box shape spongelayer fom treeid
  !! for polynomial n5 profile
  function spongelayer_box2d_roundCornerPolyn5_for_treeIDs(me, treeIDs, tree, &
    &                                                      n) result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: res_i(1), coord(3,1)
    ! --------------------------------------------------------------------------
    do i = 1,n
      !barycentric coordinate
      coord(:,1) = tem_BaryOfId( tree, treeIds(i) )
      res_i = spongeLayer_box2d_roundCornerPolyn5_for_coord(me, coord, 1)
      res(i) = res_i(1)
    end do

  end function spongelayer_box2d_roundCornerPolyn5_for_treeIDs
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the box shape spongelayer fom treeid
  !! for polynomial n5 profile
  function spongelayer_box2d_sharpCornerPolyn6_for_treeIDs(me, treeIDs, tree, &
    &                                                      n) result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: res_i(1), coord(3,1)
    ! --------------------------------------------------------------------------
    do i = 1,n
      !barycentric coordinate
      coord(:,1) = tem_BaryOfId( tree, treeIds(i) )
      res_i = spongeLayer_box2d_sharpCornerPolyn6_for_coord(me, coord, 1)
      res(i) = res_i(1)
    end do

  end function spongelayer_box2d_sharpCornerPolyn6_for_treeIDs
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the box shape spongelayer fom treeid
  !! for polynomial n5 profile
  function spongelayer_box2d_roundCornerPolyn6_for_treeIDs(me, treeIDs, tree, &
    &                                                      n) result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: res_i(1), coord(3,1)
    ! --------------------------------------------------------------------------
    do i = 1,n
      !barycentric coordinate
      coord(:,1) = tem_BaryOfId( tree, treeIds(i) )
      res_i = spongeLayer_box2d_roundCornerPolyn6_for_coord(me, coord, 1)
      res(i) = res_i(1)
    end do

  end function spongelayer_box2d_roundCornerPolyn6_for_treeIDs
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the spongelayer and fills up
  !!  the res with the target state
  function spongelayer_box2d_vector_for_treeIDs(me, ncomp, treeids, tree, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_box_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> Number of entrys in each array
    integer, intent(in) :: ncomp
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> treeIds of elements in given level
    integer(kind=long_k), intent(in) :: treeIds(n)
    !> return value
    real(kind=rk) :: res(n,ncomp)
    ! --------------------------------------------------------------------------
    integer :: i
    ! --------------------------------------------------------------------------
    res(:, 1) = spongeLayer_box2d_scalar_for_treeIDs(me, treeids, tree, n)

    if (ncomp > 1) then
      do i = 1,n
        res(i,2:) = me%targetState(:)
      end do
    end if

  end function spongelayer_box2d_vector_for_treeIDs
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !
  !                            SPONGE LAYER RADIAL
  !
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function returns the sigma for the radial viscosity spongelayer
  !! for 2D and 3D
  function spongelayer_radial_scalar_for_coord(me, coord, nDim, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_radial_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> Dimension
    integer, intent(in) :: nDim
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    select case (trim(me%dampProfile))
    case ('linear', 'exponential')
      res(:) = spongeLayer_radial_expon_for_coord(me, coord, nDim, n)
    case ('polynomial_n5')
      res(:) = spongeLayer_radial_polyn5_for_coord(me, coord, nDim, n)
    case ('polynomial_n6')
      res(:) = spongeLayer_radial_polyn6_for_coord(me, coord, nDim, n)
    end select

  end function spongelayer_radial_scalar_for_coord
  ! -------------------------------------------------------------------------- !


  ! -------------------------------------------------------------------------- !
  !> This function returns the sigma for the radial viscosity spongelayer
  !! for 2D and 3D
  function spongelayer_radial_expon_for_coord(me, coord, nDim, n)  &
    &                           result(res)
    ! --------------------------------------------------------------------------
    !> Spacetime function to evaluate
    type(tem_spongeLayer_radial_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> Dimension
    integer, intent(in) :: nDim
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, rad, origin(3), outerRadius, vec(3)
    ! --------------------------------------------------------------------------
    origin(:) = me%origin
    outerRadius = me%radius + me%thickness
    vec = 0.0_rk

    do i = 1, n
      vec(1:nDim) = coord(i, 1:nDim) - origin(1:nDim)
      rad = sqrt( vec(1)**2 + vec(2)**2 + vec(3)**2 )

      if ( rad > me%radius .and. rad < outerRadius ) then
        sigma = me%dampFactor                                        &
          &   * ( (rad-me%radius)/me%thickness )**me%dampExponent
      else if (rad > outerRadius) then
        sigma = me%dampFactor
      else
        sigma = 0.0_rk
      end if

      res(i) = sigma

    end do

  end function spongelayer_radial_expon_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the spongelayer radial from coord.
  !! Sponge profile:
  !! $$ \sigma(x) = \sigma_A*3125*(L+x0-x)*(x-x0)^4)/(256*(L)^5 $$
  !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
  !!
  !! Profile is taken from:
  !! Xu, Hui; Sagaut, Pierre (2013): Analysis of the absorbing layers for the 
  !! weakly-compressible lattice Boltzmann methods. In Journal of Computational 
  !! Physics 245, pp. 14-42. DOI: 10.1016/j.jcp.2013.02.051.
  function spongeLayer_radial_polyn5_for_coord(me, coord, nDim, n) result(res)
    ! --------------------------------------------------------------------------
    !> Spatial sponge layer to evaluate
    type(tem_spongeLayer_radial_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> Dimension
    integer, intent(in) :: nDim
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: rad, outer_radius, sigma, vec(3)
    real(kind=rk) :: const_fac
    ! --------------------------------------------------------------------------
    outer_radius = me%radius + me%thickness
    vec = 0.0_rk
    const_fac = 3125_rk/(256_rk*me%thickness**5)

    do i = 1,n
      vec(1:nDim) = coord(i,1:nDim) - me%origin(1:nDim)
      rad = sqrt( vec(1)**2 + vec(2)**2 + vec(3)**2 )
      if (rad > me%radius .and. rad < outer_radius) then
        sigma = const_fac * (outer_radius - rad) * (rad - me%radius)**4
        res(i) = sigma * me%dampFactor
      else if (rad > outer_radius) then
        res(i) = me%dampFactor
      else
        res(i) = 0.0_rk
      end if
    end do

  end function spongeLayer_radial_polyn5_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the spongelayer radial from coord 
  !! for the polynomial order n6.
  !! Sponge profile:
  !! $$ \sigma(x) = \sigma_A*729*(L+x0-x)^2*(x-x0)^4)/(16*(L)^6 $$
  !! where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
  function spongeLayer_radial_polyn6_for_coord(me, coord, nDim, n) result(res)
    ! --------------------------------------------------------------------------
    !> Spatial sponge layer to evaluate
    type(tem_spongeLayer_radial_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> Dimension
    integer, intent(in) :: nDim
    !> return value
    real(kind=rk) :: res(n)
    ! --------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: rad, outer_radius, sigma, vec(3)
    real(kind=rk) :: const_fac
    ! --------------------------------------------------------------------------
    outer_radius = me%radius + me%thickness
    vec = 0.0_rk
    const_fac = 729_rk/(16_rk*me%thickness**6)

    do i = 1,n
      vec(1:nDim) = coord(i,1:nDim) - me%origin(1:nDim)
      rad = sqrt( vec(1)**2 + vec(2)**2 + vec(3)**2 )
      if (rad > me%radius) then
        sigma = const_fac * (outer_radius - rad)**2 * (rad - me%radius)**4
        res(i) = min(1.0_rk, sigma) * me%dampFactor
      else
        res(i) = 0.0_rk
      end if
    end do

  end function spongeLayer_radial_polyn6_for_coord
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This function calculates the sigma for the radial viscosity spongelayer
  !! for 2D and 3D and fills up rest with target_state.