 SUBROUTINE fluxnum(wL,wR,vn,flux,itypflux)
!=========================== DOCUMENTATION ============================= 
!DOCDEBh*    PHY/fluxnum  *
!
!    NOM
! fluxnum(wL,wR,vn,flux,itypflux)
!    DESCRIPTION 
!*
! Calcul du flux numerique 
! Calculating the numerical flux 
!*
! cas 0: Rusanov
! cas 1 : Godunov
! Cas 2 : Flux centre
! Cas 3 : VF-Roe 
! 
! 
!    ENTREES / INPUT
! wL   : objet de type "state" Etat de la partie gauche / left state
! wR   : objet de type "state" Etat de la partie droite / right state
! vn   : objet de type "vect_3", normale à la face / normal vector to the face
! itypflux   : type de flux / type of flux
!
!    SORTIES / OUTPUT
! flux : objet de type "vect_nvar", flux numerique / numerical flux
!
!    SOURCE
! Ce logiciel est regi par la licence [CeCILL-B]
! This software is governed by the [CeCILL-B] license
!=========================== DEBUT DES DECLARATIONS ==================== 
!.1-----  Implicit, Use
  use phy_typ
  use phy ,only:  riemisot, svriemann !compflux, vpmax, vitson, energie, prim2bal, riemann_vfroe,phy_typflux
  IMPLICIT NONE
!.2-----  Declaration
  integer                     :: itypflux !! type de flux / type of flux
  type(state), intent(in)     :: wL , wR  !! Etat gauche et droit / left and right state 
  type(vect_3),intent(in)     :: vn  !! Vecteur normal a la face / Normal vector to the face
  type(vect_nvar),intent(out) :: flux !! Flux
  
  real(kind=kind(0.d0))       :: unL,unR,rL,rR,pL,pR
  real(kind=kind(0.d0))       :: hL,hR,h,zL,zR,hL_reco,hR_reco
  real(kind=kind(0.d0))       :: psi,psiR,psiL,r,un,p,xi,alf
  type(vect_3)                :: uL,uR,u
!=========================== DEBUT DU CODE EXECUTABLE ==================
!-------  Initialisations
  flux%v=0.d0
  u%v=0.d0
  uL%v=0.d0
  uR%v=0.d0

  SELECT CASE (itypflux)
  !---------------------------------------------------
      !~CASE (0) ! flux de HLL ou Rusanov
      !~    CALL compflux(wL,vn,fluxL)
      !~    CALL compflux(wR,vn,fluxR)
      !~    SELECT CASE (imodel_phy)
      !~        CASE DEFAULT
      !~            vv=vn%v(1)**2+vn%v(2)**2+vn%v(3)**2-1.d0
      !~            IF (dabs(vv).gt.1.d-10) THEN
      !~               write(*,*) vv
      !~               STOP
      !~            END IF
      !~            vp=0.d0    ! pour plus tard: maillage mobile...
      !~            uL%v(1)=wL%vprim%v(2)
      !~            uL%v(2)=wL%vprim%v(3)
      !~            uL%v(3)=wL%vprim%v(4)
      !~            CALL vitson(wL,cL)
      !~            unL=uL%v(1)*vn%v(1)+uL%v(2)*vn%v(2)+uL%v(3)*vn%v(3)
      !~            uR%v(1)=wR%vprim%v(2)
      !~            uR%v(2)=wR%vprim%v(3)
      !~            uR%v(3)=wR%vprim%v(4)
      !~            CALL vitson(wR,cR)
      !~            unR=uR%v(1)*vn%v(1)+uR%v(2)*vn%v(2)+uR%v(3)*vn%v(3)
      !~            smin=unL-vp-cL
      !~            IF (smin.gt.unR-vp-cR) smin=unR-vp-cR
      !~            smax=unL-vp+cL
      !~            IF (smax.lt.unR-vp+cR) smax=unR-vp+cR
      !~            vmax=dabs(smin)
      !~            IF (vmax.lt.dabs(smax)) vmax=dabs(smax)
      !~    END SELECT
      !~    ! pour rusanov, décommenter ces lignes : indispensable pour le bifluide !!!!!
      !~    smax=vmax
      !~    smin=-vmax
      !~    IF (smin.ge.0.d0) THEN
      !~       flux%v=fluxL%v
      !~    ELSE IF (smax.le.0.d0) THEN
      !~       flux%v=fluxR%v
      !~    ELSE
      !~       flux%v=(smax*fluxL%v-smin*fluxR%v+smin*smax*(wR%vbal%v-wL%vbal%v))/(smax-smin)
      !~    END IF
  !---------------------------------------------------
      CASE (1) ! flux de godunov
          
          ! calculs communs
          xi=0.d0

          SELECT CASE (imodel_phy)
!---------------------------------------------------  
             CASE (2) ! bifluide isotherme non conservatif
              rL=wL%vprim%v(1)
              uL%v(1)=wL%vprim%v(2)
              uL%v(2)=wL%vprim%v(3)
              uL%v(3)=wL%vprim%v(4)
              pL=wL%vprim%v(5)
              unL=uL%v(1)*vn%v(1)+uL%v(2)*vn%v(2)+uL%v(3)*vn%v(3)

              rR=wR%vprim%v(1)
              uR%v(1)=wR%vprim%v(2)
              uR%v(2)=wR%vprim%v(3)
              uR%v(3)=wR%vprim%v(4)
              pR=wR%vprim%v(5)
              unR=uR%v(1)*vn%v(1)+uR%v(2)*vn%v(2)+uR%v(3)*vn%v(3)
                psiL=wL%vbal%v(5)
                psiR=wR%vbal%v(5)
                !if(wL%chi.eq.0.d0.and.wR%chi.ne.0.d0)then
                !  psiR=psiL
                !  rR  =rL
                !endif
                !if(wR%chi.eq.0.d0.and.wL%chi.ne.0.d0)then
                !  psiL=psiR
                !  rL  =rR
                !endif
                 CALL riemisot(rL,unL,psiL,psiL,rR,unR,psiR,psiR,xi,r,un,p,psi,alf)
                 u%v=un*vn%v
                 IF (un.gt.0.d0) THEN
                    u%v=u%v+uL%v-unL*vn%v
                 ELSE
                    u%v=u%v+uR%v-unR*vn%v
                 END IF
                 flux%v(:)=0.d0
                 flux%v(1)=r*un
                 flux%v(2)=r*u%v(1)*un+p*vn%v(1)
                 flux%v(3)=r*u%v(2)*un+p*vn%v(2)
                 flux%v(4)=r*u%v(3)*un+p*vn%v(3)
                 flux%v(6)=(r*un*un/2+r*(1.d0+log(r))*c_phy*c_phy)*un  ! flux d entropie
                 vi_phy = un   ! stockage de la vitesse du contact pour le schema nc
                 flux%v(5)=0.d0   

!---------------------------------------------------
             CASE (1) ! SaintVenant
                 hL=wL%vprim%v(1)
                 uL%v(1)=wL%vprim%v(2)
                 uL%v(2)=wL%vprim%v(3)
                 unL=uL%v(1)*vn%v(1)+uL%v(2)*vn%v(2)
                 zL=wL%vprim%v(5)
        
                 hR=wR%vprim%v(1)
                 uR%v(1)=wR%vprim%v(2)
                 uR%v(2)=wR%vprim%v(3)
                 unR=uR%v(1)*vn%v(1)+uR%v(2)*vn%v(2)
                 zR=wR%vprim%v(5)

!Reconstruction Hydrostatique
                 IF(hL.LE.hzero_phy)hL=0.d0
                 IF(hR.LE.hzero_phy)hR=0.d0
                 hL_reco=max(0.d0,hL+zL-max(zL,zR))
                 hR_reco=max(0.d0,hR+zR-max(zL,zR))
!secu!
                 IF(hL_reco.LE.hzero_phy)THEN
                   hL_reco=0.d0
                   unL=0.d0
                 ENDIF
                 IF(hR_reco.LE.hzero_phy)THEN
                   hR_reco=0.d0
                   unR=0.d0
                 ENDIF

                 xi=0.d0
                 call svriemann(hL_reco,unL,hR_reco,unR,xi,h,un)
!secu!
                 if(h.LE.hzero_phy)then
                  h=0.d0
                  un=0.d0
                 endif

!projection u_riemann sur normale a la face
                 u%v=un*vn%v
!what is it? : construction du vecteur vitesse etoile a partir de la solution 
!scalaire du probleme de riemann et des normales aux faces
                 if (un.gt.0.d0) then
                  u%v(1)=u%v(1)+uL%v(1)-unL*vn%v(1)
                  u%v(2)=u%v(2)+uL%v(2)-unL*vn%v(2)
                 else
                  u%v(1)=u%v(1)+uR%v(1)-unR*vn%v(1)
                  u%v(2)=u%v(2)+uR%v(2)-unR*vn%v(2)
                 end if

                 flux%v(:)=0.d0
                 flux%v(1)=h*un
                 flux%v(2)=h*u%v(1)*un+gpes_phy*vn%v(1)/2*h*h !+ g*vn%v(1)/2.d0*(hL**2-hL_reco**2) voir fluxnoncons
                 flux%v(3)=h*u%v(2)*un+gpes_phy*vn%v(2)/2*h*h !+ g*vn%v(2)/2.d0*(hL**2-hL_reco**2)voir fluxnoncons
                 flux%v(4)=un*(5.d-1*h*(u%v(1)**2+u%v(2)**2)+gpes_phy*h**2+gpes_phy*h*max(zL,zR)) ! flux d entropie
!---------------------------------------------------
             CASE DEFAULT 
                 call print_err('FLUXNUM','Godunov is not implemented for this model')

          END SELECT

!-------------------------------------------------------------------
      !~CASE (2) ! flux centré (instable sauf si integration en temps de rk4)
      !~    wstar%vprim%v=0.5d0*(wL%vprim%v+wR%vprim%v)
      !~    rL=wL%vprim%v(1)
      !~    rR=wR%vprim%v(1)
      !~    IF (rL.lt.rR) THEN
      !~       temp=rL
      !~       rL=rR
      !~       rR=temp
      !~    END IF
      !~    co=rR/rL
      !~    pe=6.
      !~    pm=pe-1.d0
      !~    IF (rR.ne.rL) wstar%vprim%v(1) = (co+1.D0)*rL/2.D0
      !~    CALL prim2bal(wstar)
      !~    CALL compflux(wstar,vn,flux)
      !~    visnu=0.d0
      !~    flux%v=flux%v-visnu*(wR%vbal%v-wL%vbal%v)

!-------------------------------------------------------------------
      !~CASE (3) ! schéma vfroe
      !~    CALL riemann_vfroe(wL,wR,vn,0.d0,wstar)
      !~    CALL compflux(wstar,vn,flux)
  END SELECT
!===========================   FIN DE LA ROUTINE    ====================
END SUBROUTINE fluxnum
!DOCFIN