!> definition of 2D models which are simulated:
!> Eulerian model, mass diffusion
module modelEulerian

  use data_mod
  use main_data
  use f_mapping
  use mesh_oper
  use blocks_integ
  use matrix_oper_int

  implicit none

  public:: Set_R_s_Eul
  public:: Set_K_sk_Eul

contains


  !> compute viscous fluxes R_s, s=1,2 Eulerian model, mass diffusion
  !> in integ nodes
  subroutine Set_R_s_Eul(ndimL, nbDim, iRe, Qdof, w, Dw, Re_1, R_s, xi)
    integer, intent(in) :: Qdof, nbDim, iRe, ndimL
    real, dimension(1:Qdof, 1:ndimL), intent(in):: w !state  w in #Qdof nodes
    real, dimension(1:Qdof, 1:ndimL, 1:nbDim), intent(in):: Dw !state  Dw in #Qdof nodes
    real, dimension(1:iRe, 1:Qdof), intent(in) :: Re_1        ! inverse of Reynolds number
    !real, intent(in) :: Re_1                     ! inverse of Reynolds number
    real, dimension(1:Qdof, 1:nbDim, 1:ndimL), intent(inout) :: R_s
    real, dimension(1:Qdof, 1:nbDim), intent(in):: xi ! physical coordinates
    real, dimension(:), allocatable :: oRe

    allocate( oRe(1:Qdof) )

    !u(1:Qdof) = w(1:Qdof, 2)   / w(1:Qdof, 1)
    !v(1:Qdof) = w(1:Qdof, 3)   / w(1:Qdof, 1)
    !e(1:Qdof) = w(1:Qdof, 4)   / w(1:Qdof, 1)

    oRe(1:Qdof) = state%model%alpha * Re_1(1,1:Qdof) / w(1:Qdof, 1)

    !p(1:Qdof) = state%model%kappa1 * (w(1:Qdof,4)  &
    !     - (u(1:Qdof)*w(1:Qdof, 2) + v(1:Qdof)*w(1:Qdof, 3) )/2. )

    R_s(1:Qdof, 1, 1) = oRe(1:Qdof)* Dw(1:Qdof, 1, 1)
    R_s(1:Qdof, 2, 1) = oRe(1:Qdof)* Dw(1:Qdof, 1, 2)

    R_s(1:Qdof, 1, 2) = oRe(1:Qdof)* Dw(1:Qdof, 2, 1)
    R_s(1:Qdof, 2, 2) = oRe(1:Qdof)* Dw(1:Qdof, 2, 2)

    R_s(1:Qdof, 1, 3) = oRe(1:Qdof)* Dw(1:Qdof, 3, 1)
    R_s(1:Qdof, 2, 3) = oRe(1:Qdof)* Dw(1:Qdof, 3, 2)

    R_s(1:Qdof, 1, 4) = oRe(1:Qdof)* Dw(1:Qdof, 4, 1)
    R_s(1:Qdof, 2, 4) = oRe(1:Qdof)* Dw(1:Qdof, 4, 2)
    
    ! R_s(1:Qdof, 1:nbDim, 1) = 0.

    ! R_s(1:Qdof, 1, 2) = 2./3* oRe(1:Qdof) &
    !      *(2*(Dw(1:Qdof, 2, 1) - u(1:Qdof)*Dw(1:Qdof, 1, 1) )&
    !      - (Dw(1:Qdof, 3, 2) - v(1:Qdof)*Dw(1:Qdof, 1, 2) ) )

    ! R_s(1:Qdof, 1, 3) = oRe(1:Qdof) &
    !      *((Dw(1:Qdof, 3, 1) - v(1:Qdof)*Dw(1:Qdof, 1, 1) )&
    !      + (Dw(1:Qdof, 2, 2) - u(1:Qdof)*Dw(1:Qdof, 1, 2) ) )


    ! R_s(1:Qdof, 1, 4) = u(1:Qdof) * R_s(1:Qdof, 1, 2) + v(1:Qdof) * R_s(1:Qdof, 1, 3) &
    !      + state%model%kappa/state%model%Pr *  oRe(1:Qdof) &
    !      * ( Dw(1:Qdof, 4, 1) - e(1:Qdof)*Dw(1:Qdof, 1, 1) &
    !      - (u(1:Qdof) * Dw(1:Qdof, 2, 1) + v(1:Qdof) * Dw(1:Qdof, 3, 1)) &
    !      + (u(1:Qdof) * u(1:Qdof) + v(1:Qdof) *v(1:Qdof) ) *  Dw(1:Qdof, 1, 1) )



    ! R_s(1:Qdof, 2, 2) = oRe(1:Qdof) &
    !      *((Dw(1:Qdof, 3, 1) - v(1:Qdof)*Dw(1:Qdof, 1, 1) )&
    !      + (Dw(1:Qdof, 2, 2) - u(1:Qdof)*Dw(1:Qdof, 1, 2) ) )

    ! R_s(1:Qdof, 2, 3) = 2./3* oRe(1:Qdof) &
    !      *(2*(Dw(1:Qdof, 3, 2) - v(1:Qdof)*Dw(1:Qdof, 1, 2) )&
    !      - (Dw(1:Qdof, 2, 1) - u(1:Qdof)*Dw(1:Qdof, 1, 1) ) )


    ! R_s(1:Qdof, 2, 4) = u(1:Qdof) * R_s(1:Qdof, 2, 2) + v(1:Qdof) * R_s(1:Qdof, 2, 3) &
    !      + state%model%kappa/state%model%Pr *  oRe(1:Qdof) &
    !      * ( Dw(1:Qdof, 4, 2) - e(1:Qdof)*Dw(1:Qdof, 1, 2) &
    !      - (u(1:Qdof) * Dw(1:Qdof, 2, 2) + v(1:Qdof) * Dw(1:Qdof, 3, 2)) &
    !      + (u(1:Qdof) * u(1:Qdof) + v(1:Qdof) *v(1:Qdof) ) *  Dw(1:Qdof, 1, 2) )

    !write(*,*) '===== Dw(1,1,:) ======='
    !write(*,*) Dw(1,1,:)
    !write(*,*) '--------------'
    !write(*,*) R_s(1,2,:)
    !write(*,*) '================='

    deallocate( oRe)

  end subroutine Set_R_s_Eul


  
  !> compute matrices 4x4 K_sk, s,k=1,2 for Eulerian model, mass diffusion
  !> in integ nodes
  subroutine Set_K_sk_Eul(ndimL, nbDim, iRe, Qdof, w, Dw, Re_1, K_sk, xi)
    integer, intent(in) :: ndimL, nbDim, iRe, Qdof
    real, dimension(1:Qdof, 1:ndimL), intent(in):: w !state  w in #Qdof nodes
    real, dimension(1:Qdof, 1:ndimL, 1:nbDim), intent(in):: Dw !state  Dw in #Qdof nodes, not USED
    real, dimension(1:iRe, 1:Qdof), intent(in) :: Re_1        ! inverse of Reynolds number
    !real, intent(in) :: Re_1                     ! inverse of Reynolds number
    real, dimension(1:Qdof,1:nbDim,1:nbDim,1:ndimL,1:ndimL), intent(inout) :: K_sk
    real, dimension(1:Qdof, 1:nbDim), intent(in):: xi ! physical coordinates
    real :: v1, v2, E, one_over_Rew1, gamPr
    integer :: i
    real, dimension(:), allocatable :: oRe

    allocate( oRe(1:Qdof) )

    !initialization
    K_sk(1:Qdof,1:nbDim,1:nbDim,1:ndimL,1:ndimL) = 0.

    oRe(1:Qdof) = state%model%alpha * Re_1(1,1:Qdof) / w(1:Qdof, 1)

    K_sk(1:Qdof, 1 , 1, 1, 1 ) = oRe(1:Qdof)
    K_sk(1:Qdof, 2 , 2, 1, 1 ) = oRe(1:Qdof)

    K_sk(1:Qdof, 1 , 1, 2, 2 ) = oRe(1:Qdof)
    K_sk(1:Qdof, 2 , 2, 2, 2 ) = oRe(1:Qdof)

    K_sk(1:Qdof, 1 , 1, 3, 3 ) = oRe(1:Qdof)
    K_sk(1:Qdof, 2 , 2, 3, 3 ) = oRe(1:Qdof)

    K_sk(1:Qdof, 1 , 1, 4, 4 ) = oRe(1:Qdof)
    K_sk(1:Qdof, 2 , 2, 4, 4 ) = oRe(1:Qdof)

    deallocate (oRe)
    
  end subroutine Set_K_sk_Eul



end module modelEulerian
