PROGRAM INTERPOLATE_1
c
c   This program creates continuous, parabolic interpolations
c   from zone averaged data, assuming a uniform grid.


c
      PARAMETER ( NBDY = 2 )

      ALLOCATABLE :: XL(:), A_AVG(:), AL(:)

      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 (A_AVG(1-NBDY:N+NBDY))
          ALLOCATE    (AL(1-NBDY:N+1+NBDY))
c
c   Read in the rest of the data
c
      DO I=1,N
         READ (11,*) XL(I),A_AVG(I)
      ENDDO
      READ (11,*) XL(N+1)
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:
c
c
          CALL BOUNDARY( XL(1-NBDY), A_AVG(1-NBDY), IBDY_TYPE, N,NBDY )
c
c
c   Interpolate zone interface values
c
      CALL PPM98_INTERP0_INTERFACE (A_AVG(1-NBDY), AL(1-NBDY), N, NBDY)

c   AL(I) are computed on range 1<= I <= N+1

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

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

      DEALLOCATE (XL, A_AVG, AL)

      END ! PROGRAM INTERPOLATE_1