c =========================================================
      subroutine rp1(maxmx,meqn,mwaves,mbc,mx,ql,qr,auxl,auxr,
     &		       wave,s,amdq,apdq)
c =========================================================
c
c     # solve Riemann problems for the 1D shallow water equations
c     #   (h)_t + (u h)_x = 0
c     #   (h u)_t + ( u^2 h + 0.5 grav*h^2 )_x = 0
c
c     # using an exact Riemann solver with the possibility to handle
c     # the dry-bed case h=0 in the data or in the solution.
c
c     # On input, ql contains the state vector at the left edge of each cell
c     #           qr contains the state vector at the right edge of each cell
c     # On output, wave contains the waves,
c     #            s the speeds,
c     #            amdq the  left-going flux difference  A^- \Delta q
c     #            apdq the right-going flux difference  A^+ \Delta q
c
c     # Note that the i'th Riemann problem has left state qr(i-1,:)
c     #                                    and right state ql(i,:)
c     # From the basic clawpack routine step1, rp is called with ql = qr = q.
c
c     # Here meqn=mwaves=2 should be passed from the calling routine

c     # Written by James Rossmanith and Randy LeVeque
c     # Based on the procedure in Section 13.10 of 
c     #   "Finite Volume Methods for Hyperbolic Problems" by R. J. LeVeque
c
c     # NOTATION:  For p = 1 or 2:
c     #  iwavep = 1 if p-wave is a shock
c     #         = 2 if p-wave is a rarefaction
c     #  if p-wave is a shock:
c     #    sl(p) = sr(p) = shock speed 
c     #  if p-wave is a rarefaction:
c     #    sl(p) < sr(p) are the speeds of left and right edge 

c
      implicit double precision (a-h,o-z)
      dimension   ql(1-mbc:maxmx+mbc, meqn)
      dimension   qr(1-mbc:maxmx+mbc, meqn)
      dimension    s(1-mbc:maxmx+mbc, mwaves)
      dimension wave(1-mbc:maxmx+mbc, meqn, mwaves)
      dimension amdq(1-mbc:maxmx+mbc, meqn)
      dimension apdq(1-mbc:maxmx+mbc, meqn)
      dimension auxl(1-mbc:maxmx+mbc, *)
      dimension auxr(1-mbc:maxmx+mbc, *)
      dimension sl(1:2)
      dimension sr(1:2)
      dimension fl(1:2), fr(1:2), fstar(1:2)
c
c     # local storage
c     ---------------
      dimension delta(2)
      common /comrp/ grav

      logical verbose
      data verbose /.false./  
c     # set verbose to true only for debugging -- gives lots of output!
c

      MaxIters = 20  !# for Newton's method
      TOL = 1.0d-14  !# convergence criterion for Newton
      critH = 1.d-20  !# Critical value of h: set to zero below this.


      do 30 i=2-mbc,mx+mbc

c        # Define left and right states for each i
         hl = qr(i-1,1)         ! left state height
         hr = ql(i,1)           ! right sate height
         hul = qr(i-1,2)        ! left state momentum
         hur = ql(i,2)          ! right state momentum

         if (verbose) then
            write(6,*) '  '
            write(6,*) 'i = ', i
            write(6,*) 'hl, hul:', hl,hul
            write(6,*) 'hr, hur:', hr,hur
            endif

c        # check for negative depth and reset to zero:
         if (hl .lt. 0.d0) then
            write(6,*) ' *** resetting hl to zero from ',hl
            hl = 0.d0
            endif
         if (hr .lt. 0.d0) then
            write(6,*) ' *** resetting hr to zero from ',hr
            hr = 0.d0
            endif

c        # Compute ul and ur, and check for dry states in data:
c        ------------------------------------------------------

         if (hl.gt.critH) then
             ul = hul/hl
           else
c            # dry state to left
             hl = 0.d0
             ul = 0.d0
             iwave1 = 1  !# 1-wave is a zero-strength shock
             iwave2 = 2  !# 2-wave is a rarefaction
           endif
c
         if (hr.gt.critH) then
             ur = hur/hr
           else
c            # dry state to right
             hr = 0.d0
             ur = 0.d0
             iwave1 = 2  !# 1-wave is a rarefaction
             iwave2 = 1  !# 2-wave is a zero-strength shock
             if (hl.eq.0.d0) then
c               # if both states are dry, both waves are zero-strength shocks
                iwave1 = 1 
                endif
           endif
c
c        # wave speeds:
         cl = dsqrt(grav*hl)    !# wave speed in left state
         cr = dsqrt(grav*hr)    !# wave speed in right state

c
c        ---------------------------------------
c        STEP #1: Compute intermediate height hm 
c        ---------------------------------------

         if (hl.eq.0.d0 .or .hr.eq.0.d0) then
c            # if h=0 in one state or both, then hm = 0:
             hm = 0.d0
c            # um is determined by connecting to the wet state with a
c            # rarefaction, or um=0 if both states are dry:
             um = ur + ul - 2.d0*(dsqrt(grav*hr) - dsqrt(grav*hl))
             cm = 0.d0

c            # jump to Step 2:
             go to 200
             endif

c
c        # check for dry state appearing in solution:
c        --------------------------------------------
c
c        # uml is velocity at right edge of a 1-rarefaction connected to
c        # ql and extended to hm = 0:
         uml = ul + 2.d0*cl
c        # umr is velocity at left edge of a 2-rarefaction connected to
c        # qr and extended to hm = 0:
         umr = ur - 2.d0*cr

         if (uml .le. umr) then
c           # dry state appears in solution
            if (verbose) write(6,*) 'dry state appears'
            hm = 0.d0
            cm = 0.d0
c
c           # both waves are rarefactions:
            iwave1 = 2
            iwave2 = 2
            sl(1) = ul - cl
            sr(1) = uml - cm
            sl(2) = umr + cm
            sr(2) = ur + cr

c           # um is not unique in this case, but only hm*um is used later
c           # and this will be zero since hm=0, so might as well set um=0:
            um = 0.d0

c           # jump to Step 3:
            go to 300
            endif


c        # Solve Riemann problem for hm, unless hm=0 has already been found
c        ------------------------------------------------------------------

c        # initial guess for hm based on two rarefaction waves, from (13.55)
c        # in "Finite Volume Methods for Hyperbolic Problems": 
         h = (1.d0/(16.d0*grav))*(ul - ur + 2.d0*(cl + cr))**2
c
c
         if (h.le.hl .and. h.le.hr) then
c              # two rarefaction waves is correct
               if (verbose) write(6,*) 'two rarefactions'
               iwave1 = 2
               iwave2 = 2
               hm = h
               um = ul + 2.d0*(dsqrt(grav*hl) - dsqrt(grav*hm))
               cm = dsqrt(grav*hm)
               go to 200
            endif
c
c
c        # Two rarefaction waves is not correct.  Use hm computed
c        # above as initial guess for Newton itertion to find correct hm,
c        # by solving phi_l(h) - phi_r(h) = 0 as described in Section 13.10.
c        # dphil and dphir denote derivatives phi'_l(h) and phi'_r(h)
c
c
c
c           # Newton iteration:
c
            if (verbose) write(6,*) 'beginning Newton, h = ',h
            do 50 iter=1,MaxIters

c              # left wave
               if (h.gt.hl) then
                  iwave1 = 1
                  sqterm = dsqrt(0.5d0*grav*(1.d0/h + 1.d0/hl))
                  phil = ul - (h-hl)*sqterm
                  dphil = -sqterm + grav*(h-hl) / (4.d0*h*h*sqterm)
               else
                  iwave1 = 2
                  phil = ul + 2.d0*(dsqrt(grav*hl) - dsqrt(grav*h))
                  dphil = -dsqrt(grav/h)
               endif

c              # right wave
               if (h.gt.hr) then
                  iwave2 = 1
                  sqterm = dsqrt(0.5d0*grav*(1.d0/h + 1.d0/hr))
                  phir = ur + (h-hr)*sqterm
                  dphir = sqterm - grav*(h-hr) / (4.d0*h*h*sqterm)
               else
                  iwave2 = 2
                  phir = ur - 2.d0*(dsqrt(grav*hr) - dsqrt(grav*h))
                  dphir = dsqrt(grav/h)
               endif
c
c              # Newton update:
               deltah = (phil-phir) / (dphil-dphir)
               h = h - deltah
               if (verbose) write(6,*) 'iter,h:',iter,h
               if (dabs(deltah) .lt. TOL) go to 100
   50       continue
c
            write(6,599) iter
  599       format('nonconvergence of Newton iteration, iterations = ',
     &            i4)
            stop
c
  100       continue
c
            if (verbose) then
               write(6,601) iter
  601          format('Newton iteration converged in ',i4,' iterations')
               write(6,*) ' '
            endif
c
c        # After convergence of Newton, we should have phil = phir = um,
c        # the intermediate velocity in the Riemann solution,
c        # and h is the value of hm, the intermediate depth.
c
         hm = h
         um = phil
         cm = dsqrt(grav*hm)
c
         if (dabs(phil-phir) .gt. 1d-4) then
            write(6,*) 'phil and phir do not agree:',phil,phir
            write(6,*) 'dphil, dphir:',dphil,dphir
            write(6,*) 'delta h:',deltah
            stop
            endif

  200   continue
c
c        ---------------------------------------------------------
c        STEP #2: Compute wave speeds
c        ---------------------------------------------------------
c
            if (iwave1.eq.1) then
c              # 1-wave is a shock:
               if (dabs(hl-hm).le.1d-10) then
                  s1 = ul - cl
               else
                  s1 = (hl*ul - hm*um)/(hl-hm)
               endif
               sl(1) = s1
               sr(1) = s1
            else
c              # 1-wave is a rarefaction:
               sl(1) = ul - cl
               sr(1) = um - cm
            endif
c
            if (iwave2.eq.1) then
c              # 2-wave is a shock:
               if (dabs(hr-hm).le.1d-10) then
                  s2 = ur + cr
               else
                  s2 = (hr*ur - hm*um)/(hr-hm)
               endif
               sl(2) = s2
               sr(2) = s2
            else
c              # 2-wave is a rarefaction:
               sl(2) = um + cm
               sr(2) = ur + cr
            endif
c
  300    continue
c
c        ----------------------------------------------------------------
c        STEP #3: Set hstar, ustar to Godunov values along x/t = 0
c        ----------------------------------------------------------------
c
         if (sl(1).lt.0.d0 .and. sr(1).gt.0.d0) then
c           # 1-wave: transonic rarefaction
            if (verbose) write(6,*) '1-wave is transonic'
            hstar = (1.d0/(9.d0*grav))*(ul + 2.d0*dsqrt(grav*hl))**2
            ustar = ul - 2.d0*( dsqrt(grav*hstar) - dsqrt(grav*hl) )
         elseif (sl(2).lt.0.d0 .and. sr(2).gt.0.d0) then
c           # 2-wave: transonic rarefaction
            if (verbose) write(6,*) '2-wave is transonic'
            hstar = (1.d0/(9.d0*grav))*(ur - 2.d0*dsqrt(grav*hr))**2
            ustar = ur + 2.d0*( dsqrt(grav*hstar) - dsqrt(grav*hr) )
         elseif (sl(1).gt.0.d0) then
c           # fully supersonic and rightgoing
            if (verbose) write(6,*) 'fully supersonic rightgoing'
            hstar = hl
            ustar = ul
         elseif (sr(2).lt.0.d0) then
c           # fully supersonic and leftgoing
            if (verbose) write(6,*) 'fully supersonic leftgoing'
            hstar = hr
            ustar = ur
         else
c           # 1-wave is leftgoing and
c           # 2-wave is rightgoing
            if (verbose) write(6,*) 'qstar = qm'
            hstar = hm
            ustar = um
         endif
         cstar = dsqrt(grav*hstar)
c
  400    continue
c
c        ----------------------------------------------
c        STEP #4: Compute flux fstar and fluctuations
c        ----------------------------------------------
c
c        # amdq = SUM s*wave   over left-going waves
c        # apdq = SUM s*wave   over right-going waves
c
c        # Flux f at ql, qr, and qstar
         fl(1) = ul*hl
         fr(1) = ur*hr
         fstar(1) = ustar*hstar
         fl(2) = hl*ul**2 + 0.5d0*grav*hl**2
         fr(2) = hr*ur**2 + 0.5d0*grav*hr**2
         fstar(2) = hstar*ustar**2 + 0.5d0*grav*hstar**2
c
c        # fluctuations:
         apdq(i,1) = fr(1) - fstar(1)
         apdq(i,2) = fr(2) - fstar(2)
         amdq(i,1) = fstar(1) - fl(1)
         amdq(i,2) = fstar(2) - fl(2)
c
c
c        ---------------------------------------
c        STEP #5: Compute waves and their speeds
c        ---------------------------------------

         if (verbose) then
             write(6,*) 'sl(1), sr(1): ', sl(1),sr(1)
             write(6,*) 'sl(2), sr(2): ', sl(2),sr(2)
             endif

c        # 1-wave
         wave(i,1,1) = hm - hl
         wave(i,2,1) = hm*um - hl*ul
         s(i,1) = 0.5d0*( sl(1) + sr(1) )

c        # 2-wave
         wave(i,1,2) = hr - hm
         wave(i,2,2) = ur*hr - um*hm
         s(i,2) = 0.5d0*( sl(2) + sr(2) )

         if (verbose) then
            write(6,*) 'hstar, ustar:',hstar,ustar
            write(6,*) 'amdq:',amdq(i,1),amdq(i,2)
            write(6,*) 'apdq:',apdq(i,1),apdq(i,2)
            write(6,*) 's1 and 1-wave:',s(i,1),wave(i,1,1),wave(i,2,1)
            write(6,*) 's2 and 2-wave:',s(i,2),wave(i,1,2),wave(i,2,2)
            endif
c
   30    continue
c
      return
      end

