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

c General information:

c     Subroutine bcy returns the boundary condition at y = 0 for the
c     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 bcy 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.  The value of a(0)
c     should be nonzero and the value of a(lmax) should be zero or small.

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

c     The present implementation of bcy requires an exact solution
c           function exa(x,y,t,ax,ay).
c     The boundary value u(j,0) is set equal to this exact solution,
c     so bcy returns the diagonal equations:
c           a(0) = 1;  a(l) = 0 for l > 0
c           b(j) = exa(x(j),y(0),t,ax,ay).

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 boundary condition:
      double precision a(0:ldim)

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

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

c     It may be noted that if you change boundary condition bcy, 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     Subroutine bcy requires a separate function
c           exa(x,y,t,ax,ay)
c     to return the exact solution at position (x,y) and time t.
      double precision exa
      external exa

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
      parameter (zero=0.d0,one=1.d0)

c Executable statements:

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

c     The right hand sides:
      do 200 j=0,jmax
         b(j)=exa(x(j),y(0),t,ax,ay)
 200  continue

c     All done:
      goto 900

c     Exit:

 900  return
      end