PROGRAM SECOND_DER
c
c   This program is meant as an example of one way of treating
c   boundary conditions.  Given an array of zone averages A and
c   an array of lefthand zone interface locations XL which
c   define the grid, approximate the average of the 2dn derivative
c   in each zone.  Assume reflecting boundary conditions.
c
c   Zone averages and grid interface locations are read from a file
c   "INPUT.AVG" which was created earlier by either of 
c    CREATE_ZONE_AVERAGES_SINE or CREATE_ZONE_AVERAGES_TANH.
c
c   The 2nd derivative is then approximated from
c
c       (A_AVG(I-1) - 2.0*A_AVG(I) + A_AVG(I+1))/DX**2
c
c   This requires NBDY = 1, that is one fake zone on each end of the arrays.
c
c
      ALLOCATABLE :: XL(:), A_AVG(:)
          PARAMETER ( NBDY = 1 )
c
      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 programs CREATE_ZONE_AVERAGES*.  



c   Read in the number of zones in the file
            READ (11,*)  N

c   Next, allocate the appropriate space
          ALLOCATE    (XL(1-NBDY:N+1+NBDY))
          ALLOCATE (A_AVG(1-NBDY:N+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   Assume reflecting boundaries:
c
          IBDY_TYPE = 0


          CALL BOUNDARY( XL(1-NBDY), A_AVG(1-NBDY), IBDY_TYPE, N,NBDY )
c
c   We are assuming a uniform grid: DX(I) = constant
c
	  DX = XL(2)-XL(1)
	  DXSQINV = 1.0 / (DX*DX)
          PRINT *,'      I      XMID          A_AVG(I)        D2AD2X'
          DO I=1,N
            D2AD2X = (A_AVG(I-1) - 2.0*A_AVG(I) + A_AVG(I+1))* DXSQINV
            XMID = 0.5 * (XL(I+1)+XL(I))
            WRITE (6,200) I, XMID, A_AVG(I), D2AD2X
200     FORMAT (4X,I4,t10,1P,3(E12.5,3x))
          ENDDO

      DEALLOCATE (XL, A_AVG)

          END ! PROGRAM SECOND_DER