! Copyright (c) 2017, 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2019 Harald Klimach <harald.klimach@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_d3q27Macros.inc'
! **************************************************************************** !
!> author: Kannan Masilamani
!! This module provides the definition and methods for
!! MRT advection relaxation scheme for D3Q27 stencil.
!! This implementation is based on the paper
!! K. Suga et al, "A D3Q27 multiple-relaxation-time lattice Boltzmann
!! method for turbulent flows", Computers and Mathematics with Applications 69
!! (2015) 518-529 \n
!!
!! The LB equaton using MRT is
!!    f(t+dt,x+dx) = f - M^(-1) * S * ( (M*f) - m^(eq) )
!! The moments m(1:27) = M * f(1:27) are labeled as
!!  ! Density
!!  m( 1) = rho
!!  ! momentum
!!  m( 2) = jx = rho * ux
!!  m( 3) = jy = rho * uy
!!  m( 4) = jz = rho * uz
!!  ! kinetic energy
!!  m( 5) = e = rho * (ux^2 + uy^2 + uz^2)
!!  m( 6) = XX = rho * (2*ux*2 - uy^2 - uz^2)
!!  m( 7) = WW = rho * (uy^2 - uz^2)
!!  m( 8) = XY = rho * ux * uy
!!  m( 9) = YZ = rho * ux * uy
!!  m(10) = ZX = rho * uz * ux
!!  ! fluxes of the energy and square of energy
!!  m(11) = phix = 3 * rho * (ux^2 + uy^2 + uz^2) * ux
!!  m(12) = phiy = 3 * rho * (ux^2 + uy^2 + uz^2) * uy
!!  m(13) = phiz = 3 * rho * (ux^2 + uy^2 + uz^2) * uz
!!  m(14) = psix = (9/2) * rho * (ux^2 + uy^2 + uz^2)^2 * ux
!!  m(15) = psix = (9/2) * rho * (ux^2 + uy^2 + uz^2)^2 * uy
!!  m(16) = psix = (9/2) * rho * (ux^2 + uy^2 + uz^2)^2 * uz
!!  ! Square and cube of the energy
!!  m(17) = e2 = (3/2) * rho * (ux^2 + uy^2 + uz^2)^2
!!  m(18) = e3 = (9/2) * rho * (ux^2 + uy^2 + uz^2)^3
!!  ! Product of the second order tensor and the energy
!!  m(19) = XXe = rho * (2*ux*2 - uy^2 - uz^2) * (ux*2 + uy^2 + uz^2)
!!  m(20) = WWe = rho * (uy^2 - uz^2) * (ux*2 + uy^2 + uz^2)
!!  m(21) = XY = rho * ux * uy * (ux*2 + uy^2 + uz^2)
!!  m(22) = YZ = rho * uy * uz * (ux*2 + uy^2 + uz^2)
!!  m(23) = ZX = rho * uz * ux * (ux*2 + uy^2 + uz^2)
!!  ! third order psuedo-vector and totally antisymmetric tensor XYZ
!!  m(24) = taux = rho * ux * (uy^2 - uz^2)
!!  m(25) = tauy = rho * uy * (uz^2 - ux^2)
!!  m(26) = tauz = rho * uz * (ux^2 - uy^2)
!!  m(27) = XYX = rho * ux * uy * uz
!!
!! The non-zero equilibrium moments are given by
!!  meq( 1) =                                         rho
!!  meq( 2) =                                      rho*ux
!!  meq( 3) =                                      rho*uy
!!  meq( 4) =                                      rho*uz
!!  meq( 5) =        rho*ux^2 + rho*uy^2 + rho*uz^2 - rho
!!  meq( 6) =            2*rho*ux^2 - rho*uy^2 - rho*uz^2
!!  meq( 7) =                         rho*uy^2 - rho*uz^2
!!  meq( 8) =                                   rho*ux*uy
!!  meq( 9) =                                   rho*uy*uz
!!  meq(10) =                                   rho*ux*uz
!!  meq(11) =                                   -2*rho*ux
!!  meq(12) =                                   -2*rho*uy
!!  meq(13) =                                   -2*rho*uz
!!  meq(14) =                                      rho*ux
!!  meq(15) =                                      rho*uy
!!  meq(16) =                                      rho*uz
!!  meq(17) = -2*rho*ux^2 - 2*rho*uy^2 - 2*rho*uz^2 + rho
!!  meq(18) =  3*rho*ux^2 + 3*rho*uy^2 + 3*rho*uz^2 - rho
!!  meq(19) =           -2*rho*ux^2 + rho*uy^2 + rho*uz^2
!!  meq(20) =                        -rho*uy^2 + rho*uz^2
!!  meq(21) =                                  -rho*ux*uy
!!  meq(22) =                                  -rho*uy*uz
!!  meq(23) =                                  -rho*ux*uz
!!  meq(24) =                                           0
!!  meq(25) =                                           0
!!  meq(26) =                                           0
!!  meq(27) =                                           0
!!
!! Density (rho) and velocity (ux, uy, uz) are conserved during collision.
!!  i.e. m(1) = meq(1) --> mneq(1) = 0
!!       m(2) = meq(2) --> mneq(2) = 0
!!       m(3) = meq(3) --> mneq(3) = 0
!!       m(4) = meq(4) --> mneq(4) = 0
!!
!! Collision parameters are chosen as
!! s(1:4) = max(omega, 1.0)
!! s(5) = bulk_omega
!! s(6:10) = omega
!! s(11:13) = 1.5
!! s(14:16) = 1.83
!! s(17:18) = 1.61
!! s(19:23) = 1.98
!! s(23:27) = 1.74
!!
!! 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_d3q27_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, rho0, &
    &                                 cs2inv, cs4inv, t2cs2inv, t2cs4inv

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

  implicit none

  private

  public :: mrt_advRel_d3q27
  public :: mrt_advRel_d3q27_incomp
  public :: mrt_advRel_d3q27_generic
  public :: mrt_advRel_d3q27_incomp_generic

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

  integer,parameter :: QQ   = 27   !< 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 :: qNNN = 19    !<               x-,y-,z-
  integer,parameter :: qNN1 = 20    !<               x-,y-,z+
  integer,parameter :: qN1N = 21    !<               x-,y+,z-
  integer,parameter :: qN11 = 22    !<               x-,y+,z+
  integer,parameter :: q1NN = 23    !<               x+,y-,z-
  integer,parameter :: q1N1 = 24    !<               x+,y-,z+
  integer,parameter :: q11N = 25    !<               x+,y+,z-
  integer,parameter :: q111 = 26    !<               x+,y+,z+
  integer,parameter :: q000 = 27    !< rest density is last

  ! D3Q27 MRT pdf -> moment transformation matrix
  ! How to use:
  ! do iDir = 1, QQ
  !   moment(iDir) = sum( PDF(:) * MMtrD3Q27(iDir,:) )
  ! end do
  !  W      S     B     E     N     T     BS    TS   BN    TN    BW    BE    TW
  !  TE    SW    NW    SE    NE   BSW    TSW   BNW  TNW   BSE   TSE   BNE   TNE  0
  real(kind=rk), dimension(27,27), parameter, public :: MMtrD3Q27 =           &
&  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
&   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, & ! 2
&  -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, 1.0, 0.0,-1.0,-1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0,-1.0, & ! 3
&   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,-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, & ! 4
&   0.0, 0.0, 0.0,-1.0, 1.0,-1.0, 1.0,-1.0, 1.0,-1.0, 1.0, 0.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, & ! 5
&   0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,-2.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, & ! 6
&   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, 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, & ! 7
&   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, 0.0, 0.0, 0.0, 0.0, 1.0, & ! 8
&  -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, 0.0, 1.0,-1.0,-1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, & ! 9
&   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, 1.0,-1.0,-1.0, 1.0, 0.0, & ! 10
&   0.0, 0.0, 0.0, 1.0,-1.0, 1.0,-1.0,-1.0, 1.0,-1.0, 1.0, 0.0,                &
&   4.0, 0.0, 0.0,-4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,-1.0, 1.0,-1.0, 1.0, & ! 11
&   1.0,-1.0,-1.0,-2.0,-2.0,-2.0,-2.0, 2.0, 2.0, 2.0, 2.0, 0.0,                &
&   0.0, 4.0, 0.0, 0.0,-4.0, 0.0, 1.0, 1.0,-1.0,-1.0, 0.0, 0.0, 0.0, 0.0, 1.0, & ! 12
&  -1.0, 1.0,-1.0,-2.0,-2.0, 2.0, 2.0,-2.0,-2.0, 2.0, 2.0, 0.0,                &
&   0.0, 0.0, 4.0, 0.0, 0.0,-4.0, 1.0,-1.0, 1.0,-1.0, 1.0, 1.0,-1.0,-1.0, 0.0, & ! 13
&   0.0, 0.0, 0.0,-2.0, 2.0,-2.0, 2.0,-2.0, 2.0,-2.0, 2.0, 0.0,                &
&  -4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0,-2.0, 2.0,-2.0, 2.0, & ! 14
&   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,-4.0, 0.0, 0.0, 4.0, 0.0, 2.0, 2.0,-2.0,-2.0, 0.0, 0.0, 0.0, 0.0, 2.0, & ! 15
&  -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, 0.0,-4.0, 0.0, 0.0, 4.0, 2.0,-2.0, 2.0,-2.0, 2.0, 2.0,-2.0,-2.0, 0.0, & ! 16
&   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,-1.0,-1.0,-1.0,-1.0,-1.0, & ! 17
&  -1.0,-1.0,-1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 4.0,                &
&   4.0, 4.0, 4.0, 4.0, 4.0, 4.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0, & ! 18
&  -2.0,-2.0,-2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,-8.0,                &
&  -4.0, 2.0, 2.0,-4.0, 2.0, 2.0,-2.0,-2.0,-2.0,-2.0, 1.0, 1.0, 1.0, 1.0, 1.0, & ! 19
&   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,-2.0, 2.0, 0.0,-2.0, 2.0, 0.0, 0.0, 0.0, 0.0,-1.0,-1.0,-1.0,-1.0, 1.0, & ! 20
&   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, 0.0, 0.0, 0.0, 0.0,-2.0, & ! 21
&   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, 0.0, 0.0, 0.0, 0.0, 0.0,-2.0, 2.0, 2.0,-2.0, 0.0, 0.0, 0.0, 0.0, 0.0, & ! 22
&   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,-2.0, 2.0, 2.0,-2.0, 0.0, & ! 23
&   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, 1.0,-1.0, 1.0,-1.0,-1.0, & ! 24
&  -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, 1.0, & ! 25
&  -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, & ! 26
&   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.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & ! 27
&   0.0, 0.0, 0.0,-1.0, 1.0, 1.0,-1.0, 1.0,-1.0,-1.0, 1.0, 0.0                 &
&  /),(/27,27/), order=(/ 2,1 /) )

  ! 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(27,27),parameter,public  :: MMIvD3Q27 =             &
& reshape((/                                                                   &
& 1/27.0, -1/18.0, 0.0, 0.0, -1/18.0, 1/18.0, 0.0, 0.0, 0.0, 0.0, 1/18.0,      &
& 0.0, 0.0, -1/18.0, 0.0, 0.0, 0.0, 1/54.0, -1/18.0, 0.0, 0.0, 0.0,            &
& 0.0, 0.0, 0.0, 0.0, 0.0,                                                     &
& 1/27.0, 0.0, -1/18.0, 0.0, -1/18.0, -1/36.0, 1/12.0, 0.0, 0.0, 0.0, 0.0,     &
& 1/18.0, 0.0, 0.0, -1/18.0, 0.0, 0.0, 1/54.0, 1/36.0, -1/12.0, 0.0, 0.0,      &
& 0.0, 0.0, 0.0, 0.0, 0.0,                                                     &
& 1/27.0, 0.0, 0.0, -1/18.0, -1/18.0, -1/36.0, -1/12.0, 0.0, 0.0, 0.0, 0.0,    &
& 0.0, 1/18.0, 0.0, 0.0, -1/18.0, 0.0, 1/54.0, 1/36.0, 1/12.0, 0.0, 0.0,       &
& 0.0, 0.0, 0.0, 0.0, 0.0,                                                     &
& 1/27.0, 1/18.0, 0.0, 0.0, -1/18.0, 1/18.0, 0.0, 0.0, 0.0, 0.0, -1/18.0,      &
& 0.0, 0.0, 1/18.0, 0.0, 0.0, 0.0, 1/54.0, -1/18.0, 0.0, 0.0, 0.0,             &
& 0.0, 0.0, 0.0, 0.0, 0.0,                                                     &
& 1/27.0, 0.0, 1/18.0, 0.0, -1/18.0, -1/36.0, 1/12.0, 0.0, 0.0, 0.0, 0.0,      &
& -1/18.0, 0.0, 0.0, 1/18.0, 0.0, 0.0, 1/54.0, 1/36.0, -1/12.0, 0.0, 0.0,      &
& 0.0, 0.0, 0.0, 0.0, 0.0,                                                     &
& 1/27.0, 0.0, 0.0, 1/18.0, -1/18.0, -1/36.0, -1/12.0, 0.0, 0.0, 0.0, 0.0,     &
& 0.0, -1/18.0, 0.0, 0.0, 1/18.0, 0.0, 1/54.0, 1/36.0, 1/12.0, 0.0, 0.0,       &
& 0.0, 0.0, 0.0, 0.0, 0.0,                                                     &
& 1/27.0, 0.0, -1/18.0, -1/18.0, 0.0, -1/18.0, 0.0, 0.0, 1/12.0, 0.0, 0.0,     &
& 1/72.0, 1/72.0, 0.0, 1/36.0, 1/36.0, -1/36.0, -1/108.0, -1/36.0, 0.0, 0.0, -1/12.0, &
& 0.0, 0.0, -1/8.0, 1/8.0, 0.0,                                                &
& 1/27.0, 0.0, -1/18.0, 1/18.0, 0.0, -1/18.0, 0.0, 0.0, -1/12.0, 0.0, 0.0,     &
& 1/72.0, -1/72.0, 0.0, 1/36.0, -1/36.0, -1/36.0, -1/108.0, -1/36.0, 0.0, 0.0, 1/12.0, &
& 0.0, 0.0, -1/8.0, -1/8.0, 0.0,                                               &
& 1/27.0, 0.0, 1/18.0, -1/18.0, 0.0, -1/18.0, 0.0, 0.0, -1/12.0, 0.0, 0.0,     &
& -1/72.0, 1/72.0, 0.0, -1/36.0, 1/36.0, -1/36.0, -1/108.0, -1/36.0, 0.0, 0.0, 1/12.0, &
& 0.0, 0.0, 1/8.0, 1/8.0, 0.0,                                                 &
& 1/27.0, 0.0, 1/18.0, 1/18.0, 0.0, -1/18.0, 0.0, 0.0, 1/12.0, 0.0, 0.0,       &
& -1/72.0, -1/72.0, 0.0, -1/36.0, -1/36.0, -1/36.0, -1/108.0, -1/36.0, 0.0, 0.0, -1/12.0, &
& 0.0, 0.0, 1/8.0, -1/8.0, 0.0,                                                &
& 1/27.0, -1/18.0, 0.0, -1/18.0, 0.0, 1/36.0, -1/12.0, 0.0, 0.0, 1/12.0, 1/72.0, &
& 0.0, 1/72.0, 1/36.0, 0.0, 1/36.0, -1/36.0, -1/108.0, 1/72.0, -1/24.0, 0.0, 0.0, &
& -1/12.0, 1/8.0, 0.0, -1/8.0, 0.0,                                            &
& 1/27.0, 1/18.0, 0.0, -1/18.0, 0.0, 1/36.0, -1/12.0, 0.0, 0.0, -1/12.0, -1/72.0, &
& 0.0, 1/72.0, -1/36.0, 0.0, 1/36.0, -1/36.0, -1/108.0, 1/72.0, -1/24.0, 0.0, 0.0, &
& 1/12.0, -1/8.0, 0.0, -1/8.0, 0.0,                                            &
& 1/27.0, -1/18.0, 0.0, 1/18.0, 0.0, 1/36.0, -1/12.0, 0.0, 0.0, -1/12.0, 1/72.0, &
& 0.0, -1/72.0, 1/36.0, 0.0, -1/36.0, -1/36.0, -1/108.0, 1/72.0, -1/24.0, 0.0, 0.0, &
& 1/12.0, 1/8.0, 0.0, 1/8.0, 0.0,                                              &
& 1/27.0, 1/18.0, 0.0, 1/18.0, 0.0, 1/36.0, -1/12.0, 0.0, 0.0, 1/12.0, -1/72.0, &
& 0.0, -1/72.0, -1/36.0, 0.0, -1/36.0, -1/36.0, -1/108.0, 1/72.0, -1/24.0, 0.0, 0.0, &
& -1/12.0, -1/8.0, 0.0, 1/8.0, 0.0,                                            &
& 1/27.0, -1/18.0, -1/18.0, 0.0, 0.0, 1/36.0, 1/12.0, 1/12.0, 0.0, 0.0, 1/72.0, &
& 1/72.0, 0.0, 1/36.0, 1/36.0, 0.0, -1/36.0, -1/108.0, 1/72.0, 1/24.0, -1/12.0, 0.0, &
& 0.0, -1/8.0, 1/8.0, 0.0, 0.0,                                                &
& 1/27.0, -1/18.0, 1/18.0, 0.0, 0.0, 1/36.0, 1/12.0, -1/12.0, 0.0, 0.0, 1/72.0, &
& -1/72.0, 0.0, 1/36.0, -1/36.0, 0.0, -1/36.0, -1/108.0, 1/72.0, 1/24.0, 1/12.0, 0.0, &
& 0.0, -1/8.0, -1/8.0, 0.0, 0.0,                                               &
& 1/27.0, 1/18.0, -1/18.0, 0.0, 0.0, 1/36.0, 1/12.0, -1/12.0, 0.0, 0.0, -1/72.0, &
& 1/72.0, 0.0, -1/36.0, 1/36.0, 0.0, -1/36.0, -1/108.0, 1/72.0, 1/24.0, 1/12.0, 0.0, &
& 0.0, 1/8.0, 1/8.0, 0.0, 0.0,                                                 &
& 1/27.0, 1/18.0, 1/18.0, 0.0, 0.0, 1/36.0, 1/12.0, 1/12.0, 0.0, 0.0, -1/72.0, &
& -1/72.0, 0.0, -1/36.0, -1/36.0, 0.0, -1/36.0, -1/108.0, 1/72.0, 1/24.0, -1/12.0, 0.0, &
& 0.0, 1/8.0, -1/8.0, 0.0, 0.0,                                                &
& 1/27.0, -1/18.0, -1/18.0, -1/18.0, 1/18.0, 0.0, 0.0, 1/12.0, 1/12.0, 1/12.0, -1/36.0, &
& -1/36.0, -1/36.0, -1/72.0, -1/72.0, -1/72.0, 1/36.0, 1/216.0, 0.0, 0.0, 1/24.0, 1/24.0, &
& 1/24.0, 0.0, 0.0, 0.0, -1/8.0,                                               &
& 1/27.0, -1/18.0, -1/18.0, 1/18.0, 1/18.0, 0.0, 0.0, 1/12.0, -1/12.0, -1/12.0, -1/36.0, &
& -1/36.0, 1/36.0, -1/72.0, -1/72.0, 1/72.0, 1/36.0, 1/216.0, 0.0, 0.0, 1/24.0, -1/24.0, &
& -1/24.0, 0.0, 0.0, 0.0, 1/8.0,                                               &
& 1/27.0, -1/18.0, 1/18.0, -1/18.0, 1/18.0, 0.0, 0.0, -1/12.0, -1/12.0, 1/12.0, -1/36.0, &
& 1/36.0, -1/36.0, -1/72.0, 1/72.0, -1/72.0, 1/36.0, 1/216.0, 0.0, 0.0, -1/24.0, -1/24.0, &
& 1/24.0, 0.0, 0.0, 0.0, 1/8.0,                                                &
& 1/27.0, -1/18.0, 1/18.0, 1/18.0, 1/18.0, 0.0, 0.0, -1/12.0, 1/12.0, -1/12.0, -1/36.0, &
& 1/36.0, 1/36.0, -1/72.0, 1/72.0, 1/72.0, 1/36.0, 1/216.0, 0.0, 0.0, -1/24.0, 1/24.0, &
& -1/24.0, 0.0, 0.0, 0.0, -1/8.0,                                              &
& 1/27.0, 1/18.0, -1/18.0, -1/18.0, 1/18.0, 0.0, 0.0, -1/12.0, 1/12.0, -1/12.0, 1/36.0, &
& -1/36.0, -1/36.0, 1/72.0, -1/72.0, -1/72.0, 1/36.0, 1/216.0, 0.0, 0.0, -1/24.0, 1/24.0, &
& -1/24.0, 0.0, 0.0, 0.0, 1/8.0,                                               &
& 1/27.0, 1/18.0, -1/18.0, 1/18.0, 1/18.0, 0.0, 0.0, -1/12.0, -1/12.0, 1/12.0, 1/36.0, &
& -1/36.0, 1/36.0, 1/72.0, -1/72.0, 1/72.0, 1/36.0, 1/216.0, 0.0, 0.0, -1/24.0, -1/24.0, &
& 1/24.0, 0.0, 0.0, 0.0, -1/8.0,                                               &
& 1/27.0, 1/18.0, 1/18.0, -1/18.0, 1/18.0, 0.0, 0.0, 1/12.0, -1/12.0, -1/12.0, 1/36.0, &
& 1/36.0, -1/36.0, 1/72.0, 1/72.0, -1/72.0, 1/36.0, 1/216.0, 0.0, 0.0, 1/24.0, -1/24.0, &
& -1/24.0, 0.0, 0.0, 0.0, -1/8.0,                                              &
& 1/27.0, 1/18.0, 1/18.0, 1/18.0, 1/18.0, 0.0, 0.0, 1/12.0, 1/12.0, 1/12.0, 1/36.0, &
& 1/36.0, 1/36.0, 1/72.0, 1/72.0, 1/72.0, 1/36.0, 1/216.0, 0.0, 0.0, 1/24.0, 1/24.0, &
& 1/24.0, 0.0, 0.0, 0.0, 1/8.0,                                                &
& 1/27.0, 0.0, 0.0, 0.0, -1/9.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,                 &
& 0.0, 0.0, 0.0, 0.0, 0.0, 1/9.0, -1/27.0, 0.0, 0.0, 0.0, 0.0,                 &
& 0.0, 0.0, 0.0, 0.0, 0.0                                                      &
&  /),(/27,27/), order=(/ 2,1 /) )

contains

! **************************************************************************** !
  !> Semi-optimized 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_d3q27( 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
    integer :: nScalars
?? copy :: var27( 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
    real(kind=rk) :: usq
    real(kind=rk) :: sum_e1, sum_e2, sum_e3
    real(kind=rk) :: sum_ux1, sum_ux2, sum_ux3
    real(kind=rk) :: sum_uy1, sum_uy2, sum_uy3
    real(kind=rk) :: sum_uz1, sum_uz2, sum_uz3
    ! 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) :: omegaBulk
    ! ---------------------------------------------------------------------------
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)
    nScalars = varSys%nScalars

    omegaBulk = fieldProp(1)%fluid%omegaBulkLvl(level)
    ! overwrite omegaKine term in the element loop
    s_mrt = mrt_d3q27(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
      fN00 = inState(?FETCH?( qN00, 1, iElem, QQ, nScalars, nElems,neigh ))
      f0N0 = inState(?FETCH?( q0N0, 1, iElem, QQ, nScalars, nElems,neigh ))
      f00N = inState(?FETCH?( q00N, 1, iElem, QQ, nScalars, nElems,neigh ))
      f100 = inState(?FETCH?( q100, 1, iElem, QQ, nScalars, nElems,neigh ))
      f010 = inState(?FETCH?( q010, 1, iElem, QQ, nScalars, nElems,neigh ))
      f001 = inState(?FETCH?( q001, 1, iElem, QQ, nScalars, nElems,neigh ))

      f0NN = inState(?FETCH?( q0NN, 1, iElem, QQ, nScalars, nElems,neigh ))
      f0N1 = inState(?FETCH?( q0N1, 1, iElem, QQ, nScalars, nElems,neigh ))
      f01N = inState(?FETCH?( q01N, 1, iElem, QQ, nScalars, nElems,neigh ))
      f011 = inState(?FETCH?( q011, 1, iElem, QQ, nScalars, nElems,neigh ))
      fN0N = inState(?FETCH?( qN0N, 1, iElem, QQ, nScalars, nElems,neigh ))
      f10N = inState(?FETCH?( q10N, 1, iElem, QQ, nScalars, nElems,neigh ))
      fN01 = inState(?FETCH?( qN01, 1, iElem, QQ, nScalars, nElems,neigh ))
      f101 = inState(?FETCH?( q101, 1, iElem, QQ, nScalars, nElems,neigh ))
      fNN0 = inState(?FETCH?( qNN0, 1, iElem, QQ, nScalars, nElems,neigh ))
      fN10 = inState(?FETCH?( qN10, 1, iElem, QQ, nScalars, nElems,neigh ))
      f1N0 = inState(?FETCH?( q1N0, 1, iElem, QQ, nScalars, nElems,neigh ))
      f110 = inState(?FETCH?( q110, 1, iElem, QQ, nScalars, nElems,neigh ))

      fNNN = inState(?FETCH?( qNNN, 1, iElem, QQ, nScalars, nElems,neigh ))
      fNN1 = inState(?FETCH?( qNN1, 1, iElem, QQ, nScalars, nElems,neigh ))
      fN1N = inState(?FETCH?( qN1N, 1, iElem, QQ, nScalars, nElems,neigh ))
      fN11 = inState(?FETCH?( qN11, 1, iElem, QQ, nScalars, nElems,neigh ))
      f1NN = inState(?FETCH?( q1NN, 1, iElem, QQ, nScalars, nElems,neigh ))
      f1N1 = inState(?FETCH?( q1N1, 1, iElem, QQ, nScalars, nElems,neigh ))
      f11N = inState(?FETCH?( q11N, 1, iElem, QQ, nScalars, nElems,neigh ))
      f111 = inState(?FETCH?( q111, 1, iElem, QQ, nScalars, nElems,neigh ))

      f000 = inState(?FETCH?( q000, 1, iElem, QQ, nScalars, nElems,neigh ))

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

      ! Zero moment
      mom(1) = rho

      ! First moments
      mom( 2) = rho*u_x
      mom( 3) = rho*u_y
      mom( 4) = rho*u_z

      sum_e1 = f001 + f00N + f010 + f0N0 + f100 + fN00
      sum_e3 = f111 + f11N + f1N1 + f1NN + fN11 + fN1N + fNN1 + fNNN
      sum_e2 = f011 + f01N + f0N1 + f0NN + f101 + f10N &
        &    + f110 + f1N0 + fN01 + fN0N + fN10 + fNN0
      ! kinematic energy
      mom( 5) = - 2.0_rk*f000 - sum_e1 + sum_e3
      ! square of energy
      mom(17) = 4.0_rk*f000 - sum_e2 + sum_e3
      ! cube of energy
      mom(18) = - 8.0_rk*f000 + 4.0_rk * sum_e1 - 2.0_rk * sum_e2 + sum_e3

      ! Second order tensors
      mom( 6) = 2.0_rk*( f100 + fN00 - f011 - f01N - f0N1 - f0NN)   &
        &     - f001 - f00N - f010 - f0N0                           &
        &     + f101 + f10N + f110 + f1N0 + fN01 + fN0N + fN10 + fNN0
      mom( 7) = - f001 - f00N + f010 + f0N0                            &
        &       - f101 - f10N + f110 + f1N0 - fN01 - fN0N + fN10 + fNN0
      mom( 8) = f110 + f111 + f11N - f1N0 - f1N1 - f1NN - fN10 - fN11 &
        &     - fN1N + fNN0 + fNN1 + fNNN
      mom( 9) = f011 - f01N - f0N1 + f0NN + f111 - f11N - f1N1 + f1NN &
        &     + fN11 - fN1N - fNN1 + fNNN
      mom(10) = f101 - f10N + f111 - f11N + f1N1 - f1NN - fN01 + fN0N &
        &     - fN11 + fN1N - fNN1 + fNNN

      ! fluxes of the energy and square of energy
      sum_ux1 = 4.0_rk*(f100 - fN00)
      sum_ux2 = - f101 - f10N - f110 - f1N0 +  fN01 + fN0N + fN10 + fNN0
      sum_ux3 = f111 + f11N + f1N1 + f1NN - fN11 - fN1N - fNN1 - fNNN
      mom(11) = - sum_ux1 + sum_ux2 + 2.0_rk * sum_ux3
      mom(14) =   sum_ux1 + 2.0_rk * sum_ux2 + sum_ux3

      sum_uy1 = 4.0_rk*(f010 - f0N0)
      sum_uy2 = - f011 - f01N  - f110 - fN10 + f0N1 + f0NN + fNN0 + f1N0
      sum_uy3 = f111 + f11N - f1N1 - f1NN + fN11 + fN1N - fNN1 - fNNN
      mom(12) = - sum_uy1 + sum_uy2 + 2.0_rk * sum_uy3
      mom(15) =   sum_uy1 + 2.0_rk * sum_uy2 + sum_uy3

      sum_uz1 = 4.0_rk*(f001 - f00N)
      sum_uz2 = - f011 - f0N1 - f101 - fN01 + f01N + f0NN + f10N + fN0N
      sum_uz3 = f111 - f11N + f1N1 - f1NN + fN11 - fN1N + fNN1 - fNNN
      mom(13) = - sum_uz1 + sum_uz2 + 2.0_rk * sum_uz3
      mom(16) =   sum_uz1 + 2.0_rk * sum_uz2 + sum_uz3

      !product of the second order tensors and the energy
      mom(19) = - 4.0_rk*(f100 + fN00)                                         &
        &     + 2.0_rk*( f001 + f00N + f010 - f011 - f01N + f0N0 - f0N1 - f0NN ) &
        &     + f101 + f10N + f110 + f1N0 + fN01 + fN0N + fN10 + fNN0
      mom(20) = 2.0_rk*(f001 + f00N - f010 - f0N0) &
        &     - f101 - f10N + f110 + f1N0 - fN01 - fN0N + fN10 + fNN0
      mom(21) = 2.0_rk*(-f110 + f1N0 + fN10 - fNN0) &
        &     + f111 + f11N - f1N1 - f1NN - fN11 - fN1N + fNN1 + fNNN
      mom(22) = 2.0_rk*(-f011 + f01N + f0N1 - f0NN ) &
        &     + f111 - f11N - f1N1 + f1NN + fN11 - fN1N - fNN1 + fNNN
      mom(23) = 2.0_rk*(-f101 + f10N + fN01 - fN0N ) &
        &     + f111 - f11N + f1N1 - f1NN - fN11 + fN1N - fNN1 + fNNN

      ! third order pseudo-vector
      mom(24) = -f101 - f10N + f110 + f1N0 + fN01 + fN0N - fN10 - fNN0
      mom(25) = f011 + f01N - f0N1 - f0NN - f110 + f1N0 - fN10 + fNN0
      mom(26) = -f011 + f01N - f0N1 + f0NN + f101 - f10N + fN01 - fN0N
      ! totally antisymmetric tensor XYZ
      mom(27) = f111 - f11N - f1N1 + f1NN - fN11 + fN1N + fNN1 - fNNN

      ! square of velocity
      usq = u_x*u_x + u_y*u_y + u_z*u_z
      ! -------------------------------------------------------------------------
      ! Equilibrium moments
      meq(1:QQ) =  0.0_rk
      meq( 1) = rho
      meq( 2) = rho*u_x
      meq( 3) = rho*u_y
      meq( 4) = rho*u_z

      meq(14) = meq( 2)
      meq(15) = meq( 3)
      meq(16) = meq( 4)

      meq(11) = -2.0_rk*meq( 2)
      meq(12) = -2.0_rk*meq( 3)
      meq(13) = -2.0_rk*meq( 4)

      meq( 5) = rho*usq - rho
      meq( 6) = rho * (2.0_rk*u_x*u_x - u_y*u_y - u_z*u_z)
      meq(19) = -meq( 6)

      meq( 7) = rho * (u_y*u_y - u_z*u_z)
      meq( 8) = rho*u_x*u_y
      meq( 9) = rho*u_y*u_z
      meq(10) = rho*u_x*u_z

      meq(20) = -meq( 7)
      meq(21) = -meq( 8)
      meq(22) = -meq( 9)
      meq(23) = -meq(10)

      meq(17) = -2.0_rk*rho*usq + rho
      meq(18) =  3.0_rk*rho*usq - rho

      ! update kinematic omega
      s_mrt(6:10) = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem)
      ! compute neq moment
      mneq(:) = s_mrt(:) * ( mom(:) - meq(:) )

      ! compute fNeq
      do iDir = 1, QQ
        fneq(iDir) = sum( MMIvD3Q27(iDir,:) * mneq(:) )
      end do

      outState(?SAVE?( qN00, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fN00 - fneq(qN00)
      outState(?SAVE?( q0N0, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f0N0 - fneq(q0N0)
      outState(?SAVE?( q00N, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f00N - fneq(q00N)
      outState(?SAVE?( q100, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f100 - fneq(q100)
      outState(?SAVE?( q010, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f010 - fneq(q010)
      outState(?SAVE?( q001, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f001 - fneq(q001)

      outState(?SAVE?( q0NN, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f0NN - fneq(q0NN)
      outState(?SAVE?( q0N1, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f0N1 - fneq(q0N1)
      outState(?SAVE?( q01N, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f01N - fneq(q01N)
      outState(?SAVE?( q011, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f011 - fneq(q011)
      outState(?SAVE?( qN0N, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fN0N - fneq(qN0N)
      outState(?SAVE?( q10N, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f10N - fneq(q10N)
      outState(?SAVE?( qN01, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fN01 - fneq(qN01)
      outState(?SAVE?( q101, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f101 - fneq(q101)
      outState(?SAVE?( qNN0, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fNN0 - fneq(qNN0)
      outState(?SAVE?( qN10, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fN10 - fneq(qN10)
      outState(?SAVE?( q1N0, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f1N0 - fneq(q1N0)
      outState(?SAVE?( q110, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f110 - fneq(q110)

      outState(?SAVE?( qNNN, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fNNN - fneq(qNNN)
      outState(?SAVE?( qNN1, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fNN1 - fneq(qNN1)
      outState(?SAVE?( qN1N, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fN1N - fneq(qN1N)
      outState(?SAVE?( qN11, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fN11 - fneq(qN11)
      outState(?SAVE?( q1NN, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f1NN - fneq(q1NN)
      outState(?SAVE?( q1N1, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f1N1 - fneq(q1N1)
      outState(?SAVE?( q11N, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f11N - fneq(q11N)
      outState(?SAVE?( q111, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f111 - fneq(q111)

      outState(?SAVE?( q000, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f000 - fneq(q000)

    enddo nodeloop

  end subroutine mrt_advRel_d3q27
! **************************************************************************** !


! **************************************************************************** !
  !> Semi-optimized explicit implementation for incompressible model
  !!
  !! 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_d3q27_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, iDir
    integer :: nScalars
?? copy :: var27( 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
    real(kind=rk) :: usq
    real(kind=rk) :: sum_e1, sum_e2, sum_e3
    real(kind=rk) :: sum_ux1, sum_ux2, sum_ux3
    real(kind=rk) :: sum_uy1, sum_uy2, sum_uy3
    real(kind=rk) :: sum_uz1, sum_uz2, sum_uz3
    ! 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) :: omegaBulk
    ! ---------------------------------------------------------------------------
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)
    nScalars = varSys%nScalars

    omegaBulk = fieldProp(1)%fluid%omegaBulkLvl(level)
    ! overwrite omegaKine term in the element loop
    s_mrt = mrt_d3q27(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
      fN00 = inState(?FETCH?( qN00, 1, iElem, QQ, nScalars, nElems,neigh ))
      f0N0 = inState(?FETCH?( q0N0, 1, iElem, QQ, nScalars, nElems,neigh ))
      f00N = inState(?FETCH?( q00N, 1, iElem, QQ, nScalars, nElems,neigh ))
      f100 = inState(?FETCH?( q100, 1, iElem, QQ, nScalars, nElems,neigh ))
      f010 = inState(?FETCH?( q010, 1, iElem, QQ, nScalars, nElems,neigh ))
      f001 = inState(?FETCH?( q001, 1, iElem, QQ, nScalars, nElems,neigh ))

      f0NN = inState(?FETCH?( q0NN, 1, iElem, QQ, nScalars, nElems,neigh ))
      f0N1 = inState(?FETCH?( q0N1, 1, iElem, QQ, nScalars, nElems,neigh ))
      f01N = inState(?FETCH?( q01N, 1, iElem, QQ, nScalars, nElems,neigh ))
      f011 = inState(?FETCH?( q011, 1, iElem, QQ, nScalars, nElems,neigh ))
      fN0N = inState(?FETCH?( qN0N, 1, iElem, QQ, nScalars, nElems,neigh ))
      f10N = inState(?FETCH?( q10N, 1, iElem, QQ, nScalars, nElems,neigh ))
      fN01 = inState(?FETCH?( qN01, 1, iElem, QQ, nScalars, nElems,neigh ))
      f101 = inState(?FETCH?( q101, 1, iElem, QQ, nScalars, nElems,neigh ))
      fNN0 = inState(?FETCH?( qNN0, 1, iElem, QQ, nScalars, nElems,neigh ))
      fN10 = inState(?FETCH?( qN10, 1, iElem, QQ, nScalars, nElems,neigh ))
      f1N0 = inState(?FETCH?( q1N0, 1, iElem, QQ, nScalars, nElems,neigh ))
      f110 = inState(?FETCH?( q110, 1, iElem, QQ, nScalars, nElems,neigh ))

      fNNN = inState(?FETCH?( qNNN, 1, iElem, QQ, nScalars, nElems,neigh ))
      fNN1 = inState(?FETCH?( qNN1, 1, iElem, QQ, nScalars, nElems,neigh ))
      fN1N = inState(?FETCH?( qN1N, 1, iElem, QQ, nScalars, nElems,neigh ))
      fN11 = inState(?FETCH?( qN11, 1, iElem, QQ, nScalars, nElems,neigh ))
      f1NN = inState(?FETCH?( q1NN, 1, iElem, QQ, nScalars, nElems,neigh ))
      f1N1 = inState(?FETCH?( q1N1, 1, iElem, QQ, nScalars, nElems,neigh ))
      f11N = inState(?FETCH?( q11N, 1, iElem, QQ, nScalars, nElems,neigh ))
      f111 = inState(?FETCH?( q111, 1, iElem, QQ, nScalars, nElems,neigh ))

      f000 = inState(?FETCH?( q000, 1, iElem, QQ, nScalars, nElems,neigh ))

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

      ! Zero moment
      mom(1) = rho

      ! First moments
      mom( 2) = rho*u_x
      mom( 3) = rho*u_y
      mom( 4) = rho*u_z

      sum_e1 = f001 + f00N + f010 + f0N0 + f100 + fN00
      sum_e3 = f111 + f11N + f1N1 + f1NN + fN11 + fN1N + fNN1 + fNNN
      sum_e2 = f011 + f01N + f0N1 + f0NN + f101 + f10N &
        &    + f110 + f1N0 + fN01 + fN0N + fN10 + fNN0
      ! kinematic energy
      mom( 5) = - 2.0_rk*f000 - sum_e1 + sum_e3
      ! square of energy
      mom(17) = 4.0_rk*f000 - sum_e2 + sum_e3
      ! cube of energy
      mom(18) = - 8.0_rk*f000 + 4.0_rk * sum_e1 - 2.0_rk * sum_e2 + sum_e3

      ! Second order tensors
      mom( 6) = 2.0_rk*( f100 + fN00 - f011 - f01N - f0N1 - f0NN)   &
        &     - f001 - f00N - f010 - f0N0                           &
        &     + f101 + f10N + f110 + f1N0 + fN01 + fN0N + fN10 + fNN0
      mom( 7) = - f001 - f00N + f010 + f0N0                            &
        &       - f101 - f10N + f110 + f1N0 - fN01 - fN0N + fN10 + fNN0
      mom( 8) = f110 + f111 + f11N - f1N0 - f1N1 - f1NN - fN10 - fN11 &
        &     - fN1N + fNN0 + fNN1 + fNNN
      mom( 9) = f011 - f01N - f0N1 + f0NN + f111 - f11N - f1N1 + f1NN &
        &     + fN11 - fN1N - fNN1 + fNNN
      mom(10) = f101 - f10N + f111 - f11N + f1N1 - f1NN - fN01 + fN0N &
        &     - fN11 + fN1N - fNN1 + fNNN

      ! fluxes of the energy and square of energy
      sum_ux1 = 4.0_rk*(f100 - fN00)
      sum_ux2 = - f101 - f10N - f110 - f1N0 +  fN01 + fN0N + fN10 + fNN0
      sum_ux3 = f111 + f11N + f1N1 + f1NN - fN11 - fN1N - fNN1 - fNNN
      mom(11) = - sum_ux1 + sum_ux2 + 2.0_rk * sum_ux3
      mom(14) =   sum_ux1 + 2.0_rk * sum_ux2 + sum_ux3

      sum_uy1 = 4.0_rk*(f010 - f0N0)
      sum_uy2 = - f011 - f01N  - f110 - fN10 + f0N1 + f0NN + fNN0 + f1N0
      sum_uy3 = f111 + f11N - f1N1 - f1NN + fN11 + fN1N - fNN1 - fNNN
      mom(12) = - sum_uy1 + sum_uy2 + 2.0_rk * sum_uy3
      mom(15) =   sum_uy1 + 2.0_rk * sum_uy2 + sum_uy3

      sum_uz1 = 4.0_rk*(f001 - f00N)
      sum_uz2 = - f011 - f0N1 - f101 - fN01 + f01N + f0NN + f10N + fN0N
      sum_uz3 = f111 - f11N + f1N1 - f1NN + fN11 - fN1N + fNN1 - fNNN
      mom(13) = - sum_uz1 + sum_uz2 + 2.0_rk * sum_uz3
      mom(16) =   sum_uz1 + 2.0_rk * sum_uz2 + sum_uz3

      !product of the second order tensors and the energy
      mom(19) = - 4.0_rk*(f100 + fN00)                                         &
        &     + 2.0_rk*( f001 + f00N + f010 - f011 - f01N + f0N0 - f0N1 - f0NN ) &
        &     + f101 + f10N + f110 + f1N0 + fN01 + fN0N + fN10 + fNN0
      mom(20) = 2.0_rk*(f001 + f00N - f010 - f0N0) &
        &     - f101 - f10N + f110 + f1N0 - fN01 - fN0N + fN10 + fNN0
      mom(21) = 2.0_rk*(-f110 + f1N0 + fN10 - fNN0) &
        &     + f111 + f11N - f1N1 - f1NN - fN11 - fN1N + fNN1 + fNNN
      mom(22) = 2.0_rk*(-f011 + f01N + f0N1 - f0NN ) &
        &     + f111 - f11N - f1N1 + f1NN + fN11 - fN1N - fNN1 + fNNN
      mom(23) = 2.0_rk*(-f101 + f10N + fN01 - fN0N ) &
        &     + f111 - f11N + f1N1 - f1NN - fN11 + fN1N - fNN1 + fNNN

      ! third order pseudo-vector
      mom(24) = -f101 - f10N + f110 + f1N0 + fN01 + fN0N - fN10 - fNN0
      mom(25) = f011 + f01N - f0N1 - f0NN - f110 + f1N0 - fN10 + fNN0
      mom(26) = -f011 + f01N - f0N1 + f0NN + f101 - f10N + fN01 - fN0N
      ! totally antisymmetric tensor XYZ
      mom(27) = f111 - f11N - f1N1 + f1NN - fN11 + fN1N + fNN1 - fNNN

      ! square of velocity
      usq = u_x*u_x + u_y*u_y + u_z*u_z
      ! -------------------------------------------------------------------------
      ! Equilibrium moments
      meq(1:QQ) =  0.0_rk
      meq( 1) = rho
      meq( 2) = rho0*u_x
      meq( 3) = rho0*u_y
      meq( 4) = rho0*u_z

      meq(14) = meq( 2)
      meq(15) = meq( 3)
      meq(16) = meq( 4)

      meq(11) = -2.0_rk*meq( 2)
      meq(12) = -2.0_rk*meq( 3)
      meq(13) = -2.0_rk*meq( 4)

      meq( 5) = rho0*usq - rho
      meq( 6) = rho0 * (2.0_rk*u_x*u_x - u_y*u_y - u_z*u_z)
      meq(19) = -meq( 6)

      meq( 7) = rho0 * (u_y*u_y - u_z*u_z)
      meq( 8) = rho0*u_x*u_y
      meq( 9) = rho0*u_y*u_z
      meq(10) = rho0*u_x*u_z

      meq(20) = -meq( 7)
      meq(21) = -meq( 8)
      meq(22) = -meq( 9)
      meq(23) = -meq(10)

      meq(17) = -2.0_rk*rho0*usq + rho
      meq(18) =  3.0_rk*rho0*usq - rho

      ! update kinematic omega
      s_mrt(6:10) = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem)
      ! compute neq moment
      mneq(:) = s_mrt(:) * ( mom(:) - meq(:) )

      ! compute fNeq
      do iDir = 1, QQ
        fneq(iDir) = sum( MMIvD3Q27(iDir,:) * mneq(:) )
      end do

      outState(?SAVE?( qN00, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fN00 - fneq(qN00)
      outState(?SAVE?( q0N0, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f0N0 - fneq(q0N0)
      outState(?SAVE?( q00N, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f00N - fneq(q00N)
      outState(?SAVE?( q100, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f100 - fneq(q100)
      outState(?SAVE?( q010, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f010 - fneq(q010)
      outState(?SAVE?( q001, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f001 - fneq(q001)

      outState(?SAVE?( q0NN, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f0NN - fneq(q0NN)
      outState(?SAVE?( q0N1, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f0N1 - fneq(q0N1)
      outState(?SAVE?( q01N, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f01N - fneq(q01N)
      outState(?SAVE?( q011, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f011 - fneq(q011)
      outState(?SAVE?( qN0N, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fN0N - fneq(qN0N)
      outState(?SAVE?( q10N, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f10N - fneq(q10N)
      outState(?SAVE?( qN01, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fN01 - fneq(qN01)
      outState(?SAVE?( q101, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f101 - fneq(q101)
      outState(?SAVE?( qNN0, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fNN0 - fneq(qNN0)
      outState(?SAVE?( qN10, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fN10 - fneq(qN10)
      outState(?SAVE?( q1N0, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f1N0 - fneq(q1N0)
      outState(?SAVE?( q110, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f110 - fneq(q110)

      outState(?SAVE?( qNNN, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fNNN - fneq(qNNN)
      outState(?SAVE?( qNN1, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fNN1 - fneq(qNN1)
      outState(?SAVE?( qN1N, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fN1N - fneq(qN1N)
      outState(?SAVE?( qN11, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & fN11 - fneq(qN11)
      outState(?SAVE?( q1NN, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f1NN - fneq(q1NN)
      outState(?SAVE?( q1N1, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f1N1 - fneq(q1N1)
      outState(?SAVE?( q11N, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f11N - fneq(q11N)
      outState(?SAVE?( q111, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f111 - fneq(q111)

      outState(?SAVE?( q000, 1, iElem, QQ, nScalars, nElems,neigh )) = &
        & f000 - fneq(q000)

    enddo nodeloop

  end subroutine mrt_advRel_d3q27_incomp
! **************************************************************************** !


! **************************************************************************** !
  !> 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_d3q27_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
    real(kind=rk) :: usq, ucx
    real(kind=rk) :: fEq( QQ ) !< equilibrium distribution
    ! MRT Variables
    real(kind=rk) :: s_mrt( QQ )
    real(kind=rk) :: mneq( QQ )
    real(kind=rk) :: fneq( QQ )
    real(kind=rk) :: omegaBulk
    integer :: dens_pos, vel_pos(3), elemOff
    ! ---------------------------------------------------------------------------
    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)
    ! collision matrix
    s_mrt = mrt_d3q27(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

      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 * ( 1.d0 + ucx*cs2inv        &
          &         + ucx*ucx*t2cs4inv - usq*t2cs2inv )

      enddo

      ! Calculate the non-equilibrium part
      fnEq(:) = pdfTmp(:) - fEq(:)

      ! convert non-equilibrium PDF into moments
      do iDir = 1, QQ
        mneq(iDir) =  sum( fnEq(:) * MMtrD3Q27(iDir,:) )
      end do

      ! viscosity omega
      s_mrt( 6:10 ) = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem)

      ! multiply it with collision matrix
      mneq(:) = s_mrt(:) * mneq(:)

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

    enddo nodeloop

  end subroutine mrt_advRel_d3q27_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_d3q27_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
    real(kind=rk) :: usq, ucx
    real(kind=rk) :: fEq( QQ ) !< equilibrium distribution
    ! MRT Variables
    real(kind=rk) :: s_mrt( QQ )
    real(kind=rk) :: mneq( QQ )
    real(kind=rk) :: fneq( QQ )
    real(kind=rk) :: omegaBulk
    integer :: dens_pos, vel_pos(3), elemOff
    ! ---------------------------------------------------------------------------
    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)
    ! collision matrix
    s_mrt = mrt_d3q27(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

      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

      ! Calculate the non-equilibrium part
      fnEq(:) = pdfTmp(:) - fEq(:)

      ! convert non-equilibrium PDF into moments
      do iDir = 1, QQ
        mneq(iDir) =  sum( fnEq(:) * MMtrD3Q27(iDir,:) )
      end do

      ! viscosity omega
      s_mrt( 6:10 ) = fieldProp(1)%fluid%viscKine%omLvl(level)%val(iElem)

      ! multiply it with collision matrix
      mneq(:) = s_mrt(:) * mneq(:)

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

    enddo nodeloop

  end subroutine mrt_advRel_d3q27_incomp_generic
! **************************************************************************** !

end module mus_mrt_d3q27_module


