subroutine ini(jdim,ldim,io,
     & ax,ay,xmax,ymax,tmax,jmax,lmax,nmax,x,y)

c General information:

c     Subroutine ini initializes the computation of the convection
c     equation
c        u_t + ax u_x + ay u_y = 0
c     on a rectangular region of size xmax by ymax.  See the description
c     of the arguments below for information about what variables
c     are initialized.

c     Subroutine ini needs an input file 'ini.in' of the general form:
c
c        <Value of convection velocity component ax>               <comment>
c        <Value of convection velocity component ay>               <comment>
c        <Value of the horizontal domain size xmax>                <comment>
c        <Value of the vertical domain size ymax>                  <comment>
c        <Value of the last time to compute tmax>                  <comment>
c        <Value of the number of horizontal mesh-intervals jmax>   <comment>
c        <Value of the number of vertical mesh-intervals lmax>     <comment>
c        <Value of the number of time-steps to take nmax>          <comment>
c
c     Blank lines may be used in file ini.in to increase readability.

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: declared maximum mesh-point indices in the x- and y-directions:
      integer jdim,ldim

c     Input: Fortran I/O unit number of the output file:
      integer io

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

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

c     Output: the final time to compute:
      double precision tmax

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

c     Output: the number of time-steps to take:
      integer nmax

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

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

c     It may be noted that if you change the boundary conditions. for
c     example to use mirror points or staggering, you may also have to
c     change subroutine ini to change the mesh point positions.

c Local variables:

c     The Fortran I/O unit of the input file:
      integer ioin

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

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

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

c Executable statements:


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='ini.in',status='old')

c     Read the convection velocity in the x-direction:
      read(ioin,*)ax
      write(io,110)ax
 110  format('Convection velocity in the x-direction:',g15.7)
      if(ax.le.zero)then
         print*,'*** ini: '//
     &   'The x-convection velocity must be positive, not',ax
         print*,'         The problem is improperly posed! Continuing..'
      endif

c     Read the convection velocity in the y-direction:
      read(ioin,*)ay
      write(io,120)ay
 120  format('Convection velocity in the y-direction:',g15.7)
      if(ay.le.zero)then
         print*,'*** ini: '//
     &   'The y-convection velocity must be positive, not',ay
         print*,'         The problem is improperly posed! Continuing..'
      endif

c     Read the size of the domain in the x-direction:
      read(ioin,*)xmax
      write(io,130)xmax
 130  format('Size of the domain in the x-direction:',g15.7)
      if(xmax.le.zero)then
         print*,'*** ini: '//
     &   'The x-size must be positive but is',xmax
         stop
      endif

c     Read the size of the domain in the y-direction:
      read(ioin,*)ymax
      write(io,140)ymax
 140  format('Size of the domain in the y-direction:',g15.7)
      if(ymax.le.zero)then
         print*,'*** ini: '//
     &   'The y-size must be positive but is',ymax
         stop
      endif

c     Read the maximum time:
      read(ioin,*)tmax
      write(io,150)tmax
 150  format('The maximum time to compute:',g15.7)
      if(tmax.lt.zero)then
         print*,'*** ini: '//
     &   'The maximum time must be positive or zero but is',tmax
         stop
      endif

c     Read the number of mesh segments to use in the x-direction:
      read(ioin,*)jmax
      write(io,160)jmax
 160  format('Number of mesh segments in the x-direction:',i4)
      if(jmax.lt.2)then
         print*,'*** ini: '//
     &   'The number of x-segments must be at least two but is',jmax
         stop
      endif
      if(jmax.gt.jdim)then
         print*,'*** ini: '//
     &   'The requested number of x-segments,',jmax
         print*,'            '//
     &   'exceeds the maximum dimensioned limit',jdim
         stop
      endif

c     Read the number of mesh points to use in the y-direction:
      read(ioin,*)lmax
      write(io,170)lmax
 170  format('Number of mesh segments in the y-direction:',i4)
      if(lmax.lt.2)then
         print*,'*** ini: '//
     &   'The number of y-segments must be at least two but is',lmax
         stop
      endif
      if(lmax.gt.ldim)then
         print*,'*** ini: '//
     &   'The requested number of y-segments,',lmax
         print*,'            '//
     &   'exceeds the maximum dimensioned limit',ldim
         stop
      endif

c     Read the number of time levels to use:
      read(ioin,*)nmax
      write(io,180)nmax
 180  format('Number of time levels:',i5)
      if(nmax.lt.0)then
         print*,'*** ini: '//
     &   'The number of time levels cannot be negative.'
         stop
      endif

c     Close the input file:
      close(ioin)

c     Initialize the mesh points.  We will extend the mesh one mesh cell
c     outside the domain to reduce the effect of using artificial
c     boundary conditions on the computed solution.

      do 200 j=0,jmax
         x(j)=(j+zero)/(jmax-1)*xmax
 200  continue

      do 300 l=0,lmax
         y(l)=(l+zero)/(lmax-1)*ymax
 300  continue

c     Exit:

      return
      end