subroutine riemann ( lmin,lmax,gamma,prgh, urgh, vrgh,
     &                     plft, ulft, vlft, pmid, umid )
c 
c Solve the Riemann shock tube problem for the left and right input states,
C using the Newton interation procedure described in van Leer (1979).
C
      include 'sweepsize.h'
C
      integer l, lmin, lmax, n, nloop
      real gamfac1, gamfac2
      real plft (maxsweep), prgh (maxsweep), 
     &     ulft (maxsweep), urgh (maxsweep), 
     &     vlft (maxsweep), vrgh (maxsweep), 
     &     wlft (maxsweep), wrgh (maxsweep), 
     &     clft (maxsweep), crgh (maxsweep),
     &     zlft (maxsweep), zrgh (maxsweep),
     &     umidl(maxsweep), umidr(maxsweep),
     &     pmid (maxsweep), umid (maxsweep)
c
      gamfac1 = 0.5*(gamma+1.0)/gamma
      gamfac2 = gamma + 1.0
      smallp  = 1.e-35
c
C Input variables are:
C
C   lmin   = zone number of first physical zone
C   lmax   = zone number of first ghost zone on right (lmax=nmax+1)
C   gamma    = equation of state gamma
C   prgh   = pressure state on the right side of the boundary
C   plft   = pressure state on the left side of the boundary
C   urgh   = velocity state on the right side of the boundary
C   ulft   = velocity state on the left side of the boundary
C   vrgh   = density state on the right side of the boundary
C   vlft   = density state on the left side of the boundary
C             (vlft and vrgh are inverted to get the specific volume)
C
C Output variables are:
C 
C   umid   = time averaged velocity at the zone interface
C   pmid   = time averaged pressure at the zone interface
C
C----------------------------------------------------------------------
C
C Obtain first guess for Pmid by assuming Wlft, Wrgh = Clft, Crgh
C
      do l = lmin, lmax
        clft(l) = sqrt(gamma*plft(l)*vlft(l))
        crgh(l) = sqrt(gamma*prgh(l)*vrgh(l))
        vlft(l) = 1./vlft(l)
        vrgh(l) = 1./vrgh(l)
        pmid(l) = prgh(l) - plft(l) - crgh(l)*(urgh(l)-ulft(l))     
        pmid(l) = plft(l) + pmid(l) * clft(l)/(clft(l)+crgh(l))      
        pmid(l) = max(smallp,pmid(l)) 

C Iterate 4 times using Newton's method to converge on correct Pmid
C     -use previous guess for pmid to get wavespeeds: wlft, wrgh
C     -find the slope in the u-P plane for each state: zlft, zrgh
C     -use the wavespeeds and pmid to guess umid on each side: umidl, umidr
C     -project tangents from (pmid,umidl) and (pmid,umidr) to get new pmid
C     -make sure pmid does not fall below floor value for pressure

        wlft (l) = 1.0 + gamfac1*(pmid(l) - plft(l)) / plft(l)   
        wrgh (l) = 1.0 + gamfac1*(pmid(l) - prgh(l)) / prgh(l)
        wlft (l) = clft(l) * sqrt(wlft(l))      
        wrgh (l) = crgh(l) * sqrt(wrgh(l))   
        zlft (l) = 4.0 * vlft(l) * wlft(l) * wlft(l)      
        zrgh (l) = 4.0 * vrgh(l) * wrgh(l) * wrgh(l)   
        zlft (l) = -zlft(l) * wlft(l)/
     &             (zlft(l) - gamfac2*(pmid(l) - plft(l)))      
        zrgh (l) =  zrgh(l) * wrgh(l)/
     &             (zrgh(l) - gamfac2*(pmid(l) - prgh(l)))     
        umidl(l) = ulft(l) - (pmid(l) - plft(l)) / wlft(l)      
        umidr(l) = urgh(l) + (pmid(l) - prgh(l)) / wrgh(l)   
        pmid (l) = pmid(l) + (umidr(l) - umidl(l))
     &                     * (zlft (l) * zrgh (l))     
     &                     / (zrgh (l) - zlft (l))   
        pmid (l) = max(smallp,pmid(l))

        wlft (l) = 1.0 + gamfac1*(pmid(l) - plft(l)) / plft(l)   
        wrgh (l) = 1.0 + gamfac1*(pmid(l) - prgh(l)) / prgh(l)
        wlft (l) = clft(l) * sqrt(wlft(l))      
        wrgh (l) = crgh(l) * sqrt(wrgh(l))   
        zlft (l) = 4.0 * vlft(l) * wlft(l) * wlft(l)      
        zrgh (l) = 4.0 * vrgh(l) * wrgh(l) * wrgh(l)   
        zlft (l) = -zlft(l) * wlft(l)/
     &             (zlft(l) - gamfac2*(pmid(l) - plft(l)))      
        zrgh (l) =  zrgh(l) * wrgh(l)/
     &             (zrgh(l) - gamfac2*(pmid(l) - prgh(l)))     
        umidl(l) = ulft(l) - (pmid(l) - plft(l)) / wlft(l)      
        umidr(l) = urgh(l) + (pmid(l) - prgh(l)) / wrgh(l)   
        pmid (l) = pmid(l) + (umidr(l) - umidl(l))
     &                     * (zlft (l) * zrgh (l))     
     &                     / (zrgh (l) - zlft (l))   
        pmid (l) = max(smallp,pmid(l))

        wlft (l) = 1.0 + gamfac1*(pmid(l) - plft(l)) / plft(l)   
        wrgh (l) = 1.0 + gamfac1*(pmid(l) - prgh(l)) / prgh(l)
        wlft (l) = clft(l) * sqrt(wlft(l))      
        wrgh (l) = crgh(l) * sqrt(wrgh(l))   
        zlft (l) = 4.0 * vlft(l) * wlft(l) * wlft(l)      
        zrgh (l) = 4.0 * vrgh(l) * wrgh(l) * wrgh(l)   
        zlft (l) = -zlft(l) * wlft(l)/
     &             (zlft(l) - gamfac2*(pmid(l) - plft(l)))      
        zrgh (l) =  zrgh(l) * wrgh(l)/
     &             (zrgh(l) - gamfac2*(pmid(l) - prgh(l)))     
        umidl(l) = ulft(l) - (pmid(l) - plft(l)) / wlft(l)      
        umidr(l) = urgh(l) + (pmid(l) - prgh(l)) / wrgh(l)   
        pmid (l) = pmid(l) + (umidr(l) - umidl(l))
     &                     * (zlft (l) * zrgh (l))     
     &                     / (zrgh (l) - zlft (l))   
        pmid (l) = max(smallp,pmid(l))

        wlft (l) = 1.0 + gamfac1*(pmid(l) - plft(l)) / plft(l)   
        wrgh (l) = 1.0 + gamfac1*(pmid(l) - prgh(l)) / prgh(l)
        wlft (l) = clft(l) * sqrt(wlft(l))      
        wrgh (l) = crgh(l) * sqrt(wrgh(l))   
        zlft (l) = 4.0 * vlft(l) * wlft(l) * wlft(l)      
        zrgh (l) = 4.0 * vrgh(l) * wrgh(l) * wrgh(l)   
        zlft (l) = -zlft(l) * wlft(l)/
     &             (zlft(l) - gamfac2*(pmid(l) - plft(l)))      
        zrgh (l) =  zrgh(l) * wrgh(l)/
     &             (zrgh(l) - gamfac2*(pmid(l) - prgh(l)))     
        umidl(l) = ulft(l) - (pmid(l) - plft(l)) / wlft(l)      
        umidr(l) = urgh(l) + (pmid(l) - prgh(l)) / wrgh(l)   
        pmid (l) = pmid(l) + (umidr(l) - umidl(l))
     &                     * (zlft (l) * zrgh (l))     
     &                     / (zrgh (l) - zlft (l))   
        pmid (l) = max(smallp,pmid(l))

C Calculate umid by averaging umidl, umidr based on new pmid

        umidl(l) = ulft(l) - (pmid(l) - plft(l)) / wlft(l)      
        umidr(l) = urgh(l) + (pmid(l) - prgh(l)) / wrgh(l)   
        umid (l) = 0.5*(umidl(l) + umidr(l)) 

      enddo
 
      return     
      end