c
c
c
c =========================================================
      subroutine rp1(maxmx,meqn,mwaves,mbc,mx,ql,qr,auxl,auxr,
     &           fwave,s,amdq,apdq)
c =========================================================
c
c     f-wave solver for the Riemann problems for the 1D gas dynamics 
c     equations with a general equation of state (EOS) using a variant of 
c     the "Suliciu relaxation solver" described in Sec. 2.4.5 of 
c     Francois Bouchut's book "Nonlinear Stability of Finite 
c     Volume methods for Hyperbolic Conservation Laws and Well-Balanced
c     Schemes for Sources" (Birkh\"auser 2004).  This solver works well
c     near vacuum states.  The EOS should be specified in the subroutine
c     pres_psi, which returns a vector pres containing the pressure
c     and cs, the sound speed, in each cell.
c     
c     Note that in this Riemann solver it is assumed that ql=qr and only
c     ql is used.
c
c     The solver described in Bouchut has been modified by RJL and
c     Jeremiah Murphy to work in f-wave form so that source terms can also
c     be incorporated.  [cite paper, to appear].
c
c     The subroutine pres_psi should also return 
c     the source term psi*dx at the left edge of each cell:  psidx(i) is
c     used to modify dp at the interface between cell i-1 and i, where dp
c     is the pressure difference between these cells.
c
c     -------------------------------------------------------
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 fwave(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,*)
      common /comxt/ dtcom,dxcom,tcom
      double precision npoly, nstar
      common /ceos/ grav, gamma, gamma1, sK, npoly, nstar, Crho, ieos

c
c     # local storage
c     ---------------
      parameter (max2 = 10002)  !# assumes at most 10000 grid points with mbc=2
      dimension   pres(-1:max2)
      dimension   cs(-1:max2)
      dimension   psidx(-1:max2)

      if (mx+mbc .gt. max2) then
         write(6,*) '*** Increase max2 in rp1'
         stop
         endif
c
c
c     # compute the pressure pres, sound speed cs, and source term psi*dx:
      call pres_psi(ql,pres,cs,psidx,mx,maxmx,mbc,meqn)
      write(60,601) 1200,pres(1200),cs(1200),psidx(1200)
      write(60,601) 1201,pres(1201),cs(1201),psidx(1201)
 601  format(i5,3e16.6)
c
c
      do i=2-mbc,mx+mbc

c        find primative values to the left and right of the contact
c        discontinuity.
         pl = pres(i-1)
         pr = pres(i)
         rhol = ql(i-1,1)
         rhor = ql(i,1)
         ul = ql(i-1,2)/rhol
         ur = ql(i,2)/rhor

         alpha = (gamma+1.d0)*0.5d0


c        # Jump in pressure is modified by source term correction for
c        # well-balanced method:
         dp = pr - pl - psidx(i)

c        # Useful for debugging:
c        write(40,401) i,rhol,rhor,ul,ur
c        write(40,401) i,pl,pr,psidx(i),dp
c        write(40,*) ' '
c 401    format(i5,4e15.6)

c        # Formulas (2.145, 2.146) from Bouchut:
c        # ccl and ccr are the $\tilde c_l$ and $\tilde c_r$ values there,
c        # which are sound speeds multiplied by rho values

         if (dp.gt.0.d0) then
            cl = cs(i-1) + alpha*max((dp/rhor/cs(i) +ul-ur),0.d0)
            cr = cs(i) + alpha*max((-dp/rhol/cl + ul-ur),0.d0)
         else
            cr = cs(i) + alpha*max((-dp/rhol/cs(i-1) + ul-ur),0.d0)
            cl = cs(i-1) + alpha*max((dp/rhor/cr + ul-ur),0.d0)
         endif
         ccl = rhol*cl
         ccr = rhor*cr

c        # Internal energy:
         eccl = ql(i-1,3)-0.5d0*rhol*ul**2
         eccr = ql(i,3)-0.5d0*rhor*ur**2
c
c        # Enthalpies, not used:
c        hl = (eccl+pl)/rhol
c        hr = (eccr+pr)/rhor

c        # Equation (2.135) from Bouchut:
         ustar = (ul*ccl + ur*ccr - dp)/(ccr + ccl)
         pstar = (pl*ccr + pr*ccl + (ul - ur)*ccr*ccl)
     &        /(ccr+ccl)
         rstarl = 1.d0/
     &        (1.d0/rhol + (ccr*(ur-ul) - dp)/ccl/(ccl+ccr))
         rstarr = 1.d0/
     &        (1.d0/rhor + (ccl*(ur-ul) + dp)/ccr/(ccl+ccr))
c 
c        # From Jeremiah Murphy:
c        pstarl = pstar - ccl*psidx/(ccl+ccr)
c        pstarr = pstar + ccr*psidx/(ccl+ccr)

c        # Equivalent expressions:
         pstarl = pl  + (ccl*dp + (ul-ur)*ccr*ccl)/(ccl+ccr)
         pstarr = pr  - (ccr*dp - (ul-ur)*ccr*ccl)/(ccl+ccr)

         estarl = rstarl*(eccl/rhol + (pstarl**2 - pl**2)/2.d0/ccl**2)
         estarr = rstarr*(eccr/rhor + (pstarr**2 - pr**2)/2.d0/ccr**2)

c
c        # Compute the waves.
c
c        f(qstarl) - f(ql)
         fwave(i,1,1) = rstarl*ustar - rhol*ul
         fwave(i,2,1) = rstarl*ustar**2 + pstarl - rhol*ul**2 - pl
         fwave(i,3,1) = ustar*(estarl + 0.5d0*rstarl*ustar**2 + pstarl)
     &        - ul*(eccl + 0.5d0*rhol*ul**2 + pl)
         s(i,1) = ul - cl
c
c        f(qstarr) - f(qstarl)
         fwave(i,1,2) = (rstarr - rstarl)*ustar
         fwave(i,2,2) = (rstarr - rstarl)*ustar**2 
         fwave(i,3,2) = ustar*(estarr - estarl 
     &        + 0.5d0*(rstarr - rstarl)*ustar**2) 

         s(i,2) = ustar
c
c        f(qr) - f(qstarr)
         fwave(i,1,3) = rhor*ur - rstarr*ustar
         fwave(i,2,3) = rhor*ur**2 + pr - rstarr*ustar**2 - pstarr
         fwave(i,3,3) = ur*(eccr + 0.5d0*rhor*ur**2 + pr)
     &        - ustar*(estarr + 0.5d0*rstarr*ustar**2 + pstarr)
         s(i,3) = ur + cr

c        write(40,401) i,fwave(i,2,1),fwave(i,2,2),fwave(i,2,3)
c        write(40,*) ' '

      enddo
c
c
c
c     # amdq = SUM s*wave   over left-going waves
c     # apdq = SUM s*wave   over right-going waves
c

      stol = 0.d0
      do m=1,meqn
         do i=2-mbc, mx+mbc
            amdq(i,m) = 0.d0
            apdq(i,m) = 0.d0
            do mw=1,mwaves
               if (s(i,mw) .lt. stol) then
                  amdq(i,m) = amdq(i,m) + fwave(i,m,mw)
               else if (s(i,mw) .gt. stol) then
                  apdq(i,m) = apdq(i,m) + fwave(i,m,mw)
               else
                  amdq(i,m) = amdq(i,m) + 0.5d0*fwave(i,m,mw)
                  apdq(i,m) = apdq(i,m) + 0.5d0*fwave(i,m,mw)
               endif
            enddo
         enddo
      enddo

      return
      end
