SUBROUTINE BOUNDARY ( XL, A, IBDY_TYPE, NZONES, NBDY)
c
c
C  WHAT              This routine implements one of two type of common
c                    boundary conditions: reflecting or periodic.  The
c                    boundary conditions are implemented by filling fake
c                    or boundary zones appended to each end of the arrays
c                    XL and A.  It is assumed that there are NZONES real
c                    (non-boundary) zones in the problem, NBDY fake zones
c                    to set at both the left and right ends, and that 
c                    NBDY<=NZONES.
c
c                    The grid is defined by the locations of the lefthand
c                    zone interfaces XL.  Thus zone I would have its left
c                    edge at XL(I) and its right edge at XL(I+1).  Therefore
c                    a grid of N zones requires N+1 XL values.
c
C  DATA
C                  Input arrays are:
c                      XL:  the locations of the lefthand zone interfaces
c                       A:  an array of zone averages
c
C                  Input scalers are:
c                      IBDY_TYPE:  a flag indicating the type of boundary
c                                  conditions to impose:
c                             IBDY_TYPE = 0 -> reflecting
c                             IBDY_TYPE = 1 -> periodic
c                      NZONES:  the number of real (non-boundary) zones
c                                 in the problem.
c                      NBDY:    the number of non-boundary zones to append
c                                 to each end of the arrays.
c                      N:        the number of zone in the grid.
c
C                  Input arrays are:
c                      XL:  the locations of the lefthand zone interfaces, the
c                            elements (1-NBDY:0) and (NZONES+2:NZONES+NBDY+1)
c                            are filled with appropriate values for the 
c                            selected boundary conditions.
c                      AL:  the array of zone averages, the elements (1-NBDY:0) 
c                           and (NZONES+1:NZONES+NBDY)
c                            are filled with appropriate values for the 
c                            selected boundary conditions.
c  DIMENSIONS
c                   XL [1-NBDY:N+NBDY+1]
c                    A [1-NBDY:N+NBDY]
c
c
c
c
c
c
      REAL XL(1-NBDY:NZONES+NBDY+1)
      REAL  A(1-NBDY:NZONES+NBDY)

      IF ( IBDY_TYPE .NE. 0 ) GOTO 500
c
c   Reflecting boundary conditions
c
c   Left:
        DO I=1-NBDY, 0
           XL(I) = 2.0 * XL(1)  - XL(2 - I)
           A(I) = A(1-I)
        ENDDO
c
c
c   Right:
        DO I = NZONES+1, NZONES+NBDY
           XL(I+1) = 2.0 * XL(NZONES+1)  
     ^                   - XL(2*NZONES+1 - I)
           A(I) = A(2*NZONES+1-I)
        ENDDO

      RETURN
500   CONTINUE
c
      IF ( IBDY_TYPE .NE. 1 ) GOTO 900
c
c   Periodic boundary conditions
c
c   Left:
      DXGRID = XL(NZONES+1) - XL(1)
      DO I=1-NBDY,0
         XL(i) = XL(I+NZONES) - DXGRID
         A(i) = A(I+NZONES) 
      ENDDO
c
c   Right:
      DO   I = NZONES+1,NZONES+NBDY
         XL(I+1) = XL(I+1-NZONES) + DXGRID
         A(I) = A(I-NZONES) 
      ENDDO
c
      RETURN
900   CONTINUE
c
c   Error, print advisory and terminate program.
       PRINT *,' UNKNOWN BOUNDARY TYPE:',IBDY_TYPE
       CALL ABORT()
c
      END  ! SUBROUTINE BOUNDARY