! Copyright (c) 2011-2016, 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2011, 2017, 2019 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2011-2014 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) 2012 Sathish Krishnan P S <s.krishnan@grs-sim.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2018 Raphael Haupt <raphael.haupt@uni-siegen.de>
! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
?? include 'header/lbm_macros.inc'
?? include 'header/lbm_d3q19Macros.inc'
! ****************************************************************************** !
!> author: Jiaxing Qi
!! This module provides the definition and methods for
!! MRT advection relaxation scheme.
!! The LB equaton using MRT is
!!    f(t+dt,x+dx) = f - M^(-1) * S * ( (M*f) - m^(eq) )
!!
!! The moments m(1:19) = M * f(1:19) are labeled as
!!  m( 1) = rho
!!  m( 2) = e = rho * (ux^2 + uy^2 + uz^2)
!!  m( 3) = epsilon
!!  m( 4) = jx = rho * ux
!!  m( 5) = qx
!!  m( 6) = jy = rho * uy
!!  m( 7) = qy
!!  m( 8) = jz = rho * uz
!!  m( 9) = qz
!!  m(10) = 3 * pxx = rho * (2ux^2 - uy^2 - uz^2)
!!  m(11) = 3 * Pixx
!!  m(12) = pzz  = rho * (uy^2 - uz^2)
!!  m(13) = Piww
!!  m(14) = pxy  = rho * ux * uy
!!  m(15) = pyz  = rho * uy * uz
!!  m(16) = pzx  = rho * uz * ux
!!  m(17) = mx
!!  m(18) = my
!!  m(19) = mz
!!
!! The non-zero equilibirium moments are given by
!!  meq( 1) = rho
!!  meq( 2) = rho0 * ( ux^2 + uy^2 + uz^2 )
!!  meq( 4) = rho0 * ux
!!  meq( 6) = rho0 * uy
!!  meq( 8) = rho0 * uz
!!  meq(10) = rho0 * ( 2*ux^2 - uy^2 - uz^2 )
!!  meq(12) = rho0 * ( uy^2 - uz^2 )
!!  meq(14) = rho0 * ux * uy
!!  meq(15) = rho0 * uy * uz
!!  meq(16) = rho0 * ux * uz
!!
!! Density (rho) and velocity (ux, uy, uz) are conserved during collision.
!!  i.e. m(1) = meq(1) --> mneq(1) = 0
!!       m(4) = meq(4) --> mneq(4) = 0
!!       m(6) = meq(6) --> mneq(6) = 0
!!       m(8) = meq(8) --> mneq(8) = 0
!!
!! The collision parameters S correspondes to the omega in BGK model.
!!
!! The MRT implementation here is taken from:\n
!! J. Toelke, S. Freudiger, and M. Krafczyk,
!! "An adaptive scheme using hierarchical grids for lattice Boltzmann
!! multi-phase flow simulations," Comput. Fluids, vol. 35, pp. 820–830,
!! 2006. \n
!! Notice that the collision matrix S used in this papar corresponds to
!! -omega in BGK model, because it express the LB equation is slightly
!! different way.
!! In this paper, the following notions are used:\n
!!  s(a) = s(2)
!!  s(b) = s(3)
!!  s(c) = s(5) = s(7) = s(9)
!!  s(d) = s(11) = s(13
!!  s(e) = s(17) = s(18) = s(19)
!!  s(w) = s(10) = s(12) = s(14) = s(15) = s(16)
!! It is suggested that, for D3Q19,
!!  s(a) = s(b) = s(c) = s(d) = s(e) = max( s(w), -1.0 )
!!
!! SubGrid Stress model (SGS)
!! The implementation here is taken from:\n
!! M. Stiebler, M. Krafczyk, S. Freudiger, M. Geier
!! "Lattice Boltzmann large eddy simulation of subcritical flows around a sphere
!! on non-uniform grids", Computers and Mathematics with Applications, vol. 61
!! (2011), pp. 3475-3484
!! Equation 12:\n
!! tau_{total} = 3 * nu0 + dt * 0.5
!!               + 0.5 * ( sqrt( tau0*tau0  + 18 * Cs * Cs * dt * dt * Q ) - tau0 )
!!             = 0.5 * ( tau0 + sqrt( tau0 * tau0 + 18 * Cs * Cs * dt * dt * Q) )
!! Q = sqrt( 2.0 * sum( Pi^{neq} * Pi^{neq} ) )
!!
!! For single field LBM: QQ=nScalars
!!
module mus_mrt_module
  use iso_c_binding, only: c_f_pointer

  ! include treelm modules
  use env_module,               only: rk
  use tem_varSys_module,        only: tem_varSys_type, tem_varSys_op_type
  use tem_param_module,         only: div1_2, div1_3, div1_4, div1_6, div1_8, &
    &                                 div1_12, div1_16, div1_18, div1_24, div1_36,&
    &                                 div1_48, div1_72, &
    &                                 cs2inv, cs4inv, t2cs2inv, t2cs4inv, rho0

  ! include musubi modules
  use mus_field_prop_module,    only: mus_field_prop_type
  use mus_scheme_layout_module, only: mus_scheme_layout_type
  use mus_param_module,         only: mus_param_type
  use mus_varSys_module,        only: mus_varSys_data_type
  use mus_derVarPos_module,     only: mus_derVarPos_type
  use mus_mrtRelaxation_module, only: mrt_d3q19

  implicit none

  private

  public :: mrt_advRel_d3q19
  public :: mrt_advRel_d3q19_generic

  public :: mrt_advRel_d3q19_incomp
  public :: mrt_advRel_d3q19_incomp_generic


  public :: mrt_advRel_generic
  public :: mrt_advRel_incomp_generic

  !=============================================================================
  ! D3Q19 flow model
  !=============================================================================
  !> Definition of the discrete velocity set

  integer,parameter :: QQ   = 19   !< number of pdf directions

  integer,parameter :: qN00 = 1     !< west             x-
  integer,parameter :: q0N0 = 2     !< south            y-
  integer,parameter :: q00N = 3     !< bottom           z-
  integer,parameter :: q100 = 4     !< east             x+
  integer,parameter :: q010 = 5     !< north            y+
  integer,parameter :: q001 = 6     !< top              z+
  integer,parameter :: q0NN = 7     !<                  z-,y-
  integer,parameter :: q0N1 = 8     !<                  z+,y-
  integer,parameter :: q01N = 9     !<                  z-,y+
  integer,parameter :: q011 = 10    !<                  z+,y+
  integer,parameter :: qN0N = 11    !<                  x-,z-
  integer,parameter :: q10N = 12    !<                  x+,z-
  integer,parameter :: qN01 = 13    !<                  x-,z+
  integer,parameter :: q101 = 14    !<                  x+,z+
  integer,parameter :: qNN0 = 15    !<                  y-,x-
  integer,parameter :: qN10 = 16    !<                  y+,x-
  integer,parameter :: q1N0 = 17    !<                  y-,x+
  integer,parameter :: q110 = 18    !<                  y+,x+
  integer,parameter :: q000 = 19    !< rest density is last

  ! D3Q19 MRT pdf -> moment transformation matrix
  ! How to use:
  ! do iDir = 1, QQ
  !   moment(iDir) = sum( PDF(:) * MMtrD3Q19(iDir,:) )
  ! end do
  !  W      S     B     E     N     T     BS    TS   BN    TN    BW    BE    TW    TE    SW    NW    SE    NE     0
  real(kind=rk), dimension(19,19),parameter,public  :: MMtrD3Q19 = &
  reshape((/ &
   1.0,   1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0, &
   0.0,   0.0,  0.0,  0.0,  0.0,  0.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0, -1.0, &
  -2.0,  -2.0, -2.0, -2.0, -2.0, -2.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0, &
  -1.0,   0.0,  0.0,  1.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, -1.0,  1.0, -1.0,  1.0, -1.0, -1.0,  1.0,  1.0,  0.0, &
   2.0,   0.0,  0.0, -2.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, -1.0,  1.0, -1.0,  1.0, -1.0, -1.0,  1.0,  1.0,  0.0, &
   0.0,  -1.0,  0.0,  0.0,  1.0,  0.0, -1.0, -1.0,  1.0,  1.0,  0.0,  0.0,  0.0,  0.0, -1.0,  1.0, -1.0,  1.0,  0.0, &
   0.0,   2.0,  0.0,  0.0, -2.0,  0.0, -1.0, -1.0,  1.0,  1.0,  0.0,  0.0,  0.0,  0.0, -1.0,  1.0, -1.0,  1.0,  0.0, &
   0.0,   0.0, -1.0,  0.0,  0.0,  1.0, -1.0,  1.0, -1.0,  1.0, -1.0, -1.0,  1.0,  1.0,  0.0,  0.0,  0.0,  0.0,  0.0, &
   0.0,   0.0,  2.0,  0.0,  0.0, -2.0, -1.0,  1.0, -1.0,  1.0, -1.0, -1.0,  1.0,  1.0,  0.0,  0.0,  0.0,  0.0,  0.0, &
   2.0,  -1.0, -1.0,  2.0, -1.0, -1.0, -2.0, -2.0, -2.0, -2.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  0.0, &
  -2.0,   1.0,  1.0, -2.0,  1.0,  1.0, -2.0, -2.0, -2.0, -2.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  0.0, &
   0.0,   1.0, -1.0,  0.0,  1.0, -1.0,  0.0,  0.0,  0.0,  0.0, -1.0, -1.0, -1.0, -1.0,  1.0,  1.0,  1.0,  1.0,  0.0, &
   0.0,  -1.0,  1.0,  0.0, -1.0,  1.0,  0.0,  0.0,  0.0,  0.0, -1.0, -1.0, -1.0, -1.0,  1.0,  1.0,  1.0,  1.0,  0.0, &
   0.0,   0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  1.0, -1.0, -1.0,  1.0,  0.0, &
   0.0,   0.0,  0.0,  0.0,  0.0,  0.0,  1.0, -1.0, -1.0,  1.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, &
   0.0,   0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  1.0, -1.0, -1.0,  1.0,  0.0,  0.0,  0.0,  0.0,  0.0, &
   0.0,   0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  1.0, -1.0,  1.0, -1.0, -1.0, -1.0,  1.0,  1.0,  0.0, &
   0.0,   0.0,  0.0,  0.0,  0.0,  0.0, -1.0, -1.0,  1.0,  1.0,  0.0,  0.0,  0.0,  0.0,  1.0, -1.0,  1.0, -1.0,  0.0, &
   0.0,   0.0,  0.0,  0.0,  0.0,  0.0,  1.0, -1.0,  1.0, -1.0, -1.0, -1.0,  1.0,  1.0,  0.0,  0.0,  0.0,  0.0,  0.0  &
  /),(/19,19/), order=(/ 2,1 /) )

!   real(kind=rk), dimension(19,19),parameter,public  :: MMivD3Q19 = &
!   reshape((/ &
!   div1_18,        0,  -div1_18,   -div1_6,    div1_6,         0,         0,         0,         0,   div1_12,  -div1_12, &
!         0,         0,       0,      0,      0,      0,      0,      0, &
!   div1_18,        0,  -div1_18,         0,         0,   -div1_6,    div1_6,         0,         0,  -div1_24,   div1_24, &
!    div1_8,   -div1_8,       0,      0,      0,      0,      0,      0, &
!   div1_18,        0,  -div1_18,         0,         0,         0,         0,   -div1_6,    div1_6,  -div1_24,   div1_24, &
!   -div1_8,    div1_8,       0,      0,      0,      0,      0,      0, &
!   div1_18,        0,  -div1_18,    div1_6,   -div1_6,         0,         0,         0,         0,   div1_12,  -div1_12, &
!         0,         0,       0,      0,      0,      0,      0,      0, &
!   div1_18,        0,  -div1_18,         0,         0,    div1_6,   -div1_6,         0,         0,  -div1_24,   div1_24, &
!    div1_8,   -div1_8,       0,      0,      0,      0,      0,      0, &
!   div1_18,        0,  -div1_18,         0,         0,         0,         0,    div1_6,   -div1_6,  -div1_24,   div1_24, &
!   -div1_8,    div1_8,       0,      0,      0,      0,      0,      0, &
!   div1_36,  div1_24,   div1_72,         0,         0,  -div1_12,  -div1_24,  -div1_12,  -div1_24,  -div1_24,  -div1_24, &
!         0,         0,       0,    div1_4,      0,      0,   -div1_8,    div1_8, &
!   div1_36,  div1_24,   div1_72,         0,         0,  -div1_12,  -div1_24,   div1_12,   div1_24,  -div1_24,  -div1_24, &
!         0,         0,       0,   -div1_4,      0,      0,   -div1_8,   -div1_8, &
!   div1_36,  div1_24,   div1_72,         0,         0,   div1_12,   div1_24,  -div1_12,  -div1_24,  -div1_24,  -div1_24, &
!         0,         0,       0,   -div1_4,      0,      0,    div1_8,    div1_8, &
!   div1_36,  div1_24,   div1_72,         0,         0,   div1_12,   div1_24,   div1_12,   div1_24,  -div1_24,  -div1_24, &
!         0,         0,       0,    div1_4,      0,      0,    div1_8,   -div1_8, &
!   div1_36,  div1_24,   div1_72,  -div1_12,  -div1_24,         0,         0,  -div1_12,  -div1_24,   div1_48,   div1_48, &
!  -div1_16,  -div1_16,       0,      0,    div1_4,    div1_8,      0,   -div1_8, &
!   div1_36,  div1_24,   div1_72,   div1_12,   div1_24,         0,         0,  -div1_12,  -div1_24,   div1_48,   div1_48, &
!  -div1_16,  -div1_16,       0,      0,   -div1_4,   -div1_8,      0,   -div1_8, &
!   div1_36,  div1_24,   div1_72,  -div1_12,  -div1_24,         0,         0,   div1_12,   div1_24,   div1_48,   div1_48, &
!  -div1_16,  -div1_16,       0,      0,   -div1_4,    div1_8,      0,    div1_8, &
!   div1_36,  div1_24,   div1_72,   div1_12,   div1_24,         0,         0,   div1_12,   div1_24,   div1_48,   div1_48, &
!  -div1_16,  -div1_16,       0,      0,    div1_4,   -div1_8,      0,    div1_8, &
!   div1_36,  div1_24,   div1_72,  -div1_12,  -div1_24,  -div1_12,  -div1_24,         0,         0,   div1_48,   div1_48, &
!   div1_16,   div1_16,  div1_4,      0,      0,   -div1_8,    div1_8,      0, &
!   div1_36,  div1_24,   div1_72,  -div1_12,  -div1_24,   div1_12,   div1_24,         0,         0,   div1_48,   div1_48, &
!   div1_16,   div1_16, -div1_4,      0,      0,   -div1_8,   -div1_8,      0, &
!   div1_36,  div1_24,   div1_72,   div1_12,   div1_24,  -div1_12,  -div1_24,         0,         0,   div1_48,   div1_48, &
!   div1_16,   div1_16, -div1_4,      0,      0,    div1_8,    div1_8,      0, &
!   div1_36,  div1_24,   div1_72,   div1_12,   div1_24,   div1_12,   div1_24,         0,         0,   div1_48,   div1_48, &
!   div1_16,   div1_16,  div1_4,      0,      0,    div1_8,   -div1_8,      0, &
!    div1_3,  -div1_2,    div1_6,         0,         0,         0,         0,         0,         0,         0,         0, &
!         0,         0,       0,      0,      0,      0,      0,      0 &
!   /),(/19,19/))

  ! D3Q19 MRT moment --> PDF transformation matrix
  ! How to use:
  ! do iDir = 1, QQ
  !   fneq(iDir) = sum( MMIvD3Q19(iDir,:) * mneq(:) )
  ! end do
  real(kind=rk), dimension(19,19),parameter,public  :: MMivD3Q19 = &
  reshape((/ &
   1/18.0,     0.0, -1/18.0,  -1/6.0,   1/6.0,     0.0,      0.0,     0.0,     0.0,  1/12.0, -1/12.0,     0.0,     0.0,    0.0,&
   0.0,    0.0,    0.0,    0.0,    0.0,&
   1/18.0,     0.0, -1/18.0,     0.0,     0.0,  -1/6.0,    1/6.0,     0.0,     0.0, -1/24.0,  1/24.0,   1/8.0,  -1/8.0,    0.0,&
   0.0,    0.0,    0.0,    0.0,    0.0,&
   1/18.0,     0.0, -1/18.0,     0.0,     0.0,     0.0,      0.0,  -1/6.0,   1/6.0, -1/24.0,  1/24.0,  -1/8.0,   1/8.0,    0.0,&
   0.0,    0.0,    0.0,    0.0,    0.0,&
   1/18.0,     0.0, -1/18.0,   1/6.0,  -1/6.0,     0.0,      0.0,     0.0,     0.0,  1/12.0, -1/12.0,     0.0,     0.0,    0.0,&
   0.0,    0.0,    0.0,    0.0,    0.0,&
   1/18.0,     0.0, -1/18.0,     0.0,     0.0,   1/6.0,   -1/6.0,     0.0,     0.0, -1/24.0,  1/24.0,   1/8.0,  -1/8.0,    0.0,&
   0.0,    0.0,    0.0,    0.0,    0.0,&
   1/18.0,     0.0, -1/18.0,     0.0,     0.0,     0.0,      0.0,   1/6.0,  -1/6.0, -1/24.0,  1/24.0,  -1/8.0,   1/8.0,    0.0,&
   0.0,    0.0,    0.0,    0.0,    0.0,&
   1/36.0,  1/24.0,  1/72.0,     0.0,     0.0, -1/12.0,  -1/24.0, -1/12.0, -1/24.0, -1/24.0, -1/24.0,     0.0,     0.0,    0.0,&
 1/4.0,    0.0,    0.0, -1/8.0,  1/8.0,&
   1/36.0,  1/24.0,  1/72.0,     0.0,     0.0, -1/12.0,  -1/24.0,  1/12.0,  1/24.0, -1/24.0, -1/24.0,     0.0,     0.0,    0.0,&
-1/4.0,    0.0,    0.0, -1/8.0, -1/8.0,&
   1/36.0,  1/24.0,  1/72.0,     0.0,     0.0,  1/12.0,   1/24.0, -1/12.0, -1/24.0, -1/24.0, -1/24.0,     0.0,     0.0,    0.0,&
-1/4.0,    0.0,    0.0,  1/8.0,  1/8.0,&
   1/36.0,  1/24.0,  1/72.0,     0.0,     0.0,  1/12.0,   1/24.0,  1/12.0,  1/24.0, -1/24.0, -1/24.0,     0.0,     0.0,    0.0,&
 1/4.0,    0.0,    0.0,  1/8.0, -1/8.0,&
   1/36.0,  1/24.0,  1/72.0, -1/12.0, -1/24.0,     0.0,      0.0, -1/12.0, -1/24.0,  1/48.0,  1/48.0, -1/16.0, -1/16.0,    0.0,&
   0.0,  1/4.0,  1/8.0,    0.0, -1/8.0,&
   1/36.0,  1/24.0,  1/72.0,  1/12.0,  1/24.0,     0.0,      0.0, -1/12.0, -1/24.0,  1/48.0,  1/48.0, -1/16.0, -1/16.0,    0.0,&
   0.0, -1/4.0, -1/8.0,    0.0, -1/8.0,&
   1/36.0,  1/24.0,  1/72.0, -1/12.0, -1/24.0,     0.0,      0.0,  1/12.0,  1/24.0,  1/48.0,  1/48.0, -1/16.0, -1/16.0,    0.0,&
   0.0, -1/4.0,  1/8.0,    0.0,  1/8.0,&
   1/36.0,  1/24.0,  1/72.0,  1/12.0,  1/24.0,     0.0,      0.0,  1/12.0,  1/24.0,  1/48.0,  1/48.0, -1/16.0, -1/16.0,    0.0,&
   0.0,  1/4.0, -1/8.0,    0.0,  1/8.0,&
   1/36.0,  1/24.0,  1/72.0, -1/12.0, -1/24.0, -1/12.0,  -1/24.0,     0.0,     0.0,  1/48.0,  1/48.0,  1/16.0,  1/16.0,  1/4.0,&
   0.0,    0.0, -1/8.0,  1/8.0,    0.0,&
   1/36.0,  1/24.0,  1/72.0, -1/12.0, -1/24.0,  1/12.0,   1/24.0,     0.0,     0.0,  1/48.0,  1/48.0,  1/16.0,  1/16.0, -1/4.0,&
   0.0,    0.0, -1/8.0, -1/8.0,    0.0,&
   1/36.0,  1/24.0,  1/72.0,  1/12.0,  1/24.0, -1/12.0,  -1/24.0,     0.0,     0.0,  1/48.0,  1/48.0,  1/16.0,  1/16.0, -1/4.0,&
   0.0,    0.0,  1/8.0,  1/8.0,    0.0,&
   1/36.0,  1/24.0,  1/72.0,  1/12.0,  1/24.0,  1/12.0,   1/24.0,     0.0,     0.0,  1/48.0,  1/48.0,  1/16.0,  1/16.0,  1/4.0,&
   0.0,    0.0,  1/8.0, -1/8.0,    0.0,&
    1/3.0,  -1/2.0,   1/6.0,     0.0,     0.0,     0.0,      0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.0,    0.0,&
   0.0,    0.0,    0.0,    0.0,    0.0 &
  /),(/19,19/), order=(/ 2,1 /) )

contains

! ****************************************************************************** !
  !> Advection relaxation routine for the MRT model.
  !! This routine has roughly 260 FLOPS per elements.
  !!
  !! This subroutine interface must match the abstract interface definition
  !! [[kernel]] in scheme/[[mus_scheme_type_module]].f90 in order to be callable
  !! via [[mus_scheme_type:compute]] function pointer.
  subroutine mrt_advRel_d3q19( fieldProp, inState, outState, auxField, &
    &                          neigh, nElems, nSolve, level, layout,   &
    &                          params, varSys, derVarPos               )
    ! -------------------------------------------------------------------- !
    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), intent(in) :: fieldProp(:)
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> current layout
    type(mus_scheme_layout_type), intent(in) :: layout
    !> number of elements in state Array
    integer, intent(in) :: nElems
    !> input  pdf vector
    real(kind=rk), intent(in)  ::  inState(nElems * varSys%nScalars)
    !> output pdf vector
    real(kind=rk), intent(out) :: outState(nElems * varSys%nScalars)
    !> Auxiliary field computed from pre-collision state
    !! Is updated with correct velocity field for multicomponent models
    real(kind=rk), intent(inout) :: auxField(nElems * varSys%nAuxScalars)
    !> connectivity vector
    integer, intent(in) :: neigh(nElems * layout%fStencil%QQ)
    !> number of elements solved in kernel
    integer, intent(in) :: nSolve
    !> current level
    integer,intent(in) :: level
    !> global parameters
    type(mus_param_type),intent(in) :: params
    !> position of derived quantities in varsys for all fields
    type( mus_derVarPos_type ), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer       :: iElem
?? copy :: pdfTmp19( f )
    real(kind=rk) :: rho     ! local density
    real(kind=rk) :: u_x     ! local x-velocity
    real(kind=rk) :: u_y     ! local y-velocity
    real(kind=rk) :: u_z     ! local z-velocity
    ! MRT Variables
    real(kind=rk) :: meq2, meq10, meq12
    real(kind=rk) :: m2, m6, m8, m14, m15, m16
    ! m6, m8 are temporary var
    ! mout1 are temp var
?? copy :: var19( mout )
?? copy :: var19( sum )
    real(kind=rk) :: c0, c1, c2, c3, c4, c5, c6
    real(kind=rk) :: mout5_4, mout7_4, mout9_4
    real(kind=rk) :: sum_c1_c2, sub_c1_c2
    real(kind=rk) :: sum_5_17, sub_7_18, d1,d2,d3,d4
    real(kind=rk) :: sum_9_19, sub_5_17, e1,e2,e3,e4
    real(kind=rk) :: sum_7_18, sub_9_19, g1,g2,g3,g4
    integer :: dens_pos, vel_pos(3), elemOff
    real(kind=rk) :: omegaKine, omegaBulk, s_mrt(QQ)
    ! ---------------------------------------------------------------------------
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    omegaBulk = fieldProp(1)%fluid%omegaBulkLvl(level)
    ! MRT omegas
    ! overwrite omegaKine term in the element loop
    s_mrt = mrt_d3q19(omegaKine=1.0_rk, omegaBulk=omegaBulk, QQ=QQ)

    s_mrt(2)  = s_mrt( 2) * div1_24
    s_mrt(3)  = s_mrt( 3) * div1_72
    s_mrt(5)  = s_mrt( 5) * div1_24
    s_mrt(7)  = s_mrt( 7) * div1_24
    s_mrt(9)  = s_mrt( 9) * div1_24
    s_mrt(17) = s_mrt(17) * div1_8
    s_mrt(18) = s_mrt(18) * div1_8
    s_mrt(19) = s_mrt(19) * div1_8

!$omp do schedule(static)
    !NEC$ ivdep
?? IF (SOA) THEN
?? copy :: dir_vector
?? ELSE
?? copy :: dir_novec
?? END IF
    nodeloop: do iElem = 1, nSolve

      !> First load all local values into temp array
      fN00 = inState(?FETCH?( qN00, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f0N0 = inState(?FETCH?( q0N0, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f00N = inState(?FETCH?( q00N, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f100 = inState(?FETCH?( q100, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f010 = inState(?FETCH?( q010, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f001 = inState(?FETCH?( q001, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f0NN = inState(?FETCH?( q0NN, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f0N1 = inState(?FETCH?( q0N1, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f01N = inState(?FETCH?( q01N, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f011 = inState(?FETCH?( q011, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      fN0N = inState(?FETCH?( qN0N, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f10N = inState(?FETCH?( q10N, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      fN01 = inState(?FETCH?( qN01, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f101 = inState(?FETCH?( q101, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      fNN0 = inState(?FETCH?( qNN0, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      fN10 = inState(?FETCH?( qN10, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f1N0 = inState(?FETCH?( q1N0, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f110 = inState(?FETCH?( q110, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f000 = inState(?FETCH?( q000, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))

      m6     = f101 + fN0N + f10N + fN01
      m8     = f011 + f0NN + f01N + f0N1
      sum1   = f110 + fNN0 + f1N0 + fN10
      m2     =-f000 + sum1 + m6 + m8

      sum2 = f010 + f0N0
      sum3 = f001 + f00N
      sum4 = 2._rk * ( f100+ fN00 )
      sum5 = sum2 + sum3

      ! epsilon
      mout3 = ( 2._rk*(f000 - sum5) - sum4 + m2 ) * s_mrt(3)

      ! element offset for auxField array
      elemOff = (iElem-1)*varSys%nAuxScalars
      ! local density
      rho = auxField(elemOff + dens_pos)

      ! local x-, y- and z-velocity
      u_x = auxField(elemOff + vel_pos(1))
      u_y = auxField(elemOff + vel_pos(2))
      u_z = auxField(elemOff + vel_pos(3))

      omegaKine = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem)
      s_mrt(10) = omegaKine
      s_mrt(12) = omegaKine
      s_mrt(14) = div1_4 * omegaKine
      s_mrt(15) = div1_4 * omegaKine
      s_mrt(16) = div1_4 * omegaKine

      ! Equilibrium moments
      ! non zero meq are only: 2, 10, 12, 14, 15, 16
      meq2  = rho * ( u_x*u_x + u_y*u_y + u_z*u_z )
      meq10 = rho * 3.0_rk*u_x*u_x - meq2
      meq12 = rho * ( u_y*u_y - u_z*u_z )

      mout2  = s_mrt(2) * (m2 - meq2)

      !pxy
      m14 = f110 + fNN0 - f1N0 - fN10
      mout14 = s_mrt(14) * (m14 - rho*u_x*u_y)
      !pyz
      m15 = f011 + f0NN - f01N - f0N1
      mout15 = s_mrt(15) * (m15 - rho*u_y*u_z)
      !pxz
      m16 = f101 + fN0N - f10N - fN01
      mout16 = s_mrt(16) * (m16 - rho*u_x*u_z)

      sum6 = sum1 + m6 - m8 * 2.0_rk
      sum7 = sum4 - sum5
      !3pxx
      mout10 = ( sum7 + sum6 - meq10 ) * s_mrt(10)
      !3pixx
      mout11 = ( - sum7 + sum6 ) * s_mrt(11)

      sum8 = sum1 - m6
      sum9 = sum2 - sum3
      !pww
      mout12 = ( sum8 + sum9 - meq12 ) * s_mrt(12)
      !piww
      mout13 = ( sum8 - sum9 ) * s_mrt(13)

c1 = f110 - fNN0
c2 = f1N0 - fN10
c3 = f101 - fN0N
c4 = f10N - fN01

      sum10 = c1 + c2
      sum11 = c3 + c4
      !qx
      mout5  = ( sum10 + sum11 - 2.0_rk*(f100 - fN00) ) * s_mrt(5)
      !mx
      mout17 = ( sum10 - sum11 ) * s_mrt(17)

c5 = f011 - f0NN
c6 = f01N - f0N1

      sum12 = c1 - c2
      sum13 = c5 + c6
      !qy
      mout7  = ( sum12 + sum13 - 2.0_rk*(f010 - f0N0) ) * s_mrt(7)
      !my
      mout18 = ( - sum12 + sum13 ) * s_mrt(18)

      sum14 = c3 - c4
      sum15 = c5 - c6
      !qz
      mout9  = ( sum14 + sum15 - 2.0_rk*(f001 - f00N) ) * s_mrt(9)
      !mz
      mout19 = ( sum14 - sum15 ) * s_mrt(19)

      ! Transformation back to PDF
      outstate( ?SAVE?( 19,1,iElem,QQ,QQ,nElems,neigh )) = f000 + 12._rk*(mout2-mout3)

! -------------------------------------------------------------------------------
      c0 = - 4._rk*mout3 + div1_12*(mout10 - mout11)
      mout5_4 = mout5 * 4.0_rk
      outstate( ?SAVE?(  4,1,iElem,QQ,QQ,nElems,neigh )) = f100 - ( c0 - mout5_4 )
      outstate( ?SAVE?(  1,1,iElem,QQ,QQ,nElems,neigh )) = fN00 - ( c0 + mout5_4 )
! -------------------------------------------------------------------------------

! -------------------------------------------------------------------------------
      c1 = - 4._rk*mout3 - div1_24*(mout10 - mout11)
      c2 = div1_8 *(mout12 - mout13)

      sum_c1_c2 = c1 + c2
      mout7_4 = mout7 * 4.0_rk
      outstate( ?SAVE?(  5,1,iElem,QQ,QQ,nElems,neigh)) = f010 - (sum_c1_c2 - mout7_4)
      outstate( ?SAVE?(  2,1,iElem,QQ,QQ,nElems,neigh)) = f0N0 - (sum_c1_c2 + mout7_4)

      sub_c1_c2 = c1 - c2
      mout9_4 = mout9 * 4.0_rk
      outstate( ?SAVE?(  6,1,iElem,QQ,QQ,nElems,neigh)) = f001 - (sub_c1_c2 - mout9_4)
      outstate( ?SAVE?(  3,1,iElem,QQ,QQ,nElems,neigh)) = f00N - (sub_c1_c2 + mout9_4)
! -------------------------------------------------------------------------------

! -------------------------------------------------------------------------------
      mout1 = mout2 + mout3
      c3 =   mout1  + div1_48 * (mout10 + mout11)   &
       &            + div1_16 * (mout12 + mout13)
      sum_5_17 = mout5 + mout17
      sub_7_18 = mout7 - mout18

      d1 = c3 + mout14
      d2 = sum_5_17 + sub_7_18
      outstate( ?SAVE?( 18,1,iElem,QQ,QQ,nElems,neigh)) = f110-(d1+d2)
      outstate( ?SAVE?( 15,1,iElem,QQ,QQ,nElems,neigh)) = fNN0-(d1-d2)

      d3 = c3 - mout14
      d4 = sum_5_17 - sub_7_18
      outstate( ?SAVE?( 17,1,iElem,QQ,QQ,nElems,neigh)) = f1N0-(d3+d4)
      outstate( ?SAVE?( 16,1,iElem,QQ,QQ,nElems,neigh)) = fN10-(d3-d4)
! -------------------------------------------------------------------------------

! -------------------------------------------------------------------------------
      c4 = c3 - div1_8*(mout12+mout13)

      sum_9_19 = mout9 + mout19
      sub_5_17 = mout5 - mout17

      e1 = c4 + mout16
      e2 = sum_9_19 + sub_5_17
      outstate( ?SAVE?( 14,1,iElem,QQ,QQ,nElems,neigh )) = f101 - ( e1 + e2 )
      outstate( ?SAVE?( 11,1,iElem,QQ,QQ,nElems,neigh )) = fN0N - ( e1 - e2 )

      e3 = c4 - mout16
      e4 = sum_9_19 - sub_5_17
      outstate( ?SAVE?( 12,1,iElem,QQ,QQ,nElems,neigh )) = f10N - ( e3 - e4 )
      outstate( ?SAVE?( 13,1,iElem,QQ,QQ,nElems,neigh )) = fN01 - ( e3 + e4 )
! -------------------------------------------------------------------------------

! -------------------------------------------------------------------------------
      c5 = mout1 - div1_24*(mout10+mout11)
      sum_7_18 = mout7 + mout18
      sub_9_19 = mout9 - mout19

      g1 = c5 + mout15
      g2 = sum_7_18 + sub_9_19
      outstate( ?SAVE?( 10,1,iElem,QQ,QQ,nElems,neigh )) = f011 - ( g1 + g2 )
      outstate( ?SAVE?(  7,1,iElem,QQ,QQ,nElems,neigh )) = f0NN - ( g1 - g2 )

      g3 = c5 - mout15
      g4 = sum_7_18 - sub_9_19
      outstate( ?SAVE?(  9,1,iElem,QQ,QQ,nElems,neigh )) = f01N - ( g3 + g4 )
      outstate( ?SAVE?(  8,1,iElem,QQ,QQ,nElems,neigh )) = f0N1 - ( g3 - g4 )
! -------------------------------------------------------------------------------

    enddo nodeloop
!$omp end do

  end subroutine mrt_advRel_d3q19
! ****************************************************************************** !


! ****************************************************************************** !
  !> Advection relaxation routine for the MRT model.
  !! This routine has roughly 205 FLOPS per element.
  !!
  !! This subroutine interface must match the abstract interface definition
  !! [[kernel]] in scheme/[[mus_scheme_type_module]].f90 in order to be callable
  !! via [[mus_scheme_type:compute]] function pointer.
  subroutine mrt_advRel_d3q19_incomp( fieldProp, inState, outState, auxField, &
    &                                 neigh, nElems, nSolve, level, layout,   &
    &                                 params, varSys, derVarPos               )
    ! -------------------------------------------------------------------- !
    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), intent(in) :: fieldProp(:)
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> current layout
    type(mus_scheme_layout_type), intent(in) :: layout
    !> number of elements in state Array
    integer, intent(in) :: nElems
    !> input  pdf vector
    real(kind=rk), intent(in)  ::  inState(nElems * varSys%nScalars)
    !> output pdf vector
    real(kind=rk), intent(out) :: outState(nElems * varSys%nScalars)
    !> Auxiliary field computed from pre-collision state
    !! Is updated with correct velocity field for multicomponent models
    real(kind=rk), intent(inout) :: auxField(nElems * varSys%nAuxScalars)
    !> connectivity vector
    integer, intent(in) :: neigh(nElems * layout%fStencil%QQ)
    !> number of elements solved in kernel
    integer, intent(in) :: nSolve
    !> current level
    integer,intent(in) :: level
    !> global parameters
    type(mus_param_type),intent(in) :: params
    !> position of derived quantities in varsys for all fields
    type( mus_derVarPos_type ), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer       :: iElem
?? copy :: pdfTmp19( f )
    real(kind=rk) :: u_x     ! local x-velocity
    real(kind=rk) :: u_y     ! local y-velocity
    real(kind=rk) :: u_z     ! local z-velocity
    ! MRT Variables
    real(kind=rk) :: meq2, meq10, meq12
    real(kind=rk) :: m2, m6, m8, m14, m15, m16
    ! m6, m8 are temporary var
    ! mout1 are temp var
?? copy :: var19( mout )
?? copy :: var19( sum )
    real(kind=rk) :: c0, c1, c2, c3, c4, c5, c6
    real(kind=rk) :: mout5_4, mout7_4, mout9_4
    real(kind=rk) :: sum_c1_c2, sub_c1_c2
    real(kind=rk) :: sum_5_17, sub_7_18, d1,d2,d3,d4
    real(kind=rk) :: sum_9_19, sub_5_17, e1,e2,e3,e4
    real(kind=rk) :: sum_7_18, sub_9_19, g1,g2,g3,g4
    integer :: dens_pos, vel_pos(3), elemOff
    real(kind=rk) :: omegaKine, omegaBulk, s_mrt(QQ)
    ! ---------------------------------------------------------------------------
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    omegaBulk = fieldProp(1)%fluid%omegaBulkLvl(level)
    ! MRT omegas
    ! overwrite omegaKine term in the element loop
    s_mrt = mrt_d3q19(omegaKine=1.0_rk, omegaBulk=omegaBulk, QQ=QQ)

    s_mrt(2)  = s_mrt( 2) * div1_24
    s_mrt(3)  = s_mrt( 3) * div1_72
    s_mrt(5)  = s_mrt( 5) * div1_24
    s_mrt(7)  = s_mrt( 7) * div1_24
    s_mrt(9)  = s_mrt( 9) * div1_24
    s_mrt(17) = s_mrt(17) * div1_8
    s_mrt(18) = s_mrt(18) * div1_8
    s_mrt(19) = s_mrt(19) * div1_8

    !NEC$ ivdep
?? IF (SOA) THEN
?? copy :: dir_vector
?? ELSE
?? copy :: dir_novec
?? END IF
    nodeloop: do iElem = 1, nSolve

      !> First load all local values into temp array
      fN00 = inState(?FETCH?( qN00, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f0N0 = inState(?FETCH?( q0N0, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f00N = inState(?FETCH?( q00N, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f100 = inState(?FETCH?( q100, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f010 = inState(?FETCH?( q010, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f001 = inState(?FETCH?( q001, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f0NN = inState(?FETCH?( q0NN, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f0N1 = inState(?FETCH?( q0N1, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f01N = inState(?FETCH?( q01N, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f011 = inState(?FETCH?( q011, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      fN0N = inState(?FETCH?( qN0N, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f10N = inState(?FETCH?( q10N, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      fN01 = inState(?FETCH?( qN01, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f101 = inState(?FETCH?( q101, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      fNN0 = inState(?FETCH?( qNN0, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      fN10 = inState(?FETCH?( qN10, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f1N0 = inState(?FETCH?( q1N0, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f110 = inState(?FETCH?( q110, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))
      f000 = inState(?FETCH?( q000, 1, iElem, QQ, varSys%nScalars, nElems,neigh ))

      ! element offset for auxField array
      elemOff = (iElem-1)*varSys%nAuxScalars

      ! local x-, y- and z-velocity
      u_x = auxField(elemOff + vel_pos(1))
      u_y = auxField(elemOff + vel_pos(2))
      u_z = auxField(elemOff + vel_pos(3))

      ! Equilibrium moments
      ! non zero meq are only: 2, 10, 12, 14, 15, 16
      meq2  = u_x*u_x + u_y*u_y + u_z*u_z
      meq10 = 3.0_rk*u_x*u_x - meq2
      meq12 = u_y*u_y - u_z*u_z

      omegaKine = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem)
      s_mrt(10) = omegaKine
      s_mrt(12) = omegaKine
      s_mrt(14) = div1_4 * omegaKine
      s_mrt(15) = div1_4 * omegaKine
      s_mrt(16) = div1_4 * omegaKine

      !pxy
      m14 = f110 + fNN0 - f1N0 - fN10
      mout14 = s_mrt(14) * (m14 - u_x*u_y)
      !pyz
      m15 = f011 + f0NN - f01N - f0N1
      mout15 = s_mrt(15) * (m15 - u_y*u_z)
      !pxz
      m16 = f101 + fN0N - f10N - fN01
      mout16 = s_mrt(16) * (m16 - u_x*u_z)

      ! -------------------------------------------------------------------------
      ! now calculate mrt against the bgk
      m6 = f101 + fN0N + f10N + fN01
      m8 = f011 + f0NN + f01N + f0N1

      sum1   =   f110 + fNN0 + f1N0 + fN10
      m2     = - f000 + sum1 + m6 + m8
      mout2  =   s_mrt(2)  * (m2 - meq2)

      sum2 = f010 + f0N0
      sum3 = f001 + f00N
      sum4 = 2._rk * ( f100+ fN00 )
      sum5 = sum2 + sum3

      ! epsilon
      mout3 = ( 2._rk*(f000 - sum5) - sum4 + m2 ) * s_mrt(3)

      sum6 = sum1 + m6 - m8 * 2.0_rk
      sum7 = sum4 - sum5
      !3pxx
      mout10 = ( sum7 + sum6 - meq10 ) * s_mrt(10)
      !3pixx
      mout11 = ( - sum7 + sum6 ) * s_mrt(11)

      sum8 = sum1 - m6
      sum9 = sum2 - sum3
      !pww
      mout12 = ( sum8 + sum9 - meq12 ) * s_mrt(12)
      !piww
      mout13 = ( sum8 - sum9 ) * s_mrt(13)

c1 = f110 - fNN0
c2 = f1N0 - fN10
c3 = f101 - fN0N
c4 = f10N - fN01

      sum10 = c1 + c2
      sum11 = c3 + c4
      !qx
      mout5  = ( sum10 + sum11 - 2.0_rk*(f100 - fN00) ) * s_mrt(5)
      !mx
      mout17 = ( sum10 - sum11 ) * s_mrt(17)

c5 = f011 - f0NN
c6 = f01N - f0N1

      sum12 = c1 - c2
      sum13 = c5 + c6
      !qy
      mout7  = ( sum12 + sum13 - 2.0_rk*(f010 - f0N0) ) * s_mrt(7)
      !my
      mout18 = ( - sum12 + sum13 ) * s_mrt(18)

      sum14 = c3 - c4
      sum15 = c5 - c6
      !qz
      mout9  = ( sum14 + sum15 - 2.0_rk*(f001 - f00N) ) * s_mrt(9)
      !mz
      mout19 = ( sum14 - sum15 ) * s_mrt(19)

      ! Transformation back to PDF
      outstate( ?SAVE?( 19,1,iElem,QQ,QQ,nElems,neigh )) = f000 + 12._rk*(mout2-mout3)

! -------------------------------------------------------------------------------
      c0 = - 4._rk*mout3 + div1_12*(mout10 - mout11)
      mout5_4 = mout5 * 4.0_rk
      outstate( ?SAVE?(  4,1,iElem,QQ,QQ,nElems,neigh )) = f100 - ( c0 - mout5_4 )
      outstate( ?SAVE?(  1,1,iElem,QQ,QQ,nElems,neigh )) = fN00 - ( c0 + mout5_4 )
! -------------------------------------------------------------------------------

! -------------------------------------------------------------------------------
      c1 = - 4._rk*mout3 - div1_24*(mout10 - mout11)
      c2 = div1_8 *(mout12 - mout13)

      sum_c1_c2 = c1 + c2
      mout7_4 = mout7 * 4.0_rk
      outstate( ?SAVE?(  5,1,iElem,QQ,QQ,nElems,neigh)) = f010 - (sum_c1_c2 - mout7_4)
      outstate( ?SAVE?(  2,1,iElem,QQ,QQ,nElems,neigh)) = f0N0 - (sum_c1_c2 + mout7_4)

      sub_c1_c2 = c1 - c2
      mout9_4 = mout9 * 4.0_rk
      outstate( ?SAVE?(  6,1,iElem,QQ,QQ,nElems,neigh)) = f001 - (sub_c1_c2 - mout9_4)
      outstate( ?SAVE?(  3,1,iElem,QQ,QQ,nElems,neigh)) = f00N - (sub_c1_c2 + mout9_4)
! -------------------------------------------------------------------------------

! -------------------------------------------------------------------------------
      mout1 = mout2 + mout3
      c3 =   mout1  + div1_48 * (mout10 + mout11)   &
       &            + div1_16 * (mout12 + mout13)
      sum_5_17 = mout5 + mout17
      sub_7_18 = mout7 - mout18

      d1 = c3 + mout14
      d2 = sum_5_17 + sub_7_18
      outstate( ?SAVE?( 18,1,iElem,QQ,QQ,nElems,neigh)) = f110-(d1+d2)
      outstate( ?SAVE?( 15,1,iElem,QQ,QQ,nElems,neigh)) = fNN0-(d1-d2)

      d3 = c3 - mout14
      d4 = sum_5_17 - sub_7_18
      outstate( ?SAVE?( 17,1,iElem,QQ,QQ,nElems,neigh)) = f1N0-(d3+d4)
      outstate( ?SAVE?( 16,1,iElem,QQ,QQ,nElems,neigh)) = fN10-(d3-d4)
! -------------------------------------------------------------------------------

! -------------------------------------------------------------------------------
      c4 = c3 - div1_8*(mout12+mout13)

      sum_9_19 = mout9 + mout19
      sub_5_17 = mout5 - mout17

      e1 = c4 + mout16
      e2 = sum_9_19 + sub_5_17
      outstate( ?SAVE?( 14,1,iElem,QQ,QQ,nElems,neigh )) = f101 - ( e1 + e2 )
      outstate( ?SAVE?( 11,1,iElem,QQ,QQ,nElems,neigh )) = fN0N - ( e1 - e2 )

      e3 = c4 - mout16
      e4 = sum_9_19 - sub_5_17
      outstate( ?SAVE?( 12,1,iElem,QQ,QQ,nElems,neigh )) = f10N - ( e3 - e4 )
      outstate( ?SAVE?( 13,1,iElem,QQ,QQ,nElems,neigh )) = fN01 - ( e3 + e4 )
! -------------------------------------------------------------------------------

! -------------------------------------------------------------------------------
      c5 = mout1 - div1_24*(mout10+mout11)
      sum_7_18 = mout7 + mout18
      sub_9_19 = mout9 - mout19

      g1 = c5 + mout15
      g2 = sum_7_18 + sub_9_19
      outstate( ?SAVE?( 10,1,iElem,QQ,QQ,nElems,neigh )) = f011 - ( g1 + g2 )
      outstate( ?SAVE?(  7,1,iElem,QQ,QQ,nElems,neigh )) = f0NN - ( g1 - g2 )

      g3 = c5 - mout15
      g4 = sum_7_18 - sub_9_19
      outstate( ?SAVE?(  9,1,iElem,QQ,QQ,nElems,neigh )) = f01N - ( g3 + g4 )
      outstate( ?SAVE?(  8,1,iElem,QQ,QQ,nElems,neigh )) = f0N1 - ( g3 - g4 )
! -------------------------------------------------------------------------------

    enddo nodeloop

  end subroutine mrt_advRel_d3q19_incomp
! ****************************************************************************** !


! ****************************************************************************** !
  !> No comment yet!
  !!
  !! TODO add comment
  !!
  !! This subroutine interface must match the abstract interface definition
  !! [[kernel]] in scheme/[[mus_scheme_type_module]].f90 in order to be callable
  !! via [[mus_scheme_type:compute]] function pointer.
  subroutine mrt_advRel_d3q19_incomp_generic( fieldProp, inState, outState,    &
    &                                         auxField, neigh, nElems, nSolve, &
    &                                         level, layout, params, varSys,   &
    &                                         derVarPos                        )
    ! -------------------------------------------------------------------- !
    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), intent(in) :: fieldProp(:)
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> current layout
    type(mus_scheme_layout_type), intent(in) :: layout
    !> number of elements in state Array
    integer, intent(in) :: nElems
    !> input  pdf vector
    real(kind=rk), intent(in)  ::  inState(nElems * varSys%nScalars)
    !> output pdf vector
    real(kind=rk), intent(out) :: outState(nElems * varSys%nScalars)
    !> Auxiliary field computed from pre-collision state
    !! Is updated with correct velocity field for multicomponent models
    real(kind=rk), intent(inout) :: auxField(nElems * varSys%nAuxScalars)
    !> connectivity vector
    integer, intent(in) :: neigh(nElems * layout%fStencil%QQ)
    !> number of elements solved in kernel
    integer, intent(in) :: nSolve
    !> current level
    integer,intent(in) :: level
    !> global parameters
    type(mus_param_type),intent(in) :: params
    !> position of derived quantities in varsys for all fields
    type( mus_derVarPos_type ), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer       :: iElem, iDir
    real(kind=rk) :: pdfTmp( QQ ) ! temporary local pdf values
    real(kind=rk) :: rho     ! local density
    real(kind=rk) :: u_x     ! local x-velocity
    real(kind=rk) :: u_y     ! local y-velocity
    real(kind=rk) :: u_z     ! local z-velocity
    ! MRT Variables
    real(kind=rk) :: s_mrt( QQ ), meq( QQ )
    real(kind=rk) :: mneq( QQ ), mom( QQ )
    real(kind=rk) :: fneq( QQ )
    integer :: dens_pos, vel_pos(3), elemOff
    real(kind=rk) :: omegaKine, omegaBulk
    ! ---------------------------------------------------------------------------
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    omegaBulk = fieldProp(1)%fluid%omegaBulkLvl(level)
    ! MRT omegas
    ! overwrite omegaKine term in the element loop
    s_mrt = fieldProp(1)%fluid%mrtPtr(omegaKine=1.0_rk, omegaBulk=omegaBulk, &
      &                               QQ=QQ)

    !NEC$ ivdep
?? IF (SOA) THEN
?? copy :: dir_vector
?? ELSE
?? copy :: dir_novec
?? END IF
    nodeloop: do iElem = 1, nSolve

      !> First load all local values into temp array
      do iDir = 1, QQ
        pdfTmp( iDir ) = inState( ?FETCH?( iDir, 1, iElem, QQ, QQ, nElems,neigh ))
      end do

      ! element offset for auxField array
      elemOff = (iElem-1)*varSys%nAuxScalars
      ! local density
      rho = auxField(elemOff + dens_pos)

      ! local x-, y- and z-velocity
      u_x = auxField(elemOff + vel_pos(1))
      u_y = auxField(elemOff + vel_pos(2))
      u_z = auxField(elemOff + vel_pos(3))

      ! -------------------------------------------------------------------------
      ! Equilibrium moments
      meq(1:QQ) =  0.0_rk

      meq( 1) = rho
      meq( 2) = rho0 * ( u_x * u_x + u_y * u_y + u_z * u_z )
      meq( 4) = u_x*rho0
      meq( 6) = u_y*rho0
      meq( 8) = u_z*rho0
      meq(10) = rho0*(2.0_rk*u_x*u_x - ( u_y*u_y + u_z*u_z ))
      meq(12) = rho0*( u_y*u_y - u_z*u_z )
      meq(14) = rho0*u_x*u_y
      meq(15) = rho0*u_y*u_z
      meq(16) = rho0*u_x*u_z

      ! convert pdf into moment
      do iDir = 1, QQ
        mom(iDir) = sum( pdfTmp(:) * MMtrD3Q19(iDir,:) )
      end do

      ! update kinematic omega part in relaxation matrix
      omegaKine = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem)
      s_mrt(10) = omegaKine
      s_mrt(12) = omegaKine
      s_mrt(14:16) = omegaKine

      ! compute neq moment
      mneq(:) = s_mrt(:) * ( mom(:) - meq(:) )

      ! compute fNeq
      do iDir = 1, QQ
        fneq(iDir) = sum( MMIvD3Q19(iDir,:) * mneq(:) )
        outState(?SAVE?(iDir,1,iElem,QQ,QQ,nElems,neigh)) = pdfTmp(iDir) - fneq(iDir)
      end do

    enddo nodeloop

  end subroutine mrt_advRel_d3q19_incomp_generic
! ****************************************************************************** !

! ****************************************************************************** !
  !> No comment yet!
  !!
  !! TODO add comment
  !!
  !! This subroutine interface must match the abstract interface definition
  !! [[kernel]] in scheme/[[mus_scheme_type_module]].f90 in order to be callable
  !! via [[mus_scheme_type:compute]] function pointer.
  subroutine mrt_advRel_d3q19_generic( fieldProp, inState, outState, auxField, &
    &                                  neigh, nElems, nSolve, level, layout,   &
    &                                  params, varSys, derVarPos               )
    ! -------------------------------------------------------------------- !
    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), intent(in) :: fieldProp(:)
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> current layout
    type(mus_scheme_layout_type), intent(in) :: layout
    !> number of elements in state Array
    integer, intent(in) :: nElems
    !> input  pdf vector
    real(kind=rk), intent(in)  ::  inState(nElems * varSys%nScalars)
    !> output pdf vector
    real(kind=rk), intent(out) :: outState(nElems * varSys%nScalars)
    !> Auxiliary field computed from pre-collision state
    !! Is updated with correct velocity field for multicomponent models
    real(kind=rk), intent(inout) :: auxField(nElems * varSys%nAuxScalars)
    !> connectivity vector
    integer, intent(in) :: neigh(nElems * layout%fStencil%QQ)
    !> number of elements solved in kernel
    integer, intent(in) :: nSolve
    !> current level
    integer,intent(in) :: level
    !> global parameters
    type(mus_param_type),intent(in) :: params
    !> position of derived quantities in varsys for all fields
    type( mus_derVarPos_type ), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer       :: iElem, iDir
    real(kind=rk) :: pdfTmp( QQ ) ! temporary local pdf values
    real(kind=rk) :: rho     ! local density
    real(kind=rk) :: u_x     ! local x-velocity
    real(kind=rk) :: u_y     ! local y-velocity
    real(kind=rk) :: u_z     ! local z-velocity
    ! MRT Variables
    real(kind=rk) :: s_mrt( QQ ), meq( QQ )
    real(kind=rk) :: mneq( QQ ), mom( QQ )
    real(kind=rk) :: fneq( QQ )
    integer :: dens_pos, vel_pos(3), elemOff
    real(kind=rk) :: omegaKine, omegaBulk
    ! ---------------------------------------------------------------------------
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)
    omegaBulk = fieldProp(1)%fluid%omegaBulkLvl(level)
    ! MRT omegas
    ! overwrite omegaKine term in the element loop
    s_mrt = mrt_d3q19(omegaKine=1.0_rk, omegaBulk=omegaBulk, QQ=QQ)

    !NEC$ ivdep
?? IF (SOA) THEN
?? copy :: dir_vector
?? ELSE
?? copy :: dir_novec
?? END IF
    nodeloop: do iElem = 1,nSolve

      !> First load all local values into temp array
      do iDir = 1, QQ
        pdfTmp( iDir ) = inState( ?FETCH?( iDir, 1, iElem, QQ, QQ, nElems,neigh ))
      end do

      ! element offset for auxField array
      elemOff = (iElem-1)*varSys%nAuxScalars
      ! local density
      rho = auxField(elemOff + dens_pos)

      ! local x-, y- and z-velocity
      u_x = auxField(elemOff + vel_pos(1))
      u_y = auxField(elemOff + vel_pos(2))
      u_z = auxField(elemOff + vel_pos(3))

      ! -------------------------------------------------------------------------
      ! Equilibrium moments
      meq(1:QQ) =  0.0_rk

      meq( 1) = rho
      meq( 2) = rho * ( u_x * u_x + u_y * u_y + u_z * u_z )
      meq( 4) = u_x*rho
      meq( 6) = u_y*rho
      meq( 8) = u_z*rho
      meq(10) = rho*(2.0_rk*u_x*u_x - ( u_y*u_y + u_z*u_z ))
      meq(12) = rho*( u_y*u_y - u_z*u_z )
      meq(14) = rho*u_x*u_y
      meq(15) = rho*u_y*u_z
      meq(16) = rho*u_x*u_z

      ! convert pdf into moment
      do iDir = 1, QQ
        mom(iDir) = sum( pdfTmp(:) * MMtrD3Q19(iDir,:) )
      end do

      ! update kinematic omega part in relaxation matrix
      omegaKine = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem)
      s_mrt(10) = omegaKine
      s_mrt(12) = omegaKine
      s_mrt(14:16) = omegaKine

      ! compute neq moment
      mneq(:) = s_mrt(:) * ( mom(:) - meq(:) )

      ! compute fNeq
      do iDir = 1, QQ
        fneq(iDir) = sum( MMIvD3Q19(iDir,:) * mneq(:) )
        outState(?SAVE?(iDir,1,iElem,QQ,QQ,nElems,neigh)) = pdfTmp(iDir) - fneq(iDir)
      end do

    enddo nodeloop

  end subroutine mrt_advRel_d3q19_generic
! ****************************************************************************** !


! **************************************************************************** !
  !> Unoptimized explicit implementation
  !!
  !! This subroutine interface must match the abstract interface definition
  !! [[kernel]] in scheme/[[mus_scheme_type_module]].f90 in order to be callable
  !! via [[mus_scheme_type:compute]] function pointer.
  subroutine mrt_advRel_generic( fieldProp, inState, outState, auxField, &
    &                            neigh, nElems, nSolve, level, layout,   &
    &                            params, varSys, derVarPos )
    ! -------------------------------------------------------------------- !
    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), intent(in) :: fieldProp(:)
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> current layout
    type(mus_scheme_layout_type), intent(in) :: layout
    !> number of elements in state Array
    integer, intent(in) :: nElems
    !> input  pdf vector
    real(kind=rk), intent(in)  ::  inState(nElems * varSys%nScalars)
    !> output pdf vector
    real(kind=rk), intent(out) :: outState(nElems * varSys%nScalars)
    !> Auxiliary field computed from pre-collision state
    !! Is updated with correct velocity field for multicomponent models
    real(kind=rk), intent(inout) :: auxField(nElems * varSys%nAuxScalars)
    !> connectivity vector
    integer, intent(in) :: neigh(nElems * layout%fStencil%QQ)
    !> number of elements solved in kernel
    integer, intent(in) :: nSolve
    !> current level
    integer,intent(in) :: level
    !> global parameters
    type(mus_param_type),intent(in) :: params
    !> position of derived quantities in varsys for all fields
    type( mus_derVarPos_type ), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer       :: iElem, iDir
    real(kind=rk) :: pdfTmp( layout%fStencil%QQ ) ! temporary local pdf values
    real(kind=rk) :: rho     ! local density
    real(kind=rk) :: u_x     ! local x-velocity
    real(kind=rk) :: u_y     ! local y-velocity
    real(kind=rk) :: u_z     ! local z-velocity
    real(kind=rk) :: usq, ucx
    real(kind=rk) :: fEq( layout%fStencil%QQ ) !< equilibrium distribution
    ! MRT Variables
    real(kind=rk) :: fneq( layout%fStencil%QQ )
    real(kind=rk) :: mneq( layout%fStencil%QQ )
    real(kind=rk) :: s_mrt( layout%fStencil%QQ )
    real(kind=rk) :: mInvXomega(layout%fStencil%QQ, layout%fStencil%QQ)
    integer :: QQ
    integer :: dens_pos, vel_pos(3), elemOff
    real(kind=rk) :: omegaKine, omegaBulk, colli_term
    ! ---------------------------------------------------------------------------
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    QQ = layout%fStencil%QQ
    omegaBulk = fieldProp(1)%fluid%omegaBulkLvl(level)

    !NEC$ ivdep
?? IF (SOA) THEN
?? copy :: dir_vector
?? ELSE
?? copy :: dir_novec
?? END IF
    nodeloop: do iElem = 1,nSolve

      !> First load all local values into temp array
      do iDir = 1, QQ
        pdfTmp( iDir ) = inState( ?FETCH?( iDir, 1, iElem, QQ, QQ, nElems,neigh ))
      end do

      ! element offset for auxField array
      elemOff = (iElem-1)*varSys%nAuxScalars
      ! local density
      rho = auxField(elemOff + dens_pos)

      ! local x-, y- and z-velocity
      u_x = auxField(elemOff + vel_pos(1))
      u_y = auxField(elemOff + vel_pos(2))
      u_z = auxField(elemOff + vel_pos(3))

      ! velocity squared
      usq = u_x*u_x + u_y*u_y + u_z*u_z

      ! -------------------------------------------------------------------------
      ! Calculate the equilibrium distribution function
      do iDir = 1, QQ

        ! velocity times lattice unit velocity
        ucx =   layout%fStencil%cxDirRK(1,iDir) * u_x       &
          &   + layout%fStencil%cxDirRK(2,iDir) * u_y       &
          &   + layout%fStencil%cxDirRK(3,iDir) * u_z

        ! calculate equilibrium density
        fEq( iDir ) = layout%weight( iDir ) * rho * ( 1._rk + ucx*cs2inv       &
          &         + ucx*ucx*t2cs4inv - usq*t2cs2inv )

      enddo

      ! compute fNeq
      fneq = pdfTmp - fEq

      ! convert fneq to moments neq
      mneq = matmul( layout%moment%toMoments%A(1:QQ, 1:QQ), fneq(1:QQ) )

      ! MRT omegas
      ! overwrite omegaKine term in the element loop
      omegaKine = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem)
      s_mrt = fieldProp(1)%fluid%mrtPtr(omegaKine=omegaKine, &
        &                               omegaBulk=omegaBulk, &
        &                               QQ=QQ                )

      ! compute Minv times relaxation matrix
      do iDir = 1, QQ
        mInvXOmega(1:QQ,iDir) = layout%moment%toPDF%A(1:QQ, iDir) * s_mrt(iDir)
      end do

      do iDir = 1, QQ
        ! collision term: M^-1 S M (f-feq)
        colli_term = dot_product( mInvXOmega(iDir, 1:QQ), mneq(1:QQ) )

        outState(?SAVE?(iDir,1,iElem,QQ,QQ,nElems,neigh)) &
          & = pdfTmp(iDir) - colli_term
      end do

    enddo nodeloop

  end subroutine mrt_advRel_generic
! **************************************************************************** !

! **************************************************************************** !
  !> Unoptimized explicit implementation
  !!
  !! This subroutine interface must match the abstract interface definition
  !! [[kernel]] in scheme/[[mus_scheme_type_module]].f90 in order to be callable
  !! via [[mus_scheme_type:compute]] function pointer.
  subroutine mrt_advRel_incomp_generic( fieldProp, inState, outState,    &
    &                                   auxField, neigh, nElems, nSolve, &
    &                                   level, layout, params, varSys,   &
    &                                   derVarPos                        )
    ! -------------------------------------------------------------------- !
    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), intent(in) :: fieldProp(:)
    !> variable system definition
    type(tem_varSys_type), intent(in) :: varSys
    !> current layout
    type(mus_scheme_layout_type), intent(in) :: layout
    !> number of elements in state Array
    integer, intent(in) :: nElems
    !> input  pdf vector
    real(kind=rk), intent(in)  ::  inState(nElems * varSys%nScalars)
    !> output pdf vector
    real(kind=rk), intent(out) :: outState(nElems * varSys%nScalars)
    !> Auxiliary field computed from pre-collision state
    !! Is updated with correct velocity field for multicomponent models
    real(kind=rk), intent(inout) :: auxField(nElems * varSys%nAuxScalars)
    !> connectivity vector
    integer, intent(in) :: neigh(nElems * layout%fStencil%QQ)
    !> number of elements solved in kernel
    integer, intent(in) :: nSolve
    !> current level
    integer,intent(in) :: level
    !> global parameters
    type(mus_param_type),intent(in) :: params
    !> position of derived quantities in varsys for all fields
    type( mus_derVarPos_type ), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer       :: iElem, iDir
    real(kind=rk) :: pdfTmp( layout%fStencil%QQ ) ! temporary local pdf values
    real(kind=rk) :: rho     ! local density
    real(kind=rk) :: u_x     ! local x-velocity
    real(kind=rk) :: u_y     ! local y-velocity
    real(kind=rk) :: u_z     ! local z-velocity
    real(kind=rk) :: usq, ucx
    real(kind=rk) :: fEq( layout%fStencil%QQ ) !< equilibrium distribution
    ! MRT Variables
    real(kind=rk) :: fneq( layout%fStencil%QQ )
    real(kind=rk) :: mneq( layout%fStencil%QQ )
    real(kind=rk) :: s_mrt( layout%fStencil%QQ )
    real(kind=rk) :: mInvXomega(layout%fStencil%QQ, layout%fStencil%QQ)
    integer :: QQ
    integer :: dens_pos, vel_pos(3), elemOff
    real(kind=rk) :: omegaKine, omegaBulk, colli_term
    ! ---------------------------------------------------------------------------
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    QQ = layout%fStencil%QQ
    omegaBulk = fieldProp(1)%fluid%omegaBulkLvl(level)

    !NEC$ ivdep
?? IF (SOA) THEN
?? copy :: dir_vector
?? ELSE
?? copy :: dir_novec
?? END IF
    nodeloop: do iElem = 1,nSolve

      !> First load all local values into temp array
      do iDir = 1, QQ
        pdfTmp( iDir ) = inState( ?FETCH?( iDir, 1, iElem, QQ, QQ, nElems,neigh ))
      end do

      ! element offset for auxField array
      elemOff = (iElem-1)*varSys%nAuxScalars
      ! local density
      rho = auxField(elemOff + dens_pos)

      ! local x-, y- and z-velocity
      u_x = auxField(elemOff + vel_pos(1))
      u_y = auxField(elemOff + vel_pos(2))
      u_z = auxField(elemOff + vel_pos(3))

      usq = u_x*u_x + u_y*u_y + u_z*u_z

      ! -------------------------------------------------------------------------
      ! Calculate the equilibrium distribution function
      do iDir = 1, QQ

        ! velocity times lattice unit velocity
        ucx =   layout%fStencil%cxDirRK(1,iDir) * u_x       &
          &   + layout%fStencil%cxDirRK(2,iDir) * u_y       &
          &   + layout%fStencil%cxDirRK(3,iDir) * u_z

        ! calculate equilibrium density
        fEq( iDir ) = layout%weight( iDir ) * ( rho + rho0*( ucx*cs2inv &
          &         + ucx*ucx*t2cs4inv - usq*t2cs2inv ) )

      enddo

      ! compute fNeq
      fneq = pdfTmp - fEq

      ! convert fneq to moments neq
      mneq = matmul( layout%moment%toMoments%A(1:QQ, 1:QQ), fneq(1:QQ) )

      ! MRT omegas
      ! overwrite omegaKine term in the element loop
      omegaKine = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem)
      s_mrt = fieldProp(1)%fluid%mrtPtr(omegaKine=omegaKine, &
        &                               omegaBulk=omegaBulk, &
        &                               QQ=QQ                )

      ! compute Minv times relaxation matrix
      do iDir = 1, QQ
        mInvXOmega(1:QQ,iDir) = layout%moment%toPDF%A(1:QQ, iDir) * s_mrt(iDir)
      end do

      do iDir = 1, QQ
        ! collision term: M^-1 S M (f-feq)
        colli_term = dot_product( mInvXOmega(iDir, 1:QQ), mneq(1:QQ) )
        outState(?SAVE?(iDir,1,iElem,QQ,QQ,nElems,neigh)) &
          & = pdfTmp(iDir) - colli_term
      end do

    enddo nodeloop

  end subroutine mrt_advRel_incomp_generic
! **************************************************************************** !

end module mus_mrt_module
! ****************************************************************************** !
