subroutine bcya(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,io,
     & a,b)

c General information:

c     Subroutine bcya returns an artificial boundary condition for
c     the top boundary y = ymax for the convection equation
c           u_t + ax u_x + ay u_y =0
c     on a rectangular region of size xmax by ymax.

c     The boundary condition is returned in the form of the coefficients
c     of an equation.  In particular, if u(j,l) is the mesh-value of the
c     unknown u at an arbitrary mesh point (j,l), a boundary condition
c     that bcya can handle is of the general form
c        sum (l = 0 to lmax) a(l) u(j,l) = b(j)  for j = 0, 1, ..., jmax
c     where vector a is independent of x and time.

c     Subroutine bcya returns the vector a and the right hand sides
c     b(j) of the equations for all j values.

c     The present implementation of bcya finds the boundary point by
c     linear extrapolation from the interior:
c        u(j,lmax) = 2 u(j,lmax-1) - u(j,lmax-2)

c     Copyright 1997 Leon van Dommelen
c     Version 1.0 Leon van Dommelen 1/10/97

c Arguments:

c     Avoid typos:
      implicit none

c     Input: components of the convection velocity in the x and y directions:
      double precision ax,ay

c     Input: x- and y-sizes of the rectangular domain:
      double precision xmax,ymax

c     Input: declared maximum indices in the x- and y-directions:
      integer jdim,ldim

c     Input: actual maximum indices in the x- and y-directions:
      integer jmax,lmax

c     Input: the x- and y-values of the mesh points:
      double precision x(0:jdim),y(0:ldim)

c     Input: the time at which to return the boundary condition:
      double precision t

c     Input: Fortran I/O unit to do output on:
      integer io

c     Output: the coefficients of the artificial boundary condition:
      double precision a(0:ldim)

c     Output: the right hand sides of the artificial boundary condition:
      double precision b(0:jdim)

c External variables and info for compiling or changing subroutine bcya:

c     It may be noted that if you change boundary condition bcya, for
c     example to use mirror points or staggering, you may also have to
c     change subroutine ini which generates the mesh-point y-values.

c Local variables:

c     Index in the x-direction:
      integer j

c     Index in the y-direction:
      integer l

c     Constants, defined to simplify changing precision:
      double precision zero,one,two
      parameter (zero=0.d0,one=1.d0,two=2.d0)

c Executable statements:

c     Check the mesh size:
      if(lmax.lt.2)then
         print*,'*** bcya: the number of mesh intervals '//
     &   'in y-direction must be at least 2,'
         print*,'          but is:',lmax
         stop
      endif

c     Set the coefficients a(j) in the equation:
      do 100 l=0,lmax-3
         a(l)=zero
 100  continue
      a(lmax-2)=one
      a(lmax-1)=-two
      a(lmax)=one

c     The right hand sides:
      do 200 j=0,jmax
         b(j)=zero
 200  continue

c     All done:
      goto 900

c     Exit:

 900  return
      end