subroutine ppm( nleft, nright, ngeom, paracoff )
C
C Using the 1D arrays of rho, u, and P, perform a 1D lagrangian
C hydrodynamics evolution:
C   - obtain parabolic interpolations of rho, u, P
C   - compute input states from these interpolations for Riemann problem
C   - call the Riemann solver to find the time averages umid, pmid
C   - evolve the finite difference equations to get updated
C      values of rho, u, and E
C 
      include 'global.h'
      include 'sweep.h'
C
      integer nleft, nright, ngeom, iflat
      real dr  (maxsweep), du  (maxsweep), dp  (maxsweep), 
     &     r6  (maxsweep), u6  (maxsweep), p6  (maxsweep),
     &     rl  (maxsweep), ul  (maxsweep), pl  (maxsweep), 
     &     rrgh(maxsweep), urgh(maxsweep), prgh(maxsweep),
     &     rlft(maxsweep), ulft(maxsweep), plft(maxsweep), 
     &     flat(maxsweep), umid(maxsweep), pmid(maxsweep)
      real paracoff( 9, maxsweep )
C
C-----------------------------------------------------------------------
C
C Flattening is used in the interpolation procedure if iflat = 1
      iflat = 1
C
      if(iflat .eq. 1) call flatten( flat )
C
      call parabola(nmin,nmax,paracoff,p  ,dp,p6,pl,flat,iflat,0,ngeom)
      call parabola(nmin,nmax,paracoff,rho,dr,r6,rl,flat,iflat,0,ngeom)
      call parabola(nmin,nmax,paracoff,u  ,du,u6,ul,flat,iflat,0,ngeom)
C
      call states( ngeom, pl, ul, rl, p6, u6, r6, dp, du,
     &             dr, plft, ulft, rlft, prgh, urgh, rrgh )
C
C Apply boundary conditions for states at the array boundaries
C
      if (nleft .eq. 0) then
C                                        reflecting boundary condition
         plft (nmin)  =  prgh (nmin)
         ulft (nmin)  = -urgh (nmin)
         rlft (nmin)  =  rrgh (nmin)
      else if (nleft .eq. 1) then
C                                        in/out boundary condition
         plft (nmin)  =  prgh (nmin)
         ulft (nmin)  =  urgh (nmin)
         rlft (nmin)  =  rrgh (nmin)
      else if (nleft .eq. 2) then
C                                        fixed inflow boundary condition
         plft (nmin)  =  pinflo
         ulft (nmin)  =  uinflo
         rlft (nmin)  =  dinflo
      else if (nleft .eq. 3) then
C                                         periodic boundary
         plft (nmin)  =  plft (nmax + 1)
         ulft (nmin)  =  ulft (nmax + 1)
         rlft (nmin)  =  rlft (nmax + 1)
      else if (nleft .eq. 5) then
C                                         fixed inflow/transverse velocity
         plft (nmin)  =  pinflo
         ulft (nmin)  =  vinflo
         rlft (nmin)  =  dinflo
      else
         write(*,*) 'Boundary condition ',nleft,' not implemented'
         stop
      endif
      if (nright .eq. 0) then
C                                        reflecting boundary condition
         prgh (nmax+1) =  plft (nmax+1)
         urgh (nmax+1) = -ulft (nmax+1)
         rrgh (nmax+1) =  rlft (nmax+1)
      else if (nright .eq. 1) then
C                                        in/out boundary condition
         prgh (nmax+1) =  plft (nmax+1)
         urgh (nmax+1) =  ulft (nmax+1)
         rrgh (nmax+1) =  rlft (nmax+1)
      else if (nright .eq. 2) then
C                                        fixed inflow boundary condition
         prgh (nmax+1) =  pinflo
         urgh (nmax+1) =  uinflo
         rrgh (nmax+1) =  dinflo
      else if (nright .eq. 3) then
C                                         periodic boundary
         prgh (nmax+1) =  prgh (nmin)
         urgh (nmax+1) =  urgh (nmin)
         rrgh (nmax+1) =  rrgh (nmin)
      else if (nright .eq. 5) then
C                                         powerlaw boundary
         prgh (nmax+1) =  pinflo
         urgh (nmax+1) =  vinflo
         rrgh (nmax+1) =  dinflo
      else
         write(*,*) 'Boundary condition ',nright,' not implemented'
         stop
      endif
C
C Call the Riemann solver to obtain the zone face averages, umid and pmid
C
      call riemann( nmin, nmax+1, gam, prgh, urgh, rrgh, 
     &                     plft, ulft, rlft, pmid, umid )
C
C force reflecting boundary conditions on umid
C
      if (nleft  .eq. 0) umid(nmin  ) = 0.0
      if (nright .eq. 0) umid(nmax+1) = 0.0
C
C do lagrangian update 
C
      call evolve( ngeom, umid, pmid )
C
      return
      end