 
subroutine svriemann(hg,ug,hd,ud,xi,h,u)
!=========================== DOCUMENTATION ============================= 
!DOCDEBh*    PHY/svriemann  *
!
!    NOM
! svriemann(hg,ug,hd,ud,xi,h,u)
!    DESCRIPTION  
!*
! Resolution du probleme de Riemann pour Saint Venant
! Solving the Riemann problem for Saint Venant
!*
! 
!    ENTREES / INPUT
!  hg  : hauteur d'eau gauche / water level left
!  ug  : vitesse gauche / velocity left
!  hd  : hauteur d'eau droit / water level right
!  ud  : vitesse droit / velocity right
!  xi  : x/t
!
!    SORTIES / OUTPUT
!  h   : hauteur d'eau interface / water level at the interface
!  u   : vitesse interface / interface velocity
!
!    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
  implicit none
!.2-----  Declaration
  real(kind=kind(0.d0)),intent(in)  ::  hg !! hauteur d eau gauche / water level left 
  real(kind=kind(0.d0)),intent(in)  ::  hd !! vitesse gauche / velocity left 
  real(kind=kind(0.d0)),intent(in)  ::  ug !! hauteur d eau droit / water level right 
  real(kind=kind(0.d0)),intent(in)  ::  ud !! vitesse droit / velocity right 
  real(kind=kind(0.d0)),intent(in)  ::  xi !! x/t
  real(kind=kind(0.d0)),intent(out) ::  h  !! hauteur d eau a l interface / water level at the interface 
  real(kind=kind(0.d0)),intent(out) ::  u  !! vitesse a l interface / interface velocity 

  real(kind=kind(0.d0))            :: g,eps
  integer,parameter                :: maxiter=100
  real(kind=kind(0.d0))            :: cg,cd,cstar,ustar,hstar,dh,res
  real(kind=kind(0.d0))            :: xlambda1m,xlambda1p,xlambda2m,xlambda2p
  integer                          :: config,iter
!=========================== DEBUT DU CODE EXECUTABLE ==================

eps=hzero_phy
g=gpes_phy

! test type configuration
IF ((ud-ug) .GE. (2.d0*dsqrt(g)*(dsqrt(hg)+dsqrt(hd)))) THEN
! dry zone centrale
  config = 1
ELSEIF (hg .LE. eps .and. hd .GT. eps) THEN
! dry zone gauche
  config = 2
ELSEIF (hd .LE. eps .and. hg .GT. eps) THEN
! dry zone droite
  config = 3
ELSEIF (hg .LE. eps .and. hd .LE. eps) THEN
!dry zone a droite et a gauche
  config = 4
ELSE
! pas de dry zone
  config = 0
ENDIF



select case(config)

! pas de dry zone
  case(0)
!
!-------  Initialisations
!
    hstar=.5*(hg+hd)
    
! newton-raphson pour calculer h*
      iter=0
      res=1.d0
      dh=1.d0
      do while((dabs(dh).gt.eps).or.(dabs(res).gt.eps))
           iter=iter+1
          ! calcul du residu (on veut que res=0)
           res= ud-ug+(hstar-hd)*zsv(hstar,hd)+(hstar-hg)*zsv(hstar,hg)
          ! calcul de l accroissement
           dh=zsv(hstar,hd) +(hstar-hd)*dzsv(hstar,hd) &
            + zsv(hstar,hg) +(hstar-hg)*dzsv(hstar,hg) 
           dh=-res/dh
          ! mise a jour
           hstar=hstar+dh
          ! test d arret
           if (iter.ge.maxiter) then
           stop 'Newton iteration max atteinte!'
           endif
      enddo
  
! calcul de u*
      ustar=ud+(hstar-hd)*zsv(hstar,hd)
  
! celerite du "son" a gauche droite et milieu
      cg=dsqrt(g*hg)
      cd=dsqrt(g*hd)
      cstar=dsqrt(g*hstar)
  
! calcul de lambda(i,plus) et lambda(i,moins)
   ! onde 1
      if (hstar.gt.hg) then
           !     choc
           xlambda1m=(hg*ug-hstar*ustar)/(hg-hstar)
           xlambda1p=xlambda1m
      else
           !     detente
           xlambda1m=ug-cg
           xlambda1p=ustar-cstar
      endif
   ! onde 2
      if (hstar.gt.hd) then
           !     choc
           xlambda2m=(hd*ud-hstar*ustar)/(hd-hstar)
           xlambda2p=xlambda2m
      else
           !     detente
           xlambda2p=ud+cd
           xlambda2m=ustar+cstar
      endif
  
  ! h_interface et u_interface
      if (xi .lt. xlambda1m) then
           ! etat gauche
           h=hg
           u=ug
      elseif ((xi .ge. xlambda1m).and.(xi .lt. xlambda1p)) then
           !     1-detente
           h=(2.d0*cg+ug-xi)**2/(9.d0*g)
           u=(2.d0*(cg+xi)+ug)/3.d0 
      elseif ((xi .ge. xlambda1p).and.(xi .lt. xlambda2m)) then
           !     etat milieu
           h=hstar
           u=ustar
      elseif ((xi .ge. xlambda2m).and.(xi .lt. xlambda2p)) then
           !     2-detente
           h=(2.d0*cd-ud+xi)**2/(9.d0*g) 
           u=(-2.d0*cd+2.d0*xi+ud)/3.d0  
      elseif (xi .ge. xlambda2p) then
           ! etat droit
           h=hd
           u=ud
      endif
      
      
    ! dry zone centrale    
  case(1)
  ! etat milieu
      hstar=0.d0
      ustar=0.d0
  ! celerite du "son" a gauche droite et milieu
      cg=dsqrt(g*hg)
      cd=dsqrt(g*hd)
      cstar=0.d0
  ! calcul de lambda(i,plus) et lambda(i,moins)
      ! onde 1
          xlambda1m=ug-cg
          xlambda1p=ug+2.d0*cg
      ! onde 2
          xlambda2m=ud-2.d0*cd
          xlambda2p=ud+cd
  ! h_interface et u_interface
      if (xi .lt. xlambda1m) then
           ! etat gauche
           h=hg
           u=ug
      elseif ((xi .ge. xlambda1m).and.(xi .lt. xlambda1p)) then
           !     1-detente
           h=(2.d0*cg+ug-xi)**2/(9.d0*g);
           u=(2.d0*cg+ug+2*xi)/3.d0
      elseif ((xi .ge. xlambda1p).and.(xi .lt. xlambda2m)) then
           !     etat milieu
           h=hstar
           u=ustar
      elseif ((xi .ge. xlambda2m).and.(xi .lt. xlambda2p)) then
           !     2-detente
           h=(2.d0*cd-ud+xi)**2/(9.d0*g);
           u=(-2.d0*cd+ud+2*xi)/3.d0
      elseif (0.d0 .ge. xlambda2p) then
           ! etat droit
           h=hd
           u=ud
      endif
    
    
    ! dry zone gauche
  case(2)
       cd=dsqrt(g*hd)
   ! calcul de lambda(i,plus) et lambda(i,moins)
       ! onde 2
           xlambda2m=ud-2.d0*cd
           xlambda2p=ud+cd
   ! h_interface et u_interface
       if (xi .lt. xlambda2m) then
            !     etat milieu
            h=0.d0
            u=0.d0
       elseif ((xi .ge. xlambda2m).and.(xi .lt. xlambda2p)) then
            !     2-detente
            h=(2.d0*cd-ud+xi)**2/(9.d0*g);
            u=(-2.d0*cd+ud+2*xi)/3.d0
       elseif (xi .ge. xlambda2p) then
            ! etat droit
            h=hd
            u=ud
       endif            
       
       
    ! dry zone droite
  case(3)
      cg=dsqrt(g*hg)
  ! calcul de lambda(i,plus) et lambda(i,moins)
      ! onde 1
      xlambda1m=ug-cg
      xlambda1p=ug+2*cg
  ! h_interface et u_interface
      if (xi .lt. xlambda1m) then
           ! etat gauche
           h=hg
           u=ug
      elseif ((xi .ge. xlambda1m).and.(xi .lt. xlambda1p)) then
           !     1-detente
           h=(2*cg+ug-xi)**2/(9.d0*g);
           u=(2*cg+ug+2*xi)/3.d0
      elseif (xi .ge. xlambda1p) then
           !     etat milieu
           h=0.d0
           u=0.d0
      endif


    ! dry zone a droite et a gauche
  case(4)
       h=0.d0
       u=0.d0

    
end select
            
return
!
!===========================   FIN DE LA ROUTINE    ====================
contains
    function zsv(h1,h2)
      implicit none
      real(kind=kind(0.d0))     :: gg,h1,h2,zsv,zsol
    
      gg=9.81d0
      zsol=0.d0
    
      if (h1.gt.h2) then
         zsol=dsqrt(gg/2.d0)*dsqrt((h1+h2)/(h1*h2))
      else
         zsol=2.d0*dsqrt(gg)/(dsqrt(h1)+dsqrt(h2))
      endif
    
      zsv=zsol
    end function zsv

    function dzsv(h1,h2)
      implicit none
      real(kind=kind(0.d0)) :: gg,h1,h2,dzsv,dzsol
    
      gg=9.81d0
      dzsol=0.d0
    
      if (h1.gt.h2) then
         dzsol= dsqrt(2.d0)*dsqrt(gg)/4.d0*sqrt(h1*h2/(h1+h2))* &
              (1.d0/(h1*h2)-(h1+h2)/(h1**2*h2))
      else
         dzsol=-dsqrt(gg)/(dsqrt(h1)+dsqrt(h2))**2/dsqrt(h1)
      endif
    
      dzsv=dzsol
    end function dzsv
end subroutine svriemann
!DOCFIN           
