subroutine out(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,u,io)

c General information:

c     Subroutine out performs output for solutions of the convection
c     equation on a rectangular region of size xmax by ymax.

c     Subroutine out requires an input file 'out.in' of the form:
c        <Number of times to do output imax>     <comment>
c        <Output time tout(1)>                   <comment>
c        <Output time tout(2)>                   <comment>
c        ....
c        <Output time tout(imax)>                <comment>

c     Subroutine out requires an exact solution
c           function exa(x,y,t,ax,ay).
c     The numerical results are compared to this exact solution and
c     the found maximum and root-mean-square errors are printed.

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: horizontal and vertical dimensions 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 to output
      double precision t

c     Input: the values of the dependent variable u at time t to output:
      double precision u(0:jdim,0:ldim)

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

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

c     Subroutine out 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     An integer to keep track whether out has been initialized:
      integer init

c     Dimensioned maximum number of output times:
      integer idim
      parameter (idim=10)
c     Requested number of output times:
      integer imax
      save imax
c     Index over the output times:
      integer i
c     The output times
      double precision tout(idim)
      save tout

c     Indices in the x- and y-directions:
      integer j,l

c     Mesh-point error, and maximum and RMS errors in the computed solution:
      double precision err,errmax,errrms
c     Number of errors:
      double precision count

c     Fortran I/O unit to do input on:
      integer ioin

c     Whether a file is open on an I/O unit number:
      logical open

c     The previous time at which out was called
      double precision tprev
      save tprev

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

c Data statements:

      data init/0/

c Executable statements:

c Initialize the first time:

      if(init.eq.0)then

c        Find a free I/O unit, but not the reserved units 5, 6, or 7:
         ioin=io
 100     ioin=ioin+1
         if(ioin.ge.5)ioin=max(ioin,8)
         inquire(ioin,opened=open)
         if(open)goto 100

c        Open the input file:
         open(ioin,file='out.in',status='old')

c        Read the number of output times:
         read(ioin,*)imax
         write(io,210)imax
 210     format('Number of output times:',i4)
         if(imax.lt.0)then
            print*,'*** out: '//
     &     'The number of output times cannot be negative but is',imax
            stop
         endif
         if(imax.gt.idim)then
            print*,'*** out: '//
     &      'The requested number of output times,',imax
            print*,'            '//
     &      'exceeds the maximum dimensioned limit',idim
            stop
         endif

c        Read the output times:
         do 300 i=1,imax
            read(ioin,*)tout(i)
 300     continue
         write(io,310)(tout(i),i=1,imax)
 310     format(5g15.7)

c        Initialize the previous time to an arbitrary number less than t:
         tprev=t-one-abs(t)

c        All done initializing:
         init=1
      endif

c See whether we want to do output:

      do 1000 i=1,imax
         if(tout(i).gt.tprev .and. tout(i).le.t)goto 1010
 1000 continue
      goto 9000
 1010 continue

c Perform output:

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

c     Compute the errors:
      errmax=zero
      errrms=zero
      count=zero
      do 1300 j=0,jmax
         if(x(j).lt.zero .or. x(j).gt.xmax)goto 1300
         do 1290 l=0,lmax
            if(y(l).lt.zero .or. y(l).gt.ymax)goto 1290
            err=abs(u(j,l)-exa(x(j),y(l),t,ax,ay))
            if(err.gt.errmax)errmax=err
            errrms=errrms+err**2
            count=count+one
 1290    continue
 1300 continue
      errrms=sqrt(errrms/count)
      write(1,1310)errmax,errrms
 1310 format('Maximum error:',g10.2,'   RMS error:',g10.2)
      print 1320,t,errmax,errrms
 1320 format(' t =',g11.3,
     &'   Maximum error:',g10.2,'   RMS error:',g10.2)

c     All done:
      goto 9000

c     Exit:

 9000 tprev=t
      return
      end