subroutine evolve( ngeom, umid, pmid )
C
C Use umid and pmid from Riemann solver to update velocity, density,
C and total energy.
C Physical zones are from nmin to nmax.  Zone boundary numbers run from
C nmin to nmax+1
C
      include 'global.h'
      include 'sweep.h'
      real umid(maxsweep), pmid(maxsweep), amid(maxsweep), 
     &     uold(maxsweep), xa1 (maxsweep),dvol1(maxsweep),
     &     upmid(maxsweep), dmu(maxsweep),dtbdm(maxsweep),
     &     grav0(maxsweep), grav1(maxsweep), xa2(maxsweep),
     &     fict0(maxsweep), fict1(maxsweep), xa3(maxsweep)
      integer ngeom
      real third
      third = 1./3.
      rndoff = 1.e-09
C
C------------------------------------------------------------------------
C
C grid position evolution
C
      do 10 n = nmin, nmax + 1
         dmu(n)   = rho(n) * dvol(n)
         dtbdm(n) = dt / dmu(n)
         xa1(n)   = xa(n)
         dvol1(n) = dvol(n)
         delxa    = dt * umid(n) / radius
         if(ngeom.eq.5) delxa = delxa/stheta
         xa (n)   = xa(n) + delxa
         upmid(n) = umid(n) * pmid(n)
 10   continue
      do 20 n = nmin, nmax+1
         xa2(n)   = xa1(n) + 0.5*dx(n)
         dx (n)   = xa(n+1) - xa(n)
         xa3(n)   = xa (n) + 0.5*dx(n)
 20   continue
C
C Calculate forces using coordinates at t and at t+dt, note that
C   fictitious forces at t+dt depend on updated velocity, but we ignore this
C
      call forces(xa2,u,v,w,grav0,fict0)
      call forces(xa3,u,v,w,grav1,fict1)
C
C Calculate dvolume and average area based on geometry of sweep
C
      if(ngeom.eq.0) then
         do 22 n = nmin, nmax+1
            dvol(n) = dx(n)
            amid(n) = 1.
 22      continue
      else if(ngeom.eq.1) then
         do 23 n = nmin, nmax+1
            dvol(n) = dx(n)*(xa(n)+.5*dx(n))
            amid(n) = .5*(xa(n) + xa1(n))
 23      continue
      else if(ngeom.eq.2) then
         do 24 n = nmin, nmax+1
            dvol(n) = dx(n)*(xa(n)*(xa(n)+dx(n))+dx(n)*dx(n)*third)
            amid(n) = (xa(n)-xa1(n))*(third*(xa(n)-xa1(n))+xa1(n))
     &                 +xa1(n)**2
 24      continue
      else if(ngeom.eq.3) then
         do 25 n = nmin, nmax+1
            dvol(n) = dx(n)*radius
            amid(n) = 1.
 25      continue
      else if(ngeom.eq.4) then
         do 26 n = nmin, nmax+1
            dvol(n) = (cos(xa(n))-cos(xa(n+1)))*radius
            dtheta  = xa(n) - xa1(n)
            if(dtheta .eq. 0.0) then
              amid(n) = sin(xa(n))
            else
              amid(n) = (cos(xa1(n))-cos(xa(n)))/dtheta
            endif
 26      continue
      else if(ngeom.eq.5) then
         do 27 n = nmin, nmax+1
            dvol(n) = dx(n) * radius * stheta
            amid(n) = 1.
 27      continue
      else
         write(*,*) 'Geometry', ngeom, ' not implemented.'
         stop
      endif
C                 
      do 30 n = nmin, nmax
C
C density evolution. lagrangian code, so all we have to do is watch the
C change in the geometry.
C
        rho(n) = rho(n) * ( dvol1(n) / dvol(n) )
        rho(n) = max(rho(n),smallr)
C
C velocity evolution due to pressure acceleration and forces.
C
         uold (n) = u(n)
         u(n) = u(n) - dtbdm(n) * (pmid(n+1)-pmid(n))
     &                     * 0.5 * (amid(n+1)+amid(n))
     &          +.5*dt*(grav0(n)+fict0(n)+grav1(n)+fict1(n))
C
C total energy evolution
C
         enrg(n) = enrg(n) - dtbdm(n)*
     &            (amid(n+1)*upmid(n+1) - amid(n)*upmid(n))
     &          + .5*dt*(uold(n)*grav0(n) + u(n)*grav1(n))
c         pold  = p(n)
         p(n)  = rho(n)*gamm*(enrg(n)-.5*(u(n)**2+v(n)**2+w(n)**2))
c         dlogp = abs(pold-p(n))/pold
c         if(dlogp.lt.rndoff) p(n) = pold
         p(n)  = max(p(n),smallp)
C
 30   continue
C
      return
      end