!> definition of 2D models which are simulated: Euler, Navier-Stokes
module modelRANS_2e

  use main_data
  use f_mapping
  use mesh_oper
  use blocks_integ
  use model2DNS
!  use model3DNS

  implicit none

  public:: Set_f_s_RANS_2e
  public:: Set_A_s_RANS_2e
  public:: Set_Ppm_RANS_2e

  public:: Set_R_s_RANS_2e
  public:: Set_K_sk_RANS_2e

  public:: Set_S_RANS_2e
  public:: Set_DS_RANS_2e

contains

  !> compute the Euler fluxes \f$ f_s({\bf w}),\quad s=1,2\f$
  subroutine Set_f_s_RANS_2e(ndimL, nbDim, Qdof, w, f_s, xi, ie)
    integer, intent(in) :: Qdof, nbDim, ndimL
    real, dimension(1:Qdof, 1:ndimL), intent(in):: w !state  w in #Qdof nodes
    real, dimension(1:Qdof,1:nbDim,1:ndimL), intent(inout) :: f_s
                                               ! matrices A_s in  -- " --
    real, dimension(1:Qdof,1 :nbDim), intent(in) :: xi
    integer, intent(in) :: ie


    ! laminar terms
    if(nbDim == 2) then
      !might need to take out all +2 becasue for calculate the
      !f_s not need to include turbulence dimension
       call Set_f_s_Euler(nbDim+2, nbDim, Qdof, w(1:Qdof,1:nbDim+2), &
            f_s(1:Qdof, 1:nbDim, 1:nbDim+2), xi, ie )
    elseif(nbDim == 3) then
      stop 'Stopped in Set_f_s_RANS_2e - 3D not implemented'
!       call Set_f_s_Euler3D(nbDim+2, Qdof, w(1:Qdof,1:nbDim+2), &
!            f_s(1:Qdof, 1:nbDim, 1:nbDim+2) )
    endif

    !turbulence terms
    f_s(1:Qdof, 1, 2) = f_s(1:Qdof, 1, 2) + 2./3*w(1:Qdof, 5)
    f_s(1:Qdof, 1, 4) = f_s(1:Qdof, 1, 4) + (w(1:Qdof,2)/w(1:Qdof,1))*(2./3)*w(1:Qdof, 5)
    f_s(1:Qdof, 1, 5) = w(1:Qdof,2)*w(1:Qdof,5)/w(1:Qdof,1)
    f_s(1:Qdof, 1, 6) = w(1:Qdof,2)*w(1:Qdof,6)/w(1:Qdof,1)

    f_s(1:Qdof, 2, 3) = f_s(1:Qdof, 2, 3) + 2./3*w(1:Qdof, 5)
    f_s(1:Qdof, 2, 4) = f_s(1:Qdof, 2, 4) + (w(1:Qdof,3)/w(1:Qdof,1))*(2./3)*w(1:Qdof, 5)
    f_s(1:Qdof, 2, 5) = w(1:Qdof,3)*w(1:Qdof,5)/w(1:Qdof,1)
    f_s(1:Qdof, 2, 6) = w(1:Qdof,3)*w(1:Qdof,6)/w(1:Qdof,1)

  end subroutine Set_f_s_RANS_2e


  !> compute matrices \f$ A_s = \frac{D{\bf f_s}({\bf w})}{D{\bf w}},\quad s=1,2\f$
  !> for Euler fluxes
  subroutine Set_A_s_RANS_2e(ndimL, nbDim, Qdof, w, A_s, xi, ie)
    integer, intent(in) :: ndimL, nbDim, Qdof
    real, dimension(1:Qdof, 1:ndimL), intent(in):: w !state  w in #Qdof nodes
    real, dimension(1:Qdof,1:nbDim,1:ndimL,1:ndimL), intent(inout) :: A_s
                                               ! matrices A_s in  -- " --
    real, dimension(1:Qdof,1 :nbDim), intent(in) :: xi
    integer, intent(in) :: ie

    ! laminar terms
    if(nbDim == 2) then
       call Set_A_s_Euler(nbDim+2, nbDim, Qdof, w(1:Qdof,1:nbDim+2), &
            A_s(1:Qdof, 1:nbDim, 1:nbDim+2, 1:nbDim+2), xi, ie )
    elseif(nbDim == 3) then
     stop 'Stopped in Set_A_s_RANS_2e - 3D not implemented'
     !  call Set_A_s_Euler3D(nbDim+2, Qdof, w(1:Qdof,1:nbDim+2), &
     !       A_s(1:Qdof, 1:nbDim, 1:nbDim+2, 1:nbDim+2) )
    endif

    !turbulence terms
    !Changes in the laminar equations
        A_s(1:Qdof, 1, 4, 1) = A_s(1:Qdof, 1, 4, 1) - 2./3 * &
        (w(1:Qdof,5)*w(1:Qdof,2))/(w(1:Qdof,1)**2)
        A_s(1:Qdof, 1, 4, 2) = A_s(1:Qdof, 1, 4, 2) - 2./3 * w(1:Qdof,5)/w(1:Qdof, 1)

        A_s(1:Qdof, 2, 4, 1) = A_s(1:Qdof, 2, 4, 1) - 2./3 * w(1:Qdof,5)*w(1:Qdof,3)/(w(1:Qdof, 1)**2)
        A_s(1:Qdof, 2, 4, 2) = A_s(1:Qdof, 2, 4, 2) - 2./3 * w(1:Qdof,5)/w(1:Qdof, 1)

    !Rows 5 and 6 below the laminar equations
        A_s(1:Qdof, 1, 5, 1) = - w(1:Qdof,2)*w(1:Qdof,5)/(w(1:Qdof,1)**2)
        A_s(1:Qdof, 1, 5, 1) = w(1:Qdof,5)/w(1:Qdof,1)
        A_s(1:Qdof, 1, 5, 1) = 0
        A_s(1:Qdof, 1, 5, 1) = 0

        A_s(1:Qdof, 1, 6, 1) = - w(1:Qdof,2)*w(1:Qdof,6)/(w(1:Qdof,1)**2)
        A_s(1:Qdof, 1, 6, 1) = w(1:Qdof,6)/w(1:Qdof,1)
        A_s(1:Qdof, 1, 6, 1) = 0
        A_s(1:Qdof, 1, 6, 1) = 0

        A_s(1:Qdof, 2, 5, 1) = - w(1:Qdof,3)*w(1:Qdof,5)/(w(1:Qdof,1)**2)
        A_s(1:Qdof, 2, 5, 1) = 0
        A_s(1:Qdof, 2, 5, 1) = w(1:Qdof,5)/w(1:Qdof,1)
        A_s(1:Qdof, 2, 5, 1) = 0

        A_s(1:Qdof, 1, 6, 1) = - w(1:Qdof,3)*w(1:Qdof,5)/(w(1:Qdof,1)**2)
        A_s(1:Qdof, 1, 6, 1) = 0
        A_s(1:Qdof, 1, 6, 1) = w(1:Qdof,6)/w(1:Qdof,1)
        A_s(1:Qdof, 1, 6, 1) = 0
      !Colums 5 and 6
        A_s(1:Qdof, 1, 1, 5) = 0.
        A_s(1:Qdof, 1, 2, 5) = 2./3
        A_s(1:Qdof, 1, 3, 5) = 0.
        A_s(1:Qdof, 1, 4, 5) = 2./3 * w(1:Qdof,2)/w(1:Qdof,1)
        A_s(1:Qdof, 1, 5, 5) = w(1:Qdof,2)/w(1:Qdof,1)
        A_s(1:Qdof, 1, 6, 5) = 0.
        A_s(1:Qdof, 1, 1:nDimL-1, 6) = 0
        A_s(1:Qdof, 1, 6, 6) = w(1:Qdof,2)/w(1:Qdof,1)

        A_s(1:Qdof, 2, 1, 5) = 0.
        A_s(1:Qdof, 2, 2, 5) = 0.
        A_s(1:Qdof, 2, 3, 5) = 2./3
        A_s(1:Qdof, 2, 4, 5) = 2./3 * w(1:Qdof,3)/w(1:Qdof,1)
        A_s(1:Qdof, 2, 5, 5) = w(1:Qdof,3)/w(1:Qdof,1)
        A_s(1:Qdof, 2, 6, 5) = 0.
        A_s(1:Qdof, 2, 1:nDimL-1, 6) = 0
        A_s(1:Qdof, 2, 6, 6) = w(1:Qdof,3)/w(1:Qdof,1)


  end subroutine Set_A_s_RANS_2e


  !> compute matrices 4x4 K_sk, s,k=1,2 for N.-S. equations
  !> in integ nodes
  subroutine Set_K_sk_RANS_2e(ndimL, nbDim, iRe,  Qdof, w, Dw, Re_1, K_sk, 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, 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

    ! laminar terms
    if(nbDim == 2) then
       call Set_K_sk_NS(nbDim+2, nbDim, iRe, Qdof, w(1:Qdof,1:nbDim+2), &
       Dw(1:Qdof,1:nbDim+2,1:nbDim), Re_1, K_sk(1:Qdof, 1:nbDim, 1:nbDim, &
       1:nbDim+2, 1:nbDim+2),xi )
    elseif(nbDim == 3) then
       call Set_K_sk_NS(nbDim+2, nbDim, iRe, Qdof, w(1:Qdof,1:nbDim+2), Dw(1:Qdof,1:nbDim+2,1:nbDim), Re_1, &
            K_sk(1:Qdof, 1:nbDim, 1:nbDim, 1:nbDim+2, 1:nbDim+2), xi )
    endif


    !turbulence terms
    K_sk(1:Qdof, 1:nbDim, 1:nbDim, nbDim+3:ndimL, :) = 0.
    K_sk(1:Qdof, 1:nbDim, 1:nbDim, :, nbDim+3:ndimL) = 0.

  end subroutine Set_K_sk_RANS_2e


  !> compute viscous fluxes R_s, s=1,2 for N.-S. equations
  !> in integ nodes
  subroutine Set_R_s_RANS_2e(ndimL, nbDim, iRe, Qdof, w, Dw, Re_1, R_s, 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
    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:: u, v, oRe, e, mu_t, kapparans
    real :: Pr, Pr_t, sigmak, sigmaw, mu,cp

    Pr_t = state%model%Pr_t
    Pr = state%model%Pr
    sigmak = state%model%sigmak
    sigmaw = state%model%sigmaw
    mu = 1
    mu_t = w(1:Qdof,5)*w(1:Qdof,1)/w(1:Qdof,6)
    cp = state%model%cp
    kapparans = (mu/Pr +mu_t/Pr_t)*cp

    ! laminar terms
    if(nbDim == 2) then
       call Set_R_s_NS(nbDim+2, nbDim, iRe, Qdof, w(1:Qdof,1:nbDim+2), Dw(1:Qdof,1:nbDim+2, 1:nbDim), &
            Re_1, R_s(1:Qdof, 1:nbDim, 1:nbDim+2), xi )
    elseif(nbDim == 3) then
      stop 'Stopped in Set_Ppm_RANS_2e - 3D not implemented'
      ! call Set_R_s_NS3D(Qdof, w(1:Qdof,1:nbDim+2), Dw(1:Qdof,1:nbDim+2, 1:nbDim), &
      !      Re_1, R_s(1:Qdof, 1:nbDim, 1:nbDim+2) )
    endif

    !turbulence terms
    R_s(1:Qdof, 1, 2) = R_s(1:Qdof, 1, 2) + (2./3)* Pr_t/w(1:Qdof,1)  &
         *(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) = R_s(1:Qdof, 1, 3) + Pr_t/w(1:Qdof,1) &
         *((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) &
         + kapparans/w(1:Qdof,1) &
         * ( Dw(1:Qdof, 4, 1) - e(1:Qdof)*Dw(1:Qdof, 1, 1)) &
         - kapparans*w(1:Qdof,2)/w(1:Qdof,1) * (Dw(1:Qdof, 2, 1) - w(1:Qdof,2)/w(1:Qdof,1)*Dw(1:Qdof, 2, 1))&
         +(Pr + sigmak*Pr_t - state%model%kappa)*(Dw(1:Qdof, 5, 1) - w(1:Qdof,5)/w(1:Qdof,1)*Dw(1:Qdof, 5, 1))
    R_s(1:Qdof, 1, 5) = (Pr + sigmak*Pr_t) &
    *(Dw(1:Qdof, 5, 1) - w(1:Qdof,5)/w(1:Qdof,1)*Dw(1:Qdof, 5, 1))
    R_s(1:Qdof, 1, 6) = (Pr + sigmaw*Pr_t) &
    *(Dw(1:Qdof, 6, 1) - w(1:Qdof,6)/w(1:Qdof,1)*Dw(1:Qdof, 6, 1))

    R_s(1:Qdof, 2, 2) = R_s(1:Qdof, 2, 2) + Pr_t/w(1:Qdof,1) &
    *((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) = R_s(1:Qdof, 2, 3) + (2./3)* Pr_t/w(1:Qdof,1)  &
         *(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) &
         + kapparans/w(1:Qdof,1) &
         * ( Dw(1:Qdof, 4, 2) - e(1:Qdof)*Dw(1:Qdof, 1, 2)) &
         - kapparans*w(1:Qdof,3)/w(1:Qdof,1) * (Dw(1:Qdof, 3, 2) - w(1:Qdof,3)/w(1:Qdof,1)*Dw(1:Qdof, 3, 2))&
         +(Pr + sigmak*Pr_t - state%model%kappa)*(Dw(1:Qdof, 5, 2) - w(1:Qdof,5)/w(1:Qdof,1)*Dw(1:Qdof, 5, 2))
    R_s(1:Qdof, 2, 5) = (Pr + sigmak*Pr_t) &
    *(Dw(1:Qdof, 5, 2) - w(1:Qdof,5)/w(1:Qdof,1)*Dw(1:Qdof, 5, 2))
    R_s(1:Qdof, 2, 6) = (Pr + sigmaw*Pr_t) &
    *(Dw(1:Qdof, 6, 2) - w(1:Qdof,6)/w(1:Qdof,1)*Dw(1:Qdof, 6, 2))

  end subroutine Set_R_s_RANS_2e


  !> compute matrices
  !> \f$ P**{\pm} = \left(\frac{D({\bf f_1}({\bf w})n_1+{\bf f_2}({\bf w})n_2}{D{\bf w}}
  !>  \right)**{\pm}\f$
  !> for the Euler equation
  subroutine Set_Ppm_RANS_2e(ndimL, nbDim, Qdof, w, n, xi, Ppm, one_over_area, elem)
    !!class(mesh), intent(in) :: grid
    external::dgeev
    external::dgetrf
    external::dgetri
    integer, intent(in) :: Qdof, ndimL, nbDim
    real, dimension(1:Qdof, 1:ndimL), intent(in):: w !state  w in #Qdof nodes
    real, intent(in), optional  :: one_over_area
    real, dimension(1:Qdof,1:nbDim,1:nDimL,1:nDimL), intent(inout) :: Ppm
                                               ! matrices Ppm in  -- " --
    real, dimension(1:Qdof, 1:nbDim), intent(in) :: n   ! outer normal
    real, dimension(1:Qdof, 1:nbDim),intent(in) ::  xi                    ! node on the edge?
    type(element), intent(inout), optional :: elem
    double precision, dimension(1:nDimL) :: dp, dm

    double precision, dimension(1:Qdof, 1:nDimL,1:nDimL) :: P
    !double precision, dimension(nDimL) :: WR,WI
    double precision, dimension(nDimL,nDimL) :: t,t1 !Matrix,VR,VL,
    double precision, dimension(1:Qdof) :: u, v, k, wt
    integer, dimension(nDimL) :: Ipiv
    integer ::Info,i,j,ie, iwork
    real :: kappa, kappa1

    real, dimension(:),allocatable :: Work, WR, WI
    real, dimension(:,:),allocatable :: Matrix, VR, VL

    allocate(WR(1:ndimL), WI(1:ndimL))
    allocate(VR(1:ndimL, 1:ndimL), VL(1:ndimL, 1:ndimL), Matrix(1:ndimL, 1:ndimL))
    iwork =   4*ndimL


    u(1:Qdof) = w(1:Qdof,2)/w(1:Qdof,1)
    v(1:Qdof) = w(1:Qdof,3)/w(1:Qdof,1)
    k(1:Qdof) = w(1:Qdof,5)/w(1:Qdof,1)
    wt(1:Qdof) = w(1:Qdof,6)/w(1:Qdof,1)
    kappa = state%model%kappa
    kappa1 = state%model%kappa1

    P(1:Qdof,1,1) = 0
    P(1:Qdof,2,1) = (kappa1/2)*(u(1:Qdof)**2+v(1:Qdof)**2)*n(1:Qdof,1) &
    -u(1:Qdof)*(u(1:Qdof)*n(1:Qdof,1)+v(1:Qdof)*n(1:Qdof,2))
    P(1:Qdof,3,1) = (kappa1/2)*(u(1:Qdof)**2+v(1:Qdof)**2)*n(1:Qdof,2) &
    -v(1:Qdof)*(u(1:Qdof)*n(1:Qdof,1)+v(1:Qdof)*n(1:Qdof,2))
    P(1:Qdof,4,1) = (kappa1*(u(1:Qdof)**2+v(1:Qdof)**2)-kappa*w(1:Qdof,4)/w(1:Qdof,1))* &
    (u(1:Qdof)*n(1:Qdof,1)+v(1:Qdof)*n(1:Qdof,2)) - 2*k(1:Qdof)/(3*w(1:Qdof,1))* &
    (n(1:Qdof,1)+n(1:Qdof,2))
    P(1:Qdof,5,1) = -k(1:Qdof)*(u(1:Qdof)*n(1:Qdof,1)+v(1:Qdof)*n(1:Qdof,2))
    P(1:Qdof,6,1) = -wt(1:Qdof)*(u(1:Qdof)*n(1:Qdof,1)+v(1:Qdof)*n(1:Qdof,2))

    P(1:Qdof,1,2) = n(1:Qdof,1)
    P(1:Qdof,2,2) = -(kappa-2)*u(1:Qdof)*n(1:Qdof,1)+ &
    u(1:Qdof)*n(1:Qdof,1)+v(1:Qdof)*n(1:Qdof,2)
    P(1:Qdof,3,2) = v(1:Qdof)*n(1:Qdof,1) - kappa1*u(1:Qdof)*n(1:Qdof,2)
    P(1:Qdof,4,2) = (kappa*w(1:Qdof,4)/w(1:Qdof,1)- &
    kappa1/2*(u(1:Qdof)**2+v(1:Qdof)**2))*n(1:Qdof,1) - kappa1*u(1:Qdof)* &
    (u(1:Qdof)*n(1:Qdof,1)+v(1:Qdof)*n(1:Qdof,2)) +2./3*k(1:Qdof)*u(1:Qdof)
    P(1:Qdof,5,2) = k(1:Qdof)*n(1:Qdof,1)
    P(1:Qdof,6,2) = wt(1:Qdof)*n(1:Qdof,1)

    P(1:Qdof,1,3) = n(1:Qdof,2)
    P(1:Qdof,2,3) = u(1:Qdof)*n(1:Qdof,2) - kappa1*v(1:Qdof)*n(1:Qdof,1)
    P(1:Qdof,3,3) = -(kappa-2)*v(1:Qdof)*n(1:Qdof,2)+ &
    u(1:Qdof)*n(1:Qdof,1)+v(1:Qdof)*n(1:Qdof,2)
    P(1:Qdof,4,3) = (kappa*w(1:Qdof,4)/w(1:Qdof,1)- &
    kappa1/2*(u(1:Qdof)**2+v(1:Qdof)**2))*n(1:Qdof,2) - kappa1*v(1:Qdof)* &
    (u(1:Qdof)*n(1:Qdof,1)+v(1:Qdof)*n(1:Qdof,2)) +2./3*k(1:Qdof)*v(1:Qdof)
    P(1:Qdof,5,3) = k(1:Qdof)*n(1:Qdof,2)
    P(1:Qdof,6,3) = wt(1:Qdof)*n(1:Qdof,2)

    P(1:Qdof,1,4) = 0
    P(1:Qdof,2,4) = kappa1*n(1:Qdof,1)
    P(1:Qdof,3,4) = kappa1*n(1:Qdof,2)
    P(1:Qdof,4,4) = kappa*(u(1:Qdof)*n(1:Qdof,1)+v(1:Qdof)*n(1:Qdof,2))
    P(1:Qdof,5,4) = 0
    P(1:Qdof,6,4) = 0

    P(1:Qdof,1,5) = 0
    P(1:Qdof,2,5) = (2./3) * n(1:Qdof,1)
    P(1:Qdof,3,5) = (2./3) * n(1:Qdof,2)
    P(1:Qdof,4,5) = (2./3) * (u(1:Qdof)*n(1:Qdof,1)+v(1:Qdof)*n(1:Qdof,2))
    P(1:Qdof,5,5) = u(1:Qdof)*n(1:Qdof,1)+v(1:Qdof)*n(1:Qdof,2)
    P(1:Qdof,6,5) = 0

    P(1:Qdof,1,6) = 0
    P(1:Qdof,2,6) = 0
    P(1:Qdof,3,6) = 0
    P(1:Qdof,4,6) = 0
    P(1:Qdof,5,6) = 0
    P(1:Qdof,6,6) = u(1:Qdof)*n(1:Qdof,1)+v(1:Qdof)*n(1:Qdof,2)

    WR =0.
    WI = 0.
    VL = 0.
    Vr=0.
    iwork = nDimL*50
    allocate(work(1:iwork) , source = 0.)

    do ie =1,Qdof
      Matrix(1:ndimL, 1:ndimL) = P(ie,1:ndimL, 1:ndimL)

      Matrix = 0.
      write(*,*)nDimL
      do i=1,6
        Matrix(i,i) = i
        write(*,'(30es12.4)') Matrix( i, :), WR(i), WI(i)
      enddo
      do i=1,6
        write(*,'(30es12.4)') VL( i, :), VR(i,:)
      enddo
      !iwork = -1
      !call dgeev('N','V',nDimL,Matrix,nDimL,WR,WI,VL,1,VR,nDimL,Work,-1,Info)
      write(*,*)work(1),Info
      !iwork = nDimL*50!min(nDimL*200,idint(work(1)))
      call dgeev('N','N',nDimL,Matrix(1:ndimL, 1:ndimL),nDimL, &
        WR(1:ndimL),WI(1:ndimL),VL(1:ndimL, 1:ndimL),1,VR(1:ndimL, 1:ndimL),nDimL,Work,iwork,Info)
      write(*,*) info
      t = VR
      write(*,'(6es14.6)') WR
      write(*,*) VR
      call DGETRF(nDimL, nDimL, VR, nDimL, Ipiv, info)
      call DGETRI(nDimL, VR, nDimL, ipiv, work, nDimL, info)
      t1 = VR

      dp = WR
      do i=1,nDimL
        dm(i)=0.
        if(dp(i).lt.0.)then
          dm(i)=dp(i)
          dp(i)=0.
        endif
       enddo

       do i=1,ndimL
          do j=1,ndimL
             Ppm(ie,1,i,j) = sum(t(i, 1:nDimL) * dp(1:nDimL) * t1(1:nDimL, j))
             Ppm(ie,2,i,j) = sum(t(i, 1:nDimL) * dm(1:nDimL) * t1(1:nDimL, j))
          enddo
       enddo
    enddo

!    do i=1,ndimL
!       do j=1,ndimL
!          Ppm(ie,1,i,j) = sum(VR(i, 1:ndimL) * dp(1:ndimL) * t1(1:ndimL, j) )
!          Ppm(ie,2,i,j) = sum(VR(i, 1:ndimL) * dm(1:ndimL) * t1(1:ndimL, j) )
!
!       enddo
!    enddo
  deallocate(work)
  end subroutine Set_Ppm_RANS_2e

  subroutine Set_Ppm_RANS_2e_Slip(ndimL, nbDim, e_Qdof, w_ein, n_e, Ppm )
    ! compute matrix Pp on a slip boundary
    integer, intent(in) :: e_Qdof, ndimL, nbDim
    real, dimension(1:e_Qdof, 1:ndimL), intent(in):: w_ein !state  w in #Qdof nodes
    real, dimension(1:e_Qdof,2:nbDim+1,1:ndimL), intent(inout) :: Ppm
                                               ! matrices Ppm in  -- " --
    real, dimension(1:e_Qdof,1:nbDim), intent(in) :: n_e      ! outer normal
    real :: kappa, kappa1
    real :: v(2), vv, nn(2)
    integer :: ie, i, j



    ! laminar terms
    if(nbDim == 2) then
       call Set_Ppm_Euler_Slip(nbDim+2, nbDim, e_Qdof, w_ein(1:e_Qdof,1:nbDim+2), n_e,  &
            Ppm(1:e_Qdof, 2:nbDim+1, 1:nbDim+2)  )
    elseif(nbDim == 3) then
       print*,'To be done in modelRANS_2e.f90 !! '
    !   call Set_Ppm_Euler3D_Slip(nbDim+2, e_Qdof, w_ein(1:e_Qdof,1:nbDim+2), n_e,  &
    !        Ppm(1:e_Qdof, 2:nbDim+1, 1:nbDim+2)  )
    endif

    !turbulence terms
    Ppm(1:e_Qdof, 2:nbDim+1, nbDim+3:ndimL) = 0.

  end subroutine Set_Ppm_RANS_2e_Slip


  !> compute reactive terms S in integ nodes
  subroutine Set_S_RANS_2e(ndimL, nbDim, Qdof, xi, w, Dw, S)
    ! ATTENTION IN THIS CODE IS USED W BUT IT SHOULD BE USED MIN(W,W_R)
    !BUT W_R NOT SPECIFIED IN THE PAPER
    integer, intent(in) :: ndimL, nbDim, Qdof
    real, dimension(1:Qdof, 1:nbDim), intent(in) :: xi!optimal velocity,SOLUTION OF THE EIKONAL EQUATION
    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:Qdof, 1:ndimL), intent(inout) :: S
    real, dimension(1:Qdof) :: u,v
    integer :: i,j,k
    real :: V_rho, rho, rlen, betak,betaw,alphaw
    real :: mu(1:2)
    real :: rho_minimal_allowed , rho_actual, penalty

    u = w(1:Qdof,2)/w(1:Qdof,1)
    v = w(1:Qdof,3)/w(1:Qdof,1)
    betak = state%model%betak
    betaw = state%model%betaw
    alphaw = state%model%alphaw
    mu(1) = state%model%Pr
    mu(2) = state%model%Pr_t

    S(1:Qdof, :) = 0.
    S(1:Qdof,5) = 2./3* (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))-rho*w(1:Qdof,5))*(Dw(1:Qdof, 2, 1) - u(1:Qdof)*Dw(1:Qdof, 1, 1)) &
         + (Dw(1:Qdof, 2, 2) - u(1:Qdof)*Dw(1:Qdof, 1, 2))&
         + (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))&
         + (Dw(1:Qdof, 2, 2) - u(1:Qdof)*Dw(1:Qdof, 1, 2))&
         + (Dw(1:Qdof, 3, 1) - v(1:Qdof)*Dw(1:Qdof, 1, 1))*(Dw(1:Qdof, 3, 1) - v(1:Qdof)*Dw(1:Qdof, 1, 1))&
         + 2./3* (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))-rho*w(1:Qdof,5))*(Dw(1:Qdof, 3, 2) - u(1:Qdof)*Dw(1:Qdof, 1, 2))&
        -betak*rho*w(1:Qdof,5)*exp(w(1:Qdof,6))

    S(1:Qdof,6) = alphaw/w(1:Qdof,5)*&
         (2./3* (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))-rho*w(1:Qdof,5))*(Dw(1:Qdof, 2, 1) - u(1:Qdof)*Dw(1:Qdof, 1, 1)) &
         + (Dw(1:Qdof, 2, 2) - u(1:Qdof)*Dw(1:Qdof, 1, 2))&
         + (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))&
         + (Dw(1:Qdof, 2, 2) - u(1:Qdof)*Dw(1:Qdof, 1, 2))&
         + (Dw(1:Qdof, 3, 1) - v(1:Qdof)*Dw(1:Qdof, 1, 1))*(Dw(1:Qdof, 3, 1) - v(1:Qdof)*Dw(1:Qdof, 1, 1))&
         + 2./3* (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))-rho*w(1:Qdof,5))*(Dw(1:Qdof, 3, 2) - u(1:Qdof)*Dw(1:Qdof, 1, 2)))&
        -betaw*rho*exp(w(1:Qdof,6))+ (mu(1)+state%model%sigmaw*mu(2))*w(1:Qdof,6)

  end subroutine Set_S_RANS_2e

  !> compute derivative of the reactive terms S in integ nodes
  subroutine Set_DS_RANS_2e(ndimL, nbDim, Qdof, xi, w, Dw, DS)
    integer, intent(in) :: ndimL, nbDim, Qdof
    real, dimension(1:Qdof, 1:nbDim), intent(in) :: xi
    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:Qdof, 1:ndimL, 1:ndimL), intent(inout) :: DS
    integer :: i,j,k
    real :: V_rho, rho, rlen, penalty
    real :: mu(1:2)
    real :: rho_minimal_allowed, rho_actual


    DS(1:Qdof, :,:) = 0.


  end subroutine Set_DS_RANS_2e


end module modelRANS_2e
