PROGRAM ADVECT_2
     
c
c   This program advects a given distribution of zone averages through
c   a numerical grid at a fixed Courant number.  Periodic boundary
c   conditions are assumed, as is a uniform grid.
c
c
      PARAMETER ( NBDY = 4 )

      ALLOCATABLE :: XL(:),    A_AVG(:),     AL(:), AR(:), DA(:), A6(:),
     ^          SMALLDA(:), UNSMOOTH(:), DISCONT(:), DX(:),      DAL(:),
     ^          DVOLL(:),     DVOLFL(:),    XLNU(:),       A_AVGNU(:), 
     ^          DXINV(:)
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 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       (DAL(1-NBDY:N+NBDY))
          ALLOCATE      (XLNU(1-NBDY:N+1+NBDY))
          ALLOCATE   (A_AVGNU(1-NBDY:N+NBDY))
          ALLOCATE   (SMALLDA(1-NBDY:N+NBDY))
          ALLOCATE  (UNSMOOTH(1-NBDY:N+NBDY))
          ALLOCATE   (DISCONT(1-NBDY:N+NBDY))
          ALLOCATE   (DVOLL(1-NBDY:N+NBDY))
          ALLOCATE     (DVOLFL(1-NBDY:N+NBDY))

c   The following arrays will hold  various geometrical factors 
          ALLOCATE        (DX(1-NBDY:N+NBDY))
          ALLOCATE     (DXINV(1-NBDY:N+NBDY))


      DO I=1,N
         READ (11,*) XL(I),A_AVG(I)
      ENDDO
      READ (11,*) XL(N+1)
c
c
1     PRINT *,' Carefully enter the fractional offset (between -1.0 and 1.0)'
         PRINT *,' between the two grids: DXF.  per iteration '
         READ *,DXF
         IF( DXF .LT. -1.0  .OR. DXF .GT. 1.0 ) THEN
            PRINT *,' You entered',DXF, ' Please try again.'
	 GOTO 1
         ENDIF
            PRINT *,' The grids will be separated by ',DXF,' zone widths '
            PRINT *,' Thank you.'

c
c
c   Query the user for desired number of revolutions'
c   through the grid'
c
c
2     PRINT *,' Please enter the number of revolutions'
      PRINT *,'   REVS > 0.0 (fractions allowed)'
      READ *, REVS
      IF( REVS .LT. 0.0 ) GOTO 2
c
c   at a Courant number of COURNO, a grid of N zones reauires N/|COURNO|
c   iterations for one revolution thus:
c
      NITER = INT ( FLOAT(N) * REVS / ABS(DXF) )
c
c   we will always perform at least 1 iteration
c
      NITER = MAX (1,NITER)

      print *,' NITER = ',niter
c
c   Now perform the required number of iterations.
c
      DO ITER = 1,NITER
c
c   Initialize the fake zones assuming periodic boundary conditions:

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

c
c   While a uniform grid is assumed, it is convient to calculate
c   the individual zone widths.

         DO I=1-NBDY,N+NBDY
            DX(I) = XL(I+1) - XL(I)
            DXINV(I) = 1.0 / DX(I)
         ENDDO
c
c
c    Determine the locations of the lefthand zone interface of the new grid
c    (The first and last interface locations will not be needed)
c
         DO I=2-NBDY,N+NBDY
          IF (DXF .LE. 0.0)   THEN
	   XLNU(I) = XL(I) + DXF * DX(I-1)
	 ELSE
	   XLNU(I) = XL(I) + DXF * DX(I)
	 ENDIF
         ENDDO


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

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


         DO I=2-NBDY,N+NBDY

c     Advected volume:
          DVOLL(i) = XL(I) - XLNU(I)
c     Advected volume fraction:
          IF (DVOLL(I) .GE. 0.0)   THEN
               DVOLFL(I) = DVOLL(I)*DXINV(i-1)
         ELSE
               DVOLFL(I) = -DVOLL(I)*DXINV(I)
c
         ENDIF
         ENDDO


c
c   Determinve new zone-averages
c
c   PPM98_ADVECT_PASSIVE0 interplates A_AVG assuming a uniform grid,
c        with discontinuity detection and monotonicity constraints.
c        It returns the new zone averages A_AVGNU, as well as the 
c        interface fluxes DAL and the results of the unsmoothness 
c        and discontinuity detection calculations.

         CALL PPM98_ADVECT_PASSIVE0    
     ^             ( DVOLFL(1-NBDY),   DVOLL(1-NBDY),       DX(1-NBDY),
     ^                DXINV(1-NBDY),   A_AVG(1-NBDY),  SMALLDA(1-NBDY), 
     ^             UNSMOOTH(1-NBDY), DISCONT(1-NBDY),      DAL(1-NBDY), 
     ^              A_AVGNU(1-NBDY), N, NBDY )
c
c   In Because PPM98_ADVECT_PASSIVE0 does not overwrite the input
c   data, it is necessary to replace the values in the arrays for
c   the old data with the new values.  Ordinarily the grid would be
c   replaced as well, but that is not done here to ease comparisons.
c
         DO I=1,N
	 A_AVG(I) = A_AVGNU(I)
         ENDDO
      ENDDO
c
c   Display the zone averages
      PRINT *,'  I       XMID        A_AVG'
      DO I=5-NBDY,N+NBDY-4
         XMID = 0.5*(XL(I+1)+XL(I))
         WRITE(6,100) I, XMID, A_AVG(I)
100   FORMAT(3x,i4,1p,4(1x,e12.5))
      ENDDO

      DEALLOCATE (XL, A_AVG, AL, AR, DA, A6, SMALLDA, UNSMOOTH, DISCONT)
      DEALLOCATE (DX, A_AVGNU, XLNU, DAL, DVOLL, DVOLFL, DXINV)


      END ! PROGRAM ADVECT_2