c----------------------------------------------------------------------
c     Main routine for computing integral using LSODA ODE integrator.
c----------------------------------------------------------------------
      program        integral

      implicit       none

      character*2    itoc
      character*64   catsqz
      integer        iargc,         indlnb,         vsxynt
      integer        vsrc

      integer         stderr,      stdin,     stdout
      parameter     ( stderr = 0,  stdin = 5, stdout = 6 )

      real*8          s2r8
      real*8          r8arg
      integer         s2i4
      integer         i4arg
      character*128   sarg

      real*8          r8_never
      parameter     ( r8_never = -1.0d-60 )

      integer         i4_never 
      parameter     ( i4_never = -2 000 000 000 )

      integer        neq
      parameter    ( neq = 1 )

      integer        maxout
      parameter    ( maxout = 10 000 )

      real*8         y0(neq)
      real*8         vtout(maxout), vyout(maxout,neq),
     &               work(maxout)
      integer        ntout,         itout

      integer        ieq

      logical        ltrace
      parameter    ( ltrace = .false. )

      integer        maxdump
      parameter    ( maxdump = 50 )

      logical        vstrace
      parameter    ( vstrace = .false. )

      logical        lsodatrace
      parameter    ( lsodatrace = .false. )

      logical        ivs_ok

c----------------------------------------------------------------------
c     LSODA Variables.
c----------------------------------------------------------------------
      external       fcn,        jac

      real*8         y(neq)
      real*8         tbgn,       tend
      integer        itol
      real*8         rtol,       atol
      integer        itask,      istate,     iopt
      integer        lrw

      parameter    ( lrw = 22 + neq * 16 )
      real*8         rwork(lrw)

      integer        liw
      parameter    ( liw = 20 + neq )
      integer        iwork(liw)
      integer        jt

      real*8         tol
      real*8         default_tol
      parameter    ( default_tol = 1.0d-6 )

c----------------------------------------------------------------------
c     Common communication with routine 'fcn' in 'fcn.f' ...
c----------------------------------------------------------------------

      include       'fcn.inc'

c----------------------------------------------------------------------
c     Parse command line arguments (initial values) ...
c----------------------------------------------------------------------
      if( iargc() .lt. neq ) then
            write(0,*) 'integral: Insufficient initial values '//
     &                 'specified on command line.'
            go to 900
      end if

      ivs_ok = .true.
      do ieq = 1 , neq
         y0(ieq) = r8arg(ieq,r8_never)
         if( y0(ieq) .eq. r8_never ) then
            write(0,*) 'integral: Initial value ', ieq, 
     &                 ' unrecognizable as a real number.'
            ivs_ok = .false.
         end if
      end do

      if( .not. ivs_ok ) go to 900
      tol = r8arg(neq+1,default_tol)

      if( ltrace ) then
         call dvdump(y0,neq,'integral: initial values',0)
      end if

c----------------------------------------------------------------------
c     Get output times from standard input ...
c----------------------------------------------------------------------
      call dvfrom('-',vtout,ntout,maxout)
      if( ntout .le. 0 ) then
         write(0,*) 
         write(0,*) 'integral: No output times read from standard in'
         write(0,*) 'integral: Use (e.g.) dvmesh or dvgmesh to '//
     &              'generate output times.'
         write(0,*) 
         write(0,*) 'integral: Sample usages:'
         write(0,*) 
         write(0,*) 'integral: dvmesh  0.0 1.0 101 | integral ... '
         write(0,*) 'integral: dvgmesh 0.0 1.0 101 10 | integral ...'
         write(0,*) 
         stop
      end if
      if( ltrace ) then
         if( ntout .le. maxdump ) then
            call dvdump(vtout,ntout,'integral: output times',0)
         else 
            call dvdmp1(vtout,work,
     &                  int(1.0d0 * ntout / maxdump + 0.5d0),
     &                  ntout,'integral: selected output times',0)
         end if
         write(0,*)
         write(0,*) '***** INITIAL TIME: ', vtout(1)
         write(0,*)
      end if

c----------------------------------------------------------------------
c     Set LSODA parameters ...
c----------------------------------------------------------------------
      itol = 1
      rtol = tol
      atol = tol
      itask = 1 
      iopt = 0
      jt = 2

c----------------------------------------------------------------------
c     Initialize the solution ...
c-----------------------------------------------------------------------
      call dvcopy(y0,y,neq)

c----------------------------------------------------------------------
c     Do the integration ...
c----------------------------------------------------------------------
      do itout = 1 , ntout
         do ieq = 1 , neq
            vyout(itout,ieq) = y(ieq)
         end do
         if( itout .eq. ntout ) go to 500
            
         istate = 1
c----------------------------------------------------------------------
c        Need these temporaries since lsoda overwrites tend ...
c----------------------------------------------------------------------
         tbgn = vtout(itout)
         tend = vtout(itout+1)
         call lsoda(fcn,neq,y,tbgn,tend,
     &              itol,rtol,atol,itask,    
     &              istate,iopt,rwork,lrw,iwork,liw,jac,jt)
         if( lsodatrace ) then
            write(0,1000) itout, ntout, vtout(itout), vtout(itout+1)
1000        format(' integral: Step ',i4,'  t = ',1pe10.3,' .. ',
     &             1pe10.3)
            write(0,*) 'integral: lsoda reurns ', istate
         end if

         if( istate .lt. 0 ) then
            write(0,1500) istate, itout, ntout, vtout(itout),
     &                    vtout(itout+1)
1500        format(/' integral: Error return ',i2,
     &              ' from integrator LSODA.'/
     &              ' integral: At output time ',i5,' of ',i5/
     &              ' integral: Interval ',1pe11.3,' .. ',1pe11.3/)
            go to 500
         end if
      end do

500   continue
      if( vstrace ) then
         do ieq = 1 , neq
            vsrc = vsxynt(catsqz(catsqz('y(',itoc(ieq)),')'),
     &                    0.0d0,vtout,vyout(1,ieq),itout)
         end do
      end if
      write(*,*) vtout(itout), ( vyout(itout,ieq) , ieq = 1 , neq )
      stop

900   continue
         write(0,*) 'usage: integral <iv 1> <iv 2> ... <iv neq> [tol]'
         write(0,*) 
         write(0,*) '       Sode reads requested output times'
         write(0,*) '       from standard input, one number '
         write(0,*) '       per line.'
         write(0,*)
         write(0,*) '       Current default tolerance ', default_tol
      stop
      end