! Copyright (c) 2015 Nikhil Anand <nikhil.anand@uni-siegen.de>
! Copyright (c) 2015-2017, 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016-2017 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2016-2017 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.
! **************************************************************************** !


module atl_eqn_filNvrStk_var_module
  use, intrinsic :: iso_c_binding,      only: c_loc, c_f_pointer, c_ptr
  use env_module,                       only: rk, labelLen

  use tem_varSys_module,                only: tem_varSys_type,             &
    &                                         tem_varSys_init,             &
    &                                         tem_varSys_append_stateVar,  &
    &                                         tem_varSys_append_derVar,    &
    &                                         tem_varSys_proc_point,       &
    &                                         tem_varSys_proc_element,     &
    &                                         tem_varSys_proc_setparams,    &
    &                                         tem_varSys_proc_getparams,    &
    &                                         tem_varSys_proc_setupIndices, &
    &                                         tem_varSys_proc_getValOfIndex
  use tem_grow_array_module,            only: init, append, truncate
  use tem_dyn_array_module,             only: init, append, truncate, &
    &                                         PositionOfVal
  use tem_logging_module,               only: logUnit
  use tem_varMap_module,                only: tem_possible_variable_type, &
    &                                         init, append
  use tem_aux_module,                   only: tem_abort
  use tem_operation_var_module,         only: tem_divideVecByScal_forPoint,  &
    &                                         tem_divideVecByScal_fromIndex, &
    &                                         tem_opVar_setupIndices,        &
    &                                         tem_get_new_varSys_data_ptr
  use tem_varSys_module,                only: tem_varSys_type
  use tem_time_module,                  only: tem_time_type

  use ply_oversample_module,            only: ply_convert2oversample,   &
    &                                         ply_convertFromoversample
  use ply_poly_project_module,          only: ply_poly_project_type, &
    &                                         assignment(=),         &
    &                                         ply_poly_project_n2m,  &
    &                                         ply_poly_project_m2n
  use ply_leg_diff_module,              only: ply_calcDiff_leg_2d

  use atl_equation_module,              only: atl_equations_type
  use atl_source_types_module,          only: atl_eqn_sourceMap_type
  use atl_eqn_filNvrStk_derive_module,  only: atl_Rans_pressure_getElement, &
    &                                         atl_Rans_pressure_getPoint
  use atl_varSys_module,                only: atl_varSys_solverData_type,    &
    &                                         atl_varSys_data_type,          &
    &                                         atl_varSys_getStateForElement, &
    &                                         atl_varSys_getStateForPoint,   &
    &                                         atl_get_new_varSys_data_ptr,   &
    &                                         atl_varSys_setupStateIndices,  &
    &                                         atl_varSys_getStateValofIndex
  use atl_operator_module,              only: atl_op_divideVecByScal_forElement, &
    &                                         atl_opVar_setupIndices
  use atl_source_types_module,          only: atl_source_op_type
  use atl_cube_elem_module,             only: atl_cube_elem_type

  implicit none

  private

  public :: atl_init_rans_vars
  public :: atl_init_rans_2d_vars
  public :: atl_append_rans_consVars
  public :: atl_append_rans_primVars
  public :: atl_init_filNvrStk_sourceTerms
  public :: atl_get_pointwise_velocity_gradient_2D
  public :: atl_get_pointwise_visc_stress_tensor_2D
  public :: atl_get_lower_bound_turb_disscipation
  public :: atl_init_RANS_closure_coeffs

contains

  subroutine atl_init_RANS_closure_coeffs( equation )
    ! --------------------------------------------------------------------------
    !> The equation system
    type(atl_equations_type), intent(inout) :: equation
    ! --------------------------------------------------------------------------
    equation%FiltNavierStokes%rans%turb_prandtl_num =  0.9_rk
    equation%FiltNavierStokes%rans%sig_k = 0.5_rk
    equation%FiltNavierStokes%rans%beta_k = 0.09_rk
    equation%FiltNavierStokes%rans%sig_omg = 0.5_rk
    equation%FiltNavierStokes%rans%alpha_omg = 5.0/9.0
    equation%FiltNavierStokes%rans%beta_omg = 3.0/40.0
    equation%FiltNavierStokes%rans%c_mu = 1.0
    equation%FiltNavierStokes%rans%alpha = 1.0

  end subroutine atl_init_RANS_closure_coeffs

  ! *******************************************************************************


  ! *******************************************************************************
  !> Init the variable system for filtered NavierStokes equation.
  !!
  !! The variable system describes, which variables are to be used and how
  !! they are organized in the memory.
  !! The first few variables up to the sys_mark are those, describing the
  !! state, and thus describe the output for regular restart files.
  !! Here these are the conservative variables density, momentum and energy.
  !! After the mark, there are additional values described that can be deduced
  !! from the state variables.
  subroutine atl_init_Rans_vars( equation, solverData )
    ! --------------------------------------------------------------------------
    !> The equation system
    type(atl_equations_type), intent(inout) :: equation
    !> the pointer to the data required for the varsys
    type(atl_varSys_solverData_type) :: solverData
    ! --------------------------------------------------------------------------
    integer :: nDim = 3
    ! --------------------------------------------------------------------------

    ! Initialize variable system
    call tem_varSys_init( me         = equation%varSys,      &
      &                   systemName = 'rans'                )

    ! Append the required conservative variables for the RANS model
    call atl_append_rans_consVars(equation, solverData)

    ! Append primitive variables
    equation%hasPrimitiveVariables = .true.
    call atl_append_rans_primVars(equation%varSys, equation%primVar,   &
     &                            solverData, nDim )

    ! Append derived quantities (also sets derive routines for primitive vars)
    ! @todo : check if all the derived vars are relevent
    !call append_euler_derivedVars(equation%varSys, solverData)

    ! for Navier Stokes 3D
    equation%temp%overSamp = 1
    equation%temp%modal = 0
    equation%temp%nodal = 3
    equation%temp%nScal = equation%varsys%nScalars

  end subroutine atl_init_Rans_vars

  ! *******************************************************************************

  subroutine atl_init_rans_2d_vars( equation, solverData )
    ! --------------------------------------------------------------------------
    !> The equation system
    type(atl_equations_type), intent(inout) :: equation
    !> the pointer to the data required for the varsys
    type(atl_varSys_solverData_type) :: solverData
    ! --------------------------------------------------------------------------
    integer :: nDim = 2
    ! --------------------------------------------------------------------------

    ! Initialize variable system
    call tem_varSys_init( me         = equation%varSys,      &
      &                   systemName = 'rans_2d'             )

    ! Append the required conservative variables for the RANS model
    call atl_append_rans_consVars(equation, solverData)

    ! Append primitive variables
    equation%hasPrimitiveVariables = .true.
    call atl_append_rans_primVars(equation%varSys, equation%primVar,   &
     &                            solverData, nDim                     )

    ! Append derived quantities (also sets derive routines for primitive vars)
    call append_rans_2d_derivedVars(equation%varSys, solverData)

    ! for Rans 2d
    equation%temp%overSamp = 1
    equation%temp%modal = 0
    equation%temp%nodal = 3
    equation%temp%nScal = equation%varsys%nScalars

  end subroutine atl_init_rans_2d_vars

  ! *******************************************************************************
  subroutine atl_append_rans_consVars(equation, solverData)
    ! ---------------------------------------------------------------------------
    !> The equation type.
    type(atl_equations_type), intent(inout) :: equation
    !> the pointer to the data required for the varsys
    type(atl_varSys_solverData_type), target :: solverData
    ! ---------------------------------------------------------------------------
    procedure(tem_varSys_proc_point), pointer :: get_point => null()
    procedure(tem_varSys_proc_element), pointer :: get_element => null()
    procedure(tem_varSys_proc_setparams), pointer :: set_params => NULL()
    procedure(tem_varSys_proc_getparams), pointer :: get_params => NULL()
    procedure(tem_varSys_proc_setupIndices), pointer :: setup_indices => NULL()
    procedure(tem_varSys_proc_getValOfIndex), pointer :: get_valOfIndex => NULL()
    ! ---------------------------------------------------------------------------

    allocate(equation%stateVar(5))

    get_element => atl_varSys_getStateForElement
    get_point => atl_varSys_getStateForPoint
    setup_indices => atl_varSys_setupStateIndices
    get_valOfIndex => atl_varSys_getStateValofIndex

    ! Append conservative variables to varSys
    call tem_varSys_append_stateVar(                              &
      & me             = equation%varSys,                         &
      & varName        = 'density',                               &
      & nComponents    = 1,                                       &
      & method_data    = atl_get_new_varSys_data_ptr(solverData), &
      & get_point      = get_point,                               &
      & get_element    = get_element,                             &
      & set_params     = set_params,                              &
      & get_params     = get_params,                              &
      & setup_indices  = setup_indices,                           &
      & get_valOfIndex = get_valOfIndex,                          &
      & pos            = equation%stateVar(1)                     )

    call tem_varSys_append_stateVar(                              &
      & me             = equation%varSys,                         &
      & varName        = 'momentum',                              &
      & nComponents    = equation%nDimensions,                    &
      & method_data    = atl_get_new_varSys_data_ptr(solverData), &
      & get_point      = get_point,                               &
      & get_element    = get_element,                             &
      & set_params     = set_params,                              &
      & get_params     = get_params,                              &
      & setup_indices  = setup_indices,                           &
      & get_valOfIndex = get_valOfIndex,                          &
      & pos            = equation%stateVar(2)                     )

    call tem_varSys_append_stateVar(                             &
      & me             = equation%varSys,                         &
      & varName        = 'energy',                                &
      & nComponents    = 1,                                       &
      & method_data    = atl_get_new_varSys_data_ptr(solverData), &
      & get_point      = get_point,                               &
      & get_element    = get_element,                             &
      & set_params     = set_params,                              &
      & get_params     = get_params,                              &
      & setup_indices  = setup_indices,                           &
      & get_valOfIndex = get_valOfIndex,                          &
      & pos            = equation%stateVar(3)                     )

    ! Append the Turbulent Kinetic Energy
    call tem_varSys_append_stateVar(                              &
      & me             = equation%varSys,                         &
      & varName        = 'k',                                     &
      & nComponents    = 1,                                       &
      & method_data    = atl_get_new_varSys_data_ptr(solverData), &
      & get_point      = get_point,                               &
      & get_element    = get_element,                             &
      & set_params     = set_params,                              &
      & get_params     = get_params,                              &
      & setup_indices  = setup_indices,                           &
      & get_valOfIndex = get_valOfIndex,                          &
      & pos            = equation%stateVar(4)                     )

    ! Append the variable for specific dissipation rate
    call tem_varSys_append_stateVar(                              &
      & me             = equation%varSys,                         &
      & varName        = 'omega',                                 &
      & nComponents    = 1,                                       &
      & method_data    = atl_get_new_varSys_data_ptr(solverData), &
      & get_point      = get_point,                               &
      & get_element    = get_element,                             &
      & set_params     = set_params,                              &
      & get_params     = get_params,                              &
      & setup_indices  = setup_indices,                           &
      & get_valOfIndex = get_valOfIndex,                          &
      & pos            = equation%stateVar(5)                     )

  end subroutine atl_append_rans_consVars

  ! *******************************************************************************

  subroutine atl_append_rans_primVars(varSys, primVar, solverData, nDim)
    ! --------------------------------------------------------------------------
    !> The variable system to append primitive variables to
    type(tem_varSys_type), intent(inout)  :: varSys

    !> Indices of the primitive variables in the overall system.
    integer, allocatable, intent(out)     :: primVar(:)

    !> the pointer to the data required for the varsys to fulfill all operations
    !! and derivations on the variables
    type(atl_varSys_solverData_type), target :: solverData
    integer, intent(in) :: nDim
    ! --------------------------------------------------------------------------
    logical :: wasAdded
    procedure(tem_varSys_proc_point), pointer :: get_point => NULL()
    procedure(tem_varSys_proc_element), pointer :: get_element => NULL()
    procedure(tem_varSys_proc_setparams), pointer :: set_params => NULL()
    procedure(tem_varSys_proc_getparams), pointer :: get_params => NULL()
    procedure(tem_varSys_proc_setupIndices), pointer :: setup_indices => NULL()
    procedure(tem_varSys_proc_getValOfIndex), pointer :: get_valOfIndex => NULL()
    type(c_ptr) :: solver_bundle
    ! --------------------------------------------------------------------------

    ! Save var position in primVar
    allocate(primVar(5))

    ! First primitive variable is the density.
    ! It is already part of the variable system, thus we have to find it.
    primVar(1) = PositionOfVal(varSys%varname, 'density')

    ! Second primitive variable is the velocity.
    ! To compute it, density and momentum of the state are required.
    get_element => atl_op_divideVecByScal_forElement
    get_point => tem_divideVecByScal_forPoint
    setup_indices => tem_opVar_setupIndices
    get_valOfIndex => tem_divideVecByScal_fromIndex

    ! get c_ptr from solver method_Data and store it tem_varSys_op_data_type
    solver_bundle = atl_get_new_varSys_data_ptr(solverData)

    call tem_varSys_append_derVar(                                   &
      & me             = varSys,                                     &
      & varName        = 'velocity',                                 &
      & operType       = 'divide_vector_by_scalar',                  &
      & nComponents    = nDim,                                       &
      & input_varname  = ['momentum', 'density '],                   &
      & method_data    = tem_get_new_varSys_data_ptr(solver_bundle), &
      & get_point      = get_point,                                  &
      & get_element    = get_element,                                &
      & set_params     = set_params,                                 &
      & get_params     = get_params,                                 &
      & setup_indices  = setup_indices,                              &
      & get_valOfIndex = get_valOfIndex,                             &
      & wasAdded       = wasAdded,                                   &
      & pos            = primVar(2)                                  )

    if (wasAdded) then
      write(logUnit(10),*) 'Appended variable: velocity'
    else
      write(logUnit(1),*) 'Error: variable velocity not appended'
      call tem_abort()
    end if

    ! Third primitive variable is the pressure.
    ! To compute it all three state variables are required.
    get_element => atl_Rans_pressure_getElement
    get_point => atl_Rans_pressure_getPoint
    setup_indices => atl_opVar_setupIndices

    call tem_varSys_append_derVar(                                &
      & me             = varSys,                                  &
      & varName        = 'pressure',                              &
      & nComponents    = 1,                                       &
      & input_varname  = ['density ', 'momentum',                 &
      &                   'energy  ', 'k       '],                &
      & method_data    = atl_get_new_varSys_data_ptr(solverData), &
      & get_point      = get_point,                               &
      & get_element    = get_element,                             &
      & set_params     = set_params,                              &
      & get_params     = get_params,                              &
      & setup_indices  = setup_indices,                           &
      & get_valOfIndex = get_valOfIndex,                          &
      & pos            = primVar(3)                               )

    if (wasAdded) then
      write(logUnit(10),*) 'Appended variable: Pressure'
    else
      write(logUnit(1),*) 'Error: variable Pressure not appended'
      call tem_abort()
    end if

    ! Fourth primitive variable is the turbulent KE.
    ! It is already part of the variable system, thus we have to find it.
    primVar(4) = PositionOfVal(varSys%varname, 'k')


    ! Fifth primitive variable is the specific dissipation rate.
    ! It is already part of the variable system, thus we have to find it.
    primVar(5) = PositionOfVal(varSys%varname, 'omega')


    get_element => null()
    get_point => null()
    set_params => null()
    get_params => null()
    setup_Indices => null()
    get_valOfindex => null()

  end subroutine atl_append_rans_primVars

! *******************************************************************************

  subroutine append_rans_2d_derivedVars( varSys, solverData )
    ! ------------------------------------------------------------------------
    type(tem_varSys_type), intent(inout) :: varSys
    !> the pointer to the data required for the varsys to fulfill all operations
    !! and derivations on the variables
    type(atl_varSys_solverData_type), target :: solverData
    ! ------------------------------------------------------------------------
!    integer :: nDerivedVars, iVar, nComponents
!    character(len=20), allocatable :: derVarName(:)
!    character(len=labelLen), allocatable :: invar_name(:)
!    character(len=labelLen) :: varname
!    logical :: wasAdded
    procedure(tem_varSys_proc_point), pointer :: get_point => NULL()
    procedure(tem_varSys_proc_element), pointer :: get_element => NULL()
    procedure(tem_varSys_proc_setparams), pointer :: set_params => NULL()
    procedure(tem_varSys_proc_getparams), pointer :: get_params => NULL()
    procedure(tem_varSys_proc_setupIndices), pointer :: setup_indices => NULL()
    procedure(tem_varSys_proc_getValOfIndex), pointer :: get_valOfIndex => NULL()
    ! ------------------------------------------------------------------------
!    wasAdded = .False.
!
!    !@todo Write derive routines for the primitive variables
!    nDerivedVars = 1
!    allocate(derVarName(nDerivedVars))
!
!
!    derVarName    = [ 'track_source  ']
!
!    do iVar = 1, nDerivedVars
!      varname = trim(adjustl(derVarName(iVar)))
!      select case(varname)
!!
!!      ! Define the speed of sound as new possible variable to derive.
!      case ('track_source')
!        get_point => Null()
!        get_element => Null()
!        nComponents = 6
!        allocate(invar_name(5))
!        invar_name(1) = 'density'
!        invar_name(2) = 'momentum'
!        invar_name(3) = 'energy'
!        invar_name(4) = 'k'
!        invar_name(5) = 'omega'
!
!      case default
!        write(logUnit(1),*) 'WARNING: Unknown variable: '//trim(varname)
!        cycle !go to next variable
!
!      end select
!
!      ! append variable to varSys
!      call tem_varSys_append_derVar(                                &
!        & me             = varSys,                                  &
!        & varName        = varname,                                 &
!        & nComponents    = nComponents,                             &
!        & input_varname  = invar_name,                              &
!        & method_data    = atl_get_new_varSys_data_ptr(methoddata), &
!        & get_point      = get_point,                               &
!        & get_element    = get_element,                             &
!        & wasAdded       = wasAdded                                 )
!
!      call outputAppendResult( wasAdded, varname )
!
!      deallocate(invar_name)
!
!    end do
!
!    deallocate(derVarName)

  end subroutine append_rans_2d_derivedVars
 ! *******************************************************************************


  !> Init source terms for flow simulations.
  !> This routine initializes possible source variables and returns the filled
  !! up list of the poss_srcVars
  subroutine atl_init_filNvrStk_sourceTerms(possVars, eval_source, model_type)
    ! --------------------------------------------------------------------------
    type(tem_possible_variable_type), intent(inout)  :: possVars
    type(atl_eqn_sourceMap_type), allocatable, intent(out) :: eval_source(:)
    character(len=32), intent(in) :: model_type
    ! --------------------------------------------------------------------------
    integer :: pos
    ! --------------------------------------------------------------------------

    allocate(eval_source(1))
    call init(me = possVars%varName, length = 1)


    select case (trim(model_type))

    case('rans')
!NA!      ! Add the arbitrary Source term
!NA!      call append( me          = possVars,     &
!NA!        &          varName     = 'ransSource', &
!NA!        &          nComponents = 7,            &
!NA!        &          pos         = pos           )
!NA!      eval_source(pos)%do => eval_source_rans

    case('rans_2d')

      ! Add the source term in the
      call append( me          = possVars,       &
        &          varName     = 'rans2dsource', &
        &          nComponents = 6,              &
        &          pos         = pos             )
      eval_source(pos)%do => eval_source_rans2D

    end select


    call truncate(possVars%varName)
    call truncate(possVars%nComponents)

  end subroutine atl_init_filNvrStk_sourceTerms

  subroutine eval_source_rans2D( fun, varSys, time, mesh, poly_proj,       &
    &                            currentLevel, state, material, sourcedata )
    ! ---------------------------------------------------------------------------
    !> Description of method to update source
    class(atl_source_op_type), intent(in) :: fun
    !> variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> Point in time at which to evaluate the variable.
    type(tem_time_type), intent(in)  :: time
    !> Current level mesh information
    type(atl_cube_elem_type), intent(in) :: mesh
    !> Parameters for projection
    type(ply_poly_project_type), intent(inout) :: poly_proj
    !> current level
    integer, intent(in) :: currentLevel
    !> The state in modal space.
    !! This is needed for several source terms that have to be applied to the
    !! current state
    real(kind=rk), intent(in) :: state(:,:,:)
    !> Material description for the background on the current level
    real(kind=rk), intent(in) :: material(:)
    !> The source data to update. When all source terms are added to this
    !! buffer, it is applied to the state.
    real(kind=rk), intent(inout) :: sourcedata(:,:,:)
    ! ---------------------------------------------------------------------------
    !procedure(compute_source_interface) , pointer:: evaluate_source
    ! The modal coefficients of the current element in the loop.
    real(kind=rk), allocatable :: modalCoeffs(:,:), modal_grad(:,:,:)
    ! Nodal representation of the polynomial with in each cell.
    real(kind=rk), allocatable :: pointVal(:,:), pointVal_grad(:,:,:)
    real(kind=rk), allocatable :: rhsNodal(:,:), rhsModal(:,:), rhs(:,:)
    real(kind=rk) :: rey_stress_tensor(2,2)
    integer :: iElem, iDir
    integer :: iPoint
    real(kind=rk) :: beta_k, omega_r, alp_omg, beta_omg, sig_omg
    real(kind=rk) :: mu_turb, mu, c_mu, alpha,  source(2), omega, k_bar
    real(kind=rk) :: velGrad(2,2), ViscStressTensor(2,2)
    type(atl_varSys_data_type), pointer :: fPtr
    type(atl_equations_type), pointer :: eqn
    real(kind=rk) :: d_omega1, d_omega2
    ! ---------------------------------------------------------------------------
    call C_F_POINTER( varSys%method%val(fun%srcTerm_varPos)%method_Data, fPtr )
    eqn => fPtr%solverData%equationPtr

    beta_k   = eqn%FiltNavierStokes%rans%beta_k
    beta_omg = eqn%FiltNavierStokes%rans%beta_omg
    alp_omg  = eqn%FiltNavierStokes%rans%alpha_omg
    c_mu     = eqn%FiltNavierStokes%rans%c_mu
    alpha    = eqn%FiltNavierStokes%rans%alpha
    mu       = eqn%NavierStokes%mu
    sig_omg  = eqn%FiltNavierStokes%rans%sig_omg


    allocate( modalCoeffs (poly_proj%Body_2D%OverSamp_dofs ,varSys%nScalars ))
    allocate( modal_grad (poly_proj%Body_2D%OverSamp_dofs,varSys%nScalars,2 ))
    allocate( pointVal (poly_proj%Body_2D%nQuadPoints,varSys%nScalars )     )
    allocate( pointVal_grad (poly_proj%Body_2D%nQuadPoints,varSys%nScalars,2 ))
    allocate( rhsNodal (poly_proj%Body_2D%nQuadPoints, varSys%nScalars )    )
    allocate( rhsModal (poly_proj%Body_2D%OverSamp_dofs, varSys%nScalars )  )
    allocate( rhs (poly_proj%body_2d%nDofs, varSys%nScalars ))

    rhsNodal = 0.0_rk
    do iElem = 1, fun%elems(currentLevel)%nElems

      ! --> modal space
      call ply_convert2oversample(state       = state(iElem,:,:),  &
        &                         poly_proj   = poly_proj,         &
        &                         nDim        = 2,                 &
        &                         modalCoeffs = modalCoeffs        )
      ! --> oversamp modal space

      ! Calculate the gradient of modal Coeffs
      call ply_calcDiff_leg_2d( legCoeffs     = modalCoeffs,               &
        &                       legCoeffsDiff = modal_grad,                &
        &                       maxPolyDegree = poly_proj%oversamp_degree, &
        &                       nVars         = varsys%nScalars,           &
        &                       elemLength    = mesh%length                )

      ! Now, we transform the modal representation of this element to nodal
      ! space by making use of fast polynomial transformations (FPT)
      call ply_poly_project_m2n(me = poly_proj,                  &
       &                       dim = 2 ,                         &
       &                       nVars = varSys%nScalars,          &
       &                       nodal_data= pointVal,             &
       &                       modal_data= modalCoeffs           )
      ! --> oversamp nodal space

      ! Now, transform modal gradient representation of this element to
      ! nodal space by making use of fast polynomial transformations (FPT)
      do iDir = 1, 2
        call ply_poly_project_m2n(me   = poly_proj,               &
          &                 dim        = 2,                       &
          &                 nVars      = varSys%nScalars,         &
          &                 nodal_data = pointVal_grad(:,:,iDir), &
          &                 modal_data = modal_grad(:,:,iDir)     )
      end do

      do iPoint = 1, poly_proj%Body_2D%nQuadPoints


        call atl_get_pointwise_velocity_gradient_2D( &
          & pointVal_grad(iPoint,:,:),               &
          & pointVal(iPoint,:),                      &
          & velGrad                                  )

        call atl_get_pointwise_visc_stress_tensor_2D(         &
          &          velGrad = velGrad,                       &
          &          S       = ViscStressTensor               )

        ! get the value of \omega_r
        call atl_get_lower_bound_turb_disscipation(                      &
          &          S       = ViscStressTensor,                         &
          &        c_mu      = c_mu,                                     &
          &        omega     = pointVal(iPoint, 6 )/ pointVal(iPoint,1), &
          &        omega_r   = omega_r                                   )

        ! Calculate the limited turbulent eddy viscosity
        mu_turb = alpha * max(pointVal(iPoint,5),0.0)* exp(-omega_r)

        !Careful - This routine returns the reynolds tensor divided by
        ! the limited "k" (\bar{k})
        call get_rans_reynolds_tensor_2D(                      &
          &    state             = pointVal(iPoint,:),         &
          &    alpha             = alpha,                      &
          &    omega_r           = omega_r,                    &
          &    velGrad           = velGrad,                    &
          &    rey_stress_tensor = rey_stress_tensor           )

        omega = exp( pointVal(iPoint,6)/ pointVal(iPoint,1) )
        k_bar = max(pointVal(iPoint,5)/pointVal(iPoint,1),0.0)
        d_omega1 = (pointVal_grad(iPoint,6,1) - pointVal(iPoint,6)         &
          &           *pointval_grad(iPoint,1,1)/pointVal(iPoint,1)    )   &
          &              *pointVal(iPoint,1)
        d_omega2 = (pointVal_grad(iPoint,6,2) - pointVal(iPoint,6)         &
          &           *pointval_grad(iPoint,1,2)/pointVal(iPoint,1)    )   &
          &              / pointVal(iPoint,1)

        ! Note: The first four comonents of the source term = 0
        ! This is the source corresponding to the fifth var i.e turbulent KE
        ! the multiplication with k_bar as the rey_stress_tensor used here is
        ! divided by k_bar
        source(1) = rey_stress_tensor(1,1)*velGrad(1,1)*k_bar             &
          &       + rey_stress_tensor(1,2)*velGrad(1,2)*k_bar             &
          &       + rey_stress_tensor(2,1)*velGrad(2,1)*k_bar             &
          &       + rey_stress_tensor(2,2)*velGrad(2,2)*k_bar             &
          &       - beta_k*max(pointVal(iPoint,5),0.0)&
                          *exp(omega_r)

        ! This is the source corresponding to the sixth conservative var
        ! i.e specific dissipation rate
        source(2) = (                                         &
          &   rey_stress_tensor(1,1)*velGrad(1,1)             &
          & + rey_stress_tensor(1,2)*velGrad(1,2)             &
          & + rey_stress_tensor(2,1)*velGrad(2,1)             &
          & + rey_stress_tensor(2,2)*velGrad(2,2) )           &
          &   *alp_omg*omega                                  &
          & - beta_omg*pointVal(iPoint,1)*exp(omega_r)        &
          & + (mu + sig_omg*mu_turb)                          &
          & * ( d_omega1*d_omega1                             &
          & +   d_omega2*d_omega2                             )

          !write(*,*) "**>", iPoint, source, pointVal(iPoint,:)
          ! Debug statement
          rhsNodal(iPoint,5:6) = rhsNodal(iPoint,5:6) +  source(1:2)

      end do

      ! Convert everything back to modal space
      call ply_poly_project_n2m( me         = poly_proj,       &
        &                        dim        = 2,               &
        &                        nVars      = varSys%nScalars, &
        &                        nodal_data = rhsNodal(:,:),   &
        &                        modal_data = rhsModal(:,:)    )

      ! --> oversamp modal space
      call ply_convertFromoversample( modalCoeffs = rhsModal(:,:), &
        &                             poly_proj   = poly_proj,     &
        &                             nDim        = 2,             &
        &                             state       = rhs(:,:)       )
      ! --> oversampled modal space

      ! Add rhsmodal to the source representation
      sourcedata(iElem,:,:) = sourcedata(iElem,:,:) + rhs(:,:)


    end do ! Elems Loop

    deallocate( modalCoeffs )
    deallocate( modal_grad )
    deallocate( pointVal )
    deallocate( pointVal_grad )
    deallocate( rhsNodal )
    deallocate( rhsModal )
    deallocate( rhs )


  end subroutine eval_source_rans2D

! *******************************************************************************

  ! Calculates the velocity gradient from the available state gradient
  subroutine atl_get_pointwise_velocity_gradient_2D( state_grad, state, res )
    ! ---------------------------------------------------------------------------
    real(kind=rk), intent(in) :: state_grad(:,:)
    ! The gradient of the state
    real(kind=rk), intent(in) :: state(:)
    ! Reynolds_stress_tensor
    real(kind=rk), intent(inout) :: res(:,:)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: dens_sq
    ! ---------------------------------------------------------------------------
    dens_sq = state(1)*state(1)
    ! fill up the four velocity gradients.
    res(1,1) = state_grad(2,1)/ state(1) - state(2) * state_grad(1,1)/ dens_sq
    res(2,1) = state_grad(3,1)/ state(1) - state(3) * state_grad(1,1)/ dens_sq
    res(1,2) = state_grad(2,2)/ state(1) - state(2) * state_grad(1,2)/ dens_sq
    res(2,2) = state_grad(3,2)/ state(1) - state(3) * state_grad(1,2)/ dens_sq


  end subroutine atl_get_pointwise_velocity_gradient_2D

! *******************************************************************************

  ! Takes the pointwise velocity gradient as the input and evaluates
  ! the four components of viscous stress tensor
  subroutine atl_get_pointwise_visc_stress_tensor_2D( velGrad, S )
    ! ---------------------------------------------------------------------------
    ! The gradient of the state
    real(kind=rk), intent(in) :: velGrad(2,2)
    ! Viscous stress tensor - the output
    real(kind=rk), intent(inout) :: S(2,2)
    ! ---------------------------------------------------------------------------

    S(1,1) = 4.0* velGrad(1,1)/ 3.0 - 2.0* velGrad(2,2)/ 3.0
    S(1,2) = velGrad(1,2) + velGrad(2,1)
    S(2,1) = S(1,2)
    S(2,2) = 4.0 * velGrad(2,2)/ 3.0 - 2.0* velGrad(1,1)/ 3.0

  end subroutine atl_get_pointwise_visc_stress_tensor_2D

! *******************************************************************************

  ! takes pointwise input and evaluates the lower bound of the turbulent
  ! dissipation to ensure the positivity of normal turbulent stresses.
  ! The fulfilement of the inequality solved in this routine is the same as
  ! one proposed by Wilcox for the k-\omega turbulence model and also used
  ! in the Hartmann paper. Here we just solve this quadratic inequality
  subroutine atl_get_lower_bound_turb_disscipation( S, c_mu, omega, Omega_r )
    ! Viscous stress tensor
    real(kind=rk), intent(in) :: S(2,2)
    ! The constant c_mu
    real(kind=rk), intent(in) :: c_mu
    ! The state omega
    real(kind=rk), intent(in) :: omega
    ! The lower bound on omega
    real(kind=rk), intent(out) :: omega_r
    ! ---------------------------------------------------------------------------
    ! The constants needed
    real(kind=rk) :: C3,C4
    ! The omegas
    real(kind=rk) :: omega_r0
    real(kind=rk) :: omg1, omg2, omg3
    ! ---------------------------------------------------------------------------

    if (S(1,1) > 0.0_rk) then
      omg1 =   3.0 *c_mu * (S(1,1) ) / 2.0
    else
      omg1 = 0.0_rk
    end if

    if (S(2,2) > 0.0_rk) then
      omg2 =   3.0 *c_mu * (S(2,2) ) / 2.0
    else
      omg2 = 0.0_rk
    end if

    C3 = - 3.0 *c_mu * ( S(1,1) + S(2,2))/2.0
    C4 =   9.0 *c_mu * c_mu * (S(1,1)* S(2,2) - S(1,2)*S(1,2))/ 4.0

    if ( (C3 * C3 - 4.0*C4) > 0 ) then
      omg3 = - (C3 + sqrt( C3 * C3 - 4.0*C4) )/ 2.0
    else
      omg3 = 0.0_rk
    end if

    omega_r0 = max(omg1, omg2, omg3, tiny(omega))
    omega_r0 = log(omega_r0)

    omega_r = max(omega, omega_r0)


  end subroutine atl_get_lower_bound_turb_disscipation


! *******************************************************************************
  ! This routine calculates the reynolds tensor and returns the result divided
  ! by \bar{k} = \frac{\bar{\mu_t}}{\alpha \rho e^{\tilde{\omega_r}}}
  subroutine get_rans_reynolds_tensor_2D( state, alpha, omega_r, velGrad, &
    &                                     rey_stress_tensor )
    ! ---------------------------------------------------------------------------
    !! current state
    real(kind=rk), intent(in) :: state(:)
    ! \alpha for the evaluation of limited \mu_t (turbulent eddy viscosity)
    real(kind=rk), intent(in) :: alpha
    ! lower bound of \omega = \omega_r for the evaluation of limited
    ! \mu_t (turbulent eddy viscosity)
    real(kind=rk), intent(in) :: omega_r
    ! The velocity gradient
    real(kind=rk), intent(in) :: velGrad(2,2)
    ! Reynolds_stress_tensor
    real(kind=rk), intent(inout) :: rey_stress_tensor(:,:)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: fact
    ! ---------------------------------------------------------------------------

    fact = alpha*state(1)*exp(-omega_r)


    rey_stress_tensor(1,1) =  (4.0*velGrad(1,1)  - 2.0*velGrad(2,2) )  &
      &                          *fact / 3.0   - 2.0 * state(1) / 3.0

    rey_stress_tensor(1,2) = fact*(velGrad(1,2) + velGrad(2,1))

    rey_stress_tensor(2,1) = rey_stress_tensor(1,2)

    rey_stress_tensor(2,2) = fact*( 4.0*velGrad(2,2) - 2.0*velGrad(1,1)) / 3.0 &
      &                        - 2.0 * state(1) / 3.0

  end subroutine get_rans_reynolds_tensor_2D

! *******************************************************************************

  subroutine get_rans_reynolds_tensor( state, state_gradient, &
    &                                  rey_stress_tensor )
    ! ---------------------------------------------------------------------------
    !! current state
    real(kind=rk), intent(in) :: state(:)
    ! The gradient of the state
    real(kind=rk), intent(in) :: state_gradient(:,:)
    ! Reynolds_stress_tensor
    real(kind=rk), intent(inout) :: rey_stress_tensor(:,:)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: mu, lambda
        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
    ! ---------------------------------------------------------------------------
    !mu         =fun%equation%NavierStokes%mu
    !lambda     =fun%equation%NavierStokes%lambda

    ! 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 finished.
    call tem_abort( "This routine doesn't use properly initialized variables" )

    velocity(1:3) = state(2:4)/state(1)
    mu_turb = state(1)*state(6)/state(7)
    lam_turb = 2.0*mu_turb/3.0
    mu_eff = mu + mu_turb
    lam_eff = lambda + lam_turb

    rey_stress_tensor(1,1) =                                      &
      & ( (-2*mu_eff + lam_eff)* velocity(1)*state_gradient(1,1)  &
      & + lam_eff*velocity(2)*state_gradient(1,2)                 &
      & + lam_eff*velocity(3)*state_gradient(1,3)                 &
      & + (2*mu_eff - lam_eff)*state_gradient(2,1)                &
      & - lam_eff*state_gradient(3,2)                             &
      & - lam_eff*state_gradient(4,3) )/state(1)

    rey_stress_tensor(1,2) =                                      &
     &  ( - state_gradient(2,2)                                   &
     &    + velocity(1)*state_gradient(1,2)                       &
     &    - state_gradient(3,1)                                   &
     &    + velocity(2)*state_gradient(1,1) )*mu_eff / state(1)

    rey_stress_tensor(1,3) =                                      &
     &  ( - state_gradient(4,1)                                   &
     &    + velocity(3)*state_gradient(1,1)                       &
     &    - state_gradient(2,3)                                   &
     &    + velocity(1)*state_gradient(1,3) )*mu_eff / state(1)

    rey_stress_tensor(2,1) =                                      &
     &  ( - state_gradient(2,2)                                   &
     &    + velocity(1)*state_gradient(1,2)                       &
     &    - state_gradient(3,1)                                   &
     &    + velocity(2)*state_gradient(1,1) )*mu_eff / state(1)

    rey_stress_tensor(2,2) =                                      &
      & ( lam_eff*velocity(1)*state_gradient(1,1)                 &
      & + (-2*mu_eff + lam_eff)* velocity(2)*state_gradient(1,2)  &
      & + lam_eff*velocity(3)*state_gradient(1,3)                 &
      & - lam_eff*state_gradient(2,1)                             &
      & + (2*mu_eff - lam_eff)*state_gradient(3,2)                &
      & - lam_eff*state_gradient(4,3) )/state(1)

    rey_stress_tensor(2,3) =                                      &
     &  ( - state_gradient(4,2)                                   &
     &    + velocity(3)*state_gradient(1,2)                       &
     &    - state_gradient(3,3)                                   &
     &    + velocity(2)*state_gradient(1,3) )*mu_eff / state(1)

    rey_stress_tensor(3,1) =                                      &
     &  ( - state_gradient(2,3)                                   &
     &    + velocity(1)*state_gradient(1,3)                       &
     &    - state_gradient(4,1)                                   &
     &    + velocity(3)*state_gradient(1,1) )*mu_eff / state(1)

    rey_stress_tensor(3,2) =                                      &
     &  ( - state_gradient(3,3)                                   &
     &    + velocity(2)*state_gradient(1,3)                       &
     &    - state_gradient(4,2)                                   &
     &    + velocity(3)*state_gradient(1,2) )*mu_eff / state(1)

    rey_stress_tensor(3,3) =                                      &
      & ( lam_eff*velocity(1)*state_gradient(1,1)                 &
      & + lam_eff*velocity(2)*state_gradient(1,2)                 &
      & + (-2*mu_eff + lam_eff)* velocity(3)*state_gradient(1,3)  &
      & - lam_eff*state_gradient(2,1)                             &
      & - lam_eff*state_gradient(3,2)                             &
      & + (2*mu_eff - lam_eff)*state_gradient(4,3)) /state(1)

  end subroutine get_rans_reynolds_tensor

! *******************************************************************************

end module atl_eqn_filNvrStk_var_module
