subroutine out(j_dim,l_dim,j_max,l_max,c_x,c_y,x,y,t,n,u)

c     Performs output.  Writes to Fortran unit 1.

c     Avoid typos:
      implicit none

c Input variables:

c     Dimensioned mesh sizes in the x- and y-directions:
      integer j_dim,l_dim

c     Actual mesh sizes in the x- and y-directions:
      integer j_max,l_max

c     Components of the convection velocity in the x and y directions:
      double precision c_x,c_y

c     Coordinates x and y
      double precision x(0:j_dim),y(0:l_dim)

c     Time t:
      double precision t
c     Time step index:
      integer n

c     The values of the dependent variable:
      double precision u(0:j_dim,0:l_dim)

c Local variables:

c     Error in the computed solution, and its maximum and RMS values:
      double precision err,errmax,errrms

c     Indices in the x- and y-directions:
      integer j,l
c     Increments in these indices:
      integer j_inc,l_inc

c     An exact solution of the equation, used to test the method:
      double precision u_exa
      external u_exa

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

c Executable statements:

c     Write a header in the output file:
      write(1,100)t,n
 100  format(/'t =',g13.5,'   n =',i5)

c     Compute the maximum and root mean square errors:
      errmax=zero
      errrms=zero
      do 300 j=1,j_max
         do 200 l=1,l_max
            err=abs(u(j,l)-u_exa(x(j),y(l),t,c_x,c_y))
            if(err.gt.errmax)errmax=err
            errrms=errrms+err**2
 200     continue
 300  continue
      errrms=sqrt(errrms/j_max/l_max)
      write(1,310)errmax,errrms
 310  format('Maximum error:',g10.2,'   RMS error:',g10.2)
      print320,t,errmax,errrms
 320  format(' t =',g11.3,
     & '   Maximum error:',g10.2,'   RMS error:',g10.2)

c     Write abbreviated mesh point values
      j_inc=max(1,(j_max+7)/8)
      l_inc=max(1,(l_max+7)/8)
      write(1,380)(j,j=0,j_max,j_inc)
 380  format('Selected mesh point u-values:'/'   l\\j:',9i8)
      do 400 l=0,l_max,l_inc
         write(1,390)l,(u(j,l),j=0,j_max,j_inc)
 390     format(i4,'   ',9f8.4)
 400  continue

      return
      end