subroutine parabola( nmin, nmax, para, a, deltaa, a6, al, flat,
     &                     iflat, iskip, ngeom )
C
C Colella and Woodward, JCompPhys 54, 174-201 (1984) eq 1.5-1.8,1.10
C
C parabola calculates the parabolas themselves. call paraset first
C for a given grid-spacing to set the constants, which can be reused
C each time parabola is called.
C
C flatening coefficients are calculated externally in flaten. 
C
      include 'sweepsize.h'
C
      integer n, nmin, nmax, iflat, iskip, ngeom
      real  a(maxsweep), deltaa(maxsweep),     a6(maxsweep), 
     &     al(maxsweep),   flat(maxsweep), scrch(maxsweep),
     &     da(maxsweep),     ar(maxsweep),  diffa(maxsweep),
     &     almon(maxsweep),  armon(maxsweep)
      real para( 9, maxsweep )
C
C----------------------------------------------------------------------
C
C Skip to monotonicity when interpolating total Energy for remap
C since deltaa is already set up for us.
C
      if(iskip.eq.1) then
         do 50 n = nmin, nmax
            ar(n) = deltaa(n) + al(n)
 50      continue
      else
C
         do 100 n = nmin-2, nmax+1
            diffa(n) = a(n+1) - a(n)
 100     continue
C
C                                                       Equation 1.7 of C&W
C     da(j) = D1 * (a(j+1) - a(j)) + D2 * (a(j) - a(j-1))
C     zero out da(n) if a(n) is a local max/min
C
         do 200 n = nmin-1, nmax+1
            if(diffa(n-1)*diffa(n) .lt. 0.0) then
              da(n) = 0.0
            else
              da(n) = para(4,n) * diffa(n) + para(5,n) * diffa(n-1)
              da(n) = sign( min(abs(da(n)), 2.*abs(diffa(n-1)),
     &                          2.*abs(diffa(n))), da(n) )
            endif
 200     continue
C
C                                                       Equation 1.6 of C&W
C     a(j+.5) = a(j) + C1 * (a(j+1)-a(j)) + C2 * dma(j+1) + C3 * dma(j)
C MONOT: Limit ar(n) to the range defined by a(n) and a(n+1)
C      
          do 400 n = nmin-1, nmax
            ar(n) = a(n) + para(1,n)*diffa(n) + para(2,n)*da(n+1) +
     &           para(3,n)*da(n)
            al(n+1) = ar(n)
 400     continue
C
C eqn. 4.1 - flaten interpolation in zones with a shock ( flat(n)->1. )
C
         if( iflat .eq. 1 ) then
            do 450 n = nmin, nmax
               onemfl= 1.0 - flat(n)
               ar(n) = flat(n) * a(n) + onemfl * ar(n)
               al(n) = flat(n) * a(n) + onemfl * al(n)
 450        continue
         endif
      endif
C
C MONOTONICITY constraints:
C
C compute delta_a, a_6
C MONOT: if a is a local max/min, flaten zone structure ar,al -> a.
C MONOT: compute monotonized values using eq. 1.10 of C&W
C        if parabola exceeds al/ar, reset ar/al so that slope -> 0.
C Recalculate delta_a and a_6
C
      do 510 n = nmin, nmax
         deltaa(n)= ar(n) - al(n)
         a6(n)    = 6. * (a(n) - .5 * (al(n) + ar(n)))
         scrch(n) = (ar(n) - a(n)) * (a(n)-al(n)) 
         if(scrch(n) .lt. 0.0) then
            ar(n) = a(n)
            al(n) = a(n)
         endif
         almon(n) = 3. * a(n) - 2. * ar(n)
         armon(n) = 3. * a(n) - 2. * al(n)
         scrch(n) = deltaa(n) * (deltaa(n) - a6(n))
         if(scrch(n) .lt. 0.0) al(n) = almon(n)        
         scrch(n) = deltaa(n) * (a6(n) + deltaa(n))
         if(scrch(n) .lt. 0.0) ar(n) = armon(n)
         deltaa(n)= ar(n) - al(n)
         a6(n)    = 6. * (a(n) - .5 * (al(n) + ar(n)))
510   continue
C
      return
      end