      program kneedeep
c----67-------------------------------------------------------flowtec-|
c     solves the riemann problem for st.venant equations (flat bottom)
c     analytically by considering the following zones:
c     (e.g. for a rarefaction-shock configuration)
c
c      (L)     (1)        (C)         (2)        (R)
c      ------         
c            \\\\\\                       ----------
c             \\\\\                     |
c              \\\\                     |
c               \\\                     |       ^ state q
c                \\                     |       |
c                 \                     |       |
c                  ---------------------        -----> (x/t)
c----67---------------------------------------------------------------|
c     12.02.99
c     BUGS: -?-
c----67---------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      parameter(npoimx=500)
      dimension x(npoimx),u(npoimx),fh(npoimx)
      logical lvacuum,l1simple,l2simple,lldry,lrdry
      parameter (f23=2.d0/3.d0,f13=1.d0/3.d0,f19=1.d0/9.d0)
      character*20 outfile
      parameter (iunit=12,outfile='riemann.dat')
c----67---------------------------------------------------------------|      
      parameter (g=2.,time=1.d0,x0=2.5d0)
      parameter (nnwtmx=20,small=1.d-6,width=1.2d0)
c----67---------------------------------------------------------------|      

c     /* set the two initial states L/R of diaphragm    */
      write(*,10)g
      write(*,*)'hL= ?'
      read(*,*)fhL
      write(*,*)'uL= ?'
      read(*,*)uL
      write(*,*)'hR= ?'
      read(*,*)fhR
      write(*,*)'uR= ?'
      read(*,*)uR
      cL=dsqrt(g*fhL)
      cR=dsqrt(g*fhR)
      lldry=.FALSE.
      lrdry=.FALSE.
      if(fhl.ne.0.d0)then
         if(fhR.eq.0.d0)then
            write(*,*)'dry bed on the right'
            lrdry=.TRUE.
         endif
         B=fhR/fhL
         C=(uR-uL)/cL
      else
         write(*,*)'dry bed on the left'
         lldry=.TRUE.
         B=0.d0
         C=0.d0
      endif

c     /* check for vacuum (i.e. dry bed)               */
      lvacuum=.FALSE.
      if(2.d0*(1.d0+dsqrt(B)).le.C)lvacuum=.TRUE.
      write(*,*)'appearance of vacuum:',lvacuum
      if(lvacuum)write(*,*)'will set velocity arbitrarily to zero'
c     /* check of which type are the 1- and 2-families */
      l1simple=.FALSE.
      l2simple=.FALSE.
      if(B.gt.0.d0)
     $     check1=f1(dlog(B))*dsqrt(B)
      if(check1.lt.C.or.lrdry)l1simple=.TRUE.
      if(B.gt.0.d0)
     $     check2=f1(-dlog(b))
      if(check2.lt.C.or.lldry)l2simple=.TRUE.
      write(*,*)'1-family:'
      if(l1simple)write(*,*)'simple wave'
      if(.not.l1simple.and..not.lldry)write(*,*)'shock'
      write(*,*)'2-family:'
      if(l2simple)write(*,*)'simple wave'
      if(.not.l2simple.and..not.lrdry)write(*,*)'shock'
      

c     /* solve the 1-family                            */
      if(l1simple)then
         if(lvacuum.or.lrdry)then
            fhC=0.d0
            uc=0.d0
         else
            x1=-2.d0*dlog((1.d0+dsqrt(B))/2.-C/4.d0)
            fhC=fhL*dexp(-x1)
            CC=dsqrt(g*fhc)
            uC=uL+cL*2.d0*(1.d0-dexp(-x1/2.d0))
            write(*,*)'1-rarefaction: hC=',fhC,' uC=',uC
         endif
      else
         if(l2simple)then
            if(lvacuum.or.lldry)then
               fhC=0.d0
               uc=0.d0
            else
               x2=-2.d0*dlog(-C/(dsqrt(B)*4.d0)+1.d0/(2.d0*dsqrt(B))+
     $              .5d0)
               fhC=fhR*dexp(-x2)
               CC=dsqrt(g*fhc)
               uC=uR-cC*2.d0*(dexp(x2/2.d0)-1.d0)
               write(*,*)'2-rarefaction: hC=',fhC,' uC=',uC
            endif
         else
c     /* we have in fact two shocks: need to iterate              */
            z1=1.1d0
            do n=1,nnwtmx
               zold=z1
               z1=zold-funcz1(zold,B,C)/dfuncz1(zold,B)
               if(dabs(z1-zold).le.small)goto 111
            enddo
            write(*,*)'newton iteration not converged',z1,zold
            stop
 111        continue
            fhC=z1*fhL
            cc=dsqrt(fhc*g)
            uC=uL-cL*(z1-1.d0)*dsqrt((z1+1.d0)/(2.d0*z1))
         endif
      endif
c     /* we now know about the center state (c)                    */

c     /* calculate limits of zone (1) and (2)                      */
c     /* note: "positions" x** are actually given in "xi/t",i.e. velocities */
      if(l1simple)then
         xL1=uL-cL
         if(.not.lvacuum.and..not.lrdry)then
            x1C=uC-cC
         else
c     /* limit zone (1) such that vacuum is just reached through    */
c     /* the expansion */
            x1C=uL+2.d0*cL
         endif
      else
c     /* find shock position: calc shock speed from jump cond.       */
c     /* note: in this case, the zone has zero width                 */
         if(.not.lldry)then
            z=fhC/fhL
            write(*,*)'1-shock: hL/hC=',z
            sigma=(uC*z-uL)/(z-1.d0)
            xl1=sigma
            x1C=sigma
         else
c     /* leave this bound undefined */
         endif
      endif
      if(l2simple)then
         x2R=uR+cR
         if(.not.lvacuum.and..not.lldry)then
            xC2=uC+cC
         else
            xC2=uR-2.d0*cR
         endif
      else
         if(.not.lrdry)then
            z=fhR/fhC
            write(*,*)'2-shock: hR/hC=',z
            sigma=(uR*z-uC)/(z-1.d0)
            x2R=sigma
            xC2=sigma
         else
c     /* since nothing happens in the dry zone, extend it a bit from (1)(C) */
            x2R=x1C+dabs(.1d0*x1c)
            xC2=x2R
         endif
      endif
c     /* get back to the case of a dry bed on the left */
      if(lldry)then
         x1C=xC2-dabs(.1d0*xC2)
         xl1=x1c
      endif

c     /* set a discrete grid, a bit wider than interesting zone */
      flx=(x2R-xL1)
      foverlap=(width-1.d0)/2.d0
      do i=1,npoimx
         x(i)=(xL1-foverlap*flx)+
     $        dfloat(i-1)*flx*width/dfloat(npoimx-1)
      enddo

c     /* set the variables at grid points                       */
      do i=1,npoimx
         if(x(i).le.xL1)then
            u(i)=uL
            fh(i)=fhl
         elseif(x(i).le.x1C)then
            u(i)=f23*x(i)+f13*(uL+2.d0*cL)
            fh(i)=f19/g*(uL+2.d0*cL-x(i))**2
         elseif(x(i).le.xC2)then
            u(i)=uC
            fh(i)=fhC
         elseif(x(i).le.x2R)then
            u(i)=f23*x(i)+f13*(uR-2.d0*cR)
            fh(i)=f19/g*(x(i)-uR+2.d0*cR)**2
         else
            u(i)=uR
            fh(i)=fhR
         endif
      enddo

c     /* output                                                  */
      open(iunit,file=outfile)
      write(iunit,112)
      do i=1,npoimx
         if(fh(i).ne.0.d0)then
            froude=u(i)/dsqrt(g*fh(i))
         else
            froude=0.d0
         endif
         write(iunit,*)x(i),(x(i)+x0)*time,fh(i),u(i),fh(i)*u(i),
     $        dsqrt(fh(i)*g),
     $        u(i)+2.d0*dsqrt(fh(i)*g),u(i)-2.d0*dsqrt(fh(i)*g),
     $        froude
      enddo
      close(iunit)

c     /* format statements */
 10   format('Riemann problem:   (L)   |   (R) ; gravitational ',
     $     'accel. = ',e8.3)
 112  format('#1:x/t 2:x 3:h 4:u 5:h*u 6:c 7:u+2c 8:u-2c 9:Froude')

c     /* finalize                                                 */
      stop
      end
c----67---------------------------------------------------------------|      
      double precision function f1(x)
      implicit double precision (a-h,o-z)
c
      if(x.lt.0.d0)then
         f1=-(dexp(-x)-1.d0)*dsqrt((dexp(-x)+1)/(2.d0*dexp(-x)))
      elseif(x.gt.0.d0)then
         f1=2.d0*(1.d0-dexp(-x/2.d0))
      else
         f1=0.d0
      endif
      end
c----67---------------------------------------------------------------|      
      double precision function f2(x)
      implicit double precision (a-h,o-z)
c
      if(x.lt.0.d0)then
         f2=(dexp(x)-1.d0)*dsqrt((dexp(x)+1.d0)/(2.d0*dexp(x)))
      elseif(x.gt.0.d0)then
         f2=2.d0*(dexp(x/2.d0)-1.d0)
      else
         f2=0.d0
      endif
      end
c----67---------------------------------------------------------------|      
      double precision function funcz1(z,B,C)
      implicit double precision (a-h,o-z)
c
      funcz1= -(z-1.D0)*dsqrt(2.D0)*dsqrt((z+1.D0)/z)/2.D0-
     $     (z/B-1.D0)*dsqrt(2.D0)*dsqrt((z/B+1.D0)/z*B)*
     $     dsqrt(B)/2.D0-C
c      
      end
c----67---------------------------------------------------------------|      
      double precision function dfuncz1(z,B)
      implicit double precision (a-h,o-z)
c
      dfuncz1=  -dsqrt(2.D0)*dsqrt((z+1.D0)/z)/2.D0-
     $     (z-1.D0)*dsqrt(2.D0)/dsqrt((z+1.D0)/z)*
     $     (1.D0/z-(z+1.D0)/z**2)/4.D0-1.D0/dsqrt(B)*
     $     dsqrt(2.D0)*dsqrt((z/B+1.D0)/z*B)/2.D0-
     $     (z/B-1.D0)*dsqrt(2.D0)/dsqrt((z/B+1.D0)/z*B)*
     $     dsqrt(B)*(1.D0/z-(z/B+1.D0)/z**2*B)/4.D0
c      
      end
c----67---------------------------------------------------------------|      
