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

c General information:

c     Subroutine bcx returns the boundary condition at x = 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 bcx can handle is of the general form
c        sum (j = 0 to jmax) a(j) u(j,l) = b(l)  for l = 0, 1, ..., lmax
c     where vector a is independent of y and time.  The value of a(0)
c     should be nonzero and the value of a(jmax) should be zero or small.

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

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

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

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

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

c     Subroutine bcx 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(j) in the equation:
      a(0)=one
      do 100 j=1,jmax
         a(j)=zero
 100  continue

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

c     All done:
      goto 900

c     Exit:

 900  return
      end