! Copyright (c) 2015 Nikhil Anand <nikhil.anand@uni-siegen.de>
! Copyright (c) 2016-2017, 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016 Tobias Girresser <tobias.girresser@student.uni-siegen.de>
! Copyright (c) 2017 Daniel Petró <daniel.petro@student.uni-siegen.de>
!
! Permission to use, copy, modify, and distribute this software for any
! purpose with or without fee is hereby granted, provided that the above
! copyright notice and this permission notice appear in all copies.
!
! THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHORS DISCLAIM ALL WARRANTIES
! WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
! MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
! ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
! WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
! ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
! OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
! **************************************************************************** !

!> author: Jens Zudrop
!! Collects all functions related to the physical fluxes of the compressible Navier-Stokes equations.
module atl_physFluxFilNvrStk_module
  ! Treelm modules
  use env_module,            only: rk
  use tem_aux_module,        only: tem_abort
  use atl_eqn_nvrstk_module, only: atl_Navier_stokes_rans_type
  use atl_eqn_filNvrStk_var_module, &
    &                        only: atl_get_pointwise_velocity_gradient_2D,  &
    &                              atl_get_pointwise_visc_stress_tensor_2D, &
    &                              atl_get_lower_bound_turb_disscipation

  implicit none
  private


  public :: atl_viscPhysFluxRans
  public :: atl_PhysFluxRans
  public :: atl_PhysFluxRans_2d
  public :: atl_viscPhysFluxRans_2d,                      &
          & atl_mult_nu11_Rans_2d, atl_mult_nu21_Rans_2d, &
          & atl_mult_nu12_Rans_2d, atl_mult_nu22_Rans_2d

contains

  !> Physical flux calculation along x direction for the
  !  filtered Navier Stokes equation.
  function atl_physFluxRans(state, isenCoeff, penalty_char, porosity)       &
    &                       result(physFlux)
    ! ---------------------------------------------------------------------------
    !> The state in nodal space. Dimension is the number of vars, i.e. 5 for Euler
    !  and here =7 for renolds averaged navier stokes
    real(kind=rk), intent(in) :: state(:)
    !> Adiabatice index, also known as isentropic expansion factor.
    real(kind=rk), intent(in) :: isenCoeff
    !> The value of the characteristic function (stemming from penalization)
    real(kind=rk), intent(in) :: penalty_char
    !> The porosity at the current point
    real(kind=rk), intent(in) :: porosity
    !> The physical flux along the x axis for all variables
    real(kind=rk) :: physFlux(7)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: pressure, velocity(1:3), pressure_eff
    ! ---------------------------------------------------------------------------

    ! calculate pressure
    pressure = (isenCoeff-1.0_rk) * ( &
                 & state(5) - 0.5_rk*(sum(state(2:4)**2))/state(1) &
                           - state(6)           )
    pressure_eff = pressure + 2.0_rk * state(6) / 3.0_rk

    !> @todo JZ: here, we divide by a polynomial, we should be careful! We are leaving
    !! the polynomial space here!
    velocity(1:3) = state(2:4)/state(1)

    ! calculate the nonlinear term for different varibales now.
    ! ... density
    physFlux(1) = (1.0_rk + ((1.0_rk/porosity)-1.0_rk)*penalty_char) * state(2)
    ! ... x-velocity
    physFlux(2) = pressure_eff + state(2)*velocity(1)
    ! ... y-velocity
    physFlux(3) = state(2)*velocity(2)
    ! ... z-velocity
    physFlux(4) = state(2)*velocity(3)
    ! ... total energy
    physFlux(5) = velocity(1) * ( state(5) + pressure_eff )
    ! ... turbulent KE
    physFlux(6) = state(2)*state(6)/ state(1)
    ! ... spec dissipation rate
    physFlux(7) = state(2)*state(7)/ state(1)

  end function atl_physFluxRans


  function atl_physFluxRans_2d(state, isenCoeff, penalty_char, porosity)       &
    &                       result(physFlux)
    ! ---------------------------------------------------------------------------
    !> The state in nodal space. Dimension is the number of vars, i.e. 5 for Euler
    !  and here =7 for renolds averaged navier stokes
    real(kind=rk), intent(in) :: state(:)
    !> Adiabatice index, also known as isentropic expansion factor.
    real(kind=rk), intent(in) :: isenCoeff
    !> The value of the characteristic function (stemming from penalization)
    real(kind=rk), intent(in) :: penalty_char
    !> The porosity at the current point
    real(kind=rk), intent(in) :: porosity
    !> The physical flux along the x axis for all variables
    real(kind=rk) :: physFlux(6)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: pressure, velocity(1:2), pressure_eff
    ! ---------------------------------------------------------------------------

    ! calculate pressure
    pressure = (isenCoeff-1.0_rk) * ( &
                 & state(4) - 0.5_rk*(sum(state(2:3)**2))/state(1) &
                           - state(5)           )
    pressure_eff = pressure + 2.0_rk * state(5) / 3.0_rk

    !> @todo JZ: here, we divide by a polynomial, we should be careful! We are leaving
    !! the polynomial space here!
    velocity(1:2) = state(2:3)/state(1)

    ! calculate the nonlinear term for different varibales now.
    ! ... density
    physFlux(1) = (1.0_rk + ((1.0_rk/porosity)-1.0_rk)*penalty_char) * state(2)
    ! ... x-velocity
    physFlux(2) = pressure_eff + state(2)*velocity(1)
    ! ... y-velocity
    physFlux(3) = state(2)*velocity(2)
    ! ... total energy
    physFlux(4) = velocity(1) * ( state(4) + pressure_eff )
    ! ... turbulent KE
    physFlux(5) = velocity(1)*state(5)
    ! ... spec dissipation rate
    physFlux(6) = velocity(1)*state(6)

  end function atl_physFluxRans_2d


  !> Physical flux calculation along x direction for Euler equation.
  function atl_viscPhysFluxRans(state, state_gradient, isenCoeff, mu, &
                                   & lambda, thermCond, heatCap ) result(physFlux)
    ! ------------------------------------------------------------------------------------
    !> The state in nodal space. Dimension is the number of vars, i.e. 5 for Navier-Stokes.
    real(kind=rk), intent(in) :: state(:)
    !> The state in nodal space. First dimension is the number of vars, i.e. 5 for Navier-Stokes.
    !! Second dimension is the dimension, e.g. 3 in two dimensions.
    real(kind=rk), intent(in) :: state_gradient(:,:)
    !> Adiabatice index, also known as isentropic expansion factor.
    real(kind=rk), intent(in) :: isenCoeff
    !> Dynamic Viscosity
    real(kind=rk), intent(in) :: mu
    !> Viscosity
    real(kind=rk), intent(in) :: lambda
    !> The thermal cond
    real(kind=rk), intent(in) :: thermCond
    !> The specific heat capacity (per mass unit mass, at constant volume)
    real(kind=rk), intent(in) :: heatCap
    !> The physical flux along the x axis for all variables
    real(kind=rk) :: physFlux(7)
    ! -------------------------------------------------------------------------------
    real(kind=rk) :: velocity(1:3)
    ! The turbulent viscosity (dynamic)
    real(kind=rk) :: mu_turb
    ! The turbulent viscosity
    real(kind=rk) :: lam_turb
    ! The effective viscosity (dynamic)
    real(kind=rk) :: mu_eff
    ! The effective turbulent viscosity
    real(kind=rk) :: lam_eff
    real(kind=rk) :: energy_coeff
    real(kind=rk) :: sig_k
    real(kind=rk) :: sig_omg
    real(kind=rk) :: turbulent_prandtl_number
    ! -------------------------------------------------------------------------------

    velocity(1:3) = state(2:4)/state(1)
    !mu_turb = rans_params%c_mu*state(5)/exp(state(6)/state(1))

    ! The code above was deactivated without adding a proper initialization
    ! for the variables. I'm letting the code  abort when it is executed to
    ! inform the user that it is not yet usable.
    call tem_abort( "This routine doesn't use properly initialized variables" )

    lam_turb = 2.0*mu_turb/3.0
    mu_eff = mu + mu_turb
    lam_eff = lambda + lam_turb
    turbulent_prandtl_number = 0.9_rk
    sig_k = 0.5_rk
    sig_omg = 0.5_rk
    energy_coeff = thermCond/(heatCap*state(1))      &
    &                 + mu_turb*isenCoeff/(turbulent_prandtl_number*state(1))

    ! Viscous flux for density
    physFlux(1) = 0.0_rk

    ! Viscous flux for momentum in x
    physFlux(2) =                                                                 &
      ! (nu_{1,1})_{2,1} (\nabla u)_{1,1} : k = 1, i = 1
      & ((-2.0_rk * mu_eff + lam_eff) / state(1))*velocity(1)*state_gradient(1,1) &
      ! (nu_{1,2})_{2,1} (\nabla u)_{1,2} : k = 2, i = 1
      & + (lam_eff / state(1))*velocity(2)*state_gradient(1,2)                    &
      ! (nu_{1,3})_{2,1} (\nabla u)_{1,3} : k = 3, i = 1
      & + (lam_eff / state(1))*velocity(3)*state_gradient(1,3)                    &
      ! (nu_{1,1})_{2,2} (\nabla u)_{2,1} : k = 1, i = 2
      & + ((2.0_rk * mu_eff - lam_eff) / state(1))*state_gradient(2,1)            &
      ! (nu_{1,2})_{2,2} (\nabla u)_{2,2} : k = 2, i = 2
      & + 0.0_rk                                                                  &
      !(nu_{1,3})_{2,2} (\nabla u)_{2,3} : k = 3, i = 2
      & + 0.0_rk                                                                  &
      ! (nu_{1,1})_{2,3} (\nabla u)_{3,1} : k = 1, i = 3
      & + 0.0_rk                                                                  &
      ! (nu_{1,2})_{2,3} (\nabla u)_{3,2} : k = 2, i = 3
      & + (-lam_eff/state(1))*state_gradient(3,2)                                 &
      ! (nu_{1,3})_{2,3} (\nabla u)_{3,3} : k = 2, i = 3
      & + 0.0_rk                                                                  &
      ! (nu_{1,1})_{2,4} (\nabla u)_{4,1} : k = 1, i = 4
      & + 0.0_rk                                                                  &
      ! (nu_{1,2})_{2,4} (\nabla u)_{4,2} : k = 2, i = 4
      & + 0.0_rk                                                                  &
      ! (nu_{1,3})_{2,4} (\nabla u)_{4,3} : k = 3, i = 4
      & + (-lam_eff/state(1))*state_gradient(4,3)                                 &
      ! (nu_{1,1})_{2,5} (\nabla u)_{5,1} : k = 1, i = 5
      & + 0.0_rk                                                                  &
      ! (nu_{1,2})_{2,5} (\nabla u)_{5,2} : k = 2, i = 5
      & + 0.0_rk                                                                  &
      ! (nu_{1,3})_{2,5} (\nabla u)_{5,3} : k = 3, i = 5
      & + 0.0_rk                                                                  &
      ! (nu_{1,1})_{2,6} (\nabla u)_{6,1} : k = 1, i = 6
      & + 0.0_rk                                                                  &
      ! (nu_{1,2})_{2,6} (\nabla u)_{6,2} : k = 2, i = 6
      & + 0.0_rk                                                                  &
      ! (nu_{1,3})_{2,6} (\nabla u)_{6,3} : k = 3, i = 6
      & + 0.0_rk                                                                  &
      ! (nu_{1,1})_{2,7} (\nabla u)_{7,1} : k = 1, i = 7
      & + 0.0_rk                                                                  &
      ! (nu_{1,2})_{2,7} (\nabla u)_{7,2} : k = 2, i = 7
      & + 0.0_rk                                                                  &
      ! (nu_{1,3})_{2,7} (\nabla u)_{7,3} : k = 3, i = 7
      & + 0.0_rk

    ! Viscous flux for momentum in y
    physFlux(3) =                                                              &
      ! (nu_{1,1})_{3,1} (\nabla u)_{1,1} : k = 1, i = 1
      &   (-mu_eff / state(1))*velocity(2)*state_gradient(1,1)                 &
      ! (nu_{1,2})_{3,1} (\nabla u)_{1,2} : k = 2, i = 1
      & + (-mu_eff / state(1))*velocity(1)*state_gradient(1,2)                 &
      ! (nu_{1,3})_{3,1} (\nabla u)_{1,3} : k = 3, i = 1
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{3,2} (\nabla u)_{2,1} : k = 1, i = 2
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{3,2} (\nabla u)_{2,2} : k = 2, i = 2
      & + (mu_eff/state(1))*state_gradient(2,2)                                &
      ! (nu_{1,3})_{3,2} (\nabla u)_{2,3} : k = 3, i = 2
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{3,3} (\nabla u)_{3,1} : k = 1, i = 3
      & + (mu_eff/state(1))*state_gradient(3,1)                                &
      ! (nu_{1,2})_{3,3} (\nabla u)_{3,2} : k = 2, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{3,3} (\nabla u)_{3,3} : k = 3, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{3,4} (\nabla u)_{4,1} : k = 1, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{3,4} (\nabla u)_{4,2} : k = 2, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{3,4} (\nabla u)_{4,3} : k = 3, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{3,5} (\nabla u)_{5,1} : k = 1, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{3,5} (\nabla u)_{5,2} : k = 2, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{3,5} (\nabla u)_{5,3} : k = 3, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{3,6} (\nabla u)_{6,1} : k = 1, i = 6
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{3,6} (\nabla u)_{6,2} : k = 2, i = 6
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{3,6} (\nabla u)_{6,3} : k = 3, i = 6
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{3,7} (\nabla u)_{7,1} : k = 1, i = 7
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{3,7} (\nabla u)_{7,2} : k = 2, i = 7
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{3,7} (\nabla u)_{7,3} : k = 3, i = 7
      & + 0.0_rk


    ! Viscous flux for momentum in z
    physFlux(4) =                                                              &
      ! (nu_{1,1})_{4,1} (\nabla u)_{1,1} : k = 1, i = 1
      &   (-mu_eff / state(1))*velocity(3)*state_gradient(1,1)                 &
      ! (nu_{1,2})_{4,1} (\nabla u)_{1,2} : k = 2, i = 1
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{4,1} (\nabla u)_{1,3} : k = 3, i = 1
      & + (-mu_eff / state(1))*velocity(1)*state_gradient(1,3)                 &
      ! (nu_{1,1})_{4,2} (\nabla u)_{2,1} : k = 1, i = 2
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{4,2} (\nabla u)_{2,2} : k = 2, i = 2
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{4,2} (\nabla u)_{2,3} : k = 3, i = 2
      & + (mu_eff/state(1))*state_gradient(2,3)                                &
      ! (nu_{1,1})_{4,3} (\nabla u)_{3,1} : k = 1, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{4,3} (\nabla u)_{3,2} : k = 2, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{4,3} (\nabla u)_{3,3} : k = 3, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{4,4} (\nabla u)_{4,1} : k = 1, i = 4
      & + (mu_eff/state(1))*state_gradient(4,1)                                &
      ! (nu_{1,2})_{4,4} (\nabla u)_{4,2} : k = 2, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{4,4} (\nabla u)_{4,3} : k = 3, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{4,5} (\nabla u)_{5,1} : k = 1, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{4,5} (\nabla u)_{5,2} : k = 2, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{4,5} (\nabla u)_{5,3} : k = 3, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{4,6} (\nabla u)_{6,1} : k = 1, i = 6
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{4,6} (\nabla u)_{6,2} : k = 2, i = 6
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{4,6} (\nabla u)_{6,3} : k = 3, i = 6
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{4,7} (\nabla u)_{7,1} : k = 1, i = 7
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{4,7} (\nabla u)_{7,2} : k = 2, i = 7
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{4,7} (\nabla u)_{7,3} : k = 3, i = 7
      & + 0.0_rk

    ! Viscous flux for Energy
    physFlux(5) = &
      ! (nu_{1,1})_{5,1} (\nabla u)_{1,1} : k = 1, i = 1
      & ( ((-2.0_rk * mu_eff + lam_eff) / state(1))*(velocity(1)**2.0_rk)               &
      & +(-mu_eff/state(1))*(velocity(2)**2.0_rk)                                       &
      & +(-mu_eff/state(1))*(velocity(3)**2.0_rk)                                       &
      & -energy_coeff*(state(5)/state(1)-sum(velocity(:)**2.0_rk) - state(6)/state(1))) &
      &  * state_gradient(1,1)                                                          &
      ! (nu_{1,2})_{5,1} (\nabla u)_{1,2} : k = 2, i = 1
      & + ((-mu_eff+lam_eff)/state(1))*velocity(1)*velocity(2)*state_gradient(1,2)      &
      ! (nu_{1,3})_{5,1} (\nabla u)_{1,3} : k = 3, i = 1
      & + ((-mu_eff+lam_eff)/state(1))*velocity(1)*velocity(3)*state_gradient(1,3)      &
      ! (nu_{1,1})_{5,2} (\nabla u)_{2,1} : k = 1, i = 2
      & + ((2.0_rk*mu_eff-lam_eff)/state(1) - energy_coeff)*velocity(1)                 &
      &  * state_gradient(2,1)                                                          &
      ! (nu_{1,2})_{5,2} (\nabla u)_{2,2} : k = 2, i = 2
      & + (mu_eff/state(1))*velocity(2)*state_gradient(2,2)                             &
      ! (nu_{1,3})_{5,2} (\nabla u)_{2,3} : k = 3, i = 2
      & + (mu_eff/state(1))*velocity(3)*state_gradient(2,3)                             &
      ! (nu_{1,1})_{5,3} (\nabla u)_{3,1} : k = 1, i = 3
      & + ((mu_eff/state(1))- energy_coeff)*velocity(2)*state_gradient(3,1)             &
      ! (nu_{1,2})_{5,3} (\nabla u)_{3,2} : k = 2, i = 3
      & + (-lam_eff/state(1))*velocity(1)*state_gradient(3,2)                           &
      ! (nu_{1,3})_{5,3} (\nabla u)_{3,3} : k = 3, i = 3
      & + 0.0_rk                                                                        &
      ! (nu_{1,1})_{5,4} (\nabla u)_{4,1} : k = 1, i = 4
      & + ((mu_eff/state(1))- energy_coeff)*velocity(3)*state_gradient(4,1)             &
      ! (nu_{1,2})_{5,4} (\nabla u)_{4,2} : k = 2, i = 4
      & + 0.0_rk                                                                        &
      ! (nu_{1,3})_{5,4} (\nabla u)_{4,3} : k = 3, i = 4
      & + (-lam_eff/state(1))*velocity(1)*state_gradient(4,3)                           &
      ! (nu_{1,1})_{5,5} (\nabla u)_{5,1} : k = 1, i = 5
      & + energy_coeff*state_gradient(5,1)                                              &
      ! (nu_{1,2})_{5,5} (\nabla u)_{5,2} : k = 2, i = 5
      & + 0.0_rk                                                                        &
      ! (nu_{1,3})_{5,5} (\nabla u)_{5,3} : k = 3, i = 5
      & + 0.0_rk                                                                        &
      ! (nu_{1,1})_{5,6} (\nabla u)_{6,1} : k = 1, i = 6
      & + ((mu + sig_k*mu_turb)/state(1) - energy_coeff)*state_gradient(6,1)            &
      ! (nu_{1,2})_{5,6} (\nabla u)_{6,2} : k = 2, i = 6
      & + 0.0_rk                                                                        &
      ! (nu_{1,3})_{5,6} (\nabla u)_{6,3} : k = 3, i = 6
      & + 0.0_rk                                                                        &
      ! (nu_{1,1})_{5,7} (\nabla u)_{7,1} : k = 1, i = 7
      & + 0.0_rk                                                                        &
      ! (nu_{1,2})_{5,7} (\nabla u)_{7,2} : k = 2, i = 7
      & + 0.0_rk                                                                        &
      ! (nu_{1,3})_{5,7} (\nabla u)_{7,3} : k = 3, i = 7
      & + 0.0_rk


    ! Viscous flux for turbulent KE
    physFlux(6) =                                                              &
      ! (nu_{1,1})_{6,1} (\nabla u)_{1,1} : k = 1, i = 1
      & ( -(mu + sig_k*mu_turb)/state(1))*(state(6)/state(1))                  &
      &                                  *state_gradient(1,1)                  &
      ! (nu_{1,2})_{6,1} (\nabla u)_{1,2} : k = 2, i = 1
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{6,1} (\nabla u)_{1,3} : k = 3, i = 1
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{6,2} (\nabla u)_{2,1} : k = 1, i = 2
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{6,2} (\nabla u)_{2,2} : k = 2, i = 2
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{6,2} (\nabla u)_{2,3} : k = 3, i = 2
      & + (mu_eff/state(1))*state_gradient(2,3)                                &
      ! (nu_{1,1})_{6,3} (\nabla u)_{3,1} : k = 1, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{6,3} (\nabla u)_{3,2} : k = 2, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{6,3} (\nabla u)_{3,3} : k = 3, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{6,4} (\nabla u)_{4,1} : k = 1, i = 4
      & + (mu_eff/state(1))*state_gradient(4,1)                                &
      ! (nu_{1,2})_{6,4} (\nabla u)_{4,2} : k = 2, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{6,4} (\nabla u)_{4,3} : k = 3, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{6,5} (\nabla u)_{5,1} : k = 1, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{6,5} (\nabla u)_{5,2} : k = 2, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{6,5} (\nabla u)_{5,3} : k = 3, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{6,6} (\nabla u)_{6,1} : k = 1, i = 6
      & + (mu + sig_k*mu_turb)/state(1)*state_gradient(6,1)                    &
      ! (nu_{1,2})_{6,6} (\nabla u)_{6,2} : k = 2, i = 6
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{6,6} (\nabla u)_{6,3} : k = 3, i = 6
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{6,7} (\nabla u)_{7,1} : k = 1, i = 7
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{6,7} (\nabla u)_{7,2} : k = 2, i = 7
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{6,7} (\nabla u)_{7,3} : k = 3, i = 7
      & + 0.0_rk

    ! Viscous flux for specific Dissipation rate (\Omega)
    physFlux(7) =                                                              &
      ! (nu_{1,1})_{7,1} (\nabla u)_{1,1} : k = 1, i = 1
      & ( -(mu + sig_k*mu_turb)/state(1))*(state(7)/state(1))                  &
      &                                  *state_gradient(1,1)                  &
      ! (nu_{1,2})_{7,1} (\nabla u)_{1,2} : k = 2, i = 1
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{7,1} (\nabla u)_{1,3} : k = 3, i = 1
      & + 0.0_rk                &
      ! (nu_{1,1})_{7,2} (\nabla u)_{2,1} : k = 1, i = 2
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{7,2} (\nabla u)_{2,2} : k = 2, i = 2
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{7,2} (\nabla u)_{2,3} : k = 3, i = 2
      & + (mu_eff/state(1))*state_gradient(2,3)                                &
      ! (nu_{1,1})_{7,3} (\nabla u)_{3,1} : k = 1, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{7,3} (\nabla u)_{3,2} : k = 2, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{7,3} (\nabla u)_{3,3} : k = 3, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{7,4} (\nabla u)_{4,1} : k = 1, i = 4
      & + (mu_eff/state(1))*state_gradient(4,1)                                &
      ! (nu_{1,2})_{7,4} (\nabla u)_{4,2} : k = 2, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{7,4} (\nabla u)_{4,3} : k = 3, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{7,5} (\nabla u)_{5,1} : k = 1, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{7,5} (\nabla u)_{5,2} : k = 2, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{7,5} (\nabla u)_{5,3} : k = 3, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{7,6} (\nabla u)_{6,1} : k = 1, i = 6
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{7,6} (\nabla u)_{6,2} : k = 2, i = 6
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{7,6} (\nabla u)_{6,3} : k = 3, i = 6
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{7,7} (\nabla u)_{7,1} : k = 1, i = 7
      & + ((mu + sig_omg*mu_turb)/state(1))*state_gradient(7,1)                &
      ! (nu_{1,2})_{7,7} (\nabla u)_{7,2} : k = 2, i = 7
      & + 0.0_rk                                                               &
      ! (nu_{1,3})_{7,7} (\nabla u)_{7,3} : k = 3, i = 7
      & + 0.0_rk


  end function atl_viscPhysFluxRans


  function atl_viscPhysFluxRans_2d(state, state_gradient, isenCoeff, mu, &
                                   & lambda, thermCond, rans_params, heatCap ) result(physFlux)
    ! ------------------------------------------------------------------------------------
    !> The state in nodal space. Dimension is the number of vars, i.e. 4 for Navier-Stokes.
    real(kind=rk), intent(in) :: state(:)
    !> The state in nodal space. First dimension is the number of vars, i.e. 4 for Navier-Stokes.
    !! Second dimension is the dimension, e.g. 2 in two dimensions.
    real(kind=rk), intent(in) :: state_gradient(:,:)
    !> Adiabatice index, also known as isentropic expansion factor.
    real(kind=rk), intent(in) :: isenCoeff
    !> Dynamic Viscosity
    real(kind=rk), intent(in) :: mu
    !> Viscosity
    real(kind=rk), intent(in) :: lambda
    !> The thermal cond
    real(kind=rk), intent(in) :: thermCond
    !> The specific heat capacity (per mass unit mass, at constant volume)
    real(kind=rk), intent(in) :: heatCap
    !> The physical flux along the x axis for all variables
    real(kind=rk) :: physFlux(6)
    !> The constants for the Rans eqn
    type(atl_Navier_stokes_rans_type), intent(in) :: rans_params
    ! -----------------------------------------------------------------------------
    real(kind=rk) :: velocity(1:2)
    ! The turbulent viscosity (dynamic)
    real(kind=rk) :: mu_turb
    ! The turbulent viscosity
    real(kind=rk) :: lam_turb
    ! The effective viscosity (dynamic)
    real(kind=rk) :: mu_eff
    ! The effective turbulent viscosity
    real(kind=rk) :: lam_eff
    real(kind=rk) :: energy_coeff
    real(kind=rk) :: k_eff
    real(kind=rk) :: turbPrNum, limited_eddy_visc, turb_coeff1, turb_coeff2
    real(kind=rk) :: velGrad(2,2), ViscStressTensor(2,2),  omega_r
    ! -----------------------------------------------------------------------------
    mu_turb = rans_params%c_mu*state(5)/exp(state(6)/state(1))

    lam_turb = 2.0*mu_turb/3.0
    mu_eff = mu + mu_turb
    lam_eff = lambda + lam_turb
    turbPrNum = rans_params%turb_prandtl_num

    k_eff = max(state(6),0.0_rk)
    energy_coeff = thermCond/(heatCap*state(1))               &
    &                 + mu_turb*isenCoeff/(turbPrNum*state(1) )
    velocity(1:2) = state(2:3)/state(1)

    ! Now to calculate the coefficient \omega_r: We need to do the following
    ! Step-1: Get the velocity gradient from the state gradient present
    call atl_get_pointwise_velocity_gradient_2D ( state_gradient, &
     &                           state, velGrad                   )

    ! Step-2: Calculate the viscous stress tensor
    call atl_get_pointwise_visc_stress_tensor_2D(         &
      &          velGrad = velGrad,                       &
      &          S       = ViscStressTensor               )

    ! Step-3: Get the value of \omega_r
    call atl_get_lower_bound_turb_disscipation(           &
      &          S       = ViscStressTensor,              &
      &        c_mu      = rans_params%c_mu,              &
      &        omega     = state(6)/ state(1),            &
      &        omega_r   = omega_r                        )


    limited_eddy_visc = rans_params%alpha*k_eff*exp(-omega_r)
    turb_coeff1 = (mu + rans_params%sig_k*limited_eddy_visc)/ state(1)
    turb_coeff2 = (mu + rans_params%sig_omg*limited_eddy_visc)/ state(1)

    ! Viscous flux for density
    physFlux(1) = 0.0_rk

    ! Viscous flux for momentum in x
    physFlux(2) =                                                               &
      ! (nu_{1,1})_{2,1} (\nabla u)_{1,1} : k = 1, i = 1
      & ((-2.0_rk * mu_eff + lam_eff) / state(1))*velocity(1)*state_gradient(1,1) &
      ! (nu_{1,2})_{2,1} (\nabla u)_{1,2} : k = 2, i = 1
      & + (lam_eff / state(1))*velocity(2)*state_gradient(1,2)                  &
      ! (nu_{1,1})_{2,2} (\nabla u)_{2,1} : k = 1, i = 2
      & + ((2.0_rk * mu_eff - lam_eff) / state(1))*state_gradient(2,1)         &
      ! (nu_{1,2})_{2,2} (\nabla u)_{2,2} : k = 2, i = 2
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{2,3} (\nabla u)_{3,1} : k = 1, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{2,3} (\nabla u)_{3,2} : k = 2, i = 3
      & + (-lam_eff/state(1))*state_gradient(3,2)                              &
      ! (nu_{1,1})_{2,4} (\nabla u)_{4,1} : k = 1, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{2,4} (\nabla u)_{4,2} : k = 2, i = 4
      & + 0.0_rk

    ! Viscous flux for momentum in y
    physFlux(3) =                                                              &
      ! (nu_{1,1})_{3,1} (\nabla u)_{1,1} : k = 1, i = 1
      &   (-mu_eff / state(1))*velocity(2)*state_gradient(1,1)                 &
      ! (nu_{1,2})_{3,1} (\nabla u)_{1,2} : k = 2, i = 1
      & + (-mu_eff / state(1))*velocity(1)*state_gradient(1,2)                 &
      ! (nu_{1,1})_{3,2} (\nabla u)_{2,1} : k = 1, i = 2
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{3,2} (\nabla u)_{2,2} : k = 2, i = 2
      & + (mu_eff/state(1))*state_gradient(2,2)                                &
      ! (nu_{1,1})_{3,3} (\nabla u)_{3,1} : k = 1, i = 3
      & + (mu_eff/state(1))*state_gradient(3,1)                                &
      ! (nu_{1,2})_{3,3} (\nabla u)_{3,2} : k = 2, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{3,4} (\nabla u)_{4,1} : k = 1, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{3,4} (\nabla u)_{4,2} : k = 2, i = 4
      & + 0.0_rk

    ! Viscous flux for Energy
    physFlux(4) =                                                                   &
     ! (nu_{1,1})_{4,1} (\nabla u)_{1,1} : k = 1, i = 1
     & ( ((-2.0_rk * mu_eff + lam_eff) / state(1))*(velocity(1)**2.0_rk)            &
     & + (-mu_eff/state(1))*(velocity(2)**2.0_rk)                                   &
     & -energy_coeff*(state(4)/state(1)-sum(velocity(:)**2.0_rk) - state(5)/state(1)) &
     & - turb_coeff1*(state(5)/state(1)) )                                          &
     &  * state_gradient(1,1)                                                       &
     ! (nu_{1,2})_{4,1} (\nabla u)_{1,2} : k = 2, i = 1
     & + ((-mu_eff+lam_eff)/state(1))*velocity(1)*velocity(2)*state_gradient(1,2)   &
     ! (nu_{1,1})_{4,2} (\nabla u)_{2,1} : k = 1, i = 2
     & + ((2.0_rk*mu_eff-lam_eff)/state(1) - energy_coeff)*velocity(1)              &
     &  * state_gradient(2,1)                                                       &
     ! (nu_{1,2})_{4,2} (\nabla u)_{2,2} : k = 2, i = 2
     & + (mu_eff/state(1))*velocity(2)*state_gradient(2,2)                          &
     ! (nu_{1,1})_{4,3} (\nabla u)_{3,1} : k = 1, i = 3
     & + ((mu_eff/state(1)) - energy_coeff)*velocity(2)*state_gradient(3,1)    &
     ! (nu_{1,2})_{4,3} (\nabla u)_{3,2} : k = 2, i = 3
     & + (-lam_eff/state(1))*velocity(1)*state_gradient(3,2)                  &
     ! (nu_{1,1})_{4,4} (\nabla u)_{4,1} : k = 1, i = 4
     & + energy_coeff*state_gradient(4,1)                                     &
     ! (nu_{1,2})_{4,4} (\nabla u)_{4,2} : k = 2, i = 4
     & + 0.0_rk                                                               &
     ! (nu_{1,1})_{4,5} (\nabla u)_{5,1} : k = 1, i = 5
     & + (-energy_coeff + turb_coeff1)*state_gradient(5,1)                    &
     ! (nu_{1,2})_{4,5} (\nabla u)_{5,2} : k = 2, i = 5
     & + 0.0_rk                                                               &
     ! (nu_{1,1})_{4,6} (\nabla u)_{6,1} : k = 1, i = 6
     & + 0.0_rk                                                               &
     ! (nu_{1,2})_{4,6} (\nabla u)_{6,2} : k = 2, i = 6
     & + 0.0_rk

    physFlux(5) =                                                              &
      ! (nu_{1,1})_{5,1} (\nabla u)_{1,1} : k = 1, i = 1
      & -turb_coeff1*(state(5)/state(1)) *state_gradient(1,1)                  &
      ! (nu_{1,2})_{5,1} (\nabla u)_{1,2} : k = 2, i = 1
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{5,2} (\nabla u)_{2,1} : k = 1, i = 2
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{5,2} (\nabla u)_{2,2} : k = 2, i = 2
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{5,3} (\nabla u)_{3,1} : k = 1, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{5,3} (\nabla u)_{3,2} : k = 2, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{5,4} (\nabla u)_{4,1} : k = 1, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{5,4} (\nabla u)_{4,2} : k = 2, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{5,5} (\nabla u)_{5,1} : k = 1, i = 5
      & + turb_coeff1*state_gradient(5,1)                                      &
      ! (nu_{1,2})_{5,5} (\nabla u)_{5,2} : k = 2, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{5,6} (\nabla u)_{6,1} : k = 1, i = 6
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{5,6} (\nabla u)_{6,2} : k = 2, i = 6
      & + 0.0_rk

    physFlux(6) =                                                              &
      ! (nu_{1,1})_{6,1} (\nabla u)_{1,1} : k = 1, i = 1
      & ( -turb_coeff2*state(6)/state(1) )*state_gradient(1,1)                 &
      ! (nu_{1,2})_{6,1} (\nabla u)_{1,2} : k = 2, i = 1
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{6,2} (\nabla u)_{2,1} : k = 1, i = 2
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{6,2} (\nabla u)_{2,2} : k = 2, i = 2
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{6,3} (\nabla u)_{3,1} : k = 1, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{6,3} (\nabla u)_{3,2} : k = 2, i = 3
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{6,4} (\nabla u)_{4,1} : k = 1, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{6,4} (\nabla u)_{4,2} : k = 2, i = 4
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{6,5} (\nabla u)_{5,1} : k = 1, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,2})_{6,5} (\nabla u)_{5,2} : k = 2, i = 5
      & + 0.0_rk                                                               &
      ! (nu_{1,1})_{6,6} (\nabla u)_{6,1} : k = 1, i = 6
      & + turb_coeff2*state_gradient(6,1)                                      &
      ! (nu_{1,2})_{6,6} (\nabla u)_{6,2} : k = 2, i = 6
      & + 0.0_rk

  end function atl_viscPhysFluxRans_2d


  ! Multiplies the viscous flux matrux nu_11 with a given vector
  function atl_mult_nu11_Rans_2d( state, velocity, inVec,                     &
    &                isenCoeff, mu, lambda, thermCond, rans_params, heatCap)  &
    & result( outVec )
    ! -------------------------------------------------------------------------!
    !> The velocity
    real(kind=rk), intent(in) :: velocity(2)
    !> The state array
    real(kind=rk), intent(in) :: state(6)
    !> Vector to be multiplied with nu11
    real(kind=rk), intent(in) :: inVec(6)
    !> Adiabatice index, also known as isentropic expansion factor.
    real(kind=rk), intent(in) :: isenCoeff
    !> Dynamic Viscosity
    real(kind=rk), intent(in) :: mu
    !> Viscosity
    real(kind=rk), intent(in) :: lambda
    !> The thermal cond
    real(kind=rk), intent(in) :: thermCond
    !> The specific heat capacity (per mass unit mass, at constant volume)
    real(kind=rk), intent(in) :: heatCap
    !> The constants for the Rans eqn
    type(atl_Navier_stokes_rans_type), intent(in) :: rans_params
    !> The result of the matrix vector product
    real(kind=rk) :: outVec(6)
    ! -------------------------------------------------------------------------!
    !> mu_turb
    real(kind=rk) :: mu_turb, mu_eff, lam_eff, limited_eddy_visc, energy_coeff
    real(kind=rk) :: turb_coeff1, turb_coeff2
    ! -------------------------------------------------------------------------!
    mu_turb = rans_params%c_mu*state(5)/exp(state(6)/state(1))

    mu_eff = mu + mu_turb
    lam_eff = lambda + 2.0*mu_turb/3.0

    ! Stuff neeeded
    energy_coeff = thermCond/(heatCap*state(1))                               &
    &             + mu_turb*isenCoeff/(rans_params%turb_prandtl_num*state(1))

!NA!    ! Now to calculate the coefficient \omega_r: We need to do the following
!NA!    ! Step-1: Get the velocity gradient from the state gradient present
!NA!    call atl_get_pointwise_velocity_gradient_2D ( inVec, &
!NA!     &                           state, velGrad               )
!NA!
!NA!    ! Step-2: Calculate the viscous stress tensor
!NA!    call atl_get_pointwise_visc_stress_tensor_2D(             &
!NA!      &          velGrad = velGrad,                       &
!NA!      &          S       = ViscStressTensor               )
!NA!
!NA!    ! Step-3: Get the value of \omega_r
!NA!    call atl_get_lower_bound_turb_disscipation(                          &
!NA!      &          S       = ViscStressTensor,                         &
!NA!      &        c_mu      = rans_params%c_mu,                         &
!NA!      &        omega     = state(6)/ state(1),                       &
!NA!      &        omega_r   = omega_r                                   )


    ! @todo : NA : Calculate limited_eddy_visc Appropriately
!NA!    limited_eddy_visc = rans_params%alpha*max(state(5),0.0)*exp(-omega_r)
    limited_eddy_visc = rans_params%alpha*max(state(5),0.0)              &
      &                                  *exp(state(6)/state(1))
    turb_coeff1 = (mu + rans_params%sig_k*limited_eddy_visc)/ state(1)
    turb_coeff2 = (mu + rans_params%sig_omg*limited_eddy_visc)/ state(1)

    ! First row has zeros only
    outVec(1) = 0.0_rk

    ! Second row
    outVec(2) = ((-2.0_rk*mu_eff + lam_eff)/state(1))*velocity(1) * inVec(1)   &
            & + ((2.0_rk*mu_eff - lam_eff)/state(1)) * inVec(2)

    ! Third row
    outVec(3) = (-mu_eff/state(1)) * velocity(2) * inVec(1)                    &
            & + (mu_eff/state(1)) * inVec(3)

    ! Fourth row
    outVec(4) = (                                                              &
      &         ((-2.0_rk*mu_eff + lam_eff)/state(1))*velocity(1)*velocity(1)  &
      &           + (-mu_eff/state(1)) * velocity(2)*velocity(2)               &
      &           - energy_coeff*(state(4)/state(1) - sum(velocity(:)**2)      &
      &                                      - state(5)/state(1) )             &
      &           - turb_coeff1*(state(5)/state(1))                )* inVec(1) &
      & + (                                                                    &
      &     ((2.0_rk*mu_eff-lam_eff)/state(1) - energy_coeff)*velocity(1)      &
      &   ) * inVec(2)                                                         &
      & + (  (mu_eff/state(1) - energy_coeff)*velocity(2))* inVec(3)           &
      & + (  energy_coeff  ) * inVec(4)                                        &
      & + (  -energy_coeff + turb_coeff1 )* inVec(5)

    ! Fifth Row
    outVec(5) = (                                                              &
      &     -turb_coeff1*(state(5)/state(1)) )* inVec(1)                        &
      &     + turb_coeff1*inVec(5)

    ! Sixth Row
    outVec(6) = ( -turb_coeff2*state(6)/state(1)  )* inVec(1)                   &
      &         +  ( turb_coeff2 )* inVec(6)

  end function atl_mult_nu11_Rans_2d

  ! Multiplies the viscous flux matrux nu_21 with a given vector
  function atl_mult_nu21_Rans_2d( state, velocity, inVec, mu, lambda, &
    &                             rans_params)  &
    & result( outVec )
    ! -------------------------------------------------------------------------!
    !> The velocity
    real(kind=rk), intent(in) :: velocity(2)
    !> The state array
    real(kind=rk), intent(in) :: state(6)
    !> Vector to be multiplied with nu11
    real(kind=rk), intent(in) :: inVec(6)
    !> Dynamic Viscosity
    real(kind=rk), intent(in) :: mu
    !> Viscosity
    real(kind=rk), intent(in) :: lambda
    !> The constants for the Rans eqn
    type(atl_Navier_stokes_rans_type), intent(in) :: rans_params
    !> The result of the matrix vector product
    real(kind=rk) :: outVec(6)
    ! -------------------------------------------------------------------------!
    !> mu_turb
    real(kind=rk) :: mu_turb, mu_eff, lam_eff
    ! -------------------------------------------------------------------------!
    mu_turb = rans_params%c_mu*state(5)/exp(state(6)/state(1))

    mu_eff = mu + mu_turb
    lam_eff = lambda + 2.0*mu_turb/3.0

    ! First row has zeros only
    outVec(1) = 0.0_rk

    ! Second row
    outVec(2) = (-mu_eff/state(1)) * velocity(2) * inVec(1) &
            & + (mu_eff/state(1)) * inVec(3)

    ! Third row
    outVec(3) = (lam_eff/state(1)) * velocity(1) * inVec(1) &
            & + (-lam_eff/state(1)) * inVec(2)

    ! Fourth row
    outVec(4) = (                                                                           &
            &      ((-mu_eff+lam_eff)/state(1))*velocity(1)*velocity(2)                     &
            &   ) * inVec(1)                                                                &
            & + (                                                                           &
            &      (-lam_eff/state(1))*velocity(2)                                          &
            &   ) * inVec(2)                                                                &
            & + (                                                                           &
            &      (mu_eff/state(1))*velocity(1)                                            &
            &   ) * inVec(3)

    outVec(5) = 0.0_rk


    outVec(6) = 0.0_rk

  end function atl_mult_nu21_Rans_2d


  ! Multiplies the viscous flux matrux nu_12 with a given vector
  function atl_mult_nu12_Rans_2d( state, velocity, inVec, mu, lambda, &
    &                             rans_params )  &
    & result( outVec )
    ! -------------------------------------------------------------------------!
    !> The velocity
    real(kind=rk), intent(in) :: velocity(2)
    !> The state array
    real(kind=rk), intent(in) :: state(6)
    !> Vector to be multiplied with nu11
    real(kind=rk), intent(in) :: inVec(6)
    !> Dynamic Viscosity
    real(kind=rk), intent(in) :: mu
    !> Viscosity
    real(kind=rk), intent(in) :: lambda
    !> The constants for the Rans eqn
    type(atl_Navier_stokes_rans_type), intent(in) :: rans_params
    !> The result of the matrix vector product
    real(kind=rk) :: outVec(6)
    ! -------------------------------------------------------------------------!
    !> mu_turb
    real(kind=rk) :: mu_turb, mu_eff, lam_eff
    ! -------------------------------------------------------------------------!
    mu_turb = rans_params%c_mu*state(5)/exp(state(6)/state(1))

    mu_eff = mu + mu_turb
    lam_eff = lambda + 2.0*mu_turb/3.0

    ! First row has zeros only
    outVec(1) = 0.0_rk

    ! Second row
    outVec(2) = (lam_eff/state(1)) * velocity(2) * inVec(1) &
      &         + (-lam_eff/state(1)) * inVec(3)

    ! Third row
    outVec(3) = (-mu_eff/state(1)) * velocity(1) * inVec(1) &
       &        + (mu_eff/state(1)) * inVec(2)

    ! Fourth row
    outVec(4) = (                                                                           &
            &      ((-mu_eff+lam_eff)/state(1))*velocity(1)*velocity(2)                     &
            &   ) * inVec(1)                                                                &
            & + (                                                                           &
            &      (mu_eff/state(1))*velocity(2)                                            &
            &   ) * inVec(2)                                                                &
            & + (                                                                           &
            &      (-lam_eff/state(1))*velocity(1)                                          &
            &   ) * inVec(3)


     outVec(5) = 0.0_rk
     outVec(6) = 0.0_rk

  end function atl_mult_nu12_Rans_2d

  ! Multiplies the viscous flux matrux nu_22 with a given vector
  function atl_mult_nu22_Rans_2d( state, velocity, inVec,                     &
    &                isenCoeff, mu, lambda, thermCond, rans_params, heatCap)  &
    & result( outVec )
    ! -------------------------------------------------------------------------!
    !> The velocity
    real(kind=rk), intent(in) :: velocity(2)
    !> The state array
    real(kind=rk), intent(in) :: state(6)
    !> Vector to be multiplied with nu11
    real(kind=rk), intent(in) :: inVec(6)
    !> Adiabatice index, also known as isentropic expansion factor.
    real(kind=rk), intent(in) :: isenCoeff
    !> Dynamic Viscosity
    real(kind=rk), intent(in) :: mu
    !> Viscosity
    real(kind=rk), intent(in) :: lambda
    !> The thermal cond
    real(kind=rk), intent(in) :: thermCond
    !> The specific heat capacity (per mass unit mass, at constant volume)
    real(kind=rk), intent(in) :: heatCap
    !> The constants for the Rans eqn
    type(atl_Navier_stokes_rans_type), intent(in) :: rans_params
    !> The result of the matrix vector product
    real(kind=rk) :: outVec(6)
    ! -------------------------------------------------------------------------!
    !> mu_turb
    real(kind=rk) :: mu_turb, mu_eff, lam_eff, limited_eddy_visc, energy_coeff
    real(kind=rk) :: turb_coeff1, turb_coeff2
    ! -------------------------------------------------------------------------!
    mu_turb = rans_params%c_mu*state(5)/exp(state(6)/state(1))
    mu_eff = mu + mu_turb
    lam_eff = lambda + 2.0*mu_turb/3.0

    ! Stuff neeeded
    energy_coeff = thermCond/(heatCap*state(1))                               &
    &             + mu_turb*isenCoeff/(rans_params%turb_prandtl_num*state(1)  )

    ! @todo : NA : Calculate limited_eddy_visc Appropriately
    limited_eddy_visc = rans_params%alpha*max(state(5),0.0)              &
      &                                  *exp(state(6)/state(1))
    turb_coeff1 = (mu + rans_params%sig_k*limited_eddy_visc)/ state(1)
    turb_coeff2 = (mu + rans_params%sig_omg*limited_eddy_visc)/ state(1)

    ! First row has zeros only
    outVec(1) = 0.0_rk

    ! Second row
    outVec(2) = (-mu_eff/state(1)) * velocity(1) * inVec(1) &
            & + (mu_eff/state(1)) * inVec(2)

    ! Third row
    outVec(3) = ((-2.0_rk*mu_eff + lam_eff)/state(1))*velocity(2) * inVec(1) &
            & + ((2.0_rk*mu_eff - lam_eff)/state(1)) * inVec(3)

    ! Fourth row
    outVec(4) = (                                                                &
            &    (-mu_eff/state(1)) * velocity(1)*velocity(1)                    &
            &    + ((-2.0_rk*mu_eff + lam_eff)/state(1))*velocity(2)*velocity(2) &
            &    - energy_coeff*(state(4)/state(1) - sum(velocity(:)**2)         &
            &                    -state(5)/state(1)                    )         &
            &    - turb_coeff1*state(5)/state(1)        ) * inVec(1)             &
            & + (                                                                &
            &     (mu_eff/state(1) - energy_coeff)*velocity(1)                   &
            &   ) * inVec(2)                                                     &
            & + (                                                                &
            &     ((2.0_rk*mu_eff-lam_eff)/state(1) - energy_coeff )*velocity(2) &
            &   ) * inVec(3)                                                     &
            & + (                                                                &
            &     energy_coeff                                                   &
            &   ) * inVec(4)                                                     &
            & + (-energy_coeff + turb_coeff1 )* invec(5)

    ! Fifth row
    outVec(5) = ( -turb_coeff1*(state(5)/state(1)) )* inVec(1)                  &
      &           + ( turb_coeff1 )* inVec(5)

    ! Sixt Row
    outVec(6) = ( -turb_coeff2*state(6)/state(1)  ) * inVec(1)                   &
      &           + ( turb_coeff2 )* inVec(6)

  end function atl_mult_nu22_Rans_2d

end module atl_physFluxFilNvrStk_module

