mus_auxFieldVar_module.f90 Source File


This file depends on

sourcefile~~mus_auxfieldvar_module.f90~~EfferentGraph sourcefile~mus_auxfieldvar_module.f90 mus_auxFieldVar_module.f90 sourcefile~mus_auxfield_module.f90 mus_auxField_module.f90 sourcefile~mus_auxfieldvar_module.f90->sourcefile~mus_auxfield_module.f90 sourcefile~mus_enrtl_module.f90 mus_eNRTL_module.f90 sourcefile~mus_auxfieldvar_module.f90->sourcefile~mus_enrtl_module.f90 sourcefile~mus_connectivity_module.f90 mus_connectivity_module.f90 sourcefile~mus_auxfieldvar_module.f90->sourcefile~mus_connectivity_module.f90 sourcefile~mus_scheme_type_module.f90 mus_scheme_type_module.f90 sourcefile~mus_auxfieldvar_module.f90->sourcefile~mus_scheme_type_module.f90 sourcefile~mus_dervarpos_module.f90 mus_derVarPos_module.f90 sourcefile~mus_auxfieldvar_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_varsys_module.f90 mus_varSys_module.f90 sourcefile~mus_auxfieldvar_module.f90->sourcefile~mus_varsys_module.f90 sourcefile~mus_physics_module.f90 mus_physics_module.f90 sourcefile~mus_auxfieldvar_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_source_var_module.f90 mus_source_var_module.f90 sourcefile~mus_auxfieldvar_module.f90->sourcefile~mus_source_var_module.f90 sourcefile~mus_scheme_header_module.f90 mus_scheme_header_module.f90 sourcefile~mus_auxfieldvar_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_auxfield_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_auxfield_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_scheme_layout_module.f90 mus_scheme_layout_module.f90 sourcefile~mus_connectivity_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_bc_header_module.f90 mus_bc_header_module.f90 sourcefile~mus_connectivity_module.f90->sourcefile~mus_bc_header_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_auxfield_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_source_var_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_field_module.f90 mus_field_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_field_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_interpolate_header_module.f90 mus_interpolate_header_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_interpolate_header_module.f90 sourcefile~mus_nernstplanck_module.f90 mus_nernstPlanck_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_nernstplanck_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_bc_header_module.f90 sourcefile~mus_graddata_module.f90 mus_gradData_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_graddata_module.f90 sourcefile~mus_transport_var_module.f90 mus_transport_var_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_transport_var_module.f90 sourcefile~mus_mixture_module.f90 mus_mixture_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_mixture_module.f90 sourcefile~mus_pdf_module.f90 mus_pdf_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_pdf_module.f90 sourcefile~mus_param_module.f90 mus_param_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_param_module.f90 sourcefile~mus_field_prop_module.f90 mus_field_prop_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_field_prop_module.f90 sourcefile~mus_dervarpos_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_varsys_module.f90->sourcefile~mus_scheme_type_module.f90 sourcefile~mus_varsys_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_geom_module.f90 mus_geom_module.f90 sourcefile~mus_varsys_module.f90->sourcefile~mus_geom_module.f90 sourcefile~mus_source_var_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_source_var_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_source_var_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_field_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_field_module.f90->sourcefile~mus_source_var_module.f90 sourcefile~mus_field_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_field_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_field_module.f90->sourcefile~mus_nernstplanck_module.f90 sourcefile~mus_field_module.f90->sourcefile~mus_bc_header_module.f90 sourcefile~mus_field_module.f90->sourcefile~mus_mixture_module.f90 sourcefile~mus_field_module.f90->sourcefile~mus_field_prop_module.f90 sourcefile~mus_fluid_module.f90 mus_fluid_module.f90 sourcefile~mus_field_module.f90->sourcefile~mus_fluid_module.f90 sourcefile~mus_species_module.f90 mus_species_module.f90 sourcefile~mus_field_module.f90->sourcefile~mus_species_module.f90 sourcefile~mus_moments_type_module.f90 mus_moments_type_module.f90 sourcefile~mus_scheme_layout_module.f90->sourcefile~mus_moments_type_module.f90 sourcefile~mus_interpolate_header_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_interpolate_header_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_interpolate_header_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_interpolate_header_module.f90->sourcefile~mus_field_prop_module.f90 sourcefile~mus_directions_module.f90 mus_directions_module.f90 sourcefile~mus_interpolate_header_module.f90->sourcefile~mus_directions_module.f90 sourcefile~mus_nernstplanck_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_mixture_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_pdf_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_field_prop_module.f90 sourcefile~mus_ibm_module.f90 mus_IBM_module.f90 sourcefile~mus_geom_module.f90->sourcefile~mus_ibm_module.f90 sourcefile~mus_geomincrhead_module.f90 mus_geomIncrHead_module.f90 sourcefile~mus_geom_module.f90->sourcefile~mus_geomincrhead_module.f90 sourcefile~mus_transport_var_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_mixture_module.f90->sourcefile~mus_enrtl_module.f90 sourcefile~mus_mixture_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_mixture_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_param_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_abortcriteria_module.f90 mus_abortCriteria_module.f90 sourcefile~mus_param_module.f90->sourcefile~mus_abortcriteria_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_fluid_module.f90 sourcefile~mus_poisson_module.f90 mus_poisson_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_poisson_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_species_module.f90

Files dependent on this one

sourcefile~~mus_auxfieldvar_module.f90~~AfferentGraph sourcefile~mus_auxfieldvar_module.f90 mus_auxFieldVar_module.f90 sourcefile~mus_variable_module.f90 mus_variable_module.f90 sourcefile~mus_variable_module.f90->sourcefile~mus_auxfieldvar_module.f90 sourcefile~mus_derquanpoisson_module.f90 mus_derQuanPoisson_module.f90 sourcefile~mus_variable_module.f90->sourcefile~mus_derquanpoisson_module.f90 sourcefile~mus_derquanpoisson_module.f90->sourcefile~mus_auxfieldvar_module.f90 sourcefile~mus_scheme_module.f90 mus_scheme_module.f90 sourcefile~mus_scheme_module.f90->sourcefile~mus_auxfieldvar_module.f90 sourcefile~mus_scheme_module.f90->sourcefile~mus_variable_module.f90 sourcefile~mus_program_module.f90 mus_program_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_scheme_module.f90 sourcefile~mus_tools_module.f90 mus_tools_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_tools_module.f90 sourcefile~mus_dynloadbal_module.f90 mus_dynLoadBal_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_dynloadbal_module.f90 sourcefile~mus_aux_module.f90 mus_aux_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_aux_module.f90 sourcefile~mus_hvs_config_module.f90 mus_hvs_config_module.f90 sourcefile~mus_hvs_config_module.f90->sourcefile~mus_scheme_module.f90 sourcefile~mus_config_module.f90 mus_config_module.f90 sourcefile~mus_hvs_config_module.f90->sourcefile~mus_config_module.f90 sourcefile~mus_tools_module.f90->sourcefile~mus_scheme_module.f90 sourcefile~mus_harvesting.f90 mus_harvesting.f90 sourcefile~mus_harvesting.f90->sourcefile~mus_scheme_module.f90 sourcefile~mus_harvesting.f90->sourcefile~mus_hvs_config_module.f90 sourcefile~mus_hvs_aux_module.f90 mus_hvs_aux_module.f90 sourcefile~mus_harvesting.f90->sourcefile~mus_hvs_aux_module.f90 sourcefile~mus_dynloadbal_module.f90->sourcefile~mus_scheme_module.f90 sourcefile~mus_dynloadbal_module.f90->sourcefile~mus_tools_module.f90 sourcefile~mus_tracking_module.f90 mus_tracking_module.f90 sourcefile~mus_dynloadbal_module.f90->sourcefile~mus_tracking_module.f90 sourcefile~mus_config_module.f90->sourcefile~mus_scheme_module.f90 sourcefile~mus_config_module.f90->sourcefile~mus_tools_module.f90 sourcefile~mus_hvs_aux_module.f90->sourcefile~mus_tools_module.f90 sourcefile~mus_hvs_aux_module.f90->sourcefile~mus_tracking_module.f90 sourcefile~mus_tracking_module.f90->sourcefile~mus_tools_module.f90 sourcefile~musubi.f90 musubi.f90 sourcefile~musubi.f90->sourcefile~mus_program_module.f90 sourcefile~musubi.f90->sourcefile~mus_config_module.f90 sourcefile~musubi.f90->sourcefile~mus_aux_module.f90 sourcefile~mus_interpolate_verify_module.f90 mus_interpolate_verify_module.f90 sourcefile~mus_interpolate_verify_module.f90->sourcefile~mus_config_module.f90 sourcefile~mus_aux_module.f90->sourcefile~mus_tools_module.f90 sourcefile~mus_aux_module.f90->sourcefile~mus_tracking_module.f90 sourcefile~mus_flow_module.f90 mus_flow_module.f90 sourcefile~mus_flow_module.f90->sourcefile~mus_interpolate_verify_module.f90 sourcefile~mus_control_module.f90 mus_control_module.f90 sourcefile~mus_control_module.f90->sourcefile~mus_aux_module.f90

Contents


Source Code

! Copyright (c) 2019-2021 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2019-2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2021 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 UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! **************************************************************************** !
!> author: Kannan Masilamani
!! This module contains routine to retrieve auxiliary field variables for
!! getElement, getPoint, setupIndices and getValOfIndex.
!! Auxilary field variables are:
!!    * density and velocity for fluid
!!    * species desity and velocity for multispecies
!!    * potential for poisson
!!
! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2011 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2011 Konstantin Kleinheinz <k.kleinheinz@grs-sim.de>
! Copyright (c) 2011-2012 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012, 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2012 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2013-2015, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! Copyright (c) 2014-2015, 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2015-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.







! Copyright (c) 2015 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

! Copyright (c) 2015-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.



! Copyright (c) 2015 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

! Copyright (c) 2013 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2013-2014 Nikhil Anand <nikhil.anand@uni-siegen.de>
! Copyright (c) 2014, 2016 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2015, 2018, 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.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.
!--------------------------------------------
!    A O S - Array of structures layout new
!-------------------------------------------
! Access to get_point value output
! Access to get_element value output
module mus_auxFieldVar_module
  use iso_c_binding, only: c_loc, c_ptr, c_f_pointer

  ! include treelm modules
  use env_module,          only: rk
  use tem_aux_module,      only: tem_abort
  use tem_param_module,    only: rho0, rho0Inv, cs2, cs2inv
  use tem_varSys_module,   only: tem_varSys_type, tem_varSys_op_type
  use tem_time_module,     only: tem_time_type
  use treelmesh_module,    only: treelmesh_type
  use tem_geometry_module, only: tem_CoordOfReal, &
    &                            tem_PosofId, tem_BaryOfId
  use tem_topology_module, only: tem_IdOfCoord
  use tem_topology_module, only: tem_levelOf
  use tem_stencil_module,  only: tem_stencilHeader_type
  use tem_comm_module,     only: tem_commpattern_type, tem_communication_type
  use tem_construction_module, only: tem_levelDesc_type
  use tem_logging_module,  only: logUnit
  use tem_math_module,     only: invert_matrix

  use mus_varSys_module,         only: mus_varSys_data_type
  use mus_scheme_type_module,    only: mus_scheme_type
  use mus_scheme_header_module,  only: mus_scheme_header_type
  use mus_derVarPos_module,      only: mus_derVarPos_type
  use mus_physics_module,        only: mus_convertFac_type
  use mus_auxField_module,       only: mus_proc_calcAuxField
  use mus_source_var_module,     only: mus_source_op_type
  use mus_connectivity_module,   only: mus_getSrcElemPosForIntp
  use mus_eNRTL_module,          only: mus_calc_thermFactor

  implicit none
  private

  public :: mus_assign_calcAuxField_ptr
  public :: mus_access_auxFieldVar_forElement
  public :: mus_auxFieldVar_forPoint
  public :: mus_auxFieldVar_fromIndex
  public :: mus_addForceToAuxField_fluid
  public :: mus_addForceToAuxField_fluidIncomp
  public :: mus_addForceToAuxField_MSL
  public :: mus_addForceToAuxField_MSL_WTDF
  public :: mus_addElectricToAuxField_MSL
  public :: mus_addElectricToAuxField_MSL_WTDF
  public :: mus_addSrcToAuxField_poisson
  public :: mus_addSponFldToAuxField_fluid
  public :: mus_addDynSponFldToAuxField_fluid
  public :: mus_addSrcToAuxField_dummy


contains


  ! ************************************************************************* !
  !> This routine assign function pointer to compute auxField var
  subroutine mus_assign_calcAuxField_ptr(schemeHeader, calcAuxField)
    ! --------------------------------------------------------------------- !
    !> scheme defnition
    type(mus_scheme_header_type), intent(in) :: schemeHeader
    !> function pointer to assign
    procedure(mus_proc_calcAuxField), pointer, intent(out) :: calcAuxField
    ! --------------------------------------------------------------------- !
    ! --------------------------------------------------------------------- !
    calcAuxField => mus_calcAuxField_dummy
    select case (trim(schemeHeader%kind))
    case ('fluid')
      select case (trim(schemeHeader%layout))
      case ('d2q9')
        calcAuxField => mus_calcAuxField_fluid_d2q9
      case ('d3q19')
        calcAuxField => mus_calcAuxField_fluid_d3q19
      case ('d3q27')
        calcAuxField => mus_calcAuxField_fluid_d3q27
      case default
        calcAuxField => mus_calcAuxField_fluid
      end select
    case ('fluid_incompressible', 'isotherm_acEq')
      select case (trim(schemeHeader%layout))
      case ('d2q9')
        calcAuxField => mus_calcAuxField_fluidIncomp_d2q9
      case ('d3q19')
        calcAuxField => mus_calcAuxField_fluidIncomp_d3q19
      case ('d3q27')
        calcAuxField => mus_calcAuxField_fluidIncomp_d3q27
      case default
        calcAuxField => mus_calcAuxField_fluidIncomp
      end select

    case ('passive_scalar','poisson', 'poisson_boltzmann_linear', &
      &   'poisson_boltzmann_nonlinear')
      calcAuxField => mus_calcAuxField_zerothMoment
    case ('nernst_planck')
      calcAuxField => mus_calcAuxField_nernst_planck
    case ('multispecies_gas', 'multispecies_liquid')
      calcAuxField => mus_calcAuxField_MS
    case default
      calcAuxField => mus_calcAuxField_dummy
    end select
  end subroutine mus_assign_calcAuxField_ptr
  ! ************************************************************************* !


  ! ************************************************************************* !
  !> Return the solver aux variable for a given set of elements
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
  recursive subroutine mus_access_auxFieldVar_forElement(fun, varsys, elempos, &
    &                                                    time, tree, nElems,   &
    &                                                    nDofs, res            )
    ! -------------------------------------------------------------------- !
    !> Description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    class(tem_varSys_op_type), intent(in) :: fun

    !> The variable system to obtain the variable from.
    type(tem_varSys_type), intent(in) :: varSys

    !> Position of the TreeID of the element to get the variable for in the
    !! global treeID list.
    integer, intent(in) :: elempos(:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

    !> Number of values to obtain for this variable (vectorized access).
    integer, intent(in) :: nElems

    !> Number of degrees of freedom within an element.
    integer, intent(in) :: nDofs

    !> Resulting values for the requested variable.
    !!
    !! Linearized array dimension:
    !! (n requested entries) x (nComponents of this variable)
    !! x (nDegrees of freedom)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------- !
    integer :: statePos, iElem, iComp, iLevel
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )
    scheme => fPtr%solverData%scheme

    ! res is always AOS layout
    res = 0.0_rk
    do iElem = 1, nElems
      ! if state array is defined level wise then use levelPointer(pos)
      ! to access state array
      statePos = fPtr%solverData%geometry%levelPointer( elemPos(iElem) )
      iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )
      do iComp = 1, fun%nComponents
        res( (iElem-1)*fun%nComponents+iComp ) =                         &
          & scheme%auxField( iLevel )%val(                               &
          ! position of this aux variable in the aux array
          & (statePos-1)*varSys%nAuxScalars + fun%auxField_varPos(iComp) )
      end do !iComp
    end do !iElem

  end subroutine mus_access_auxFieldVar_forElement
  ! ************************************************************************* !

  ! ************************************************************************* !
  !> Auxilary field variable for a given set of points using linear
  !! interpolation. Unlike mus_deriveVar_forPoint which does not consider
  !! ghost and halo elements, this routine considers them because
  !! auxField vars on ghost elements are interpolated and halo elements are
  !! exchanged.
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_point]].
  !!
  recursive subroutine mus_auxFieldVar_forPoint(fun, varsys, point, time, &
    &                                           tree, nPnts, res          )
    ! -------------------------------------------------------------------- !
    !> Description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    class(tem_varSys_op_type), intent(in) :: fun

    !> The variable system to obtain the variable from.
    type(tem_varSys_type), intent(in) :: varSys

    !> Three-dimensional coordinates at which the variable should be
    !! evaluated. Only useful for variables provided as space-time functions.
    real(kind=rk), intent(in) :: point(:,:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

    !> Number of values to obtain for this variable (vectorized access).
    integer, intent(in) :: nPnts

    !> Resulting values for the requested variable.
    !!
    !! Dimension: n requested entries x nComponents of this variable
    !! Access: (iElem-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    integer :: statePos, iPnt, iLevel, elemPos
    integer :: coord(4)
    integer, allocatable :: srcElemPos(:)
    real(kind=rk), allocatable :: weights(:), srcRes(:), pntVal(:)
    integer :: iSrc, iComp, nSrcElems
    ! -------------------------------------------------------------------- !
!write(dbgUnit(1),*) 'Derive for point :'//trim(varSys%varname%val(fun%myPos))
    call C_F_POINTER( fun%method_Data, fPtr )
    scheme => fPtr%solverData%scheme
    allocate(srcElemPos(scheme%layout%fStencil%QQ))
    allocate(weights(scheme%layout%fStencil%QQ))
    allocate(srcRes(scheme%layout%fStencil%QQ*fun%nComponents))
    allocate(pntVal(fun%nComponents))

    res = 0.0_rk
    do iPnt = 1, nPnts
      weights = 0.0_rk
      srcRes = 0.0_rk
      srcElemPos = 0
      coord =  tem_CoordOfReal(tree, point(iPnt,:), tree%global%maxLevel )
      ! returns position of existing element in tree which contains point
      ! 0 if no corresponding node is found,
      ! or the negative of the found ID, if it is a virtual node.
      elemPos = abs(tem_PosofId( tem_IdOfCoord(coord), tree%treeID ))

      ! skip this point if no corresponding element is found
      if (elemPos == 0) cycle

      ! level of existing element
      iLevel = tem_levelOf( tree%treeID( elemPos ) )

      ! position of element in levelDesc total list
      statePos = fPtr%solverData%geometry%levelPointer( elemPos )

      ! get position of source element position in state array for interpolation
      ! using neighbor connectivity array.
      ! Also calculate weights for interpolation using distance between the
      ! point and source element barycenter
      call mus_getSrcElemPosForIntp(                              &
        & srcElemPos  = srcElemPos,                               &
        & weights     = weights,                                  &
        & nSrcElems   = nSrcElems,                                &
        & point       = point(iPnt,:),                            &
        & elemPos     = statePos,                                 &
        & neigh       = fPtr%solverData%scheme%pdf(iLevel)%neigh, &
        & baryOfTotal = scheme%levelDesc(iLevel)%baryOfTotal,     &
        & nElems      = scheme%pdf(iLevel)%nSize,                 &
        & nSolve      = scheme%pdf(iLevel)%nElems_solve,          &
        & stencil     = scheme%layout%fStencil,                   &
        & excludeHalo = .false.,                                  &
        & nScalars    = varSys%nScalars                           )

      ! get source element values
      do iSrc = 1, nSrcElems
        ! position of element in levelDesc total list
        statePos = srcElemPos(iSrc)
        do iComp = 1, fun%nComponents
          srcRes( (iSrc-1)*fun%nComponents + iComp )                       &
            & = scheme%auxField( iLevel )%val(                             &
            & (statePos-1)*varSys%nAuxScalars + fun%auxField_varPos(iComp) )
        end do !iComp
      end do !iSrc

      ! Linear interpolation res = sum(weight_i*phi_i)
      pntVal = 0.0_rk
      do iSrc = 1, nSrcElems
        pntVal(:) = pntVal(:) + weights(iSrc) &
          & * srcRes( (iSrc-1)*fun%nComponents+1 : iSrc*fun%nComponents )
      end do
      res( (iPnt-1)*fun%nComponents+1 : iPnt*fun%nComponents ) = pntVal

    end do !iPnt
  end subroutine mus_auxFieldVar_forPoint
  ! ************************************************************************* !

  ! ************************************************************************* !
  !> Routine to get the actual value for a given array of indices.
  !! The indices belong to the grwarray of points storing levelwise in
  !! Pointdata%pntLvl(iLevel).
  !! Hence this routines takes the indeices as input, can refer to the pointData
  !! and evaluate the variable and returns the values
  subroutine mus_auxFieldVar_fromIndex( fun, varSys, time, iLevel, idx, &
    &                                   idxLen,  nVals,  res            )
    ! -------------------------------------------------------------------------- !
    !> Description of the method to obtain the variables,
    class(tem_varSys_op_type), intent(in) :: fun

    !> The variable system to obtain the variable from.
    type(tem_varSys_type), intent(in) :: varSys

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

    !> Level on which values are requested
    integer, intent(in) :: iLevel

    !> Index of points in the growing array and variable val array to
    !! return.
    !! Size: n
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: dependes on number of first index for contiguous array,
    !! but the sum of all idxLen is equal to nVals
    integer, optional, intent(in) :: idxLen(:)

    !> Number of values to obtain for this variable (vectorized access).
    integer, intent(in) :: nVals

    !> Resulting values for the requested variable.
    !!
    !! Dimension: n requested entries x nComponents of this variable
    !! Access: (iElem-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    integer :: iComp, iSrc, iVal, first, last, loc_level
    integer :: iSrcElems, statePos, nSrcElems
    real(kind=rk), allocatable :: srcRes(:), pntVal(:)
    real(kind=rk) :: weight
    ! -------------------------------------------------------------------------- !

!    write(dbgUnit(4),*) 'get the values of indices for variable ',  &
!      &                  trim(varSys%varname%val(fun%myPos))

    call C_F_POINTER( fun%method_Data, fPtr )
    scheme => fPtr%solverData%scheme
    allocate(srcRes(scheme%layout%fStencil%QQ*fun%nComponents))
    allocate(pntVal(fun%nComponents))

    res = 0.0_rk

    ! distinguish if we have an array of index or we have contingous memory
    ! access where index are always first entries!
    if (present(idxLen)) then
      call tem_abort('Error: idxLen is not supported in get_valOfIndex for ' &
        &          //'state variable')
    end if

    ! Get the state value at the specific point
    do iVal = 1, nVals
      if (idx(iVal)>0) then

        ! elemPos in tree
        ! elemPos = fPtr%pointData%pntLvl(iLevel)%elemPos%val(idx(iVal))

        first = fPtr%pointData%pntLvl(iLevel)%srcElem%first%val(idx(iVal))
        last = fPtr%pointData%pntLvl(iLevel)%srcElem%last%val(idx(iVal))

        ! get pdf's of source elements
        srcRes = 0.0_rk
        iSrcElems = 0
        do iSrc = first, last
          iSrcElems = iSrcElems + 1

          statePos = fPtr%pointData%pntLvl(iLevel)%srcElem%elemPos%val(iSrc)
          loc_level = fPtr%pointData%pntLvl(iLevel)%pntLevel%val(idx(iVal))

          do iComp = 1, fun%nComponents
            srcRes( (iSrcElems-1)*fun%nComponents + iComp )                  &
              & = scheme%auxField(loc_level)%val(                            &
              & (statePos-1)*varSys%nAuxScalars + fun%auxField_varPos(iComp) )
          end do !iComp
        end do !iSrc

        ! Linear interpolation res = sum(weight_i*phi_i)
        pntVal = 0.0_rk
        nSrcElems = 0
        do iSrc = first, last
          weight = fPtr%pointData%pntLvl(iLevel)%srcElem%weight%val(iSrc)
          nSrcElems = nSrcElems + 1
          pntVal(:) = pntVal(:) + weight &
            & * srcRes( (nSrcElems-1)*fun%nComponents+1 &
            &         : nSrcElems*fun%nComponents )
        end do

        ! get the state value for each component of this
        res( (iVal-1)*fun%nComponents+1: iVal*fun%nComponents ) = pntVal
      end if !idx>0
    end do !iVal

  end subroutine mus_auxFieldVar_fromIndex
  ! ************************************************************************* !

  ! ************************************************************************* !
  !> This routine compute auxFields density and velocity for compressible model
  !! for fluid and nGhostFromCoarser elements
  subroutine mus_calcAuxField_fluid(auxField, state, neigh, nSize, nSolve, &
    & iLevel, stencil, varSys, derVarPos)
    ! -------------------------------------------------------------------------- !
    !> output auxField array
    real(kind=rk), intent(inout) :: auxField(:)
    !> input state array
    real(kind=rk), intent(in) :: state(:)
    !> connectivity array
    integer, intent(in) :: neigh(:)
    !> number of elements in the state array
    integer, intent(in) :: nSize
    !> number of elements excluding halos
    integer, intent(in) :: nSolve
    !> current level
    integer, intent(in) :: iLevel
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------------- !
    integer :: dens_pos, vel_pos(3), elemOff
    real(kind=rk) :: rho, pdf( stencil%QQ ), inv_rho
    integer :: iElem, iDir, QQ, nScalars
    ! -------------------------------------------------------------------------- !
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    QQ = stencil%QQ
    nScalars = varSys%nScalars

    !NEC$ ivdep
    do iElem = 1, nSolve
      !NEC$ shortloop
      do iDir = 1, QQ
        pdf(iDir) = state( neigh((idir-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      end do

      ! element offset
      elemoff = (iElem-1)*varSys%nAuxScalars

      ! density
      rho = sum( pdf )
      ! store density
      auxField(elemOff+dens_pos) = rho

      ! store velocity
      inv_rho = 1.0_rk/rho
      auxField(elemOff+vel_pos(1)) = sum( pdf * stencil%cxDirRK(1,:) ) * inv_rho
      auxField(elemOff+vel_pos(2)) = sum( pdf * stencil%cxDirRK(2,:) ) * inv_rho
      auxField(elemOff+vel_pos(3)) = sum( pdf * stencil%cxDirRK(3,:) ) * inv_rho

    end do
  end subroutine mus_calcAuxField_fluid
  ! ************************************************************************* !


  ! ************************************************************************* !
  !> This routine compute auxFields density and velocity for compressible
  !! d2q9 model for fluid and nGhostFromCoarser elements

  subroutine mus_calcAuxField_fluid_d2q9(auxField, state, neigh, nSize, nSolve, &
    & iLevel, stencil, varSys, derVarPos)
    ! -------------------------------------------------------------------------- !
    !> output auxField array
    real(kind=rk), intent(inout) :: auxField(:)
    !> input state array
    real(kind=rk), intent(in) :: state(:)
    !> connectivity array
    integer, intent(in) :: neigh(:)
    !> number of elements in the state array
    integer, intent(in) :: nSize
    !> number of elements excluding halos
    integer, intent(in) :: nSolve
    !> current level
    integer, intent(in) :: iLevel
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------------- !
    integer :: dens_pos, vel_pos(3), elemOff
    real(kind=rk) :: rho, pdf(9), inv_rho, u_x, u_y
    integer :: iElem, QQ, nScalars
    ! -------------------------------------------------------------------------- !
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    QQ = 9
    nScalars = varSys%nScalars

    !NEC$ ivdep
    do iElem = 1, nSolve

      pdf(1) = state( neigh((1-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(2) = state( neigh((2-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(3) = state( neigh((3-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(4) = state( neigh((4-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(5) = state( neigh((5-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(6) = state( neigh((6-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(7) = state( neigh((7-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(8) = state( neigh((8-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(9) = state( neigh((9-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)

      ! element offset
      elemoff = (iElem-1)*varSys%nAuxScalars
      ! density
      rho = sum(pdf)
      ! store density
      auxField(elemOff+dens_pos) = rho

      ! store velocity
      inv_rho = 1.0_rk/rho
  u_x = pdf(3) + pdf(7) + pdf(8) - pdf(1) - pdf(5) - pdf(6)
  u_y = pdf(4) + pdf(6) + pdf(8) - pdf(2) - pdf(5) - pdf(7)

  u_x = u_x * inv_rho
  u_y = u_y * inv_rho
      auxField(elemOff+vel_pos(1)) = u_x
      auxField(elemOff+vel_pos(2)) = u_y
      auxField(elemOff+vel_pos(3)) = 0.0_rk

    end do
  end subroutine mus_calcAuxField_fluid_d2q9
  ! ************************************************************************* !


  ! ************************************************************************* !
  !> This routine compute auxFields density and velocity for compressible
  !! d3q19 model for fluid and nGhostFromCoarser elements
  subroutine mus_calcAuxField_fluid_d3q19(auxField, state, neigh, nSize, &
    &  nSolve, iLevel, stencil, varSys, derVarPos)
    ! ----------------------------------------------------------------------- !
    !> output auxField array
    real(kind=rk), intent(inout) :: auxField(:)
    !> input state array
    real(kind=rk), intent(in) :: state(:)
    !> connectivity array
    integer, intent(in) :: neigh(:)
    !> number of elements in the state array
    integer, intent(in) :: nSize
    !> number of elements excluding halos
    integer, intent(in) :: nSolve
    !> current level
    integer, intent(in) :: iLevel
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ----------------------------------------------------------------------- !
    integer :: dens_pos, vel_pos(3), elemOff
    real(kind=rk) :: rho, pdf(19), inv_rho, u_x, u_y, u_z
    integer :: iElem, QQ, nScalars
    ! ----------------------------------------------------------------------- !
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    QQ = 19
    nScalars = varSys%nScalars

    !$omp parallel do                                              &
    !$omp   & private( rho, pdf, inv_rho, u_x, u_y, u_z, elemOff ) &
    !$omp   & firstprivate( dens_pos, vel_pos )
    !NEC$ ivdep
    do iElem = 1, nSolve

      pdf( 1) = state( neigh(( 1-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf( 2) = state( neigh(( 2-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf( 3) = state( neigh(( 3-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf( 4) = state( neigh(( 4-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf( 5) = state( neigh(( 5-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf( 6) = state( neigh(( 6-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf( 7) = state( neigh(( 7-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf( 8) = state( neigh(( 8-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf( 9) = state( neigh(( 9-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(10) = state( neigh((10-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(11) = state( neigh((11-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(12) = state( neigh((12-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(13) = state( neigh((13-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(14) = state( neigh((14-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(15) = state( neigh((15-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(16) = state( neigh((16-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(17) = state( neigh((17-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(18) = state( neigh((18-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(19) = state( neigh((19-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)

      ! element offset
      elemoff = (iElem-1)*varSys%nAuxScalars

      ! density
      rho = sum(pdf)
      ! store density
      auxField(elemOff+dens_pos) = rho

      ! store velocity
      inv_rho = 1.0_rk/rho
    u_x =   pdf(18) + pdf(17) + pdf( 4) &
      &     + pdf(14) + pdf(12)           &
      &     - pdf(16) - pdf( 1) - pdf(15) &
      &     - pdf(13) - pdf(11)

    u_y =   pdf(18) + pdf( 5) + pdf(16) &
      &     + pdf(10) + pdf( 9)           &
      &     - pdf(15) - pdf( 2) - pdf(17) &
      &     - pdf( 8) - pdf( 7)

    u_z =   pdf( 6) + pdf(14) + pdf(10) &
      &     + pdf(13) + pdf( 8)           &
      &     - pdf( 3) - pdf(12) - pdf( 9) &
      &     - pdf(11) - pdf( 7)

    u_x = u_x * inv_rho
    u_y = u_y * inv_rho
    u_z = u_z * inv_rho
      auxField(elemOff+vel_pos(1)) = u_x
      auxField(elemOff+vel_pos(2)) = u_y
      auxField(elemOff+vel_pos(3)) = u_z

    end do
    !$omp end parallel do

  end subroutine mus_calcAuxField_fluid_d3q19
  ! ************************************************************************* !


  ! ************************************************************************* !
  !> This routine compute auxFields density and velocity for compressible
  !! d3q27 model for fluid and nGhostFromCoarser elements
  subroutine mus_calcAuxField_fluid_d3q27(auxField, state, neigh, nSize, &
    &  nSolve, iLevel, stencil, varSys, derVarPos)
    ! ----------------------------------------------------------------------- !
    !> output auxField array
    real(kind=rk), intent(inout) :: auxField(:)
    !> input state array
    real(kind=rk), intent(in) :: state(:)
    !> connectivity array
    integer, intent(in) :: neigh(:)
    !> number of elements in the state array
    integer, intent(in) :: nSize
    !> number of elements excluding halos
    integer, intent(in) :: nSolve
    !> current level
    integer, intent(in) :: iLevel
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ----------------------------------------------------------------------- !
    integer :: dens_pos, vel_pos(3), elemOff
    real(kind=rk) :: rho, pdf(27), inv_rho, u_x, u_y, u_z
    integer :: iElem, QQ, nScalars
    ! ----------------------------------------------------------------------- !
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    QQ = 27
    nScalars = varSys%nScalars

    !$omp parallel do                                              &
    !$omp   & private( rho, pdf, inv_rho, u_x, u_y, u_z, elemOff ) &
    !$omp   & firstprivate( dens_pos, vel_pos )
    !NEC$ ivdep
    do iElem = 1, nSolve

      pdf(1) = state( neigh((1-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(2) = state( neigh((2-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(3) = state( neigh((3-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(4) = state( neigh((4-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(5) = state( neigh((5-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(6) = state( neigh((6-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(7) = state( neigh((7-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(8) = state( neigh((8-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(9) = state( neigh((9-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(10) = state( neigh((10-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(11) = state( neigh((11-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(12) = state( neigh((12-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(13) = state( neigh((13-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(14) = state( neigh((14-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(15) = state( neigh((15-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(16) = state( neigh((16-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(17) = state( neigh((17-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(18) = state( neigh((18-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(19) = state( neigh((19-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(20) = state( neigh((20-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(21) = state( neigh((21-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(22) = state( neigh((22-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(23) = state( neigh((23-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(24) = state( neigh((24-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(25) = state( neigh((25-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(26) = state( neigh((26-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(27) = state( neigh((27-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)

      ! element offset
      elemoff = (iElem-1)*varSys%nAuxScalars

      ! density
      rho = sum(pdf)
      ! store density
      auxField(elemOff+dens_pos) = rho

      ! store velocity
      inv_rho = 1.0_rk/rho
    u_x =   pdf(18) + pdf(17) + pdf( 4) &
      &     + pdf(14) + pdf(12)           &
      &     + pdf(23) + pdf(24)           &
      &     + pdf(25) + pdf(26)           &
      &     - pdf(16) - pdf( 1) - pdf(15) &
      &     - pdf(13) - pdf(11)              &
      &     - pdf(19) - pdf(20)             &
      &     - pdf(21) - pdf(22)

    u_y =   pdf(18) + pdf( 5) + pdf(16) &
      &     + pdf(10) + pdf( 9)           &
      &     + pdf(21) + pdf(22)           &
      &     + pdf(25) + pdf(26)           &
      &     - pdf(15) - pdf( 2) - pdf(17) &
      &     - pdf( 8) - pdf( 7)               &
      &     - pdf(19) - pdf(20)             &
      &     - pdf(23) - pdf(24)

    u_z =   pdf( 6) + pdf(14) + pdf(10) &
      &     + pdf(13) + pdf( 8)           &
      &     + pdf(20) + pdf(22)           &
      &     + pdf(24) + pdf(26)           &
      &     - pdf( 3) - pdf(12) - pdf( 9) &
      &     - pdf(11) - pdf( 7)               &
      &     - pdf(19) - pdf(21)             &
      &     - pdf(23) - pdf(25)

    u_x = u_x * inv_rho
    u_y = u_y * inv_rho
    u_z = u_z * inv_rho
      auxField(elemOff+vel_pos(1)) = u_x
      auxField(elemOff+vel_pos(2)) = u_y
      auxField(elemOff+vel_pos(3)) = u_z

    end do
    !$omp end parallel do

  end subroutine mus_calcAuxField_fluid_d3q27
 ! ************************************************************************* !



 ! ************************************************************************* !
 !> This routine compute auxFields density and velocity for incompressible model
 !! for fluid and nGhostFromCoarser elements
  subroutine mus_calcAuxField_fluidIncomp(auxField, state, neigh, nSize, &
    & nSolve, iLevel, stencil, varSys, derVarPos)
 ! ------------------------------------------------------------------------- !
    !> output auxField array
    real(kind=rk), intent(inout) :: auxField(:)
    !> input state array
    real(kind=rk), intent(in) :: state(:)
    !> connectivity array
    integer, intent(in) :: neigh(:)
    !> number of elements excluding halos
    integer, intent(in) :: nSolve
    !> number of elements in the state array
    integer, intent(in) :: nSize
    !> current level
    integer, intent(in) :: iLevel
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------------- !
    integer :: dens_pos, vel_pos(3), elemOff
    real(kind=rk) :: pdf( stencil%QQ )
    integer :: iElem, iDir, QQ, nScalars
    ! -------------------------------------------------------------------------- !
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    QQ = stencil%QQ
    nScalars = varSys%nScalars

    !NEC$ ivdep
    do iElem = 1, nSolve
      !NEC$ shortloop
      do iDir = 1, QQ
        pdf(iDir) = state( neigh((idir-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      end do

      ! element offset
      elemoff = (iElem-1)*varSys%nAuxScalars

      ! store density
      auxField(elemOff+dens_pos) = sum(pdf)

      ! store velocity
      auxField(elemOff+vel_pos(1)) = sum( pdf * stencil%cxDirRK(1,:) ) * rho0Inv
      auxField(elemOff+vel_pos(2)) = sum( pdf * stencil%cxDirRK(2,:) ) * rho0Inv
      auxField(elemOff+vel_pos(3)) = sum( pdf * stencil%cxDirRK(3,:) ) * rho0Inv

    end do
  end subroutine mus_calcAuxField_fluidIncomp
  ! ************************************************************************* !

  ! ************************************************************************* !
  !> This routine compute auxFields density and velocity for incompressible
  !! d2q9 model for fluid and nGhostFromCoarser elements

  subroutine mus_calcAuxField_fluidIncomp_d2q9(auxField, state, neigh, nSize, nSolve, &
    & iLevel, stencil, varSys, derVarPos)
    ! -------------------------------------------------------------------------- !
    !> output auxField array
    real(kind=rk), intent(inout) :: auxField(:)
    !> input state array
    real(kind=rk), intent(in) :: state(:)
    !> connectivity array
    integer, intent(in) :: neigh(:)
    !> number of elements in the state array
    integer, intent(in) :: nSize
    !> number of elements excluding halos
    integer, intent(in) :: nSolve
    !> current level
    integer, intent(in) :: iLevel
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------------- !
    integer :: dens_pos, vel_pos(3), elemOff
    real(kind=rk) :: rho, pdf(9), u_x, u_y
    integer :: iElem, QQ, nScalars
    ! -------------------------------------------------------------------------- !
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    QQ = 9
    nScalars = varSys%nScalars

    !NEC$ ivdep
    do iElem = 1, nSolve

      pdf(1) = state( neigh((1-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(2) = state( neigh((2-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(3) = state( neigh((3-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(4) = state( neigh((4-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(5) = state( neigh((5-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(6) = state( neigh((6-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(7) = state( neigh((7-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(8) = state( neigh((8-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(9) = state( neigh((9-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)

      ! element offset
      elemoff = (iElem-1)*varSys%nAuxScalars
      ! density
      rho = sum(pdf)
      ! store density
      auxField(elemOff+dens_pos) = rho

      ! store velocity
  u_x = pdf(3) + pdf(7) + pdf(8) - pdf(1) - pdf(5) - pdf(6)
  u_y = pdf(4) + pdf(6) + pdf(8) - pdf(2) - pdf(5) - pdf(7)

  u_x = u_x * rho0inv
  u_y = u_y * rho0inv
      auxField(elemOff+vel_pos(1)) = u_x
      auxField(elemOff+vel_pos(2)) = u_y
      auxField(elemOff+vel_pos(3)) = 0.0_rk

    end do
  end subroutine mus_calcAuxField_fluidIncomp_d2q9
  ! ************************************************************************* !


  ! ************************************************************************* !
  !> This routine compute auxFields density and velocity for incompressible
  !! d3q19 model for fluid and nGhostFromCoarser elements
  subroutine mus_calcAuxField_fluidIncomp_d3q19(auxField, state, neigh, nSize, &
    &  nSolve, iLevel, stencil, varSys, derVarPos)
    ! ----------------------------------------------------------------------- !
    !> output auxField array
    real(kind=rk), intent(inout) :: auxField(:)
    !> input state array
    real(kind=rk), intent(in) :: state(:)
    !> connectivity array
    integer, intent(in) :: neigh(:)
    !> number of elements in the state array
    integer, intent(in) :: nSize
    !> number of elements excluding halos
    integer, intent(in) :: nSolve
    !> current level
    integer, intent(in) :: iLevel
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ----------------------------------------------------------------------- !
    integer :: dens_pos, vel_pos(3), elemOff
    real(kind=rk) :: rho, pdf(19), u_x, u_y, u_z
    integer :: iElem, QQ, nScalars
    ! ----------------------------------------------------------------------- !
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    QQ = 19
    nScalars = varSys%nScalars

    !NEC$ ivdep
    do iElem = 1, nSolve

      pdf(1) = state( neigh((1-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(2) = state( neigh((2-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(3) = state( neigh((3-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(4) = state( neigh((4-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(5) = state( neigh((5-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(6) = state( neigh((6-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(7) = state( neigh((7-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(8) = state( neigh((8-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(9) = state( neigh((9-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(10) = state( neigh((10-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(11) = state( neigh((11-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(12) = state( neigh((12-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(13) = state( neigh((13-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(14) = state( neigh((14-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(15) = state( neigh((15-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(16) = state( neigh((16-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(17) = state( neigh((17-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(18) = state( neigh((18-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(19) = state( neigh((19-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)

      ! element offset
      elemoff = (iElem-1)*varSys%nAuxScalars

      ! density
      rho = sum(pdf)
      ! store density
      auxField(elemOff+dens_pos) = rho

      ! store velocity
    u_x =   pdf(18) + pdf(17) + pdf( 4) &
      &     + pdf(14) + pdf(12)           &
      &     - pdf(16) - pdf( 1) - pdf(15) &
      &     - pdf(13) - pdf(11)

    u_y =   pdf(18) + pdf( 5) + pdf(16) &
      &     + pdf(10) + pdf( 9)           &
      &     - pdf(15) - pdf( 2) - pdf(17) &
      &     - pdf( 8) - pdf( 7)

    u_z =   pdf( 6) + pdf(14) + pdf(10) &
      &     + pdf(13) + pdf( 8)           &
      &     - pdf( 3) - pdf(12) - pdf( 9) &
      &     - pdf(11) - pdf( 7)

    u_x = u_x * rho0inv
    u_y = u_y * rho0inv
    u_z = u_z * rho0inv
      auxField(elemOff+vel_pos(1)) = u_x
      auxField(elemOff+vel_pos(2)) = u_y
      auxField(elemOff+vel_pos(3)) = u_z


    end do
  end subroutine mus_calcAuxField_fluidIncomp_d3q19
  ! ************************************************************************* !


  ! ************************************************************************* !
  !> This routine compute auxFields density and velocity for incompressible
  !! d3q27 model for fluid and nGhostFromCoarser elements
  subroutine mus_calcAuxField_fluidIncomp_d3q27(auxField, state, neigh, nSize, &
    &  nSolve, iLevel, stencil, varSys, derVarPos)
    ! ----------------------------------------------------------------------- !
    !> output auxField array
    real(kind=rk), intent(inout) :: auxField(:)
    !> input state array
    real(kind=rk), intent(in) :: state(:)
    !> connectivity array
    integer, intent(in) :: neigh(:)
    !> number of elements in the state array
    integer, intent(in) :: nSize
    !> number of elements excluding halos
    integer, intent(in) :: nSolve
    !> current level
    integer, intent(in) :: iLevel
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ----------------------------------------------------------------------- !
    integer :: dens_pos, vel_pos(3), elemOff
    real(kind=rk) :: rho, pdf(27), u_x, u_y, u_z
    integer :: iElem, QQ, nScalars
    ! ----------------------------------------------------------------------- !
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    QQ = 27
    nScalars = varSys%nScalars

    !NEC$ ivdep
    do iElem = 1, nSolve

      pdf(1) = state( neigh((1-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(2) = state( neigh((2-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(3) = state( neigh((3-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(4) = state( neigh((4-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(5) = state( neigh((5-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(6) = state( neigh((6-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(7) = state( neigh((7-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(8) = state( neigh((8-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(9) = state( neigh((9-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(10) = state( neigh((10-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(11) = state( neigh((11-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(12) = state( neigh((12-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(13) = state( neigh((13-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(14) = state( neigh((14-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(15) = state( neigh((15-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(16) = state( neigh((16-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(17) = state( neigh((17-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(18) = state( neigh((18-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(19) = state( neigh((19-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(20) = state( neigh((20-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(21) = state( neigh((21-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(22) = state( neigh((22-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(23) = state( neigh((23-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(24) = state( neigh((24-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(25) = state( neigh((25-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(26) = state( neigh((26-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)
      pdf(27) = state( neigh((27-1)* nsize+ ielem)+( 1-1)* qq+ nscalars*0)

      ! element offset
      elemoff = (iElem-1)*varSys%nAuxScalars

      ! density
      rho = sum(pdf)
      ! store density
      auxField(elemOff+dens_pos) = rho

      ! store velocity
    u_x =   pdf(18) + pdf(17) + pdf( 4) &
      &     + pdf(14) + pdf(12)           &
      &     + pdf(23) + pdf(24)           &
      &     + pdf(25) + pdf(26)           &
      &     - pdf(16) - pdf( 1) - pdf(15) &
      &     - pdf(13) - pdf(11)              &
      &     - pdf(19) - pdf(20)             &
      &     - pdf(21) - pdf(22)

    u_y =   pdf(18) + pdf( 5) + pdf(16) &
      &     + pdf(10) + pdf( 9)           &
      &     + pdf(21) + pdf(22)           &
      &     + pdf(25) + pdf(26)           &
      &     - pdf(15) - pdf( 2) - pdf(17) &
      &     - pdf( 8) - pdf( 7)               &
      &     - pdf(19) - pdf(20)             &
      &     - pdf(23) - pdf(24)

    u_z =   pdf( 6) + pdf(14) + pdf(10) &
      &     + pdf(13) + pdf( 8)           &
      &     + pdf(20) + pdf(22)           &
      &     + pdf(24) + pdf(26)           &
      &     - pdf( 3) - pdf(12) - pdf( 9) &
      &     - pdf(11) - pdf( 7)               &
      &     - pdf(19) - pdf(21)             &
      &     - pdf(23) - pdf(25)

    u_x = u_x * rho0inv
    u_y = u_y * rho0inv
    u_z = u_z * rho0inv
      auxField(elemOff+vel_pos(1)) = u_x
      auxField(elemOff+vel_pos(2)) = u_y
      auxField(elemOff+vel_pos(3)) = u_z


    end do
  end subroutine mus_calcAuxField_fluidIncomp_d3q27
  ! ************************************************************************* !


  ! ************************************************************************* !
  !> This routine compute zeroth moment from state and store in auxField.
  !! use this routine only for models which requires only zeroth-order
  !! moment as auxField
  subroutine mus_calcAuxField_zerothMoment(auxField, state, neigh, nSize, &
    & nSolve, iLevel, stencil, varSys, derVarPos)
    ! -------------------------------------------------------------------------- !
    !> output auxField array
    real(kind=rk), intent(inout) :: auxField(:)
    !> input state array
    real(kind=rk), intent(in) :: state(:)
    !> connectivity array
    integer, intent(in) :: neigh(:)
    !> number of elements excluding halos
    integer, intent(in) :: nSolve
    !> number of elements in the state array
    integer, intent(in) :: nSize
    !> current level
    integer, intent(in) :: iLevel
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------------- !
    real(kind=rk) :: pdf( stencil%QQ )
    integer :: iElem, iDir
    ! -------------------------------------------------------------------------- !
    ! safety check
    if (varSys%nAuxScalars /= 1) then
      call tem_abort('CalcAuxField_zeroth moment can be used only when' &
        &          //'nAuxScalars = 1')
    end if

    !NEC$ ivdep
    do iElem = 1, nSolve
      !NEC$ shortloop
      do iDir = 1, stencil%QQ
        pdf(iDir) = state(                                                    &
          &  neigh((idir-1)* nsize+ ielem)+( 1-1)* stencil%qq+ varsys%nscalars*0)
      end do

      ! store scalar zeroth moment
      auxField(iElem) = sum(pdf)
    end do
  end subroutine mus_calcAuxField_zerothMoment
  ! ************************************************************************* !

  ! ************************************************************************* !
  !> This routine compute zeroth moment (mole density) from state for each
  !! species and store in auxField.
  subroutine mus_calcAuxField_nernst_planck(auxField, state, neigh, nSize, &
    & nSolve, iLevel, stencil, varSys, derVarPos)
    ! -------------------------------------------------------------------------- !
    !> output auxField array
    real(kind=rk), intent(inout) :: auxField(:)
    !> input state array
    real(kind=rk), intent(in) :: state(:)
    !> connectivity array
    integer, intent(in) :: neigh(:)
    !> number of elements excluding halos
    integer, intent(in) :: nSolve
    !> number of elements in the state array
    integer, intent(in) :: nSize
    !> current level
    integer, intent(in) :: iLevel
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------------- !
    real(kind=rk) :: pdf( stencil%QQ )
    integer :: iElem, iDir, iFld, nScalars, elemOff
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    ! -------------------------------------------------------------------------- !
    call C_F_POINTER( varSys%method%val(1)%method_Data, fPtr )
    scheme => fPtr%solverData%scheme
    nScalars = varSys%nScalars

    !NEC$ ivdep
    do iElem = 1, nSolve
      ! element offset
      elemOff = (iElem-1)*varSys%nAuxScalars
      ! store scalar mole density
      !NEC$ shortloop
      do iFld = 1, scheme%nFields
        do iDir = 1, stencil%QQ
          pdf(iDir) = state(                                                &
            &  neigh((idir-1)* nsize+ ielem)+( ifld-1)* stencil%qq+ nscalars*0)
        end do
        ! Since nernst_planck has only mole density in auxField array.
        ! So, instead of varSys->auxField_varPos, the array can be accessed
        ! using iFld itself
        auxField(elemOff + iFld) = sum(pdf)
      end do !iField
    end do !iElem
  end subroutine mus_calcAuxField_nernst_planck
  ! ************************************************************************* !


  ! ************************************************************************* !
  !> This routine compute auxField density and momentum for each species
  !! for multicomponent models. The momentum computed here is only momentum
  !! of transformed PDF. The momentum of original PDF is computed by solving
  !! linear equation system in compute kernel and the momentum in auxField is
  !! updated there.
  subroutine mus_calcAuxField_MS(auxField, state, neigh, nSize, &
    & nSolve, iLevel, stencil, varSys, derVarPos)
    ! ------------------------------------------------------------------------ !
    !> output auxField array
    real(kind=rk), intent(inout) :: auxField(:)
    !> input state array
    real(kind=rk), intent(in) :: state(:)
    !> connectivity array
    integer, intent(in) :: neigh(:)
    !> number of elements excluding halos
    integer, intent(in) :: nSolve
    !> number of elements in the state array
    integer, intent(in) :: nSize
    !> current level
    integer, intent(in) :: iLevel
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ------------------------------------------------------------------------ !
    integer :: iElem, iFld, iDir, elemOff, nScalars
    integer :: QQ, nFields, dens_pos, mom_pos(3)
    real(kind=rk) :: pdf( stencil%QQ )
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    ! ------------------------------------------------------------------------ !
    call C_F_POINTER( varSys%method%val(1)%method_Data, fPtr )
    scheme => fPtr%solverData%scheme

    nFields = scheme%nFields
    QQ = stencil%QQ
    nScalars = varSys%nScalars

    !NEC$ ivdep
    do iElem = 1, nSolve
      ! element offset
      elemoff = (iElem-1)*varSys%nAuxScalars

      !NEC$ shortloop
      do iFld = 1, scheme%nFields
        do iDir = 1, QQ
          pdf(iDir) = state( neigh((idir-1)* nsize+ ielem)+( ifld-1)* qq+ nscalars*0)
        end do
        ! position of density and velocity of current field in auxField array
        dens_pos = varSys%method%val(derVarPos(iFld)%density)%auxField_varPos(1)
        mom_pos = varSys%method%val(derVarPos(iFld)%momentum)%auxField_varPos(:)

        ! store mass density of species
        auxField(elemOff+dens_pos) = sum(pdf)

        ! store momentum
        auxField(elemOff+mom_pos(1)) = sum(pdf * stencil%cxDirRK(1, :))
        auxField(elemOff+mom_pos(2)) = sum(pdf * stencil%cxDirRK(2, :))
        auxField(elemOff+mom_pos(3)) = sum(pdf * stencil%cxDirRK(3, :))
      end do !iField
    end do

  end subroutine mus_calcAuxField_MS
  ! ************************************************************************* !

  ! ************************************************************************* !
  !> Dummy routine for calcAuxField
  subroutine mus_calcAuxField_dummy(auxField, state, neigh, nSize, nSolve, &
    & iLevel, stencil, varSys, derVarPos)
    ! -------------------------------------------------------------------------- !
    !> output auxField array
    real(kind=rk), intent(inout) :: auxField(:)
    !> input state array
    real(kind=rk), intent(in) :: state(:)
    !> connectivity array
    integer, intent(in) :: neigh(:)
    !> number of elements excluding halos
    integer, intent(in) :: nSolve
    !> number of elements in the state array
    integer, intent(in) :: nSize
    !> current level
    integer, intent(in) :: iLevel
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------------- !
    call tem_abort('Dummy routine for calcAuxField')
  end subroutine mus_calcAuxField_dummy
  ! ************************************************************************* !

  ! ************************************************************************** !
  !> This routine add body force to velocity in auxField for weakly-compressible
  !! model.
  subroutine mus_addForceToAuxField_fluid(fun, auxField, iLevel, time, varSys, &
    &                                   phyConvFac, derVarPos)
    ! ------------------------------------------------------------------------ !
    !> Description of method to update source
    class(mus_source_op_type), intent(inout) :: fun
    !> output auxField array
    real(kind=rk), intent(inout)         :: auxField(:)
    !> current level
    integer, intent(in)                :: iLevel
    !> current timing information
    type(tem_time_type), intent(in)    :: time
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> Physics conversion factor for current level
    type(mus_convertFac_type), intent(in) :: phyConvFac
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ------------------------------------------------------------------------ !
    integer :: dens_pos, vel_pos(3)
    real(kind=rk) :: forceTerm(3)
    integer :: iElem, nElems, posInTotal, elemOff
    real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3), inv_rho
    ! ------------------------------------------------------------------------ !
    ! position of density and velocity field in auxField
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)
    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems
    ! Get force which is refered in config file either its
    ! spacetime variable or operation variable
    call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
      & varSys  = varSys,                                   &
      & time    = time,                                     &
      & iLevel  = iLevel,                                   &
      & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
      & nVals   = nElems,                                   &
      & res     = forceField                                )

    ! convert physical to lattice
    forceField = forceField / phyConvFac%body_force

!$omp parallel do schedule(static), private( posInTotal, forceTerm, inv_rho, elemOff )
    !NEC$ ivdep
    do iElem = 1, nElems
      posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem)
      ! element offset
      elemoff = (posInTotal-1)*varSys%nAuxScalars
      ! inverse of density
      inv_rho = 1.0_rk/auxField(elemOff+dens_pos)
      ! forceterm to add to velocity: F/(2*rho)
      forceTerm = forceField((iElem-1)*3+1 : iElem*3)*0.5_rk*inv_rho
      ! add force to velocity
      auxField(elemOff+vel_pos(1)) = auxField(elemOff+vel_pos(1)) + forceTerm(1)
      auxField(elemOff+vel_pos(2)) = auxField(elemOff+vel_pos(2)) + forceTerm(2)
      auxField(elemOff+vel_pos(3)) = auxField(elemOff+vel_pos(3)) + forceTerm(3)
    end do

  end subroutine mus_addForceToAuxField_fluid
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This routine add force to velocity in auxField for incompressible model
  subroutine mus_addForceToAuxField_fluidIncomp(fun, auxField, iLevel, time, &
    & varSys, phyConvFac, derVarPos)
    ! ------------------------------------------------------------------------ !
    !> Description of method to update source
    class(mus_source_op_type), intent(inout) :: fun
    !> output auxField array
    real(kind=rk), intent(inout)         :: auxField(:)
    !> current level
    integer, intent(in)                :: iLevel
    !> current timing information
    type(tem_time_type), intent(in)    :: time
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> Physics conversion factor for current level
    type(mus_convertFac_type), intent(in) :: phyConvFac
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ------------------------------------------------------------------------ !
    integer :: vel_pos(3), elemOff
    real(kind=rk) :: forceTerm(3)
    integer :: iElem, nElems, posInTotal
    real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3)
    ! ------------------------------------------------------------------------ !
    ! position of velocity field in auxField
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)
    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems
    ! Get force which is refered in config file either its
    ! spacetime variable or operation variable
    call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
      & varSys  = varSys,                                   &
      & time    = time,                                     &
      & iLevel  = iLevel,                                   &
      & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
      & nVals   = nElems,                                   &
      & res     = forceField                                )

    ! convert physical to lattice
    forceField = forceField / phyConvFac%body_force

!$omp parallel do schedule(static), private( posInTotal, forceTerm, elemOff )
    !NEC$ ivdep
    do iElem = 1, nElems
      posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem)
      ! element offset
      elemoff = (posInTotal-1)*varSys%nAuxScalars
      ! forceterm to add to velocity: F/(2*rho)
      forceTerm = forceField((iElem-1)*3+1 : iElem*3)*0.5_rk*rho0Inv
      ! add force to velocity
      auxField(elemOff+vel_pos(1)) = auxField(elemOff+vel_pos(1)) + forceTerm(1)
      auxField(elemOff+vel_pos(2)) = auxField(elemOff+vel_pos(2)) + forceTerm(2)
      auxField(elemOff+vel_pos(3)) = auxField(elemOff+vel_pos(3)) + forceTerm(3)
    end do

  end subroutine mus_addForceToAuxField_fluidIncomp
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This routine add body force to momentum in auxField for multispecies
  !! liquid model
  !! Refer to Appendix in PhD Thesis of K. Masilamani
  !! "Coupled Simulation Framework to Simulate Electrodialysis Process for
  !! Seawater Desalination"
  subroutine mus_addForceToAuxField_MSL(fun, auxField, iLevel, time, varSys, &
    &                                   phyConvFac, derVarPos)
    ! ------------------------------------------------------------------------ !
    !> Description of method to update source
    class(mus_source_op_type), intent(inout) :: fun
    !> output auxField array
    real(kind=rk), intent(inout)         :: auxField(:)
    !> current level
    integer, intent(in)                :: iLevel
    !> current timing information
    type(tem_time_type), intent(in)    :: time
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> Physics conversion factor for current level
    type(mus_convertFac_type), intent(in) :: phyConvFac
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ------------------------------------------------------------------------ !
    integer :: dens_pos, mom_pos(3), depField
    real(kind=rk) :: forceTerm(3)
    real(kind=rk), dimension(varSys%nStateVars) :: mass_dens, massFrac
    integer :: iElem, nElems, elemOff, nInputStates, iField
    real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3)
    type(mus_varSys_data_type), pointer :: fPtr
    ! ------------------------------------------------------------------------ !
    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    ! number of pdf states this source depends on
    ! last input is spacetime function so it is neglected
    nInputStates = varSys%method%val(fun%srcTerm_varPos)%nInputs - 1

    ! Get force which is refered in config file either its
    ! spacetime variable or operation variable
    call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
      & varSys  = varSys,                                   &
      & time    = time,                                     &
      & iLevel  = iLevel,                                   &
      & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
      & nVals   = nElems,                                   &
      & res     = forceField                                )

    ! convert physical to lattice
    forceField = forceField / phyConvFac%body_force

    call c_f_pointer( varSys%method%val(1)%method_data, fPtr )
    associate( scheme => fPtr%solverData%scheme,            &
      &        posInTotal => fun%elemLvl(iLevel)%posInTotal )

!$omp parallel do schedule(static), private( forceTerm, elemOff )
      !NEC$ ivdep
      do iElem = 1, nElems
        ! element offset
        elemoff = (posInTotal(iElem)-1)*varSys%nAuxScalars

        ! forceterm to add to momentum: F/2
        forceTerm = forceField((iElem-1)*3+1 : iElem*3)*0.5_rk

        !NEC$ shortloop
        do iField = 1, scheme%nFields
          ! position of density and momentum field in auxField
          dens_pos = varSys%method%val(derVarPos(iField)%density) &
            &                     %auxField_varPos(1)
          ! mass density of species
          mass_dens(iField) = auxField( elemOff + dens_pos )
        end do

        !mass fraction
        massFrac(:) = mass_dens(:)/sum(mass_dens)

        ! Update auxField depends on nInputStates
        ! if nInputStates = 1, it is field source else it is global source
        !NEC$ shortloop
        do iField = 1, nInputStates
          depField = varSys%method%val(fun%srcTerm_varPos)%input_varPos(iField)
          ! Position of depField momentum in auxField array
          mom_pos = varSys%method%val(derVarPos(depField)%momentum) &
            &                    %auxField_varPos(1:3)

          ! add force to momentum
          auxField(elemOff+mom_pos(1)) = auxField(elemOff+mom_pos(1))    &
            &                          + massFrac(depField) * forceTerm(1)
          auxField(elemOff+mom_pos(2)) = auxField(elemOff+mom_pos(2))    &
            &                          + massFrac(depField) * forceTerm(2)
          auxField(elemOff+mom_pos(3)) = auxField(elemOff+mom_pos(3))    &
            &                          + massFrac(depField) * forceTerm(3)
        end do !iField
      end do !iElem
    end associate
  end subroutine mus_addForceToAuxField_MSL
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This routine add electric force to momentum in auxField for multispecies
  !! liquid model
  !! Refer to Appendix in PhD Thesis of K. Masilamani
  !! "Coupled Simulation Framework to Simulate Electrodialysis Process for
  !! Seawater Desalination"
  subroutine mus_addElectricToAuxField_MSL(fun, auxField, iLevel, time, &
    &                                      varSys, phyConvFac, derVarPos)
    ! ------------------------------------------------------------------------ !
    !> Description of method to update source
    class(mus_source_op_type), intent(inout) :: fun
    !> output auxField array
    real(kind=rk), intent(inout)         :: auxField(:)
    !> current level
    integer, intent(in)                :: iLevel
    !> current timing information
    type(tem_time_type), intent(in)    :: time
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> Physics conversion factor for current level
    type(mus_convertFac_type), intent(in) :: phyConvFac
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ------------------------------------------------------------------------ !
    integer :: dens_pos, mom_pos(3), depField
    real(kind=rk), dimension(varSys%nStateVars) :: mass_dens, massFrac
    integer :: iElem, nElems, elemOff, nInputStates, iField
    real(kind=rk) :: electricField(fun%elemLvl(iLevel)%nElems*3)
    real(kind=rk) :: EF_elem(3)
    real(kind=rk), dimension(3, varSys%nStateVars ) :: spcForce
    real(kind=rk) :: charge_dens, diffForce_cs2inv, minMolWeight
    real(kind=rk), dimension(varSys%nStateVars) :: chargeTerm
    type(mus_varSys_data_type), pointer :: fPtr
    ! ------------------------------------------------------------------------ !
    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    ! number of pdf states this source depends on
    ! last input is spacetime function so it is neglected
    nInputStates = varSys%method%val(fun%srcTerm_varPos)%nInputs - 1

    call c_f_pointer( varSys%method%val(1)%method_data, fPtr )
    associate( scheme => fPtr%solverData%scheme,                            &
      &        posInTotal => fun%elemLvl(iLevel)%posInTotal,                &
      &        physics => fPtr%solverData%physics,                          &
      &        mixture => fPtr%solverData%scheme%mixture,                   &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

      ! Get electrical force which is refered in config file either its
      ! spacetime variable or operation variable
      call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
        & varSys  = varSys,                                   &
        & time    = time,                                     &
        & iLevel  = iLevel,                                   &
        & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
        & nVals   = nElems,                                   &
        & res     = electricField                             )

      ! convert physical to lattice
      electricField = electricField * physics%coulomb0 &
        &           / physics%fac(iLevel)%force

      ! minimum molecular weight
      minMolWeight = minval(species(:)%molWeight)

      ! constant term to multiply forcing term
      diffForce_cs2inv = minMolWeight / ( mixture%gasConst_R_LB &
        &              * mixture%temp0LB )

!$omp parallel do schedule(static), private( elemOff )
      !NEC$ ivdep
      do iElem = 1, nElems
        ! element offset
        elemoff = (posInTotal(iElem)-1)*varSys%nAuxScalars

        !NEC$ shortloop
        do iField = 1, scheme%nFields
          ! position of density and momentum field in auxField
          dens_pos = varSys%method%val(derVarPos(iField)%density) &
            &                     %auxField_varPos(1)

          ! mass density of species
          mass_dens(iField) = auxField( elemOff + dens_pos )

          ! chargeTerm for each species: \rho_k z_k Faraday / M_k
          chargeTerm(iField) = mass_dens(iField)            &
            &                * species(iField)%molWeightInv &
            &                * species(iField)%chargeNr     &
            &                * mixture%faradayLB
        end do

        !mass fraction
        massFrac(:) = mass_dens(:)/sum(mass_dens)

        ! compute charge density: \sum_k \rho_k z_k Faraday / M_k
        charge_dens = sum(chargeTerm)

        ! electric field for current element
        EF_elem = electricField((iElem-1)*3+1 : iElem*3) * 0.5_rk

        ! compute force on each species
        ! F_k = cs2 min(M) (\rho_k z_k/M_k - y_k \sum_l \rho_l z_l / M_l)
        !      F E / (RT)
        ! forceterm to add to momentum: F_k/2
        !NEC$ shortloop
        do iField = 1, scheme%nFields
          spcForce(:, iField) = EF_elem * cs2 * diffForce_cs2inv      &
            & * ( chargeTerm(iField) - massFrac(iField) * charge_dens )
        end do

        ! Update auxField depends on nInputStates
        ! if nInputStates = 1, it is field source else it is global source
        !NEC$ shortloop
        do iField = 1, nInputStates
          depField = varSys%method%val(fun%srcTerm_varPos)%input_varPos(iField)
          ! Position of depField momentum in auxField array
          mom_pos = varSys%method%val(derVarPos(depField)%momentum) &
            &                    %auxField_varPos(1:3)

          ! add force to momentum
          auxField(elemOff+mom_pos(1)) = auxField(elemOff+mom_pos(1))    &
            &                          + spcForce(1, depField)
          auxField(elemOff+mom_pos(2)) = auxField(elemOff+mom_pos(2))    &
            &                          + spcForce(2, depField)
          auxField(elemOff+mom_pos(3)) = auxField(elemOff+mom_pos(3))    &
            &                          + spcForce(3, depField)
        end do !iField

      end do !iElem
    end associate

  end subroutine mus_addElectricToAuxField_MSL
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This routine add body force to momentum in auxField for multispecies
  !! liquid model with thermodynamic factor
  !! Refer to Appendix in PhD Thesis of K. Masilamani
  !! "Coupled Simulation Framework to Simulate Electrodialysis Process for
  !! Seawater Desalination"
  subroutine mus_addForceToAuxField_MSL_WTDF(fun, auxField, iLevel, time, &
    &                                        varSys, phyConvFac, derVarPos)
    ! ------------------------------------------------------------------------ !
    !> Description of method to update source
    class(mus_source_op_type), intent(inout) :: fun
    !> output auxField array
    real(kind=rk), intent(inout)         :: auxField(:)
    !> current level
    integer, intent(in)                :: iLevel
    !> current timing information
    type(tem_time_type), intent(in)    :: time
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> Physics conversion factor for current level
    type(mus_convertFac_type), intent(in) :: phyConvFac
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ------------------------------------------------------------------------ !
    integer :: dens_pos, mom_pos(3), depField
    real(kind=rk) :: forceTerm(3)
    real(kind=rk), dimension(varSys%nStateVars) :: mass_dens, massFrac
    real(kind=rk), dimension(varSys%nStateVars) :: num_dens, moleFrac
    integer :: iElem, nElems, elemOff, nInputStates, iField, iField_2
    real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3)
    real(kind=rk), dimension(varSys%nStateVars, varSys%nStateVars) :: &
      & thermodynamic_fac, inv_thermodyn_fac
    real(kind=rk), dimension(3, varSys%nStateVars ) :: spcForce, spcForce_WTDF
    type(mus_varSys_data_type), pointer :: fPtr
    ! ------------------------------------------------------------------------ !
    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    ! number of pdf states this source depends on
    ! last input is spacetime function so it is neglected
    nInputStates = varSys%method%val(fun%srcTerm_varPos)%nInputs - 1

    ! Get force which is refered in config file either its
    ! spacetime variable or operation variable
    call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
      & varSys  = varSys,                                   &
      & time    = time,                                     &
      & iLevel  = iLevel,                                   &
      & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
      & nVals   = nElems,                                   &
      & res     = forceField                                )

    ! convert physical to lattice
    forceField = forceField / phyConvFac%body_force

    call c_f_pointer( varSys%method%val(1)%method_data, fPtr )
    associate( scheme => fPtr%solverData%scheme,                            &
      &        posInTotal => fun%elemLvl(iLevel)%posInTotal,                &
      &        mixture => fPtr%solverData%scheme%mixture,                   &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

!$omp parallel do schedule(static), private( forceTerm, elemOff )
      !NEC$ ivdep
      do iElem = 1, nElems
        ! element offset
        elemoff = (posInTotal(iElem)-1)*varSys%nAuxScalars

        !NEC$ shortloop
        do iField = 1, scheme%nFields
          ! position of density and momentum field in auxField
          dens_pos = varSys%method%val(derVarPos(iField)%density) &
            &                     %auxField_varPos(1)
          ! mass density of species
          mass_dens(iField) = auxField( elemOff + dens_pos )

          ! number density of species
          num_dens(iField) = mass_dens(iField) * species(iField)%molWeightInv
        end do

        !mass fraction
        massFrac(:) = mass_dens(:)/sum(mass_dens)

        ! mole fraction
        moleFrac(:) = num_dens(:)/sum(num_dens)

        ! Thermodynamic factor from C++ code
        call mus_calc_thermFactor( nFields       = scheme%nFields,    &
          &                        temp          = mixture%temp0,     &
          &                        press         = mixture%atm_press, &
          &                        mole_frac     = moleFrac,          &
          &                        therm_factors = thermodynamic_fac  )

        ! invert thermodynamic factor
        inv_thermodyn_fac = invert_matrix( thermodynamic_fac )

        ! forceterm to add to momentum: F/2
        forceTerm = forceField((iElem-1)*3+1 : iElem*3)*0.5_rk

        ! compute force on each species
        ! F_k =  y_k F
        do iField = 1, scheme%nFields
          spcForce(:, iField) = massFrac(iField) * forceTerm(:)
        end do

        ! compute external forcing term
        ! d^m_k = \gamma^{-1}_{k,l} F_k
        ! F_k is diffusive forcing term
        spcForce_WTDF = 0.0_rk
        do iField = 1, scheme%nFields
          do iField_2 = 1, scheme%nFields
            spcForce_WTDF(:, iField ) = spcForce_WTDF(:, iField)           &
              &                      + inv_thermodyn_fac(iField, iField_2) &
              &                      * spcForce(:, iField_2)
          end do
        end do

        ! Update auxField depends on nInputStates
        ! if nInputStates = 1, it is field source else it is global source
        !NEC$ shortloop
        do iField = 1, nInputStates
          depField = varSys%method%val(fun%srcTerm_varPos)%input_varPos(iField)
          ! Position of depField momentum in auxField array
          mom_pos = varSys%method%val(derVarPos(depField)%momentum) &
            &                    %auxField_varPos(1:3)

          ! add force to momentum
          auxField(elemOff+mom_pos(1)) = auxField(elemOff+mom_pos(1)) &
            &                          + spcForce_WTDF(1, depField)
          auxField(elemOff+mom_pos(2)) = auxField(elemOff+mom_pos(2)) &
            &                          + spcForce_WTDF(2, depField)
          auxField(elemOff+mom_pos(3)) = auxField(elemOff+mom_pos(3)) &
            &                          + spcForce_WTDF(3, depField)
        end do !iField
      end do !iElem
    end associate
  end subroutine mus_addForceToAuxField_MSL_WTDF
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This routine add electric force to momentum in auxField for multispecies
  !! liquid model with thermodynamic factor
  !! Refer to Appendix in PhD Thesis of K. Masilamani
  !! "Coupled Simulation Framework to Simulate Electrodialysis Process for
  !! Seawater Desalination"
  subroutine mus_addElectricToAuxField_MSL_WTDF(fun, auxField, iLevel, time, &
    &                                           varSys, phyConvFac, derVarPos)
    ! ------------------------------------------------------------------------ !
    !> Description of method to update source
    class(mus_source_op_type), intent(inout) :: fun
    !> output auxField array
    real(kind=rk), intent(inout)         :: auxField(:)
    !> current level
    integer, intent(in)                :: iLevel
    !> current timing information
    type(tem_time_type), intent(in)    :: time
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> Physics conversion factor for current level
    type(mus_convertFac_type), intent(in) :: phyConvFac
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ------------------------------------------------------------------------ !
    integer :: dens_pos, mom_pos(3), depField
    real(kind=rk), dimension(varSys%nStateVars) :: mass_dens, massFrac
    real(kind=rk), dimension(varSys%nStateVars) :: num_dens, moleFrac
    integer :: iElem, nElems, elemOff, nInputStates, iField, iField_2
    real(kind=rk) :: electricField(fun%elemLvl(iLevel)%nElems*3)
    real(kind=rk), dimension(varSys%nStateVars, varSys%nStateVars) :: &
      & thermodynamic_fac, inv_thermodyn_fac
    real(kind=rk), dimension(3, varSys%nStateVars ) :: spcForce, spcForce_WTDF
    real(kind=rk) :: charge_dens, diffForce_cs2inv, minMolWeight
    real(kind=rk), dimension(varSys%nStateVars) :: chargeTerm
    type(mus_varSys_data_type), pointer :: fPtr
    ! ------------------------------------------------------------------------ !
    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    ! number of pdf states this source depends on
    ! last input is spacetime function so it is neglected
    nInputStates = varSys%method%val(fun%srcTerm_varPos)%nInputs - 1

    call c_f_pointer( varSys%method%val(1)%method_data, fPtr )
    associate( scheme => fPtr%solverData%scheme,                            &
      &        posInTotal => fun%elemLvl(iLevel)%posInTotal,                &
      &        physics => fPtr%solverData%physics,                          &
      &        mixture => fPtr%solverData%scheme%mixture,                   &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

      ! Get electrical force which is refered in config file either its
      ! spacetime variable or operation variable
      call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
        & varSys  = varSys,                                   &
        & time    = time,                                     &
        & iLevel  = iLevel,                                   &
        & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
        & nVals   = nElems,                                   &
        & res     = electricField                             )

      ! convert physical to lattice
      electricField = electricField * physics%coulomb0 &
        &           / physics%fac(iLevel)%force

      ! minimum molecular weight
      minMolWeight = minval(species(:)%molWeight)

      ! constant term to multiply forcing term
      diffForce_cs2inv = minMolWeight / ( mixture%gasConst_R_LB &
        &           * mixture%temp0LB )

!$omp parallel do schedule(static), private( elemOff )
      !NEC$ ivdep
      do iElem = 1, nElems
        ! element offset
        elemoff = (posInTotal(iElem)-1)*varSys%nAuxScalars

        !NEC$ shortloop
        do iField = 1, scheme%nFields
          ! position of density and momentum field in auxField
          dens_pos = varSys%method%val(derVarPos(iField)%density) &
            &                     %auxField_varPos(1)

          ! mass density of species
          mass_dens(iField) = auxField( elemOff + dens_pos )

          ! number density of species
          num_dens(iField) = mass_dens(iField) * species(iField)%molWeightInv

          ! chargeTerm for each species: \rho_k z_k Faraday / M_k
          chargeTerm(iField) = num_dens(iField) * species(iField)%chargeNr &
            &                * mixture%faradayLB
        end do

        !mass fraction
        massFrac(:) = mass_dens(:)/sum(mass_dens)

        ! compute charge density: \sum_k \rho_k z_k Faraday / M_k
        charge_dens = sum(chargeTerm)

        ! compute force on each species
        ! F_k = cs2 min(M) (\rho_k z_k/M_k - y_k \sum_l \rho_l z_l / M_l)
        !      F E / (RT)
        ! forceterm to add to momentum: F_k/2
        !NEC$ shortloop
        do iField = 1, scheme%nFields
          spcForce(:, iField) = electricField((iElem-1)*3+1 : iElem*3)  &
            & * ( chargeTerm(iField) - massFrac(iField) * charge_dens ) &
            & * diffForce_cs2inv * 0.5_rk * cs2
        end do

        ! mole fraction
        moleFrac(:) = num_dens(:)/sum(num_dens)

        ! Thermodynamic factor from C++ code
        call mus_calc_thermFactor( nFields       = scheme%nFields,    &
          &                        temp          = mixture%temp0,     &
          &                        press         = mixture%atm_press, &
          &                        mole_frac     = moleFrac,          &
          &                        therm_factors = thermodynamic_fac  )

        ! invert thermodynamic factor
        inv_thermodyn_fac = invert_matrix( thermodynamic_fac )

        ! compute external forcing term
        ! d^m_k = \gamma^{-1}_{k,l} F_k
        ! F_k is diffusive forcing term
        spcForce_WTDF = 0.0_rk
        do iField = 1, scheme%nFields
          do iField_2 = 1, scheme%nFields
            spcForce_WTDF(:, iField ) = spcForce_WTDF(:, iField)           &
              &                      + inv_thermodyn_fac(iField, iField_2) &
              &                      * spcForce(:, iField_2)
          end do
        end do

        ! Update auxField depends on nInputStates
        ! if nInputStates = 1, it is field source else it is global source
        !NEC$ shortloop
        do iField = 1, nInputStates
          depField = varSys%method%val(fun%srcTerm_varPos)%input_varPos(iField)
          ! Position of depField momentum in auxField array
          mom_pos = varSys%method%val(derVarPos(depField)%momentum) &
            &                    %auxField_varPos(1:3)

          ! add force to momentum
          auxField(elemOff+mom_pos(1)) = auxField(elemOff+mom_pos(1))    &
            &                          + spcForce_WTDF(1, depField)
          auxField(elemOff+mom_pos(2)) = auxField(elemOff+mom_pos(2))    &
            &                          + spcForce_WTDF(2, depField)
          auxField(elemOff+mom_pos(3)) = auxField(elemOff+mom_pos(3))    &
            &                          + spcForce_WTDF(3, depField)
        end do !iField

      end do !iElem
    end associate

  end subroutine mus_addElectricToAuxField_MSL_WTDF
  ! ************************************************************************** !


! ************************************************************************** !
  !> This routine add source term with charge density in the Poisson equation
  !! to the potential.
  !! Refer to Appendix in PhD Thesis of K. Masilamani
  !! "Coupled Simulation Framework to Simulate Electrodialysis Process for
  !! Seawater Desalination"
  subroutine mus_addSrcToAuxField_poisson(fun, auxField, iLevel, time, &
    &                                            varSys, phyConvFac, derVarPos)
    ! ------------------------------------------------------------------------ !
    !> Description of method to update source
    class(mus_source_op_type), intent(inout) :: fun
    !> output auxField array
    real(kind=rk), intent(inout)         :: auxField(:)
    !> current level
    integer, intent(in)                :: iLevel
    !> current timing information
    type(tem_time_type), intent(in)    :: time
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> Physics conversion factor for current level
    type(mus_convertFac_type), intent(in) :: phyConvFac
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ------------------------------------------------------------------------ !
    integer :: iElem, nElems
    real(kind=rk) :: rhs(fun%elemLvl(iLevel)%nElems), rhs_Fac
    type(mus_varSys_data_type), pointer :: fPtr
    ! ------------------------------------------------------------------------ !
    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    call c_f_pointer( varSys%method%val(fun%srcTerm_varPos)%method_data, fPtr )
    associate( posInTotal => fun%elemLvl(iLevel)%posInTotal,                 &
      &        poisson => fPtr%solverData%scheme%field(1)%fieldProp%poisson, &
      &        physics => fPtr%solverData%physics                            )

      ! factor to multiply rhs
      rhs_Fac = 0.5_rk * poisson%pot_diff / poisson%permittivity

      ! Get charge density which is refered in config file either its
      ! spacetime variable or operation variable
      call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
        & varSys  = varSys,                                   &
        & time    = time,                                     &
        & iLevel  = iLevel,                                   &
        & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
        & nVals   = nElems,                                   &
        & res     = rhs                                       )

      ! convert physical to lattice and multiply with rhs factor
      rhs = (rhs / physics%fac(iLevel)%chargeDens) * rhs_Fac

!$omp parallel do schedule(static)
      !NEC$ ivdep
      do iElem = 1, nElems
        ! add source term to potential
        auxField(posInTotal(iElem)) = auxField(posInTotal(iElem)) + rhs(iElem)
      end do !iElem
    end associate

  end subroutine mus_addSrcToAuxField_poisson
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This routine add sponge density and velocity field to density and velocity
  !! in auxField for weakly-compressible model.
  !! Reference:
  !! Jacob, J.; Sagaut, P. (2019): Solid wall and open boundary conditions in
  !! hybrid recursive regularized lattice Boltzmann method for compressible
  !! flows. In Physics of Fluids 31 (12), p. 126103.
  subroutine mus_addSponFldToAuxField_fluid(fun, auxField, iLevel, time, &
    &                                       varSys, phyConvFac, derVarPos)
    ! ------------------------------------------------------------------------ !
    !> Description of method to update source
    class(mus_source_op_type), intent(inout) :: fun
    !> output auxField array
    real(kind=rk), intent(inout)         :: auxField(:)
    !> current level
    integer, intent(in)                :: iLevel
    !> current timing information
    type(tem_time_type), intent(in)    :: time
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> Physics conversion factor for current level
    type(mus_convertFac_type), intent(in) :: phyConvFac
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ------------------------------------------------------------------------ !
    integer :: dens_pos, vel_pos(3)
    real(kind=rk) :: sigma
    real(kind=rk) :: inv_rho_phy, inv_vel_phy
    integer :: iElem, nElems, elemOff
    real(kind=rk) :: spongeField(fun%elemLvl(iLevel)%nElems)
    real(kind=rk) :: dens, vel(3), spongeForce(27)
    real(kind=rk) :: dens_ref, vel_ref(3), sponDens, sponVel(3)
    ! ------------------------------------------------------------------------ !
    ! position of density and velocity field in auxField
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)
    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems
    ! Get force which is refered in config file either its
    ! spacetime variable or operation variable
    call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
      & varSys  = varSys,                                   &
      & time    = time,                                     &
      & iLevel  = iLevel,                                   &
      & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
      & nVals   = nElems,                                   &
      & res     = spongeField                               )

    inv_rho_phy = 1.0_rk / phyConvFac%press * cs2inv
    inv_vel_phy = 1.0_rk / phyConvFac%vel

    ! target pressure and velocity in lattice unit
    dens_ref = fun%absLayer%config%target_pressure * inv_rho_phy + rho0
    vel_ref(1:3) = fun%absLayer%config%target_velocity(1:3) * inv_vel_phy

    associate(posInTotal => fun%elemLvl(iLevel)%posInTotal)
      !NEC$ ivdep
      do iElem = 1, nElems
        ! element offset
        elemoff = (posInTotal(iElem)-1)*varSys%nAuxScalars
        ! Local density and velocity
        dens = auxField(elemOff+dens_pos)
        vel(1) = auxField(elemOff+vel_pos(1))
        vel(2) = auxField(elemOff+vel_pos(2))
        vel(3) = auxField(elemOff+vel_pos(3))

        ! SpongeField contains: spongeStrength
        sigma = spongeField(iElem)

        ! Sponge factor for density and velocity field
        sponDens = -sigma*(dens - dens_ref)*0.5_rk
        sponVel(:) = -sigma*(vel - vel_ref)*0.5_rk

        ! add force to velocity
        auxField(elemOff+dens_pos) = dens + sponDens
        auxField(elemOff+vel_pos(1)) = vel(1) + sponVel(1)
        auxField(elemOff+vel_pos(2)) = vel(2) + sponVel(2)
        auxField(elemOff+vel_pos(3)) = vel(3) + sponVel(3)
      end do
    end associate

  end subroutine mus_addSponFldToAuxField_fluid
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This routine add sponge density and velocity field to density and velocity
  !! in auxField. Density and velocity in far field are computed by time
  !! average.
  !!
  !! Reference:
  !! Jacob, J.; Sagaut, P. (2019): Solid wall and open boundary conditions in
  !! hybrid recursive regularized lattice Boltzmann method for compressible
  !! flows. In Physics of Fluids 31 (12), p. 126103.
  subroutine mus_addDynSponFldToAuxField_fluid(fun, auxField, iLevel, time, &
    &                                          varSys, phyConvFac, derVarPos)
    ! ------------------------------------------------------------------------ !
    !> Description of method to update source
    class(mus_source_op_type), intent(inout) :: fun
    !> output auxField array
    real(kind=rk), intent(inout)         :: auxField(:)
    !> current level
    integer, intent(in)                :: iLevel
    !> current timing information
    type(tem_time_type), intent(in)    :: time
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> Physics conversion factor for current level
    type(mus_convertFac_type), intent(in) :: phyConvFac
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! ------------------------------------------------------------------------ !
    integer :: dens_pos, vel_pos(3)
    real(kind=rk) :: sigma
    integer :: iElem, nElems,  elemOff
    real(kind=rk) :: spongeField(fun%elemLvl(iLevel)%nElems)
    real(kind=rk) :: dens, vel(3)
    real(kind=rk) :: densAvgNew, velAvgNew(3)
    real(kind=rk) :: sponDens, sponVel(3)
    real(kind=rk) :: smoothFac
    real(kind=rk) :: inv_rho_phy, inv_vel_phy
    real(kind=rk) :: dens_ref, vel_ref(3)
    ! ------------------------------------------------------------------------ !
    ! position of density and velocity field in auxField
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)
    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems
    ! Get force which is refered in config file either its
    ! spacetime variable or operation variable
    call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
      & varSys  = varSys,                                   &
      & time    = time,                                     &
      & iLevel  = iLevel,                                   &
      & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
      & nVals   = nElems,                                   &
      & res     = spongeField                               )

    inv_rho_phy = 1.0_rk / phyConvFac%press * cs2inv
    inv_vel_phy = 1.0_rk / phyConvFac%vel

    ! Calculate time average value
    ! Set initial density and velocity if no initialized
    associate( dynAvg => fun%elemLvl(iLevel)%dynAvg,         &
      &        posInTotal => fun%elemLvl(iLevel)%posInTotal  )
      ! Do Exponential moving average over nRecord window
      ! P_n = P_{n-1} + smoothFac * (P - P_{n-1})
      smoothFac = fun%absLayer%smoothFac
      ! DO time average for density
      if (fun%absLayer%config%isPressDyn) then
        if (dynAvg%isInitDens) then
          dynAvg%isInitDens = .false.
          ! Initialize dynAvg%dens with initialize density
          do iElem = 1, nElems
            ! element offset
            elemoff = (posInTotal(iElem)-1)*varSys%nAuxScalars
            dynAvg%dens(iElem) = auxField(elemOff+dens_pos)
          end do
        else
          do iElem = 1, nElems
            ! element offset
            elemoff = (posInTotal(iElem)-1)*varSys%nAuxScalars
            ! Local density and velocity
            dens = auxField(elemOff+dens_pos)
            ! New average
            dynAvg%dens(iElem) = dynAvg%dens(iElem)                     &
              &                + smoothFac * (dens - dynAvg%dens(iElem) )
          end do
        end if
      else
        ! target pressure in lattice unit
        dens_ref = fun%absLayer%config%target_pressure * inv_rho_phy + rho0
        dynAvg%dens(:) = dens_ref
      end if

      ! DO time average for velocity
      if (fun%absLayer%config%isVelDyn) then
        if (dynAvg%isInitVel) then
          dynAvg%isInitVel = .false.
          ! Initialize dynAvg%vel with initialize velocity
          do iElem = 1, nElems
            ! element offset
            elemoff = (posInTotal(iElem)-1)*varSys%nAuxScalars
            dynAvg%velX(iElem) = auxField(elemOff+vel_pos(1))
            dynAvg%velY(iElem) = auxField(elemOff+vel_pos(2))
            dynAvg%velZ(iElem) = auxField(elemOff+vel_pos(3))
          end do
        else
          do iElem = 1, nElems
            ! element offset
            elemoff = (posInTotal(iElem)-1)*varSys%nAuxScalars
            ! Local velocity
            vel(1) = auxField(elemOff+vel_pos(1))
            vel(2) = auxField(elemOff+vel_pos(2))
            vel(3) = auxField(elemOff+vel_pos(3))
            ! New average
            dynAvg%velX(iElem) = dynAvg%velX(iElem)                       &
              &                + smoothFac * (vel(1) - dynAvg%velX(iElem) )
            dynAvg%velY(iElem) = dynAvg%velY(iElem)                       &
              &                + smoothFac * (vel(2) - dynAvg%velY(iElem) )
            dynAvg%velZ(iElem) = dynAvg%velZ(iElem)                       &
              &                + smoothFac * (vel(3) - dynAvg%velZ(iElem) )
          end do
        end if
      else
        ! target velocity in lattice unit
        vel_ref(1:3) = fun%absLayer%config%target_velocity(1:3) * inv_vel_phy
        dynAvg%velX(:) = vel_ref(1)
        dynAvg%velY(:) = vel_ref(2)
        dynAvg%velZ(:) = vel_ref(3)
      end if

!$omp parallel do schedule(static), private( elemOff )
      !NEC$ ivdep
      do iElem = 1, nElems
        ! element offset
        elemoff = (posInTotal(iElem)-1)*varSys%nAuxScalars
        ! Local density and velocity
        dens = auxField(elemOff+dens_pos)
        vel(1) = auxField(elemOff+vel_pos(1))
        vel(2) = auxField(elemOff+vel_pos(2))
        vel(3) = auxField(elemOff+vel_pos(3))

        !! New avg
        densAvgNew = dynAvg%dens(iElem)
        velAvgNew(1) = dynAvg%velX(iElem)
        velAvgNew(2) = dynAvg%velY(iElem)
        velAvgNew(3) = dynAVg%velZ(iElem)

        ! SpongeField contains: spongeStrength
        sigma = spongeField(iElem)

        ! Sponge factor for density and velocity field
        sponDens = -sigma*(dens - densAvgNew)*0.5_rk
        sponVel(:) = -sigma*(vel - velAvgNew)*0.5_rk

        ! add force to velocity
        auxField(elemOff+dens_pos) = dens + sponDens
        auxField(elemOff+vel_pos(1)) = vel(1) + sponVel(1)
        auxField(elemOff+vel_pos(2)) = vel(2) + sponVel(2)
        auxField(elemOff+vel_pos(3)) = vel(3) + sponVel(3)

      end do
    end associate

  end subroutine mus_addDynSponFldToAuxField_fluid
  ! ************************************************************************** !

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


end module mus_auxFieldVar_module