      program godunov
c----67-------------------------------------------------------flowtec-|
c     solves the flat-bottom st.venant equations numerically in
c     one dimension (dam-break problem) by a finite-volume 
c     method using:
c     ischeme           fluxes
c     1                 godunov  (exact riemann solver)
c     3                 roe      (lin. riemann solver, flux-diff-split)
c     2                 gallouet (lin. riemann solver, var-diff-split)
c----67-------------------------------------------------------flowtec-|
c     12.05.99
c     BUGS: -?-
c----67---------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      parameter(nmax=400,g=2.d0)
      dimension q(nmax,2),rhs(nmax,2),flux(2)
      parameter(cfl=.4d0)
c----67---------------------------------------------------------------|
      parameter(hll=1.d0,ull=0.d0,hrr=5.d-1,urr=0.d0)
c     /* flx:   length of domain; tfin: final integration time */
c     /* nitmx: max. number of time steps                      */
      parameter(flx=1.d0,tfin=.1d0,nitmx=50000)
      parameter(ischeme=1)
c----67---------------------------------------------------------------|
      time=0.d0
      nit=0
      dx=flx/dfloat(nmax)
      do i=1,nmax
         rhs(i,1)=0.d0
         rhs(i,2)=0.d0
      enddo
c     /* set the initial state (diaphragm at center)*/
      do i=1,nmax/2
         q(i,1)=hll
         q(i,2)=ull*hll
      enddo
      do i=nmax/2+1,nmax
         q(i,1)=hrr
         q(i,2)=urr*hrr
      enddo

c     /* main time loop */
 11   continue
      nit=nit+1

         dt=1.d35

c     /* flux balance for each cell-interface */
         do 20 icell=1,nmax-1

c     /* solve interfacial riemann problem analytically */
            if(ischeme.eq.1)then

               call riemann(q(icell,1),q(icell,2)/q(icell,1),
     $              q(icell+1,1),q(icell+1,2)/q(icell+1,1),
     $              g,
     $              fhc,uc,flamx)
               
            elseif(ischeme.eq.2)then
               
               call gallouet(q(icell,1),q(icell,2),
     $              q(icell+1,1),q(icell+1,2),
     $              g,
     $              fhc,uc,flamx)

            elseif(ischeme.eq.3)then

               call roeflux(q(icell,1),q(icell,2),
     $              q(icell+1,1),q(icell+1,2),
     $              g,
     $              flux,flamx)
            else
               write(*,*)'small problem here...',ischeme
            endif

c     /* calculate maximum GLOBAL time step */
            dt=dmin1(dt,dx*cfl/flamx)
            if(dt.eq.0.d0.or.flamx.eq.0.d0)then
               write(*,*)'icell =',icell,dx,cfl,flamx
               stop
            endif
            
c     /* accumulate rhs fluxes */
            if (ischeme.eq.1.or.ischeme.eq.2)then
               do ic=1,2
                  add=dflux(fhc,uc,g,ic)
                  rhs(icell,ic)=rhs(icell,ic)+add
                  rhs(icell+1,ic)=rhs(icell+1,ic)-add
               enddo
            else
               do ic=1,2
                  rhs(icell,ic)=rhs(icell,ic)+flux(ic)
                  rhs(icell+1,ic)=rhs(icell+1,ic)-flux(ic)
               enddo
            endif

 20         continue

            time=time+dt
            write(*,*)'time = ',time, ' dt= ',dt, ' tfin = ',tfin
            if(time.gt.tfin.or.nit.gt.nitmx)then
               call output(q,time,flx,nmax,g)
               goto 10
            endif

c     /* update the variables */
            do ival=2,nmax-1
               do ic=1,2
                  q(ival,ic)=q(ival,ic)-dt/dx*rhs(ival,ic)
                  rhs(ival,ic)=0.d0
               enddo
            enddo

c     /* check for non-physical values */
            do i=1,nmax
               if(q(i,1).lt.0d0)then
                  WRITE(*,*)'negative height: ',i,q(i,1),q(i,2)
                  stop
               endif
            enddo
            
            goto 11
 10      continue
      stop 
      end
c----67-------------------------------------------------------flowtec-|
      subroutine output(q,time,flx,ni,g)
      implicit double precision (a-h,o-z)
      dimension q(ni,2)
      open(12,file='sol.dat')
      write(12,100)time,ni
      do i=1,ni
         write(12,101)dfloat(i-1)/dfloat(ni-1)*flx,
     $        q(i,1),q(i,2)/q(i,1),q(i,2),sqrt(g*q(i,1))
      enddo
      close(12)
 100  format('# 1:x 2:h 3:u 4:h*u 5:c',x,2(E12.5,x))
 101  format(1x,5(E12.5,x))
      return
      end
c----67-------------------------------------------------------flowtec-|
      subroutine riemann(fhl,ul,fhr,ur,g,fhc,uc,flamx)
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)
      logical lvacuum,l1simple,l2simple,lldry,lrdry
      parameter (f23=2.d0/3.d0,f13=1.d0/3.d0,f19=1.d0/9.d0)
c----67---------------------------------------------------------------|      
c      parameter (g=2.,time=1.d0,x0=2.5d0)
      parameter (nnwtmx=20,small=1.d-6,width=1.2d0)
c----67---------------------------------------------------------------|      
      lrdry=.FALSE.
      lldry=.FALSE.
      cL=dsqrt(g*fhL)
      cR=dsqrt(g*fhR)
      if(fhl.ne.0.d0)then
         if(fhR.eq.0.d0)then
c            write(*,*)'dry bed on the right'
            lrdry=.TRUE.
         endif
         B=fhR/fhL
         C=(uR-uL)/cL
      else
c         write(*,*)'dry bed on the left'
         lldry=.TRUE.
         B=0.d0
         C=0.d0
      endif
c     /* check for equal states on left and right */
      if(fhr.eq.fhl.and.ur.eq.ul)then
         fhc=fhr
         uc=ur
         flamx=dmax1(dabs(ur+sqrt(g*fhr)),dabs(ur-sqrt(g*fhr)),
     $        dabs(ul+sqrt(g*fhl)),dabs(ul-sqrt(g*fhl)))
         goto 112
      endif

c     /* check for vacuum (i.e. dry bed)               */
      lvacuum=.FALSE.
      if(2.d0*(1.d0+dsqrt(B)).le.C)lvacuum=.TRUE.
c      write(*,*)'appearance of vacuum:',lvacuum
c      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.
c      write(*,*)fhl,ul,fhr,ur
c      write(*,*)'1-family:'
c      if(l1simple)write(*,*)'simple wave'
c      if(.not.l1simple.and..not.lldry)write(*,*)'shock'
c      write(*,*)'2-family:'
c      if(l2simple)write(*,*)'simple wave'
c      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))
c            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)
c               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
c            write(*,*)'1-shock: hL/hC=',z
            if(dabs(z-1.d0).gt.small)then
               sigma=(uC*z-uL)/(z-1.d0)
            else
               sigma=ul-sqrt(g*fhl)
            endif
            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
c            write(*,*)'2-shock: hR/hC=',z
            if(dabs(z-1.d0).gt.small)then
               sigma=(uR*z-uC)/(z-1.d0)
            else
               sigma=ur+sqrt(g*fhr)
            endif
            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     /* compute maximum wave speed from limiting velocities xL1,x2R */
      flamx=dmax1(dabs(xL1),dabs(x2R))

c     /* finalize                                                 */
 112  continue 
      return
      end
c----67-------------------------------------------------------flowtec-|
      subroutine gallouet(fhl,ql,fhr,qr,g,fhc,uc,flamx)
c----67-------------------------------------------------------flowtec-|
c     solves the riemann problem for st.venant equations (flat bottom)
c     in its linearized form, cf Buffard, Gallouet, Herard, CRAcadSci,
c     t.326,Serie I,p.385-390, 1998
c----67---------------------------------------------------------------|
c     12.02.99
c     BUGS: -?-
c----67---------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      logical lvacuum
      parameter (f23=2.d0/3.d0,f13=1.d0/3.d0,f19=1.d0/9.d0)
      parameter (small=1.d-6)
c----67---------------------------------------------------------------|      
      if(fhl.ne.0.d0)then
         ul=ql/fhl
      else
         ul=0.d0
      endif
      if(fhr.ne.0.d0)then
         ur=qr/fhr
      else
         ur=0.d0
      endif
      cl=dsqrt(g*fhl)
      cr=dsqrt(g*fhr)
      um=(ul+ur)/2.d0
      fhm=(fhr+fhl)/2.d0
      cm=(cl+cr)/2.d0

c     /* check for existence of dry bottom */
      if((ul-ur).ge.2.d0*(cl+cr))then
         lvacuum=.true.
         uc=0.d0
         cc=0.d0
         flamx=small
         goto 111
      endif

c     /* intermediate state in non-vacuum case */
      uc=(um-(cr-cl))
      cc=(cr+cl)/4.d0*(2.d0-(ur-ul)/(cr+cl))
      flamx=dmax1(dabs(um+cm),dabs(um-cm))

 111  continue
      fhc=cc*cc/g
c
      return
      end
c----67---------------------------------------------------------------|      
      subroutine roeflux(fhl,ul,fhr,ur,g,rflux,flamx)
c----67-------------------------------------------------------flowtec-|
c     solves the linearized riemann problem for st.venant equations 
c     (flat bottom) by using Roe's method (J.Comp.Physics, vol.43,1981)
c     * fluxes are passed upon exit in vector rflux(1:2)
c----67---------------------------------------------------------------|
c     21.05.99
c     BUGS: -?-
c----67---------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      dimension a(2,2),cflux(2),uflux(2),diffv(2),rflux(2)
      logical lvacuum
      parameter (f23=2.d0/3.d0,f13=1.d0/3.d0,f19=1.d0/9.d0)
      parameter (small=1.d-6)
c----67---------------------------------------------------------------|      
      cl=dsqrt(g*fhl)
      cr=dsqrt(g*fhr)

c     /* roe's average for shallow water */
      fhc=(fhl+fhr)/2.d0
      uc=(dsqrt(fhr)*ur+dsqrt(fhl)*ul)/(dsqrt(fhr)+dsqrt(fhl))
      cc=dsqrt(g*fhc)

c     /* diagonalization |A|=R*|lambda|*R1 */
      call aroe(fhc,uc,cc,a)

c     /* variable difference vector */
      diffv(1)=fhr-fhl
      diffv(2)=fhr*ur-fhl*ul

c     /* centered flux */
      do i=1,2
         cflux(i)=dflux(fhr,ur,g,i)+dflux(fhl,ul,g,i)
      enddo

c     /* evaluate roe's flux: f=(fl+fr-|a|*(qr-ql))/2 */
      do i=1,2
         rflux(i)=cflux(i)
         do j=1,2
            rflux(i)=rflux(i)-a(i,j)*diffv(j)
         enddo
         rflux(i)=rflux(i)/2.d0
      enddo

c     /* estimate the maximum signal velocity */
      flamx=dmax1(dabs(uc+cc),dabs(uc-cc))
c
      return
      end
c----67---------------------------------------------------------------|      
      subroutine aroe(hc,uc,cc,a)
c----67-------------------------------------------------------flowtec-|
c     calculates roe's matrix
c----67---------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      dimension a(2,2),r(2,2),r1(2,2),flambda(2)
c----67---------------------------------------------------------------|      
c     /* diagonalization matrix r */
      r(1,1)=1.d0
      r(1,2)=1.d0
      r(2,1)=uc-cc
      r(2,2)=uc+cc

c     /* inverse r1 */
      r1(1,1)=(uc+cc)/(2.d0*cc)
      r1(1,2)=-1.d0/(2.d0*cc)
      r1(2,1)=(cc-uc)/(2.d0*cc)
      r1(2,2)=1.d0/(2.d0*cc)

c     /* eigenvalues */
      flambda(1)=uc-cc
      flambda(2)=uc+cc

c     /* matrix*diag*matrix multiply */
      do i=1,2
         do j=1,2
            a(i,j)=0.d0
            do k=1,2
               a(i,j)=a(i,j)+r(i,k)*dabs(flambda(k))*r1(k,j)
            enddo
         enddo
      enddo
c
      return
      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---------------------------------------------------------------|      
      double precision function dflux(h,u,g,i)
      implicit double precision (a-h,o-z)
c
      if (i.eq.1)then
         dflux=h*u
      elseif(i.eq.2)then
         dflux=u*u*h+g*h*h/2.d0
      else
         write(*,*)'problem here:..',i
         stop
      endif
c      
      end
c----67---------------------------------------------------------------|      
