mus_derQuanMSLiquid_module.f90 Source File


This file depends on

sourcefile~~mus_derquanmsliquid_module.f90~~EfferentGraph sourcefile~mus_derquanmsliquid_module.f90 mus_derQuanMSLiquid_module.f90 sourcefile~mus_enrtl_module.f90 mus_eNRTL_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_enrtl_module.f90 sourcefile~mus_scheme_layout_module.f90 mus_scheme_layout_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_source_var_module.f90 mus_source_var_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_source_var_module.f90 sourcefile~mus_varsys_module.f90 mus_varSys_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_varsys_module.f90 sourcefile~mus_statevar_module.f90 mus_stateVar_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_statevar_module.f90 sourcefile~mus_scheme_type_module.f90 mus_scheme_type_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_scheme_type_module.f90 sourcefile~mus_scheme_header_module.f90 mus_scheme_header_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_dervarpos_module.f90 mus_derVarPos_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_operation_var_module.f90 mus_operation_var_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_operation_var_module.f90 sourcefile~mus_physics_module.f90 mus_physics_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_pdf_module.f90 mus_pdf_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_pdf_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_source_var_module.f90->sourcefile~mus_scheme_header_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_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_statevar_module.f90->sourcefile~mus_varsys_module.f90 sourcefile~mus_statevar_module.f90->sourcefile~mus_scheme_type_module.f90 sourcefile~mus_connectivity_module.f90 mus_connectivity_module.f90 sourcefile~mus_statevar_module.f90->sourcefile~mus_connectivity_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_scheme_layout_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_scheme_type_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_pdf_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_bc_header_module.f90 mus_bc_header_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_bc_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_interpolate_header_module.f90 mus_interpolate_header_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_interpolate_header_module.f90 sourcefile~mus_param_module.f90 mus_param_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_param_module.f90 sourcefile~mus_nernstplanck_module.f90 mus_nernstPlanck_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_nernstplanck_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_graddata_module.f90 mus_gradData_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_graddata_module.f90 sourcefile~mus_auxfield_module.f90 mus_auxField_module.f90 sourcefile~mus_scheme_type_module.f90->sourcefile~mus_auxfield_module.f90 sourcefile~mus_dervarpos_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_operation_var_module.f90->sourcefile~mus_varsys_module.f90 sourcefile~mus_derivedquantities_module.f90 mus_derivedQuantities_module.f90 sourcefile~mus_operation_var_module.f90->sourcefile~mus_derivedquantities_module.f90

Files dependent on this one

sourcefile~~mus_derquanmsliquid_module.f90~~AfferentGraph sourcefile~mus_derquanmsliquid_module.f90 mus_derQuanMSLiquid_module.f90 sourcefile~mus_bc_species_module.f90 mus_bc_species_module.f90 sourcefile~mus_bc_species_module.f90->sourcefile~mus_derquanmsliquid_module.f90 sourcefile~mus_derquanmsgas_module.f90 mus_derQuanMSGas_module.f90 sourcefile~mus_derquanmsgas_module.f90->sourcefile~mus_derquanmsliquid_module.f90 sourcefile~mus_variable_module.f90 mus_variable_module.f90 sourcefile~mus_variable_module.f90->sourcefile~mus_derquanmsliquid_module.f90 sourcefile~mus_variable_module.f90->sourcefile~mus_derquanmsgas_module.f90 sourcefile~mus_bc_general_module.f90 mus_bc_general_module.f90 sourcefile~mus_bc_general_module.f90->sourcefile~mus_bc_species_module.f90 sourcefile~mus_scheme_module.f90 mus_scheme_module.f90 sourcefile~mus_scheme_module.f90->sourcefile~mus_variable_module.f90 sourcefile~mus_harvesting.f90 mus_harvesting.f90 sourcefile~mus_harvesting.f90->sourcefile~mus_scheme_module.f90 sourcefile~mus_hvs_config_module.f90 mus_hvs_config_module.f90 sourcefile~mus_harvesting.f90->sourcefile~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_dynloadbal_module.f90 mus_dynLoadBal_module.f90 sourcefile~mus_dynloadbal_module.f90->sourcefile~mus_bc_general_module.f90 sourcefile~mus_dynloadbal_module.f90->sourcefile~mus_scheme_module.f90 sourcefile~mus_tools_module.f90 mus_tools_module.f90 sourcefile~mus_dynloadbal_module.f90->sourcefile~mus_tools_module.f90 sourcefile~mus_control_module.f90 mus_control_module.f90 sourcefile~mus_control_module.f90->sourcefile~mus_bc_general_module.f90 sourcefile~mus_debug_module.f90 mus_debug_module.f90 sourcefile~mus_control_module.f90->sourcefile~mus_debug_module.f90 sourcefile~mus_debug_module.f90->sourcefile~mus_bc_general_module.f90 sourcefile~mus_tools_module.f90->sourcefile~mus_scheme_module.f90 sourcefile~mus_config_module.f90->sourcefile~mus_scheme_module.f90 sourcefile~mus_config_module.f90->sourcefile~mus_tools_module.f90 sourcefile~mus_program_module.f90 mus_program_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_bc_general_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_scheme_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_dynloadbal_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_control_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_tools_module.f90 sourcefile~mus_hvs_aux_module.f90 mus_hvs_aux_module.f90 sourcefile~mus_hvs_aux_module.f90->sourcefile~mus_tools_module.f90 sourcefile~mus_tracking_module.f90 mus_tracking_module.f90 sourcefile~mus_tracking_module.f90->sourcefile~mus_tools_module.f90 sourcefile~musubi.f90 musubi.f90 sourcefile~musubi.f90->sourcefile~mus_control_module.f90 sourcefile~musubi.f90->sourcefile~mus_config_module.f90 sourcefile~musubi.f90->sourcefile~mus_program_module.f90 sourcefile~mus_aux_module.f90 mus_aux_module.f90 sourcefile~mus_aux_module.f90->sourcefile~mus_tools_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_construction_module.f90 mus_construction_module.f90 sourcefile~mus_construction_module.f90->sourcefile~mus_debug_module.f90

Contents


Source Code

! Copyright (c) 2013-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2013-2014 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2013, 2016 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2014 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2015-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016-2017 Raphael Haupt <raphael.haupt@uni-siegen.de>
! Copyright (c) 2019-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.
! **************************************************************************** !
!> author: Kannan Masilamani
!! This module provides the MUSUBI specific functions for calculating
!! macroscopic quantities from the state variables for multispecies liquid model.
!!
!! The depending common interface between MUSUBI and ATELES is defined in the
!! [[tem_derived_module]]. The functionality for accessing a variable from the
!! state
!! and evaluating a lua function are also provided in the
!! [[tem_derived_module]].
!!
!! Do not use get_Element or get_Point routines to update the state !
!!
! 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) 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_derQuanMSLiquid_module
  use iso_c_binding, only: c_loc, c_ptr, c_f_pointer

  ! include treelm modules
  use env_module,               only: rk, long_k, globalMaxLevels, labelLen
  use tem_param_module,         only: div1_2, div1_3, div1_54, div1_9, div3_4, &
    &                                 sqrt3, cs2inv, cs2, t2cs2inv, t2cs4inv,  &
    &                                 cs4inv
  use tem_aux_module,           only: tem_abort
  use tem_varSys_module,        only: tem_varSys_type, tem_varSys_op_type,     &
    &                                 tem_varSys_append_derVar,                &
    &                                 tem_varSys_proc_point,                   &
    &                                 tem_varSys_proc_element,                 &
    &                                 tem_varSys_proc_setParams,               &
    &                                 tem_varSys_proc_getParams,               &
    &                                 tem_varSys_proc_setupIndices,            &
    &                                 tem_varSys_proc_getValOfIndex,           &
    &                                 tem_varSys_getPoint_dummy,               &
    &                                 tem_varSys_getElement_dummy,             &
    &                                 tem_varSys_setupIndices_dummy,           &
    &                                 tem_varSys_getValOfIndex_dummy,          &
    &                                 tem_varSys_setParams_dummy,              &
    &                                 tem_varSys_getParams_dummy
  use tem_variable_module,      only: tem_variable_type
  use tem_stencil_module,       only: tem_stencilHeader_type
  use tem_logging_module,       only: logUnit
  use tem_topology_module,      only: tem_levelOf
  use tem_time_module,          only: tem_time_type
  use treelmesh_module,         only: treelmesh_type
  use tem_math_module,          only: invert_matrix
  use tem_spatial_module,       only: tem_spatial_for
  use tem_subTree_type_module,  only: tem_subTree_type, tem_treeIDfrom_subTree
  use tem_debug_module,         only: dbgUnit
  use tem_operation_var_module, only: tem_evalMag_forElement,                &
    &                                 tem_evalMag_forPoint,                  &
    &                                 tem_evalMag_fromIndex,                 &
    &                                 tem_opVar_setupIndices,                &
    &                                 tem_get_new_varSys_data_ptr,           &
    &                                 tem_evalAdd_forElement,                &
    &                                 tem_evalAdd_fromIndex,                 &
    &                                 tem_opVar_setParams, tem_opVar_getParams
  use tem_spacetime_fun_module, only: tem_spacetime_for,       &
    &                                 tem_st_fun_listElem_type
  use tem_tools_module,         only: tem_PositionInSorted
  use tem_dyn_array_module,     only: PositionOfVal
  use tem_grow_array_module,    only: grw_labelarray_type, append

  ! include musubi modules
  use mus_source_var_module,     only: mus_source_op_type
  use mus_pdf_module,            only: pdf_data_type
  use mus_scheme_header_module,  only: mus_scheme_header_type
  use mus_scheme_type_module,    only: mus_scheme_type
  use mus_eNRTL_module,          only: mus_calc_thermFactor, &
    &                                  mus_calc_MS_DiffMatrix
  use mus_scheme_layout_module,  only: mus_scheme_layout_type
  use mus_varSys_module,         only: mus_varSys_data_type,        &
    &                                  mus_varSys_solverData_type,  &
    &                                  mus_get_new_solver_ptr,      &
    &                                  mus_deriveVar_ForPoint
  use mus_stateVar_module,       only: mus_accessVar_setupIndices
  use mus_operation_var_module,  only: mus_opVar_setupIndices
  use mus_derVarPos_module,      only: mus_derVarPos_type
  use mus_physics_module,            only: mus_convertFac_type

  ! include aotus modules
  use aotus_module, only: flu_State

  implicit none

  private

  public :: mus_append_derVar_MSLiquid
  public :: mus_append_derMixVar_MS
  public :: deriveMoleDensityMS
  public :: deriveMoleDensityMS_fromIndex
  public :: derivePressureMS
  public :: deriveMoleFluxMS
  public :: deriveMoleFluxMS_fromIndex
  public :: deriveMoleFracMS
  public :: deriveMassFracMS
  public :: deriveVelocityMS, deriveVelocityMS_fromIndex

  public :: deriveEquilMSLiquid_FromMacro
  public :: deriveEqMSLiquid_FromState
  public :: deriveVelMSLiquid_FromState
  public :: deriveMomMSLiquid_FromState
  public :: deriveMomentaMSLiquid_FromState
  public :: deriveVelocitiesMSLiquid_FromState
  public :: deriveAuxMSLiquid_fromState
  public :: deriveAuxMSLiquid_fromState_WTDF
  public :: deriveEquilMSLiquid_fromAux

  public :: momentumFromMacroLSE, momentumFromMacroLSE_WTDF
  ! routines which uses thermodynamic factor and variable diffusivities
  ! public :: deriveVelocityWTDF_MSLiquid_FromState
  !@todo TO IMPLEMENT
!  public :: deriveEquilWTDF_MSLiquid_FromState

  ! source variables
  public :: applySrc_electricMSLiquid_2ndOrd
  public :: applySrc_electricMSLiquid_2ndOrd_WTDF
  public :: applySrc_electricMSLiquid_1stOrd
  public :: applySrc_electricMSLiquid_1stOrd_WTDF
  public :: applySrc_forceMSLiquid_2ndOrd
  public :: applySrc_forceMSLiquid_2ndOrd_WTDF
  public :: applySrc_forceMSLiquid_1stOrd
  public :: applySrc_forceMSLiquid_1stOrd_WTDF

contains

  ! ************************************************************************ !
  !> subroutine to add derive variables for multispecies-liquid
  !! (schemekind = 'multispecies_liquid') to the varsys.
  !! @todo KM: set proper velocity and equilbrium pointer for bgk_thermodynfac
  subroutine mus_append_derVar_MSLiquid( varSys, solverData, schemeHeader, &
    &                                    stencil, nFields, fldLabel,       &
    &                                    derVarName )
    ! -------------------------------------------------------------------- !
    !> global variable system
    type(tem_varSys_type), intent(inout)  :: varSys

    !> Contains pointer to solver data types
    type(mus_varSys_solverData_type), target, intent(in) :: solverData

    !> identifier of the scheme
    type(mus_scheme_header_type), intent(in)  :: schemeHeader

    !> compute stencil defintion
    type(tem_stencilHeader_type), intent(in)       :: stencil

    !> number of fields
    integer, intent(in)                      :: nFields

    !> array of field label prefix. Size=nFields
    character(len=*), intent(in)              :: fldLabel(:)

    !> array of derive physical variables
    type(grw_labelarray_type), intent(inout) :: derVarName
    ! -------------------------------------------------------------------- !
    ! number of derive variables
    integer :: nDerVars, iVar, nComponents, addedPos, iIn
    integer :: iField
    logical :: wasAdded
    character(len=labelLen), allocatable ::  input_varname(:)
    character(len=labelLen)  ::  varName
    procedure(tem_varSys_proc_point), pointer :: get_point => NULL()
    procedure(tem_varSys_proc_element), pointer :: get_element => NULL()
    procedure(tem_varSys_proc_setParams), pointer :: set_params => null()
    procedure(tem_varSys_proc_getParams), pointer :: get_params => null()
    procedure(tem_varSys_proc_setupIndices), pointer :: &
      &                                      setup_indices => null()
    procedure(tem_varSys_proc_getValOfIndex), pointer :: &
      &                                       get_valOfIndex => null()
    type(c_ptr) :: method_data
    character(len=labelLen), allocatable :: derVarName_loc(:)
    ! -------------------------------------------------------------------- !
    nullify(get_point, get_element, set_params, get_params, setup_indices, &
      &     get_valOfIndex)

    nDerVars = 9
    allocate(derVarName_loc(nDerVars))
    derVarName_loc    = [ 'pressure          ', 'velocity          ', &
      &                   'equilibrium       ', 'equilibrium_vel   ', &
      &                   'mole_density      ', 'mole_fraction     ', &
      &                   'mass_fraction     ', 'mole_flux         ', &
      &                   'vel_mag           '   ]

    ! append local derVarName to growing array.
    ! should be done only once for all fields
    do iVar = 1, nDerVars
      call append(derVarName, derVarName_loc(iVar))
    end do

    ! add variable for each field and add mixture variable later
    do iField = 1, nFields
      do iVar = 1, nDerVars

        ! set default pointers, overwrite if neccessary
        get_element => tem_varSys_getElement_dummy
        get_point => mus_deriveVar_ForPoint
        setup_indices => mus_opVar_setupIndices
        get_valOfIndex => tem_varSys_getValOfIndex_dummy
        method_data  = mus_get_new_solver_ptr(solverData)
        set_params => tem_varSys_setParams_dummy
        get_params => tem_varSys_getParams_dummy

        select case(trim(adjustl(derVarName_loc(iVar))))
        case ('mole_density')
          get_element => deriveMoleDensityMS
          get_valOfIndex => deriveMoleDensityMS_fromIndex
          nComponents = 1
          allocate(input_varname(1))
          input_varname(1) = 'density'

        case ('mole_fraction')
          get_element => deriveMoleFracMS
          nComponents = 1
          allocate(input_varname(1))
          input_varname(1) = 'density'

        case ('mass_fraction')
          get_element => deriveMassFracMS
          nComponents = 1
          allocate(input_varname(1))
          input_varname(1) = 'density'

        case ('velocity')
          get_element => deriveVelocityMS
          get_valOfIndex => deriveVelocityMS_fromIndex
          nComponents = 3
          allocate(input_varname(2))
          input_varname(1) = 'density'
          input_varname(2) = 'momentum'

        case ('mole_flux')
          ! momentum / molecularWeight
          get_element => deriveMoleFluxMS
          get_valOfIndex => deriveMoleFluxMS_fromIndex
          nComponents = 3
          allocate(input_varname(1))
          input_varname(1) = 'momentum'

        case ('pressure')
          get_element => derivePressureMS
          nComponents = 1
          allocate(input_varname(1))
          input_varname(1) = 'density'

        case ('equilibrium_vel')
          if (trim(schemeHeader%relaxation) == 'bgk_withthermodynfac' .or. &
            & trim(schemeHeader%relaxation) == 'mrt_withthermodynfac') then
            get_element => deriveEquilVelWTDF_MSLiquid
          else
            get_element => deriveEquilVelMSLiquid
          end if
          nComponents = 3
          allocate(input_varname(2))
          input_varname(1) = 'density'
          input_varname(2) = 'momentum'

        case ('equilibrium')
          if (trim(schemeHeader%relaxation) == 'bgk_withthermodynfac' .or. &
            & trim(schemeHeader%relaxation) == 'mrt_withthermodynfac') then
            get_element => deriveEquilWTDF_MSLiquid
          else
            get_element => deriveEquilMSLiquid
          end if
          nComponents = stencil%QQ
          allocate(input_varname(2))
          input_varname(1) = 'density'
          input_varname(2) = 'momentum'

        case ('vel_mag')
          get_element => tem_evalMag_forElement
          get_point => tem_evalMag_forPoint
          get_valOfIndex => tem_evalMag_fromIndex
          setup_indices => tem_opVar_setupIndices
          method_data = tem_get_new_varSys_data_ptr(method_data)
          nComponents = 1
          allocate(input_varname(1))
          input_varname(1) = 'velocity'

        case default
          write(logUnit(1),*) 'WARNING: Unknown variable: '//&
            &                 trim(derVarName_loc(iVar))
          cycle !go to next variable
        end select

        ! update variable names with field label
        varname = trim(fldLabel(iField))//trim(adjustl(derVarName_loc(iVar)))
        do iIn = 1, size(input_varname)
          input_varname(iIn) = trim(fldLabel(iField))//trim(input_varname(iIn))
        end do

        ! append variable to varSys
        call tem_varSys_append_derVar( me             = varSys,         &
          &                            varName        = trim(varname),  &
          &                            nComponents    = nComponents,    &
          &                            input_varname  = input_varname,  &
          &                            method_data    = method_data,    &
          &                            get_point      = get_point,      &
          &                            get_element    = get_element,    &
          &                            set_params     = set_params,     &
          &                            get_params     = get_params,     &
          &                            setup_indices  = setup_indices,  &
          &                            get_valOfIndex = get_valOfIndex, &
          &                            pos            = addedPos,       &
          &                            wasAdded       = wasAdded        )

        if (wasAdded) then
          write(logUnit(10),*) ' Appended variable: '//trim(varname)
        else if (addedpos < 1) then
          write(logUnit(1),*) 'Error: variable '//trim(varname)// &
            &                 ' is not added to variable system'
        end if

        deallocate(input_varname)
      end do !iVar
    end do !iField

    ! append mixture variable for every derived species variable
    call mus_append_derMixVar_MS( varSys     = varSys,     &
      &                           solverData = solverData, &
      &                           nFields    = nFields,    &
      &                           fldLabel   = fldLabel,   &
      &                           derVarName = derVarName  )

    ! append liquid mixture variables
    call mus_append_derLiquidMixVar( varSys     = varSys,     &
      &                              solverData = solverData, &
      &                              nFields    = nFields,    &
      &                              fldLabel   = fldLabel,   &
      &                              derVarName = derVarName  )

  end subroutine mus_append_derVar_MSLiquid
  ! ************************************************************************ !

  ! ************************************************************************** !
  !> Append mixture variables for multicomponent models
  subroutine mus_append_derMixVar_MS( varSys, solverData, nFields, fldLabel, &
    &                                 derVarName )
    ! -------------------------------------------------------------------- !
    !> global variable system
    type(tem_varSys_type), intent(inout)                 :: varSys
    !> Contains pointer to solver data types
    type(mus_varSys_solverData_type), target, intent(in) :: solverData
    !> number of fields
    integer, intent(in)                                  :: nFields
    !> array of field label prefix. Size=nFields
    character(len=*), intent(in)                         :: fldLabel(:)
    !> array of derive physical variable
    type(grw_labelarray_type), intent(in)                :: derVarName
    ! -------------------------------------------------------------------- !
    ! number of derive variables
    integer :: iVar, nComponents, addedPos
    integer, allocatable :: inPos(:)
    integer :: iField
    logical :: wasAdded
    character(len=labelLen), allocatable ::  input_varname(:)
    character(len=labelLen)  ::  varName
    procedure(tem_varSys_proc_point), pointer :: get_point => NULL()
    procedure(tem_varSys_proc_element), pointer :: get_element => NULL()
    procedure(tem_varSys_proc_setParams), pointer :: set_params => null()
    procedure(tem_varSys_proc_getParams), pointer :: get_params => null()
    procedure(tem_varSys_proc_setupIndices), pointer :: &
      &                                      setup_indices => null()
    procedure(tem_varSys_proc_getValOfIndex), pointer :: &
      &                                       get_valOfIndex => null()
    type(c_ptr) :: method_data
    ! -------------------------------------------------------------------- !

    ! append mixture variable for every derived species variable
    do iVar = 1, derVarName%nVals
      ! set pointers for mixture. mixture quantity is obtained by summing up
      ! species quantities
      get_element => tem_evalAdd_forElement
      get_point => mus_deriveVar_ForPoint
      setup_indices => tem_opVar_setupIndices
      get_valOfIndex => tem_evalAdd_fromIndex
      method_data = tem_get_new_varSys_data_ptr(method_data)
      set_params => tem_opVar_setParams
      get_params => tem_opVar_getParams

      varname = trim(adjustl(derVarName%val(iVar)))

      ! overwrite get_element and get_valOfIndex for certain mixture variable
      select case (trim(varname))
      case ('velocity')
        get_element => deriveMixVelMS
        get_valOfIndex => deriveMixVelMS_fromIndex
        method_data  = mus_get_new_solver_ptr(solverData)
        allocate(input_varname(nFields*2))
        allocate(inPos(nFields*2))
        do iField = 1, nFields
          input_varname(iField) = trim(fldlabel(iField))//'density'
          ! position of input variable in varSys
          inPos(iField) = PositionofVal(varSys%varname, input_varname(iField))
        end do
        do iField = nFields+1, nFields*2
          input_varname(iField) = trim(fldlabel(iField-nFields))//'momentum'
          ! position of input variable in varSys
          inPos(iField) = PositionofVal(varSys%varname, input_varname(iField))
        end do

        ! nComponents same as species momentum/velocity
        nComponents = 3

      case ('vel_mag')
        get_element => tem_evalMag_forElement
        get_point => tem_evalMag_forPoint
        get_valOfIndex => tem_evalMag_fromIndex
        setup_indices => tem_opVar_setupIndices
        method_data = tem_get_new_varSys_data_ptr(method_data)
        nComponents = 1
        allocate(input_varname(1))
        allocate(inPos(1))
        input_varname(1) = 'velocity'
        inPos(1) = PositionofVal(varSys%varname, input_varname(1))

      case ( 'density', 'pressure', 'mole_density', 'mole_fraction', &
        &    'mass_fraction', 'momentum', 'mole_flux'                )
        allocate(input_varname(nFields))
        allocate(inPos(nFields))
        do iField = 1, nFields
          input_varname(iField) = trim(fldlabel(iField))//trim(varname)
          ! position of input variable in varSys
          inPos(iField) = PositionofVal(varSys%varname, input_varname(iField))
        end do

        ! get nComponents from dependent variable
        if (inPos(1) > 0) then
          nComponents = varSys%method%val(inPos(1))%nComponents
        else
          write(logUnit(1),*) 'Input variable for mixture variable ' &
            &                //trim(varname)//' not found in varSys'
          call tem_abort()
        end if
      case default
        ! Not supported as a mixture variable
        cycle
      end select


      ! input variable not found in varSys. Goto next variable
      if (any(inPos <= 0)) then
        deallocate(input_varname)
        deallocate(inPos)
        cycle
      end if

      ! append variable to varSys
      call tem_varSys_append_derVar(  me             = varSys,         &
        &                             varName        = trim(varname),  &
        &                             nComponents    = nComponents,    &
        &                             input_varname  = input_varname,  &
        &                             method_data    = method_data,    &
        &                             get_point      = get_point,      &
        &                             get_element    = get_element,    &
        &                             set_params     = set_params,     &
        &                             get_params     = get_params,     &
        &                             setup_indices  = setup_indices,  &
        &                             get_valOfIndex = get_valOfIndex, &
        &                             pos            = addedPos,       &
        &                             wasAdded       = wasAdded        )

      if (wasAdded) then
        write(logUnit(10),*) ' Appended variable: '//trim(varname)
      else if (addedpos < 1) then
        write(logUnit(1),*) 'Error: variable '//trim(varname)// &
          &                 ' is not added to variable system'
      end if

      deallocate(input_varname)
      deallocate(inPos)
    end do !iVar

  end subroutine mus_append_derMixVar_MS
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Append mixture variables for multicomponent liquid models
  subroutine mus_append_derLiquidMixVar(varSys, solverData, nFields, &
    &                                   fldLabel, derVarName)
    ! -------------------------------------------------------------------- !
    !> global variable system
    type(tem_varSys_type), intent(inout)                 :: varSys
    !> Contains pointer to solver data types
    type(mus_varSys_solverData_type), target, intent(in) :: solverData
    !> number of fields
    integer, intent(in)                                  :: nFields
    !> array of field label prefix. Size=nFields
    character(len=*), intent(in)                         :: fldLabel(:)
    !> array of derive physical variable
    type(grw_labelarray_type), intent(in)                :: derVarName
    ! -------------------------------------------------------------------- !
    ! number of derive variables
    integer :: iVar, nComponents, addedPos, iIn, nInputs
    integer, allocatable :: inPos(:)
    integer :: iField
    logical :: wasAdded
    ! input varname with field label
    character(len=labelLen), allocatable ::  input_varname(:)
    character(len=labelLen)  ::  varName
    procedure(tem_varSys_proc_point), pointer :: get_point => NULL()
    procedure(tem_varSys_proc_element), pointer :: get_element => NULL()
    procedure(tem_varSys_proc_setParams), pointer :: set_params => null()
    procedure(tem_varSys_proc_getParams), pointer :: get_params => null()
    procedure(tem_varSys_proc_setupIndices), pointer :: &
      &                                      setup_indices => null()
    procedure(tem_varSys_proc_getValOfIndex), pointer :: &
      &                                       get_valOfIndex => null()
    type(c_ptr) :: method_data
    integer :: nDerVars
    character(len=labelLen), allocatable :: derVarName_loc(:)
    ! -------------------------------------------------------------------- !

    nDerVars = 3
    allocate(derVarName_loc(nDerVars))
    derVarName_loc    = [ 'kinematic_pressure', 'charge_density    ', &
      &                   'current_density   '                        ]

    ! append mixture variable dedicated to liquid model
    do iVar = 1, nDerVars
      ! append local derVarName to growing array.
      call append(derVarName, derVarName_loc(iVar))
      ! set pointers for mixture. mixture quantity is obtained by summing up
      ! species quantities
      get_element => tem_varSys_getElement_dummy
      get_point => mus_deriveVar_ForPoint
      get_valOfIndex => tem_varSys_getValOfIndex_dummy
      setup_indices => mus_opVar_setupIndices
      method_data  = mus_get_new_solver_ptr(solverData)
      set_params => tem_opVar_setParams
      get_params => tem_opVar_getParams

      varname = trim(adjustl(derVarName_loc(iVar)))

      select case (trim(varname))
      case ('kinematic_pressure')
        get_element => deriveKinePressMSLiquid
        nComponents = 1
        nInputs = 1
        allocate(input_varname(nInputs))
        input_varname(1) = 'pressure'

      case ('charge_density')
        get_element => deriveChargeDensity
        get_valOfIndex => deriveChargeDensity_fromIndex
        nComponents = 1
        nInputs = nFields
        allocate(input_varname(nInputs))
        do iField = 1, nFields
          input_varname(iField) = trim(fldlabel(iField))//'density'
        end do

      case ('current_density')
        get_element => deriveCurrentDensity
        get_valOfIndex => deriveCurrentDensity_fromIndex
        nComponents = 3
        nInputs = nFields
        allocate(input_varname(nInputs))
        do iField = 1, nFields
          input_varname(iField) = trim(fldlabel(iField))//'momentum'
        end do

      case default
        write(logUnit(1),*) 'WARNING: Unknown variable: '//trim(varname)
        cycle !go to next variable
      end select

      ! position of input variable in varSys
      allocate(inPos(nInputs))
      do iIn = 1, nInputs
        inPos(iIn) = PositionofVal(varSys%varname, input_varname(iIn))
      end do

      ! input variable not found in varSys. Goto next variable
      if (any(inPos <= 0)) then
        deallocate(input_varname)
        deallocate(inPos)
        cycle
      end if

      ! append variable to varSys
      call tem_varSys_append_derVar(  me             = varSys,         &
        &                             varName        = trim(varname),  &
        &                             nComponents    = nComponents,    &
        &                             input_varname  = input_varname,  &
        &                             method_data    = method_data,    &
        &                             get_point      = get_point,      &
        &                             get_element    = get_element,    &
        &                             set_params     = set_params,     &
        &                             get_params     = get_params,     &
        &                             setup_indices  = setup_indices,  &
        &                             get_valOfIndex = get_valOfIndex, &
        &                             pos            = addedPos,       &
        &                             wasAdded       = wasAdded        )

      if (wasAdded) then
        write(logUnit(10),*) ' Appended variable: '//trim(varname)
      else if (addedpos < 1) then
        write(logUnit(1),*) 'Error: variable '//trim(varname)// &
          &                 ' is not added to variable system'
      end if

      deallocate(input_varname)
      deallocate(inPos)

    end do !iVar

  end subroutine mus_append_derLiquidMixVar
  ! ************************************************************************** !

  ! ************************************************************************ !
  !      Subroutines with common interface for the function pointers         !
  ! ************************************************************************ !

  ! ************************************************************************ !
  !> Calculate the number density of a given element for single species
  !! from the density stored in auxField array.
  !! Mixture number density is computed by summing species number density
  !! using tem_evalAdd_forElement
  !! Number density = density/molecular weight
  !! mixture number density = sum(number_density)
  recursive subroutine deriveMoleDensityMS(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, iLevel, dens_pos, elemOff
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: iField, depField
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )

    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                      &
      &        levelPointer => fPtr%solverData%geometry%levelPointer, &
      &        auxField => fPtr%solverData%scheme%auxField,           &
      &        field => fPtr%solverData%scheme%field                  )

      ! find dependent field of current variable by comparing the position
      ! of this variable against derVarPos of this variable for all fields in
      ! scheme
      depField = 0
      do iField = 1, scheme%nFields
        if ( fun%myPos == scheme%derVarpos(iField)%moleDensity ) &
          & depField = iField
      end do

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state and auxField array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! element offset for auxField
        elemoff = (statePos-1)*varSys%nAuxScalars

        ! position of density field in auxField array
        dens_pos = varSys%method%val( fun%input_varPos(1) ) &
          &                     %auxField_varPos(1)

        ! number density of species: mass_dens / MolWeight
        res(iElem) = auxField(iLevel)%val( elemOff + dens_pos ) &
          &        * field( depField )%fieldProp%species%molWeightInv

      end do ! iElem
    end associate

  end subroutine deriveMoleDensityMS
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Calculate charge density of a given element for mixture.
  !! Charge density is mixture quantity to it returns same value for all
  !! species
  !! charge_density = Faraday * \sum_i z_i*density_i/molecularWeight_i
  recursive subroutine deriveChargeDensity(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, iLevel, iField, dens_pos, elemOff
    type(mus_varSys_data_type), pointer :: fPtr
    real(kind=rk) :: mass_dens, charge_dens
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )

    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                            &
      &        levelPointer => fPtr%solverData%geometry%levelPointer,       &
      &        auxField => fPtr%solverData%scheme%auxField,                 &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state and auxField array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! element offset for auxField
        elemoff = (statePos-1)*varSys%nAuxScalars

        charge_dens = 0.0_rk
        do iField = 1, scheme%nFields
          ! position of density field in auxField array
          dens_pos = varSys%method%val( fun%input_varPos(iField) ) &
            &                     %auxField_varPos(1)

          ! mass density of species
          mass_dens = auxField(iLevel)%val( elemOff + dens_pos )

          ! charge density
          charge_dens = charge_dens + mass_dens * species(iField)%molWeightInv &
            &         * species(iField)%chargeNr
        end do !iField

        res(iElem) = charge_dens * scheme%mixture%faradayLB

      end do ! iElem
    end associate

  end subroutine deriveChargeDensity
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Current density is computed from species momentum stored in auxField array.
  !! Current density, J = charge_density*velocity = \rho_e * u
  !! = \sum_k z_k F p_k / M_k
  !! where p_k is the species momentum
  recursive subroutine deriveCurrentDensity(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, iField, mom_pos(3)
    type(mus_varSys_data_type), pointer :: fPtr
    ! current density
    real(kind=rk) :: curr_dens(3)
    ! species momentum
    real(kind=rk) :: momentum(3)
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )

    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                            &
      &        levelPointer => fPtr%solverData%geometry%levelPointer,       &
      &        auxField => fPtr%solverData%scheme%auxField,                 &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        curr_dens = 0.0_rk
        do iField = 1, scheme%nFields
          ! position of current field momentum in auxField array
          mom_pos = varSys%method%val( fun%input_varPos(iField) ) &
            &                     %auxField_varPos(1)

          ! species momentum
          do iComp = 1, 3
            momentum(iComp) = auxField(iLevel)%val(                            &
              &               (statePos-1)*varSys%nAuxScalars + mom_pos(iComp) )
          end do

          curr_dens = curr_dens + momentum(:) * species(iField)%molWeightInv &
            &       * species(iField)%chargeNr
        end do !iField

        ! copy the results to the res
        res( (iElem-1)*3 + 1 : iElem*3) = curr_dens * scheme%mixture%faradayLB

      end do ! iElem
    end associate

  end subroutine deriveCurrentDensity
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Calculate mixture kinematic pressure.
  !! This routine requires initial total mole density which is defined in
  !! mixture table.
  !! Formula to compute kinematic pressure
  !! \( p = c^2_s (\sum_k \rho_k \phi_k - min_l (m_l) n_0)/\rho_0 \)
  !! here, \( \rho_k \) - species density, \\
  !! \( \phi_k \) - species molecular weight ratio, \\
  !! \( n_0 \) - reference mixture number density,\\
  !! \( \rho_0 \) - reference density.
  !! In tracking,
  !!```lua
  !! variable = {{"kinematic_pressure"}}
  !!```
  recursive subroutine deriveKinePressMSLiquid(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 :: press_pos
    type(mus_varSys_data_type), pointer :: fPtr
    real(kind=rk) :: min_molWeight, const_press
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )

    associate( mixture => fPtr%solverData%scheme%mixture,                   &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

      ! position of mixture pressue in glob system
      press_pos = fun%input_varPos(1)

      ! derive dependent variable, mixture pressue
      call varSys%method%val(press_pos)%get_element( varSys  = varSys,  &
        &                                            elemPos = elemPos, &
        &                                            time    = time,    &
        &                                            tree    = tree,    &
        &                                            nElems  = nElems,  &
        &                                            nDofs   = nDofs,   &
        &                                            res     = res      )

      min_molWeight = minval(species(:)%molWeight)

      const_press = min_molWeight * mixture%moleDens0 * cs2 / mixture%rho0

      ! p = cs2 (sum_k( rho_k*phi_k ) - min_l (m_l) n_0 ) / rho_0
      res = ( res / mixture%rho0 - const_press )
    end associate

  end subroutine deriveKinePressMSLiquid
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Calculate species  pressure both for gas and liquid model
  !! In case of gas mixture, it is partial pressure where as in
  !! liquid mixture this is not valid.
  !! However, it is used to compute mixture pressure and then the
  !! kinematic_pressure from the mixture pressure.
  !! Formula to calculate pressure:
  !! \( p_k = c^2_s ( \rho_k \phi_k ) \)
  !! here, \( \rho_k \) - species density, \\
  !! \( \phi_k \) - species molecular weight ratio, \\
  recursive subroutine derivePressureMS(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, iLevel, dens_Pos, iField, depField, elemOff
    type(mus_varSys_data_type), pointer :: fPtr
    real(kind=rk) :: mass_dens
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )

    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                      &
      &        levelPointer => fPtr%solverData%geometry%levelPointer, &
      &        auxField => fPtr%solverData%scheme%auxField,           &
      &        field => fPtr%solverData%scheme%field                  )

      ! find dependent field of current variable by comparing the position
      ! of this variable against derVarPos of this variable for all fields in
      ! scheme
      depField = 0
      do iField = 1, scheme%nFields
        if ( fun%myPos == scheme%derVarpos(iField)%pressure ) &
          & depField = iField
      end do

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state and auxField array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! element offset for auxField
        elemoff = (statePos-1)*varSys%nAuxScalars

        ! position of density field in auxField array
        dens_pos = varSys%method%val( fun%input_varPos(1) ) &
          &                     %auxField_varPos(1)

        ! mass density of species
        mass_dens = auxField(iLevel)%val( elemOff + dens_pos )

        ! pressue p_k = cs2 * rho_k * phi_k
        ! multiply factor cs2 after dep since for mixture pressue
        ! p = cs2 * sum_k(rho_k*phi_k)
        res(iElem) = cs2 * mass_dens                                &
          &        * field( depField )%fieldProp%species%molWeigRatio

      end do ! iElem
    end associate

  end subroutine derivePressureMS
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Calculate the velocity of a given element for single species.
  !! from the momentum and density stored in auxField array for liquid mixture.
  !! auxField was updated with momentum of untransformed PDF which was computed
  !! by solving LSE in compute kernel.
  recursive subroutine deriveVelocityMS(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, elemOff
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: dens_pos, mom_pos(3)
    ! mass density of species
    real(kind=rk) :: mass_dens
    ! species equation
    real(kind=rk) :: momentum(3)
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                      &
      &        levelPointer => fPtr%solverData%geometry%levelPointer, &
      &        auxField => fPtr%solverData%scheme%auxField            )

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! element offset for auxField
        elemoff = (statePos-1)*varSys%nAuxScalars
        ! position of current field density in auxField array
        dens_pos = varSys%method%val( fun%input_varPos(1) ) &
          &                     %auxField_varPos(1)

        ! position of current field momentum in auxField array
        mom_pos = varSys%method%val( fun%input_varPos(2) ) &
          &                     %auxField_varPos(1:3)

        ! mass density of species
        mass_dens = auxField(iLevel)%val( elemOff + dens_pos )

        ! species momentum
        do iComp = 1, 3
          momentum(iComp) = auxField(iLevel)%val( elemOff + mom_pos(iComp) )
        end do

        ! compute and store velocity
        res( (iElem-1)*3 + 1 : iElem*3) = momentum / mass_dens

      end do ! iElem
    end associate

  end subroutine deriveVelocityMS
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Equilibrium velocity from density and momentum in auxField.
  !!
  recursive subroutine deriveEquilVelMSLiquid(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, iField, depField, elemOff
    integer :: dens_pos, mom_pos(3)
    type(mus_varSys_data_type), pointer :: fPtr
    !mass density of nSpecies
    real(kind=rk) :: mass_dens( varSys%nStateVars )
    !number density of nSpecies
    real(kind=rk) :: num_dens( varSys%nStateVars )
    !mole fraction
    real(kind=rk) :: moleFraction( varSys%nStateVars )
    real(kind=rk) :: totNum_densInv
    !momentum from auxField
    real(kind=rk) :: momentum( 3 )
    real(kind=rk) :: vel( 3, varSys%nStateVars )
    real(kind=rk) :: eqVel(3)
    !mixture info
    real(kind=rk) :: paramBInv
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )

    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                            &
      &        levelPointer => fPtr%solverData%geometry%levelPointer,       &
      &        auxField => fPtr%solverData%scheme%auxField,                 &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

      ! find dependent field of current variable by comparing the position
      ! of this variable against derVarPos of this variable for all fields in
      ! scheme
      depField = 0
      do iField = 1, scheme%nFields
        if ( fun%myPos == scheme%derVarpos(iField)%equilibriumVel ) &
          & depField = iField
      end do

      paramBInv = 1.0_rk / scheme%mixture%paramB

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! element offset for auxField
        elemoff = (statePos-1)*varSys%nAuxScalars
        do iField = 1, scheme%nFields
          ! position of current field density in auxField array
          dens_pos = varSys%method%val( fun%input_varPos(iField) ) &
            &                     %auxField_varPos(1)

          ! position of current field momentum in auxField array
          mom_pos = varSys%method%val( fun%input_varPos(iField) ) &
            &                     %auxField_varPos(1:3)

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

          ! species momentum
          do iComp = 1, 3
            momentum(iComp) = auxField(iLevel)%val( elemOff + mom_pos(iComp) )
          end do

          ! species velocity
          vel(:, iField) = momentum / mass_dens(iField)

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

        end do !iField

        ! total number density Inv
        totNum_densInv = 1.0_rk/sum(num_dens(:))

        ! molefraction
        moleFraction(:) = num_dens(:)*totNum_densInv

        eqVel = equilVelFromMacro( iField       = depField,             &
          &                        moleFraction = moleFraction,         &
          &                        velocity     = vel,                  &
          &                        nFields      = scheme%nFields,       &
          &                        paramBInv    = paramBInv,            &
          &                        phi          = species(depField)     &
          &                                       %molWeigRatio,        &
          &                        resi_coeff   = species(depField)     &
          &                                       %resi_coeff(:)        )

        ! copy the results to the res
        res( (iElem-1)*3 + 1 : iElem*3) = eqVel

      end do ! iElem
    end associate

  end subroutine deriveEquilVelMSLiquid
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Equilibrium velocity from density and momentum in auxField
  !! with thermodynamic factor
  !!
  recursive subroutine deriveEquilVelWTDF_MSLiquid(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, iField, depField, elemOff
    integer :: dens_pos, mom_pos(3)
    type(mus_varSys_data_type), pointer :: fPtr
    !mass density of nSpecies
    real(kind=rk) :: mass_dens( varSys%nStateVars )
    !number density of nSpecies
    real(kind=rk) :: num_dens( varSys%nStateVars )
    !mole fraction
    real(kind=rk) :: moleFraction( varSys%nStateVars )
    real(kind=rk) :: totNum_densInv
    !momentum from auxField
    real(kind=rk) :: momentum( 3 )
    real(kind=rk) :: vel( 3, varSys%nStateVars )
    real(kind=rk) :: eqVel(3)
    real(kind=rk), dimension(varSys%nStateVars, varSys%nStateVars) :: &
      & resi_coeff, thermodynamic_fac, &
      & inv_thermodyn_fac, diff_coeff
    !mixture info
    real(kind=rk) :: paramBInv
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )

    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                            &
      &        levelPointer => fPtr%solverData%geometry%levelPointer,       &
      &        auxField => fPtr%solverData%scheme%auxField,                 &
      &        mixture => fPtr%solverData%scheme%mixture,                   &
      &        physics => fPtr%solverData%physics,                          &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

      ! find dependent field of current variable by comparing the position
      ! of this variable against derVarPos of this variable for all fields in
      ! scheme
      depField = 0
      do iField = 1, scheme%nFields
        if ( fun%myPos == scheme%derVarpos(iField)%equilibriumVel ) &
          & depField = iField
      end do

      paramBInv = 1.0_rk / mixture%paramB

      do iField = 1, scheme%nFields
        ! diffusivity coefficients
        diff_coeff(iField, :) = species(iField)%diff_coeff(:)
      end do

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! element offset for auxField
        elemoff = (statePos-1)*varSys%nAuxScalars
        do iField = 1, scheme%nFields
          ! position of current field density in auxField array
          dens_pos = varSys%method%val( fun%input_varPos(iField) ) &
            &                     %auxField_varPos(1)

          ! position of current field momentum in auxField array
          mom_pos = varSys%method%val( fun%input_varPos(iField) ) &
            &                     %auxField_varPos(1:3)

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

          ! species momentum
          do iComp = 1, 3
            momentum(iComp) = auxField(iLevel)%val( elemOff + mom_pos(iComp) )
          end do

          ! species velocity
          vel(:, iField) = momentum / mass_dens(iField)

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

        end do !iField

        ! total number density Inv
        totNum_densInv = 1.0_rk/sum(num_dens(:))

        ! molefraction
        moleFraction(:) = num_dens(:)*totNum_densInv

        ! MS-Diff coeff matrix from C++ code
        call mus_calc_MS_DiffMatrix( nFields   = scheme%nFields,             &
          &                          temp      = mixture%temp0,              &
          &                          press     = mixture%atm_press,          &
          &                          mole_dens = num_dens*physics%moleDens0, &
          &                          D_ij_out  = diff_coeff                  )

        ! Convert to lattice unit
        resi_coeff = physics%fac(iLevel)%diffusivity/diff_coeff

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

        ! invert thermodynamic factor
        inv_thermodyn_fac = invert_matrix( thermodynamic_fac )

        ! compute equilibrium velocity with thermodynamic factor
        eqVel = equilVelFromMacroWTDF( iField            = depField,          &
          &                            mass_dens         = mass_dens,         &
          &                            moleFraction      = moleFraction,      &
          &                            velocity          = vel,               &
          &                            nFields           = scheme%nFields,    &
          &                            inv_thermodyn_fac = inv_thermodyn_fac, &
          &                            paramBInv         = paramBInv,         &
          &                            phi               = species(:)         &
          &                                                %molWeigRatio,     &
          &                            resi_coeff        = resi_coeff         )

        ! copy the results to the res
        res( (iElem-1)*3 + 1 : iElem*3) = eqVel

      end do ! iElem
    end associate

  end subroutine deriveEquilVelWTDF_MSLiquid
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Equilibrium from density and momentum stored in auxField
  recursive subroutine deriveEquilMSLiquid(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, iField, depField, elemOff
    integer :: dens_pos, mom_pos(3)
    type(mus_varSys_data_type), pointer :: fPtr
    !mass density of nSpecies
    real(kind=rk) :: mass_dens( varSys%nStateVars )
    !number density of nSpecies
    real(kind=rk) :: num_dens( varSys%nStateVars )
    !mole fraction
    real(kind=rk) :: moleFraction( varSys%nStateVars )
    real(kind=rk) :: totNum_densInv
    !momentum from auxField
    real(kind=rk) :: momentum( 3 )
    real(kind=rk) :: vel( 3, varSys%nStateVars )
    !mixture info
    real(kind=rk) :: paramBInv
    real(kind=rk) :: fEq(fun%nComponents)
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )

    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                            &
      &        levelPointer => fPtr%solverData%geometry%levelPointer,       &
      &        auxField => fPtr%solverData%scheme%auxField,                 &
      &        mixture => fPtr%solverData%scheme%mixture,                   &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

      ! find dependent field of current variable by comparing the position
      ! of this variable against derVarPos of this variable for all fields in
      ! scheme
      depField = 0
      do iField = 1, scheme%nFields
        if ( fun%myPos == scheme%derVarpos(iField)%equilibrium ) &
          & depField = iField
      end do

      paramBInv = 1.0_rk / scheme%mixture%paramB

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! element offset for auxField
        elemoff = (statePos-1)*varSys%nAuxScalars
        do iField = 1, scheme%nFields
          ! position of current field density in auxField array
          dens_pos = varSys%method%val( fun%input_varPos(iField) ) &
            &                     %auxField_varPos(1)

          ! position of current field momentum in auxField array
          mom_pos = varSys%method%val( fun%input_varPos(iField) ) &
            &                     %auxField_varPos(1:3)

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

          ! species momentum
          do iComp = 1, 3
            momentum(iComp) = auxField(iLevel)%val( elemOff + mom_pos(iComp) )
          end do

          ! species velocity
          vel(:, iField) = momentum / mass_dens(iField)

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

        end do !iField

        ! total number density Inv
        totNum_densInv = 1.0_rk/sum(num_dens(:))

        ! molefraction
        moleFraction(:) = num_dens(:)*totNum_densInv

        ! compute equilibrium from macroscopic quantities
        fEq = equilFromMacro( iField       = depField,          &
          &                   mass_dens    = mass_dens,         &
          &                   moleFraction = moleFraction,      &
          &                   velocity     = vel,               &
          &                   layout       = scheme%layout,     &
          &                   nFields      = scheme%nFields,    &
          &                   paramBInv    = paramBInv,         &
          &                   phi          = species(depField)  &
          &                                  %molWeigRatio,     &
          &                   resi_coeff   = species(depField)  &
          &                                  %resi_coeff(:),    &
          &                   theta_eq     = mixture%theta_eq   )

        ! copy the results to the res
        res( (iElem-1)*fun%nComponents + 1 : iElem*fun%nComponents) = fEq

      end do ! iElem
    end associate

   end subroutine deriveEquilMSLiquid
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Equilibrium from density and momentum in auxField with thermodynamic factor
  recursive subroutine deriveEquilWTDF_MSLiquid(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, iField, depField, elemOff
    integer :: dens_pos, mom_pos(3)
    type(mus_varSys_data_type), pointer :: fPtr
    !mass density of nSpecies
    real(kind=rk) :: mass_dens( varSys%nStateVars )
    !number density of nSpecies
    real(kind=rk) :: num_dens( varSys%nStateVars )
    !mole fraction
    real(kind=rk) :: moleFraction( varSys%nStateVars )
    real(kind=rk) :: totNum_densInv
    !momentum from auxField
    real(kind=rk) :: momentum( 3 )
    real(kind=rk) :: vel( 3, varSys%nStateVars )
    !parameters from solver specific conf
    !field specific info from field table
    real(kind=rk), dimension(varSys%nStateVars, varSys%nStateVars) :: &
      & resi_coeff, thermodynamic_fac, &
      & inv_thermodyn_fac, diff_coeff
    !mixture info
    real(kind=rk) :: paramBInv
    real(kind=rk) :: fEq(fun%nComponents)
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )

    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                            &
      &        levelPointer => fPtr%solverData%geometry%levelPointer,       &
      &        auxField => fPtr%solverData%scheme%auxField,                 &
      &        mixture => fPtr%solverData%scheme%mixture,                   &
      &        physics => fPtr%solverData%physics,                          &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

      ! find dependent field of current variable by comparing the position
      ! of this variable against derVarPos of this variable for all fields in
      ! scheme
      depField = 0
      do iField = 1, scheme%nFields
        if ( fun%myPos == scheme%derVarpos(iField)%equilibrium ) &
          & depField = iField
      end do

      paramBInv = 1.0_rk / mixture%paramB

      do iField = 1, scheme%nFields
        ! diffusivity coefficients
        diff_coeff(iField, :) = species(iField)%diff_coeff(:)
      end do

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! element offset for auxField
        elemoff = (statePos-1)*varSys%nAuxScalars
        do iField = 1, scheme%nFields
          ! position of current field density in auxField array
          dens_pos = varSys%method%val( fun%input_varPos(iField) ) &
            &                     %auxField_varPos(1)

          ! position of current field momentum in auxField array
          mom_pos = varSys%method%val( fun%input_varPos(iField) ) &
            &                     %auxField_varPos(1:3)

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

          ! species momentum
          do iComp = 1, 3
            momentum(iComp) = auxField(iLevel)%val( elemOff + mom_pos(iComp) )
          end do

          ! species velocity
          vel(:, iField) = momentum / mass_dens(iField)

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

        end do !iField

        ! total number density Inv
        totNum_densInv = 1.0_rk/sum(num_dens(:))

        ! molefraction
        moleFraction(:) = num_dens(:)*totNum_densInv

        ! MS-Diff coeff matrix from C++ code
        call mus_calc_MS_DiffMatrix( nFields   = scheme%nFields,             &
          &                          temp      = mixture%temp0,              &
          &                          press     = mixture%atm_press,          &
          &                          mole_dens = num_dens*physics%moleDens0, &
          &                          D_ij_out  = diff_coeff                  )

        ! Convert to lattice unit
        resi_coeff = physics%fac(iLevel)%diffusivity/diff_coeff

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

        ! invert thermodynamic factor
        inv_thermodyn_fac = invert_matrix( thermodynamic_fac )

        ! compute equilibrium from macroscopic quantities
        fEq = equilFromMacroWTDF( iField            = depField,          &
          &                       mass_dens         = mass_dens,         &
          &                       moleFraction      = moleFraction,      &
          &                       velocity          = vel,               &
          &                       inv_thermodyn_fac = inv_thermodyn_fac, &
          &                       layout            = scheme%layout,     &
          &                       nFields           = scheme%nFields,    &
          &                       paramBInv         = paramBInv,         &
          &                       phi               = species(:)         &
          &                                           %molWeigRatio,     &
          &                       resi_coeff        = resi_coeff,        &
          &                       theta_eq          = mixture%theta_eq   )

        ! copy the results to the res
        res( (iElem-1)*fun%nComponents + 1 : iElem*fun%nComponents) = fEq

      end do ! iElem
    end associate

   end subroutine deriveEquilWTDF_MSLiquid
  ! ************************************************************************ !

  ! ************************************************************************ !
  !> mole fraction from density stored in auxField
  recursive subroutine deriveMoleFracMS(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, iLevel, depField, iField, elemOff
    integer :: dens_pos
    type(mus_varSys_data_type), pointer :: fPtr
    !number density of nSpecies
    real(kind=rk) :: num_dens( varSys%nStateVars )
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                      &
      &        levelPointer => fPtr%solverData%geometry%levelPointer, &
      &        auxField => fPtr%solverData%scheme%auxField,           &
      &        field => fPtr%solverData%scheme%field                  )

      ! find dependent field of current variable by comparing the position
      ! of this variable against derVarPos of this variable for all fields in
      ! scheme
      depField = 0
      do iField = 1, scheme%nFields
        if ( fun%myPos == scheme%derVarpos(iField)%moleFrac ) &
          & depField = iField
      end do

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state and auxField array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! element offset for auxField
        elemoff = (statePos-1)*varSys%nAuxScalars
        do iField = 1, scheme%nFields
          ! position of density field in auxField array
          dens_pos = varSys%method%val( scheme%derVarPos(iField)%density ) &
            &                     %auxField_varPos(1)

          ! number density of species: mass density / MolWeight
          num_dens(iField) = auxField(iLevel)%val( elemOff + dens_pos ) &
            &              * field( iField )%fieldProp%species%molWeightInv
        end do !iField

        ! Mole fraction for current field
        res(iElem) = num_dens(depField) / sum( num_dens(:) )

      end do ! iElem
    end associate

  end subroutine deriveMoleFracMS
  ! ************************************************************************ !

  ! ************************************************************************ !
  !> mass fraction from density stored in auxField
  recursive subroutine deriveMassFracMS(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, iLevel, depField, iField, elemOff
    integer :: dens_pos
    type(mus_varSys_data_type), pointer :: fPtr
    !density of nSpecies
    real(kind=rk) :: mass_dens( varSys%nStateVars )
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                      &
      &        levelPointer => fPtr%solverData%geometry%levelPointer, &
      &        auxField => fPtr%solverData%scheme%auxField,           &
      &        field => fPtr%solverData%scheme%field                  )

      ! find dependent field of current variable by comparing the position
      ! of this variable against derVarPos of this variable for all fields in
      ! scheme
      depField = 0
      do iField = 1, scheme%nFields
        if ( fun%myPos == scheme%derVarpos(iField)%massFrac ) &
          & depField = iField
      end do

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state and auxField array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! element offset for auxField
        elemoff = (statePos-1)*varSys%nAuxScalars
        do iField = 1, scheme%nFields
          ! position of density field in auxField array
          dens_pos = varSys%method%val( scheme%derVarPos(iField)%density ) &
            &                     %auxField_varPos(1)

          ! density of species
          mass_dens(iField) = auxField(iLevel)%val( elemOff + dens_pos )
        end do !iField

        ! mass fraction for current field
        res(iElem) = mass_dens(depField) / sum( mass_dens(:) )
      end do ! iElem
    end associate

  end subroutine deriveMassFracMS
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Compute mole flux from momentum stored in auxField.
  !! mole flux = numDens_i*velocity_i = momentum / molWeight
  recursive subroutine deriveMoleFluxMS(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, iField, depField
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: mom_pos(3)
    ! species equation
    real(kind=rk) :: momentum(3)
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                      &
      &        levelPointer => fPtr%solverData%geometry%levelPointer, &
      &        auxField => fPtr%solverData%scheme%auxField,           &
      &        field => fPtr%solverData%scheme%field                  )

      ! find dependent field of current variable by comparing the position
      ! of this variable against derVarPos of this variable for all fields in
      ! scheme
      depField = 0
      do iField = 1, scheme%nFields
        if ( fun%myPos == scheme%derVarpos(iField)%moleFlux ) &
          & depField = iField
      end do

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! position of current field momentum in auxField array
        mom_pos = varSys%method%val( fun%input_varPos(1) ) &
          &                     %auxField_varPos(1:3)

        ! species momentum
        do iComp = 1, 3
          momentum(iComp) = auxField(iLevel)%val(                            &
            &               (statePos-1)*varSys%nAuxScalars + mom_pos(iComp) )
        end do

        ! mole flux = momentum / molWeight
        res( (iElem-1)*3 + 1 : iElem*3) = momentum                  &
          &        * field( depField )%fieldProp%species%molWeightInv

      end do ! iElem
    end associate

  end subroutine deriveMoleFluxMS
  ! ************************************************************************ !

  ! ************************************************************************ !
  !> Calculate mixture velocity of a given element
  !! from the momentum and density stored in auxField array for liquid mixture.
  !! auxField was updated with momentum of untransformed PDF which was computed
  !! by solving LSE in compute kernel.
  recursive subroutine deriveMixVelMS(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, iLevel, elemOff, iField
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: dens_pos, mom_pos(3)
    ! species equation
    real(kind=rk) :: vel(3,varSys%nStateVars), mixVel(3), inv_rho
    real(kind=rk), dimension(varSys%nStateVars) :: mass_dens, massFrac
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                      &
      &        levelPointer => fPtr%solverData%geometry%levelPointer, &
      &        auxField => fPtr%solverData%scheme%auxField            )

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! element offset for auxField
        elemoff = (statePos-1)*varSys%nAuxScalars

        do iField = 1, scheme%nFields
          ! position of current field density in auxField array
          dens_pos = varSys%method%val( scheme%derVarPos(iField)%density ) &
            &                     %auxField_varPos(1)

          ! position of current field momentum in auxField array
          mom_pos = varSys%method%val( scheme%derVarPos(iField)%momentum ) &
            &                     %auxField_varPos(1:3)

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

          ! species momentum
          inv_rho = 1.0_rk / mass_dens(iField)
          vel(1, iField) = auxField(iLevel)%val( elemOff + mom_pos(1) )*inv_rho
          vel(2, iField) = auxField(iLevel)%val( elemOff + mom_pos(2) )*inv_rho
          vel(3, iField) = auxField(iLevel)%val( elemOff + mom_pos(3) )*inv_rho
        end do

        ! mass fraction
        massFrac = mass_dens / sum(mass_dens)

        ! mixture velocity: \sum_k y_k v_k
        mixVel(1) = dot_product(massFrac, vel(1, :))
        mixVel(2) = dot_product(massFrac, vel(2, :))
        mixVel(3) = dot_product(massFrac, vel(3, :))

        ! compute and store velocity
        res( (iElem-1)*3 + 1 : iElem*3) = mixVel(1:3)

      end do ! iElem
    end associate

  end subroutine deriveMixVelMS
  ! ************************************************************************ !


  ! ************************************************************************* !
  !         Subroutines with common interface for values from index           !
  ! ************************************************************************* !

  ! ************************************************************************** !
  !> Calculate mole density from species concentration for getValOfIndex
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]].
  !!
  recursive subroutine deriveMoleDensityMS_fromIndex(fun, varSys, time,   &
    &                                                iLevel, idx, idxLen, &
    &                                                nVals, 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

    !> 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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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
    integer :: dens_pos
    integer :: iField, depField
    ! -------------------------------------------------------------------- !


    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                            &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

      ! find dependent field of current variable by comparing the position
      ! of this variable against derVarPos of this variable for all fields in
      ! scheme
      depField = 0
      do iField = 1, scheme%nFields
        if ( fun%myPos == scheme%derVarpos(iField)%moleDensity ) &
          & depField = iField
      end do

      dens_pos = fun%input_varPos(1)
      ! get mass density values for IDX for iField
      call varSys%method%val( dens_pos )%get_ValOfIndex( &
        &     varSys = varSys,                           &
        &     time   = time,                             &
        &     iLevel = iLevel,                           &
        &     idx    = fPtr%opData%input_pntIndex(1)     &
        &              %indexLvl(iLevel)%val( idx(:) ),  &
        &     nVals  = nVals,                            &
        &     res    = res(:)                            )

      ! convert mass density to mole density
      res(1:nVals) = res(1:nVals) * species(depField)%molWeightInv

    end associate

  end subroutine deriveMoleDensityMS_fromIndex
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Calculate velocity from species density and momentum in auxField
  !! for getValOfIndex
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]].
  !!
  recursive subroutine deriveVelocityMS_fromIndex(fun, varSys, time, iLevel, &
    &                                             idx, idxLen, nVals, 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

    !> 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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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
    integer :: dens_pos, mom_pos, iVal
    real(kind=rk) :: mass_dens(nVals)
    ! -------------------------------------------------------------------- !


    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk

    ! get mass density values for IDX for iField
    dens_pos = fun%input_varPos(1)
    call varSys%method%val( dens_pos )%get_ValOfIndex( &
      &     varSys = varSys,                           &
      &     time   = time,                             &
      &     iLevel = iLevel,                           &
      &     idx    = fPtr%opData%input_pntIndex(1)     &
      &              %indexLvl(iLevel)%val( idx(:) ),  &
      &     nVals  = nVals,                            &
      &     res    = mass_dens                         )

    ! get species momentum values for IDX for iField
    mom_pos = fun%input_varPos(2)
    call varSys%method%val( mom_pos )%get_ValOfIndex( &
      &     varSys = varSys,                          &
      &     time   = time,                            &
      &     iLevel = iLevel,                          &
      &     idx    = fPtr%opData%input_pntIndex(2)    &
      &              %indexLvl(iLevel)%val( idx(:) ), &
      &     nVals  = nVals,                           &
      &     res    = res(:)                           )

    ! convert momentum to velocity
    do iVal = 1, nVals
      res( (iVal-1)*3 + 1 : iVal*3 ) =  res( (iVal-1)*3 + 1 : iVal*3 ) &
        &                            / mass_dens(iVal)
    end do

  end subroutine deriveVelocityMS_fromIndex
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Calculate mole flux from species momentum for getValOfIndex
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]].
  !!
  recursive subroutine deriveMoleFluxMS_fromIndex(fun, varSys, time, iLevel, &
    &                                             idx, idxLen, nVals, 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

    !> 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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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
    integer :: mom_pos
    integer :: iField, depField
    ! -------------------------------------------------------------------- !


    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                            &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

      ! find dependent field of current variable by comparing the position
      ! of this variable against derVarPos of this variable for all fields in
      ! scheme
      depField = 0
      do iField = 1, scheme%nFields
        if ( fun%myPos == scheme%derVarpos(iField)%moleDensity ) &
          & depField = iField
      end do

      mom_pos = fun%input_varPos(1)
      ! get mass density values for IDX for iField
      call varSys%method%val( mom_pos )%get_ValOfIndex(  &
        &     varSys = varSys,                           &
        &     time   = time,                             &
        &     iLevel = iLevel,                           &
        &     idx    = fPtr%opData%input_pntIndex(1)     &
        &              %indexLvl(iLevel)%val( idx(:) ),  &
        &     nVals  = nVals,                            &
        &     res    = res(:)                            )

      ! convert mass flux to mole flux
      res(:) = res(:) * species(depField)%molWeightInv

    end associate

  end subroutine deriveMoleFluxMS_fromIndex
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Calculate charge density from species concentration
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]].
  !!
  recursive subroutine deriveChargeDensity_fromIndex(fun, varSys, time,   &
    &                                                iLevel, idx, idxLen, &
    &                                                nVals, 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

    !> 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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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
    integer :: dens_pos
    integer :: iField
    real(kind=rk) :: mass_dens(nVals)
    ! -------------------------------------------------------------------- !


    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk
    associate( mixture => fPtr%solverData%scheme%mixture,                   &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

      do iField = 1, fun%nInputs
        dens_pos = fun%input_varPos(iField)
        ! get mass density values for IDX for iField
        call varSys%method%val( dens_pos )%get_ValOfIndex(  &
          &     varSys = varSys,                            &
          &     time   = time,                              &
          &     iLevel = iLevel,                            &
          &     idx    = fPtr%opData%input_pntIndex(iField) &
          &              %indexLvl(iLevel)%val( idx(:) ),   &
          &     nVals  = nVals,                             &
          &     res    = mass_dens(:)                       )

        ! compute charge density
        res(1:nVals) = res(1:nVals) + mass_dens(1:nVals) &
          &          * species(iField)%molWeightInv      &
          &          * species(iField)%chargeNr
      end do

      ! Multiply faraday constant
      res(1:nVals) = res(1:nVals) * mixture%faradayLB

    end associate

  end subroutine deriveChargeDensity_fromIndex
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Calculate current density from species momentum
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]].
  !!
  recursive subroutine deriveCurrentDensity_fromIndex(fun, varSys, time,   &
    &                                                 iLevel, idx, idxLen, &
    &                                                 nVals, 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

    !> 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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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
    integer :: mom_pos
    integer :: iField, iVal
    real(kind=rk) :: momentum(nVals*3)
    ! -------------------------------------------------------------------- !


    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk
    associate( mixture => fPtr%solverData%scheme%mixture,                   &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

      do iField = 1, fun%nInputs
        mom_pos = fun%input_varPos(iField)
        ! get mass density values for IDX for iField
        call varSys%method%val( mom_pos )%get_ValOfIndex(   &
          &     varSys = varSys,                            &
          &     time   = time,                              &
          &     iLevel = iLevel,                            &
          &     idx    = fPtr%opData%input_pntIndex(iField) &
          &              %indexLvl(iLevel)%val( idx(:) ),   &
          &     nVals  = nVals,                             &
          &     res    = momentum(:)                        )

        ! compute current density
        do iVal = 1, nVals
          res(:) = res(:) + momentum(:) * species(iField)%molWeightInv &
            &    * species(iField)%chargeNr
        end do ! iVal
      end do ! iField

      ! Multiply faraday constant
      res(:) = res(:) * mixture%faradayLB

    end associate

  end subroutine deriveCurrentDensity_fromIndex
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Calculate mixture velocity from density from species momentum
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]].
  !!
  recursive subroutine deriveMixVelMS_fromIndex(fun, varSys, time, iLevel, &
    &                                           idx, idxLen, nVals, 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

    !> 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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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(:)
    ! -------------------------------------------------------------------- !

    call tem_abort('deriveMixVelMS_fromIndex is not implemented yet!')

  end subroutine deriveMixVelMS_fromIndex
  ! ************************************************************************** !

  ! ************************************************************************* !
  !         Subroutines with common interface for apply source                !
  ! ************************************************************************* !

  ! ************************************************************************ !
  !> Update state with source variable "electric_field" with 2nd order force
  !! integration.
  !! Refer to Appendix in PhD Thesis of K. Masilamani
  !! "Coupled Simulation Framework to Simulate Electrodialysis Process for
  !! Seawater Desalination"
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[proc_apply_source]] in derived/[[mus_source_var_module]].f90 in order to
  !! be callable via [[mus_source_op_type:applySrc]] function pointer.
  subroutine applySrc_electricMSLiquid_2ndOrd( fun, inState, outState, neigh, &
    &                                          auxField, nPdfSize, iLevel,    &
    &                                          varSys, time, phyConvFac,      &
    &                                          derVarPos                      )
    ! -------------------------------------------------------------------- !
    !> Description of method to apply source terms
    class(mus_source_op_type), intent(in) :: fun

    !> input  pdf vector
    real(kind=rk), intent(in) :: inState(:)

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

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

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

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

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

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

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

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

    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    real(kind=rk) :: electricField(fun%elemLvl(iLevel)%nElems*3)
    real(kind=rk) :: EF_elem(3)
    integer :: iElem, nElems, iDir, posInTotal
    integer :: iField, depField, nScalars, QQ, nInputStates
    !mass density of nSpecies
    real(kind=rk) :: mass_dens( varSys%nStateVars )
    real(kind=rk) :: massFrac( varSys%nStateVars )
    real(kind=rk) :: charge_dens, diffForce_cs2inv, diffForce_cs2inv_sqr
    real(kind=rk) :: omegaTerm, mixVel(3), inv_rho, ucx, uMinusCX(3)
    real(kind=rk), dimension(varSys%nStateVars) :: chargeTerm
    real(kind=rk) :: minMolWeight, forceTerm
    real(kind=rk), dimension(3, varSys%nStateVars ) :: spcForce, vel
    integer :: statePos, dens_pos, mom_pos(3), elemOff
    ! -------------------------------------------------------------------- !
!write(dbgUnit(1),*) 'source variable: ', trim(varSys%varname%val(fun%srcTerm_varPos))
    ! convert c pointer to solver type fortran pointer
    call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, &
      &               fPtr )

    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    associate( scheme => fPtr%solverData%scheme,                            &
      &        auxField => fPtr%solverData%scheme%auxField,                 &
      &        mixture => fPtr%solverData%scheme%mixture,                   &
      &        physics => fPtr%solverData%physics,                          &
      &        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 )
      diffForce_cs2inv_sqr = diffForce_cs2inv * diffForce_cs2inv

      ! omega term to multiply forceTerm
      omegaTerm = 1.0_rk                                                  &
        &       / ( 1.0_rk + mixture%relaxLvl(iLevel)%omega_diff * 0.5_rk )

      ! 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

      QQ = scheme%layout%fStencil%QQ
      nScalars = varSys%nScalars

      ! update source for each element
      do iElem = 1, nElems

        ! to access level wise state array
        posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem)
        ! element offset for auxField
        elemoff = (posInTotal-1)*varSys%nAuxScalars

        ! get mass density from auxField
        do iField = 1, scheme%nFields
          ! position of current field density in auxField array
          dens_pos = varSys%method%val( scheme%derVarPos(iField)%density ) &
            &                     %auxField_varPos(1)

          ! mass density of species
          mass_dens(iField) = auxField(iLevel)%val( 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

          ! position of current field momentum in auxField array
          mom_pos = varSys%method%val( scheme%derVarPos(iField)%momentum ) &
            &                    %auxField_varPos(1:3)

          inv_rho = 1.0_rk / mass_dens(iField)
          ! species velocity
          vel(1, iField) = auxField(iLevel)%val(elemOff + mom_pos(1)) * inv_rho
          vel(2, iField) = auxField(iLevel)%val(elemOff + mom_pos(2)) * inv_rho
          vel(3, iField) = auxField(iLevel)%val(elemOff + mom_pos(3)) * inv_rho

        end do !iField

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

        ! Mixture velocity
        mixVel(1) = dot_product(massFrac(:), vel(1,:))
        mixVel(2) = dot_product(massFrac(:), vel(2,:))
        mixVel(3) = dot_product(massFrac(:), vel(3,:))

        ! 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)
        ! compute electrical migrating force each species
        ! F_k = (\rho_k z_k/M_k - y_k \sum_l \rho_l z_l / M_l) F E / (RT)
        ! Above term is multiplied by minMolWeight which comes from lattice
        ! force term
        do iField = 1, scheme%nFields
          spcForce(:, iField) = EF_elem(:) * diffForce_cs2inv * cs2 &
            & * (chargeTerm(iField) - massFrac(iField) * charge_dens )
        end do

        ! compute external forcing term
        ! d^m_k = w_m*c_m*( min_a(m_a)*F_k  )
        ! F_k is diffusive forcing term
        !
        ! Update souce depends on nInputStates
        ! if nInputStates = 1, it is field source else it is global source
        do iField = 1, nInputStates
          depField = varSys%method%val(fun%srcTerm_varPos)%input_varPos(iField)

          do iDir = 1, QQ
            ucx = dot_product( scheme%layout%fStencil%cxDirRK(:, iDir), &
              &                mixVel(:) )
            uMinusCX = scheme%layout%fStencil%cxDirRK(:, iDir) - mixVel(:)

            forceTerm = dot_product( uMinusCx * cs2inv               &
              &       + ucx * scheme%layout%fStencil%cxDirRK(:,iDir) &
              &       * cs4inv, spcForce(:, depField)                )
            !forceTerm = dot_product( scheme%layout%fStencil%cxDirRK(:, iDir), &
            !  &                      spcForce(:, depField) )

            statePos                                                        &
              & = (posintotal-1)*nscalars+idir+(depfield-1)*qq

            outState( statePos ) = outState( statePos )              &
              & + omegaTerm * scheme%layout%weight( iDir ) * forceTerm
          end do ! iDir

        end do !iField
      end do !iElem
    end associate

  end subroutine applySrc_electricMSLiquid_2ndOrd
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Update state with source variable "electric_field"
  !! Simuilar to derive routine but it updates the state whereas derive
  !! is used for tracking
  !! @todo species electricField
  !! Refer to Appendix in PhD Thesis of K. Masilamani
  !! "Coupled Simulation Framework to Simulate Electrodialysis Process for
  !! Seawater Desalination"
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[proc_apply_source]] in derived/[[mus_source_var_module]].f90 in order to
  !! be callable via [[mus_source_op_type:applySrc]] function pointer.
  subroutine applySrc_electricMSLiquid_1stOrd( fun, inState, outState, neigh, &
    &                                          auxField, nPdfSize, iLevel,    &
    &                                          varSys, time, phyConvFac,      &
    &                                          derVarPos                      )
    ! -------------------------------------------------------------------- !
    !> Description of method to apply source terms
    class(mus_source_op_type), intent(in) :: fun

    !> input  pdf vector
    real(kind=rk), intent(in) :: inState(:)

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

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

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

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

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

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

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

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

    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    real(kind=rk) :: electricField(fun%elemLvl(iLevel)%nElems*3)
    real(kind=rk) :: EF_elem(3)
    integer :: iElem, nElems, iDir, posInTotal
    integer :: iField, depField, nScalars, QQ, nInputStates
    !mass density of nSpecies
    real(kind=rk) :: mass_dens( varSys%nStateVars )
    real(kind=rk) :: massFrac( varSys%nStateVars )
    real(kind=rk) :: charge_dens, diffForce_cs2inv
    real(kind=rk), dimension(varSys%nStateVars) :: chargeTerm
    real(kind=rk) :: minMolWeight, forceTerm
    real(kind=rk), dimension(3, varSys%nStateVars ) :: spcForce
    integer :: dens_pos, elemOff
    ! -------------------------------------------------------------------- !
!write(dbgUnit(1),*) 'source variable: ', trim(varSys%varname%val(fun%srcTerm_varPos))
    ! convert c pointer to solver type fortran pointer
    call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, &
      &               fPtr )

    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    associate( scheme => fPtr%solverData%scheme,                            &
      &        auxField => fPtr%solverData%scheme%auxField,                 &
      &        mixture => fPtr%solverData%scheme%mixture,                   &
      &        physics => fPtr%solverData%physics,                          &
      &        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 )

      ! 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

      QQ = scheme%layout%fStencil%QQ
      nScalars = varSys%nScalars

      ! update source for each element
      do iElem = 1, nElems

        ! to access level wise state array
        posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem)
        ! element offset for auxField
        elemoff = (posInTotal-1)*varSys%nAuxScalars

        ! get mass density from auxField
        do iField = 1, scheme%nFields
          ! position of current field density in auxField array
          dens_pos = varSys%method%val( scheme%derVarPos(iField)%density ) &
            &                     %auxField_varPos(1)

          ! mass density of species
          mass_dens(iField) = auxField(iLevel)%val( 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 !iField

        !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)

        ! compute electrical migrating force each species
        ! F_k = (\rho_k z_k/M_k - y_k \sum_l \rho_l z_l / M_l) F E / (RT)
        ! Above term is multiplied by minMolWeight which comes from lattice
        ! force term
        do iField = 1, scheme%nFields
          spcForce(:, iField) = EF_elem(:) * diffForce_cs2inv        &
            & * (chargeTerm(iField) - massFrac(iField) * charge_dens )
        end do

        ! compute external forcing term
        ! d^m_k = w_m*c_m*( min_a(m_a)*F_k  )
        ! F_k is diffusive forcing term
        !
        ! Update souce depends on nInputStates
        ! if nInputStates = 1, it is field source else it is global source
        do iField = 1, nInputStates
          depField = varSys%method%val(fun%srcTerm_varPos)%input_varPos(iField)

          do iDir = 1, QQ
            forceTerm = scheme%layout%fStencil%cxDirRK( 1, iDir ) &
              &         * spcForce(1, depField)                   &
              &       + scheme%layout%fStencil%cxDirRK( 2, iDir ) &
              &         * spcForce(2, depField)                   &
              &       + scheme%layout%fStencil%cxDirRK( 3, iDir ) &
              &         * spcForce(3, depField)

            outState(                                                         &
              & (posintotal-1)*nscalars+idir+(depfield-1)*qq ) &
              & = outState(                                                   &
              & (posintotal-1)*nscalars+idir+(depfield-1)*qq ) &
              & + scheme%layout%weight( iDir ) * forceTerm
          end do ! iDir

        end do !iField
      end do !iElem
    end associate

  end subroutine applySrc_electricMSLiquid_1stOrd
  ! ************************************************************************ !

  ! ************************************************************************ !
  !> Update state with source variable "electric_field" with 2nd order
  !! integration of force term in LBE with thermodynamic factor
  !! Refer to Appendix in PhD Thesis of K. Masilamani
  !! "Coupled Simulation Framework to Simulate Electrodialysis Process for
  !! Seawater Desalination"
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[proc_apply_source]] in derived/[[mus_source_var_module]].f90 in order to
  !! be callable via [[mus_source_op_type:applySrc]] function pointer.
  subroutine applySrc_electricMSLiquid_2ndOrd_WTDF( fun, inState, outState,    &
    &                                               neigh, auxField, nPdfSize, &
    &                                               iLevel, varSys, time,      &
    &                                               phyConvFac, derVarPos      )
    ! -------------------------------------------------------------------- !
    !> Description of method to apply source terms
    class(mus_source_op_type), intent(in) :: fun

    !> input  pdf vector
    real(kind=rk), intent(in) :: inState(:)

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

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

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

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

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

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

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

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

    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    real(kind=rk) :: electricField(fun%elemLvl(iLevel)%nElems*3)
    real(kind=rk) :: EF_elem(3)
    integer :: iElem, nElems, iDir, posInTotal
    integer :: iField, iField_2, depField, nScalars, QQ, nInputStates
    !mass density of nSpecies
    real(kind=rk) :: mass_dens( varSys%nStateVars )
    real(kind=rk) :: massFrac( varSys%nStateVars )
    real(kind=rk) :: num_dens( varSys%nStateVars )
    real(kind=rk) :: moleFrac( varSys%nStateVars )
    real(kind=rk) :: charge_dens, diffForce_cs2inv, diffForce_cs2inv_sqr
    real(kind=rk) :: omegaTerm, mixVel(3), inv_rho, ucx, uMinusCX(3)
    real(kind=rk), dimension(varSys%nStateVars) :: chargeTerm
    real(kind=rk) :: minMolWeight, forceTerm
    real(kind=rk), dimension(3, varSys%nStateVars ) :: spcForce, vel, &
      & spcForce_WTDF
    real(kind=rk), dimension(varSys%nStateVars, varSys%nStateVars) :: &
      & thermodynamic_fac, inv_thermodyn_fac
    integer :: statePos, dens_pos, mom_pos(3), elemOff
    ! -------------------------------------------------------------------- !
!write(dbgUnit(1),*) 'source variable: ', trim(varSys%varname%val(fun%srcTerm_varPos))
    ! convert c pointer to solver type fortran pointer
    call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, &
      &               fPtr )

    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    associate( scheme => fPtr%solverData%scheme,                            &
      &        auxField => fPtr%solverData%scheme%auxField,                 &
      &        mixture => fPtr%solverData%scheme%mixture,                   &
      &        physics => fPtr%solverData%physics,                          &
      &        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 )
      diffForce_cs2inv_sqr = diffForce_cs2inv * diffForce_cs2inv

      ! omega term to multiply forceTerm
      omegaTerm = 1.0_rk/(1.0_rk + mixture%relaxLvl(iLevel)%omega_diff * 0.5_rk)

      ! 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

      QQ = scheme%layout%fStencil%QQ
      nScalars = varSys%nScalars

      ! update source for each element
      do iElem = 1, nElems

        ! to access level wise state array
        posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem)
        ! element offset for auxField
        elemoff = (posInTotal-1)*varSys%nAuxScalars

        ! get mass density from auxField
        do iField = 1, scheme%nFields
          ! position of current field density in auxField array
          dens_pos = varSys%method%val( scheme%derVarPos(iField)%density ) &
            &                     %auxField_varPos(1)

          ! mass density of species
          mass_dens(iField) = auxField(iLevel)%val( 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

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

          ! position of current field momentum in auxField array
          mom_pos = varSys%method%val( scheme%derVarPos(iField)%momentum ) &
            &                    %auxField_varPos(1:3)

          inv_rho = 1.0_rk / mass_dens(iField)
          ! species velocity
          vel(1, iField) = auxField(iLevel)%val(elemOff + mom_pos(1)) * inv_rho
          vel(2, iField) = auxField(iLevel)%val(elemOff + mom_pos(2)) * inv_rho
          vel(3, iField) = auxField(iLevel)%val(elemOff + mom_pos(3)) * inv_rho

        end do !iField

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

        ! Mixture velocity
        mixVel(1) = dot_product(massFrac(:), vel(1,:))
        mixVel(2) = dot_product(massFrac(:), vel(2,:))
        mixVel(3) = dot_product(massFrac(:), vel(3,:))

        ! 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)
        ! compute electrical migrating force each species
        ! F_k = (\rho_k z_k/M_k - y_k \sum_l \rho_l z_l / M_l) F E / (RT)
        ! Above term is multiplied by minMolWeight which comes from lattice
        ! force term
        do iField = 1, scheme%nFields
          spcForce(:, iField) = EF_elem(:)                           &
            & * (chargeTerm(iField) - massFrac(iField) * charge_dens )
        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 = w_m*c_m*( \sum_l \gamma^{-1}_{k,l} min_a(m_a)*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

        ! compute external forcing term
        ! d^m_k = w_m*c_m*( min_a(m_a)*F_k  )
        ! F_k is diffusive forcing term
        !
        ! Update souce depends on nInputStates
        ! if nInputStates = 1, it is field source else it is global source
        do iField = 1, nInputStates
          depField = varSys%method%val(fun%srcTerm_varPos)%input_varPos(iField)

          do iDir = 1, QQ
            ucx = dot_product( scheme%layout%fStencil%cxDirRK(:, iDir), &
              &                mixVel )
            uMinusCX = scheme%layout%fStencil%cxDirRK(:, iDir) - mixVel

            forceTerm = dot_product( uMinusCx * diffForce_cs2inv         &
              &       + ucx * scheme%layout%fStencil%cxDirRK(:,iDir)     &
              &       * diffForce_cs2inv_sqr, spcForce_WTDF(:, depField) )

            statePos                                                        &
              & = (posintotal-1)*nscalars+idir+(depfield-1)*qq

            outState( statePos ) = outState( statePos )              &
              & + omegaTerm * scheme%layout%weight( iDir ) * forceTerm
          end do ! iDir

        end do !iField
      end do !iElem
    end associate

  end subroutine applySrc_electricMSLiquid_2ndOrd_WTDF
  ! ************************************************************************ !

  ! ************************************************************************ !
  !> Update state with source variable "electric_field" with thermodynamic
  !! factor.
  !! Simuilar to derive routine but it updates the state whereas derive
  !! is used for tracking
  !! @todo species electricField
  !! Refer to Appendix in PhD Thesis of K. Masilamani
  !! "Coupled Simulation Framework to Simulate Electrodialysis Process for
  !! Seawater Desalination"
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[proc_apply_source]] in derived/[[mus_source_var_module]].f90 in order to
  !! be callable via [[mus_source_op_type:applySrc]] function pointer.
  subroutine applySrc_electricMSLiquid_1stOrd_WTDF( fun, inState, outState,    &
    &                                               neigh, auxField, nPdfSize, &
    &                                               iLevel, varSys, time,      &
    &                                               phyConvFac, derVarPos      )
    ! -------------------------------------------------------------------- !
    !> Description of method to apply source terms
    class(mus_source_op_type), intent(in) :: fun

    !> input  pdf vector
    real(kind=rk), intent(in) :: inState(:)

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

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

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

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

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

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

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

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

    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    real(kind=rk) :: electricField(fun%elemLvl(iLevel)%nElems*3)
    real(kind=rk) :: EF_elem(3)
    integer :: iElem, nElems, iDir, posInTotal
    integer :: iField, iField_2, depField, nScalars, QQ, nInputStates
    !mass density of nSpecies
    real(kind=rk) :: mass_dens( varSys%nStateVars )
    real(kind=rk) :: num_dens( varSys%nStateVars )
    real(kind=rk) :: massFrac( varSys%nStateVars )
    real(kind=rk) :: moleFrac( varSys%nStateVars )
    real(kind=rk) :: charge_dens, diffForce_cs2inv
    real(kind=rk), dimension(varSys%nStateVars) :: chargeTerm
    real(kind=rk) :: minMolWeight, forceTerm
    real(kind=rk), dimension(3, varSys%nStateVars ) :: spcForce, spcForce_WTDF
    real(kind=rk), dimension(varSys%nStateVars, varSys%nStateVars) :: &
      & thermodynamic_fac, inv_thermodyn_fac
    integer :: dens_pos, elemOff
    ! -------------------------------------------------------------------- !
!write(dbgUnit(1),*) 'source variable: ', trim(varSys%varname%val(fun%srcTerm_varPos))
    ! convert c pointer to solver type fortran pointer
    call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, &
      &               fPtr )

    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    associate( scheme => fPtr%solverData%scheme,                            &
      &        auxField => fPtr%solverData%scheme%auxField,                 &
      &        mixture => fPtr%solverData%scheme%mixture,                   &
      &        physics => fPtr%solverData%physics,                          &
      &        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 )

      ! 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

      QQ = scheme%layout%fStencil%QQ
      nScalars = varSys%nScalars

      ! update source for each element
      do iElem = 1, nElems

        ! to access level wise state array
        posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem)
        ! element offset for auxField
        elemoff = (posInTotal-1)*varSys%nAuxScalars

        ! get mass density from auxField
        do iField = 1, scheme%nFields
          ! position of current field density in auxField array
          dens_pos = varSys%method%val( scheme%derVarPos(iField)%density ) &
            &                     %auxField_varPos(1)

          ! mass density of species
          mass_dens(iField) = auxField(iLevel)%val( 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 !iField

        !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)

        ! compute electrical migrating force each species
        ! F_k = (\rho_k z_k/M_k - y_k \sum_l \rho_l z_l / M_l) F E / (RT)
        ! Above term is multiplied by minMolWeight which comes from lattice
        ! force term
        do iField = 1, scheme%nFields
          spcForce(:, iField) = EF_elem(:) * diffForce_cs2inv        &
            & * (chargeTerm(iField) - massFrac(iField) * charge_dens )
        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 = w_m*c_m*( \sum_l \gamma^{-1}_{k,l} min_a(m_a)*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 souce depends on nInputStates
        ! if nInputStates = 1, it is field source else it is global source
        do iField = 1, nInputStates
          depField = varSys%method%val(fun%srcTerm_varPos)%input_varPos(iField)

          do iDir = 1, QQ
            forceTerm = scheme%layout%fStencil%cxDirRK( 1, iDir ) &
              &         * spcForce_WTDF(1, depField)              &
              &       + scheme%layout%fStencil%cxDirRK( 2, iDir ) &
              &         * spcForce_WTDF(2, depField)              &
              &       + scheme%layout%fStencil%cxDirRK( 3, iDir ) &
              &         * spcForce_WTDF(3, depField)

            outState(                                                         &
              & (posintotal-1)*nscalars+idir+(depfield-1)*qq ) &
              & = outState(                                                   &
              & (posintotal-1)*nscalars+idir+(depfield-1)*qq ) &
              & + scheme%layout%weight( iDir ) * forceTerm
          end do ! iDir

        end do !iField
      end do !iElem
    end associate

  end subroutine applySrc_electricMSLiquid_1stOrd_WTDF
  ! ************************************************************************ !



  ! ************************************************************************ !
  !> Update state with source variable "force" with 2nd order integration
  !! of force in lattice Boltzmann equation.
  !! Simuilar to derive routine but it updates the state whereas derive
  !! is used for tracking
  !! Refer to Appendix in PhD Thesis of K. Masilamani
  !! "Coupled Simulation Framework to Simulate Electrodialysis Process for
  !! Seawater Desalination"
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[proc_apply_source]] in derived/[[mus_source_var_module]].f90 in order to
  !! be callable via [[mus_source_op_type:applySrc]] function pointer.
  subroutine applySrc_forceMSLiquid_2ndOrd( fun, inState, outState, neigh, &
    &                                       auxField, nPdfSize, iLevel,    &
    &                                       varSys, time, phyConvFac,      &
    &                                       derVarPos                      )
    ! -------------------------------------------------------------------- !
    !> Description of method to apply source terms
    class(mus_source_op_type), intent(in) :: fun

    !> input  pdf vector
    real(kind=rk), intent(in) :: inState(:)

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

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

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

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

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

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

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

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

    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3)
    integer :: iElem, nElems, iDir, posInTotal
    integer :: iField, depField, nScalars, QQ, nInputStates
    !mass density of nSpecies
    real(kind=rk) :: mass_dens( varSys%nStateVars )
    real(kind=rk) :: massFrac( varSys%nStateVars )
    real(kind=rk) :: forceTerm, force_elem(3), ucx, uMinusCX(3)
    real(kind=rk), dimension(3, varSys%nStateVars ) :: spcForce, vel
    real(kind=rk) :: omegaTerm, mixVel(3), inv_rho
    integer :: statePos, dens_pos, mom_pos(3), elemOff
    ! -------------------------------------------------------------------- !
    !write(*,*) 'source variable: ', trim(varSys%varname%val(fun%srcTerm_varPos))
    ! convert c pointer to solver type fortran pointer
    call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, &
      &               fPtr )

    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    associate( scheme => fPtr%solverData%scheme,                            &
      &        auxField => fPtr%solverData%scheme%auxField,                 &
      &        mixture => fPtr%solverData%scheme%mixture,                   &
      &        physics => fPtr%solverData%physics,                          &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species )

      ! Get body 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 / physics%fac(iLevel)%body_force

      ! omega term to multiply forceTerm
      omegaTerm = 1.0_rk/(1.0_rk + mixture%relaxLvl(iLevel)%omega_kine * 0.5_rk)

      ! 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

      QQ = scheme%layout%fStencil%QQ
      nScalars = varSys%nScalars

      ! update source for each element
      do iElem = 1, nElems

        ! to access level wise state array
        posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem)

        ! element offset for auxField
        elemoff = (posInTotal-1)*varSys%nAuxScalars

        ! get mass density from auxField
        do iField = 1, scheme%nFields
          ! position of current field density in auxField array
          dens_pos = varSys%method%val( scheme%derVarPos(iField)%density ) &
            &                     %auxField_varPos(1)

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

          ! position of current field momentum in auxField array
          mom_pos = varSys%method%val( scheme%derVarPos(iField)%momentum ) &
            &                    %auxField_varPos(1:3)

          inv_rho = 1.0_rk / mass_dens(iField)
          ! species velocity
          vel(1, iField) = auxField(iLevel)%val(elemOff + mom_pos(1)) * inv_rho
          vel(2, iField) = auxField(iLevel)%val(elemOff + mom_pos(2)) * inv_rho
          vel(3, iField) = auxField(iLevel)%val(elemOff + mom_pos(3)) * inv_rho

        end do !iField

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

        ! Mixture velocity
        mixVel(1) = dot_product(massFrac(:), vel(1,:))
        mixVel(2) = dot_product(massFrac(:), vel(2,:))
        mixVel(3) = dot_product(massFrac(:), vel(3,:))

        ! force field for current element
        force_elem = forceField((iElem-1)*3+1 : iElem*3)

        ! compute external force for each species
        ! F_k = y_k F, F - body force per unit volume of form
        ! \rho g or \rho_e E.
        ! Above term is multiplied by cs2inv which comes from lattice
        ! force term
        do iField = 1, scheme%nFields
          spcForce(:, iField) =  massFrac(iField) * force_elem(:)
        end do

        ! compute external forcing term
        ! d^m_k = w_m*c_m*( F_k / cs2  )
        ! F_k is the external forcing term
        !
        ! Update souce depends on nInputStates
        ! if nInputStates = 1, it is field source else it is global source
        do iField = 1, nInputStates
          depField = varSys%method%val(fun%srcTerm_varPos)%input_varPos(iField)

          do iDir = 1, QQ
            ucx = dot_product( scheme%layout%fStencil%cxDirRK(:, iDir), &
              &                mixVel )
            uMinusCX = scheme%layout%fStencil%cxDirRK(:, iDir) - mixVel

            forceTerm = dot_product( uMinusCx * cs2inv               &
              &       + ucx * scheme%layout%fStencil%cxDirRK(:,iDir) &
              &       * cs4inv, spcForce(:, depField)                )

            statePos                                                        &
              & = (posintotal-1)*nscalars+idir+(depfield-1)*qq
            outState(statePos) = outState(statePos)                &
              & + omegaTerm * scheme%layout%weight(iDir) * forceTerm
          end do ! iDir

        end do !iField
      end do !iElem
    end associate

  end subroutine applySrc_forceMSLiquid_2ndOrd
  ! ************************************************************************ !

  ! ************************************************************************ !
  !> Update state with source variable "force" with 1st order integration
  !! of force in lattice Boltzmann equation.
  !! Simuilar to derive routine but it updates the state whereas derive
  !! is used for tracking
  !! Refer to Appendix in PhD Thesis of K. Masilamani
  !! "Coupled Simulation Framework to Simulate Electrodialysis Process for
  !! Seawater Desalination"
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[proc_apply_source]] in derived/[[mus_source_var_module]].f90 in order to
  !! be callable via [[mus_source_op_type:applySrc]] function pointer.
  subroutine applySrc_forceMSLiquid_1stOrd( fun, inState, outState, neigh, &
    &                                       auxField, nPdfSize, iLevel,    &
    &                                       varSys, time, phyConvFac,      &
    &                                       derVarPos                      )
    ! -------------------------------------------------------------------- !
    !> Description of method to apply source terms
    class(mus_source_op_type), intent(in) :: fun

    !> input  pdf vector
    real(kind=rk), intent(in) :: inState(:)

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

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

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

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

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

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

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

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

    !> position of derived quantities in varsys
    type(mus_derVarPos_type