svriemann.f90

### TABLE OF CONTENTS

- 1. PHY/svriemann

## PHY/svriemann [ Modules ]

[ Top ] [ Modules ]

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