PROGRAM INTERPOLATE_2
c
c   This program creates continuous, parabolic interpolations
c   from zone averaged data, without assuming a uniform grid.
c
      PARAMETER ( NBDY = 2 )

      ALLOCATABLE :: XL(:), A_AVG(:), AL(:), DX(:), DALFAC(:)
      ALLOCATABLE :: DARFAC(:), FACDA(:), FCDAZL(:), FACDAL(:)
 
      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   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+1+NBDY))
          ALLOCATE    (DALFAC(1-NBDY:N+1+NBDY))
          ALLOCATE    (DARFAC(1-NBDY:N+1+NBDY))
          ALLOCATE     (FACDA(1-NBDY:N+1+NBDY))
          ALLOCATE    (FCDAZL(1-NBDY:N+1+NBDY))
          ALLOCATE    (FACDAL(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
          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_INTERFACE_FACTORS ( DX(1-NBDY), DALFAC(1-NBDY), 
     ^            DARFAC(1-NBDY), FACDA(1-NBDY), FCDAZL(1-NBDY),
     ^           FACDAL(1-NBDY), N, NBDY)

c
c   Interpolate zone interface values
c
      CALL PPM98_INTERP_INTERFACE ( DALFAC(1-NBDY), DARFAC(1-NBDY), 
     &               FACDA(1-NBDY), FCDAZL(1-NBDY), FACDAL(1-NBDY),
     &               A_AVG(1-NBDY),     AL(1-NBDY), N, NBDY )

c   AL(I) are defined for 1<= I <= N+1

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

      DO I=1,N
         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(I), A_INTERP
100         FORMAT( I4,1P, 2(1X,e12.5))
         ENDDO
      ENDDO

      DEALLOCATE ( XL, A_AVG, AL, DX, DALFAC )
      DEALLOCATE ( DARFAC, FACDA, FCDAZL, FACDAL )

      END ! PROGRAM INTERPOLATE_2