subroutine states( ngeom, pl, ul, rl, p6, u6, r6, dp, du, 
     &                   dr, plft, ulft, rlft, prgh, urgh, rrgh)
C
C This subroutine takes the values of rho, u, and P at the left hand
C side of the zone, the change accross the zone, and the parabolic 
C coefficients, p6, u6, and rho6, and computes the left and right states
C (integrated over the charachteristics) for each variable for input 
C to the Riemann solver.
C
      include 'global.h'
      include 'sweep.h'
C
      integer n, ngeom
      real plft(maxsweep), ulft(maxsweep), rlft(maxsweep), 
     &     prgh(maxsweep), urgh(maxsweep), rrgh(maxsweep),
     &     dp  (maxsweep), du  (maxsweep), dr  (maxsweep),
     &     pl  (maxsweep), ul  (maxsweep), rl  (maxsweep),
     &     p6  (maxsweep), u6  (maxsweep), r6  (maxsweep),
     &     areal(maxsweep), darea(maxsweep), area6(maxsweep),
     &     alft(maxsweep), argh(maxsweep), ur(maxsweep),
     &     grav(maxsweep), fict(maxsweep)
C
      real Cdtdx(maxsweep), fCdtdx(maxsweep), fothd, hdt
C
C Input variables are:
C
C    ngeom =  geometry of sweep
C    pl    =  pressure at left edge of zone
C    dp    =  pressure slope
C    p6    =  pressure parabola coeffecient
C    ul    =  velocity at left edge of zone
C    du    =  velocity slope
C    u6    =  velocity parabola coeffecient
C    rl    =  density at left edge of zone
C    dr    =  density slope
C    r6    =  density parabola coeffecient
C
C Output variables are:
C
C    plft =  average pressure in the + wave state (left of interface)
C    prgh =  average pressure in the - wave state (right of interface)
C    ulft =  average velocity in the + wave state (left of interface)
C    urgh =  average velocity in the - wave state (right of interface)
C    rlft =  average density in the  + wave state (left of interface)
C    rrgh =  average density in the  - wave state (right of interface)
C
C--------------------------------------------------------------------------
      fothd = 4./3.
      hdt   = 0.5*dt
C
C Calculate the domain of dependence along space coordinate,
C C*dt, by multiplying the Eulerian sound speed in 
C each zone by dt.  Divide by two to save flops later.
C
C If angular coordinate, then divide out by radius to get physical dimensions
C
      if(ngeom.lt.3) then
       do 10 n = nmin, nmax
         Cdtdx (n) = sqrt(gam*p(n)/rho(n))/dx(n)*hdt
         fCdtdx(n) = 1. - fothd*Cdtdx(n)
 10    continue
      else if(ngeom.eq.3 .or. ngeom.eq.4) then
       do 11 n = nmin, nmax
         Cdtdx (n) = sqrt(gam*p(n)/rho(n))/(dx(n)*radius)*hdt
         fCdtdx(n) = 1. - fothd*Cdtdx(n)
 11    continue
      else if(ngeom.eq.5) then
       do 12 n = nmin, nmax
         Cdtdx (n) = sqrt(gam*p(n)/rho(n))*hdt/(dx(n)*radius*stheta)
         fCdtdx(n) = 1. - fothd*Cdtdx(n)
 12    continue
      endif
C
C Obtain averages of rho, u, and P over the domain (+/-)Cdt
C                    lft is the + wave on the left  side of the boundary
C                    rgh is the - wave on the right side of the boundary
C
C Include gravitational and ficticious forces
C    (NOTE: we are approximating ulft,urgh with ul for the ficticious forces)
C
      ul(nmax+1) = ul(nmax) + du(nmax)
      call forces(xa,ul,v,w,grav,fict)
      do 20 n = nmin, nmax
         k = n+1
         plft(k) = pl(n)+dp(n)-Cdtdx(n)*(dp(n)-fCdtdx(n)*p6(n))
         ulft(k) = ul(n)+du(n)-Cdtdx(n)*(du(n)-fCdtdx(n)*u6(n))
         rlft(k) = rl(n)+dr(n)-Cdtdx(n)*(dr(n)-fCdtdx(n)*r6(n))
         plft(k) = max(smallp,plft(k))
         ulft(k) = ulft(k) + hdt*(grav(k)+fict(k))
 20   continue
      do 25 n = nmin, nmax
         prgh(n) = pl(n) + Cdtdx(n)*(dp(n)+fCdtdx(n)*p6(n))
         urgh(n) = ul(n) + Cdtdx(n)*(du(n)+fCdtdx(n)*u6(n))
         rrgh(n) = rl(n) + Cdtdx(n)*(dr(n)+fCdtdx(n)*r6(n))
         prgh(n) = max(smallp,prgh(n))
         urgh(n) = urgh(n) + hdt*(grav(n)+fict(n))
 25   continue
C
C Add geometry corrections to input states if diverging geometry
C
C###  NOT WORKING CORRECTLY  ###
C      if(ngeom.eq.1) then
C         do 30 n = nmin, nmax+1
C            areal(n) = xa(n)
C            darea(n) = xa(n+1) - xa(n)
C            area6(n) = 0.0
C 30      continue
C      else if(ngeom.eq.2) then
C         do 31 n = nmin, nmax+1
C            areal(n) = xa(n)**2
C            darea(n) = dx(n)*(2*xa(n)+dx(n))
C            area6(n) = -dx(n)**2
C 31      continue
C      endif
CC
C      if(ngeom.eq.1 .or. ngeom.eq.2) then
C      do 40 n = nmin, nmax
C         k = n+1
C         alft(k) = areal(k)-Cdtdx(n)*(darea(n)-fCdtdx(n)*area6(n))
C         argh(n) = areal(n)+Cdtdx(n)*(darea(n)+fCdtdx(n)*area6(n))
C         ur(k) = ul(n) + du(n)
C         ulft(k) = ur(k) + 2*(alft(k)*ulft(k) - areal(k)*ur(k))
C     &                      /(alft(k)         + areal(k) )
C         urgh(n) = ul(n) + 2*(argh(n)*urgh(n) - areal(n)*ul(n))
C     &                      /(argh(n)         + areal(n) )
C 40   continue
C      endif
C
      return
      end