subroutine parabola( nmin, nmax, nz, para, a, deltaa, a6, al, 
     &                     flat, iflat, iskip )
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
C nmin/nmax are indicies over which the parabolae are calculated

      include 'sweepsize.h'

      integer n, nmin, nmax, iflat, iskip, nz
      real  a(maxsweep), deltaa(maxsweep),     a6(maxsweep), 
     &     al(maxsweep),   flat(maxsweep), 
     &     da(maxsweep),     ar(maxsweep),  diffa(maxsweep),
     &     scrch1(maxsweep), scrch2(maxsweep)
      real para( 5, nz )

C----------------------------------------------------------------------
C Skip to monotonicity when interpolating total Energy for remap
C since deltaa is already set up for us.

      if(iskip.eq.1) then
         do n = nmin, nmax
            ar(n) = deltaa(n) + al(n)
         enddo
      else
C
         do n = nmin-2, nmax+1
            diffa(n) = a(n+1) - a(n)
         enddo

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

         do 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
         enddo

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)

          do 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)
          enddo

C eqn. 4.1 - flaten interpolation in zones with a shock ( flat(n)->1. )

         if( iflat .eq. 1 ) then
            do 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)
            enddo
         endif
      endif

C MONOTONICITY constraints:

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

      do n = nmin, nmax

         deltaa(n)= ar(n) - al(n)
         a6(n)    = 6. * (a(n) - .5 * (al(n) + ar(n)))

         scrch1(n) = (ar(n) - a(n)) * (a(n)-al(n)) 
         if(scrch1(n) .lt. 0.0) then
            ar(n) = a(n)
            al(n) = a(n)
         endif

         scrch1(n) = deltaa(n) * deltaa(n)
         scrch2(n) = deltaa(n) * a6(n)
         if(scrch1(n) .lt. +scrch2(n)) al(n) = 3. * a(n) - 2. * ar(n)       
         if(scrch1(n) .lt. -scrch2(n)) ar(n) = 3. * a(n) - 2. * al(n)

         deltaa(n)= ar(n) - al(n)
         a6(n)    = 6. * (a(n) - .5 * (al(n) + ar(n)))

      enddo

      return
      end