PROGRAM INTERPOLATE_3
c
c   This program creates piecewise parabolic interpolations
c   from zone averaged data, without assuming a uniform grid.
c   Monotonicity constraints are applied.
c
      PARAMETER ( NBDY = 3 )

      ALLOCATABLE :: XL(:), A_AVG(:), AL(:), AR(:), SMALLDA(:), 
     ^         UNSMOOTH(:)
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   (A_AVG(1-NBDY:N+NBDY))
          ALLOCATE (SMALLDA(1-NBDY:N+NBDY))
          ALLOCATE      (AL(1-NBDY:N+1+NBDY))
          ALLOCATE      (AR(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+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))
          ALLOCATE  (UNSMOOTH(1-NBDY:N+NBDY))




c   First, read in the data, which consists of N+1 lines.  Each
c   line has two numbers, a value for the location of the left edge
c   of the zone ZL, and the zone average for that zone A_AVG.  The
c   lines are assumed to be in sequential order according to XL.  The final
c   line has only one value, the final left zone edge.  A grid of N zones
c   requires N+1 left zone edge locations. An appropriate data set
c   can be created from the program CREATE_ZONE_AVERAGES.

      OPEN (unit = 11, file ='INPUT.AVG')

      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:

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_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
c   Fill an array, SMALLDA, which contains what is to be considered a trivial 
c   difference of A_AVG.  In this Program, A_AVG is not positive definite, so 
c  a small value is appropriate.

      DO I=1-NBDY,N+NBDY

         SMALLDA(I) = 0.0005
      ENDDO

c   Interpolate zone interface values
c
      CALL PPM98_INTERP_MONO ( 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), UNSMOOTH(1-NBDY),   N, NBDY)

c   AL(I) is defined for 1<= I <= N
c   AL(I), and AR(I) are defined for 1<= I <= N

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


      DO I=1,N
         DA = AR(I) - AL(I)
         A6 = 6.0 * (A_AVG(I) - 0.5*(AL(I) + AR(I)))
         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, UNSMOOTH(I)
100         FORMAT( I4,1P, 3(1X,e12.5))
         ENDDO
      ENDDO

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


      END ! PROGRAM INTERPOLATE_3