PROGRAM INTERPOLATE_4
c
c   This program creates piecewise parabolic interpolations
c   from zone averaged data, without assuming a uniform grid.
c   Monotonicity constraints are applied. For those zones
c   considered to be in a discontinuity, an interpolation more
c   appropriate for discontinuous functions is used, producing
c   a better representation of unresolved shocks.
c
      PARAMETER ( NBDY = 3 )

      ALLOCATABLE :: XL(:),    A_AVG(:),     AL(:), AR(:), DA(:), A6(:),
     ^          SMALLDA(:), UNSMOOTH(:), DISCONT(:)
c
c   The following arrays will hold  various geometrical factors which can be
c   precomputed for efficiency when several quantities are interpolated on the
c   same grid:

      ALLOCATABLE :: DX(:), DALFAC(:), DARFAC(:),  FACDA(:), FCDAZL(:),
     ^           FACDAL(:), FA6DAL(:), FA6DAR(:), FRRDAL(:), FRRDAR(:),
     ^           FLLDAL(:), FLLDAR(:)


      OPEN (unit = 11, file ='INPUT.AVG')
c  The data consists of N+2 lines.  The first line
c  contains the integer N, the number to zone averages in the data set.
c  This is followed by N lines, each line having two numbers, a value for
c  the location of the left edge of the zone XL, and the zone average for
c  that zone A_AVG.  The lines are assumed to be in sequential order according
c  to XL.  An additional line has only one value, the final left zone edge.  A
c  grid of N zones requires N+1 left zone edge locations. An appropriate data
c  set can be created from the program CREATE_ZONE_AVERAGES.  The XL values
c   are used for plotting, not in the interpolation process.


c   Read in the number of zones in the file
            READ (11,*)  N
c
c   Next, allocate the appropriate space
          ALLOCATE        (XL(1-NBDY:N+1+NBDY))
          ALLOCATE        (AL(1-NBDY:N+1+NBDY))
          ALLOCATE        (AR(1-NBDY:N+1+NBDY))
          ALLOCATE        (DA(1-NBDY:N+NBDY))
          ALLOCATE        (A6(1-NBDY:N+NBDY))
          ALLOCATE     (A_AVG(1-NBDY:N+NBDY))
          ALLOCATE   (SMALLDA(1-NBDY:N+NBDY))
          ALLOCATE  (UNSMOOTH(1-NBDY:N+NBDY))
          ALLOCATE   (DISCONT(1-NBDY:N+NBDY))

c   The following arrays will hold  various geometrical factors which can be
c   precomputed for efficiency when several quantities are interpolated on the
c   same grid:
          ALLOCATE        (DX(1-NBDY:N+NBDY))
          ALLOCATE    (DALFAC(1-NBDY:N+NBDY))
          ALLOCATE    (DARFAC(1-NBDY:N+NBDY))
          ALLOCATE     (FACDA(1-NBDY:N+NBDY))
          ALLOCATE    (FCDAZL(1-NBDY:N+NBDY))
          ALLOCATE    (FACDAL(1-NBDY:N+NBDY))
          ALLOCATE    (FA6DAL(1-NBDY:N+NBDY))
          ALLOCATE    (FA6DAR(1-NBDY:N+NBDY))
          ALLOCATE    (FRRDAL(1-NBDY:N+NBDY))
          ALLOCATE    (FRRDAR(1-NBDY:N+NBDY))
          ALLOCATE    (FLLDAL(1-NBDY:N+NBDY))
          ALLOCATE    (FLLDAR(1-NBDY:N+NBDY))


      DO I=1,N
         READ (11,*) XL(I),A_AVG(I)
      ENDDO
      READ (11,*) XL(N+1)
c
c
c   Query the user for desired boundary conditions
c
      PRINT *,' Please enter 0 for reflecting boundary conditions'
      PRINT *,'        or    1 for periodic boundary conditions'
      READ *, IBDY_TYPE
      IF ( IBDY_TYPE .EQ. 0 ) THEN
         PRINT *,' Selecting reflecting boundary conditions'
      ELSE
         IF ( IBDY_TYPE .EQ. 1 ) THEN
            PRINT *,' Selecting periodic boundary conditions'
         ELSE
            PRINT *,' Unknown boundary type:',IBDY_TYPE
            call abort()
         ENDIF

      ENDIF
c
c   Initialize the fake zones:


          CALL BOUNDARY( XL(1-NBDY), A_AVG(1-NBDY), IBDY_TYPE, N,NBDY )
c

c
c   A uniform grid is NOT assumed.  Therefore, various
c   geometrical factors must be precomputed by a call to
c   _INTERFACE_FACTORS, input to  which consists of DX(I),
c   the individual zone widths.

      DO I=1-NBDY,N+NBDY
         DX(I) = XL(I+1) - XL(I)
      ENDDO


      CALL PPM98_INTERPOLATION_FACTORS (  DX(1-NBDY), DALFAC(1-NBDY),
     ^                DARFAC(1-NBDY),  FACDA(1-NBDY), FCDAZL(1-NBDY),
     ^                FACDAL(1-NBDY), FA6DAL(1-NBDY), FA6DAR(1-NBDY),
     ^           FRRDAL(1-NBDY), FRRDAR(1-NBDY), FLLDAL(1-NBDY),
     ^           FLLDAR(1-NBDY), N, NBDY)



c   Fill an array, SMALLDA, which contains what is to be considered a trivial difference of A_AVG.
c   In this Program, A_AVG is not positive definite, so a small value is appropriate.

      DO I=1-NBDY,N+NBDY
         SMALLDA(I) = 0.0005
      ENDDO

c   Interpolate zone interface values and determine ppm parabola coefficients
c
      CALL PPM98_INTERP_DISCONTINUITY (  DALFAC(1-NBDY), DARFAC(1-NBDY),  
     ^                 FACDA(1-NBDY),    FCDAZL(1-NBDY), FACDAL(1-NBDY),   
     ^                FA6DAL(1-NBDY),    FA6DAR(1-NBDY), FLLDAL(1-NBDY), 
     ^                FLLDAR(1-NBDY),    FRRDAL(1-NBDY), FRRDAR(1-NBDY),  
     ^               SMALLDA(1-NBDY),     A_AVG(1-NBDY),     AL(1-NBDY),     
     ^                    AR(1-NBDY),        DA(1-NBDY),     A6(1-NBDY), 
     ^              UNSMOOTH(1-NBDY),   DISCONT(1-NBDY),      N, NBDY )

c   AL(I), AR(I), DA(I), A6(I), are defined for 1<= I <= N

c   Write out the interpolations for plotting [11 samples per zone]


      DO I=1,N
         DO J=0,10
            XX = 0.1*FLOAT(J)
            A_INTERP = AL(I) + XX * (DA(I) + A6(I) * (1.0 - XX))
            WRITE(6,100) I, XL(I)+XX*DX(I), A_INTERP, UNSMOOTH(I), DISCONT(I)
100         FORMAT( I4,1P, 4(1X,e12.5))
         ENDDO
      ENDDO


      DEALLOCATE (XL, A_AVG, AL, AR, DA, A6, SMALLDA, UNSMOOTH, DISCONT)
      DEALLOCATE (DX, DALFAC, DARFAC,  FACDA, FCDAZL,
     ^        FACDAL, FA6DAL, FA6DAR, FRRDAL, FRRDAR,
     ^        FLLDAL, FLLDAR)


      END ! PROGRAM INTERPOLATE_4