c===========================================================
c     History: sode.f
c
c     Driver routine which integrates ODEs defining 
c     model for deuteron.
c 
c     See class notes and Arfken, Math. Methods for 
c     Physicists, 2nd Edition, section 9.1.2
c     for more details.
c===========================================================
      program        deut

      implicit       none

      character*4    cdnm
      parameter    ( cdnm = 'deut' )

      integer        iargc,         indlnb,         i4arg
      real*8         r8arg

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

c-----------------------------------------------------------
c     Order of system.
c-----------------------------------------------------------
      integer        neq
      parameter    ( neq = 2 )

c-----------------------------------------------------------
c     Storage for solution at requested output radii.
c-----------------------------------------------------------
      integer        maxout
      parameter    ( maxout = 10 000 )

      real*8         y0(neq)
      real*8         vxout(maxout), vyout(maxout,neq),
     &               work(maxout)
      integer        nxout,         ixout,       nxout_succ

      integer        ieq

      logical        ltrace
      parameter    ( ltrace = .true. )

      integer        maxdump
      parameter    ( maxdump = 50 )

      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-8 )

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

c-----------------------------------------------------------
c     Parse command line arguments.  Deviation from 
c     'sode'.  Usage is 
c
c     deut <x0> <E> [<tol>]
c-----------------------------------------------------------
      if( iargc() .lt. 2 ) go to 900

      x0      = r8arg(1,r8_never)
      E       = r8arg(2,r8_never)
      tol     = r8arg(3,default_tol)
      if( x0 .eq. r8_never   .or.  E .eq. r8_never )
     &    go to 900

c-----------------------------------------------------------
c     Both boundary conditions are fixed in this case.
c-----------------------------------------------------------
      y0(1) = 0.0d0
      y0(2) = 1.0d0

c-----------------------------------------------------------
c     Echo command line arguments if "local tracing" 
c     enabled ...
c-----------------------------------------------------------
      if( ltrace ) then
         write(0,*) 'x0:  ', x0
         write(0,*) 'E:   ', E
         write(0,*) 'tol: ', tol
      end if

c-----------------------------------------------------------
c     Get output radii from standard input ...
c-----------------------------------------------------------
      call dvfrom('-',vxout,nxout,maxout)
      if( nxout .le. 0 ) then
         write(0,*) 
         write(0,*) cdnm//': No output radii read from standard input'
         write(0,*) cdnm//': Use (e.g.) dvmesh or dvgmesh to '//
     &              'generate output radii.'
         write(0,*) 
         write(0,*) cdnm//': Sample usages:'
         write(0,*) 
         write(0,*) cdnm//': dvmesh  0.0 1.0 101'
         write(0,*) cdnm//': dvgmesh 0.0 1.0 101 10'
         write(0,*) 
         stop
      end if
      if( ltrace ) then
         if( nxout .le. maxdump ) then
            call dvdump(vxout,nxout,cdnm//': output radii',0)
         else 
            call dvdmp1(vxout,work,
     &                  int(1.0d0 * nxout / maxdump + 0.5d0),
     &                  nxout,cdnm//': selected output radii',0)
         end if
         write(0,*)
         write(0,*) cdnm//': Initial time: ', vxout(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-----------------------------------------------------------
      do ieq = 1 , neq
         y(ieq)       = y0(ieq)
         vyout(1,ieq) = y0(ieq)
      end do

c-----------------------------------------------------------
c     Do the integration ...
c-----------------------------------------------------------
      do ixout = 2 , nxout
         istate = 1
c-----------------------------------------------------------
c        Need these temporaries since lsoda overwrites 
c        tend ... 
c-----------------------------------------------------------
         tbgn = vxout(ixout-1)
         tend = vxout(ixout)
         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) cdnm, ixout, nxout, vxout(ixout), 
     &      vxout(ixout+1)
1000        format(' ',a,': Step ',i4,'  t = ',1pe10.3,
     &      ' .. ',1pe10.3)
            write(0,*) cdnm//': lsoda reurns ', istate
         end if

         if( istate .lt. 0 ) then
            write(0,1500) cdnm, istate, ixout, nxout, 
     &                    vxout(ixout-1), vxout(ixout)
1500        format(/' ',a,': Error return ',i2,
     &              ' from integrator LSODA.'/
     &              '       At output time ',i5,' of ',i5/
     &              '       Interval ',1pe11.3,' .. ',
     &             1pe11.3/)
             nxout_succ = ixout - 1
            go to 500
         end if
         do ieq = 1 , neq
            vyout(ixout,ieq) = y(ieq)
         end do
      end do
      nxout_succ = nxout

500   continue

      do ixout = 1 , nxout_succ
         write(*,2000)   vxout(ixout),
     &                 ( vyout(ixout,ieq) , ieq = 1 , neq )
c-----------------------------------------------------------
c        NOTE: This format will not dump all values on a 
c        single line if neq > 10 !!
c-----------------------------------------------------------
2000     format(1P,10E25.16)
      end do

      stop

900   continue
         write(0,*) 'usage: '//cdnm//
     &              ' <x0> <E> [<tol>]'
         write(0,*) ' '
         write(0,*) '       Output radii read from standard input'
      stop

      end

c===========================================================
c     Some double precision vector utility routines. 
c
c     Originally from ~matt/vutil/dveclib.f
c===========================================================
 
c-----------------------------------------------------------
c     Dumps vector labelled with LABEL on UNIT.
c-----------------------------------------------------------
      subroutine dvdump(v,n,label,unit)
         implicit       none

         real*8         v(*)
         character*(*)  label
         integer        i, n, st, unit
 
         if( n .lt. 1 ) go to 130
            write(unit,100) label
 100        format(/' <<< ',A,' >>>'/)
            st = 1
 110        continue
               write(unit,120) ( v(i) , i = st , min(st+3,n))
 120           FORMAT(' ',4(1PE19.10))
               st = st + 4
            if( st .le. n ) go to 110
 
 130     continue
 
         return
      end
 
c-----------------------------------------------------------
c     Extension of dvdump which dumps every 'inc'th element.
c-----------------------------------------------------------
      subroutine dvdmp1(v,w,inc,n,label,unit)
         implicit       none

         real*8         v(*),    w(*)
         character*(*)  label
         integer        inc,     n,      unit
 
         call dvinj(v,w,inc,n)
         call dvdump(w,1+(n-1)/inc,label,unit)

         return

      end
 
c---------------------------------------------------------
c     Injects every inc'th element of v1 into v2
c---------------------------------------------------------
      subroutine dvinj(v1,v2,inc,n)
         implicit    none
   
         real*8      v1(*),     v2(*)
         integer     i,         inc,     j,      n
 
         j = 1
         do i = 1 , n , inc
            v2(j) = v1(i)
            j = j + 1
         end do
 
         return
 
      end
