c-----------------------------------------------------------------------
c     driver for GAM  code
c-----------------------------------------------------------------------
      program GAMDO
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter(nd=500,lwork=2*nd*nd+40*nd+18,liwork=nd+24)
      dimension y(nd),yref(nd),work(lwork),iwork(liwork)
      external fevalr, jevalr, solout
      character driver*5, problm*8, solver*5
      parameter (driver = 'GAMDO', solver='GAM')
c-----------------------------------------------------------------------
c     get the problem dependent parameters
c-----------------------------------------------------------------------
      call prob(problm,neqn,tbegin,tend,ijac,mljac,mujac)
c-----------------------------------------------------------------------
c     get the initial values
c-----------------------------------------------------------------------
      call init(neqn,y)
      itol = 0
      iout = 0
c-----------------------------------------------------------------------
c     read the tolerances and initial stepsize
c-----------------------------------------------------------------------
 
      write(6,*) 'give relative error tolerance'
      read(5,*)   rtol
      write(6,*) 'give absolute error tolerance'
      read(5,*)   atol
      write(6,*) 'give initial  stepsize '
      read(5,*)   h0
      h = h0
      t0   = tbegin
      tf   = tend
      do 10 i=1, 18 
         iwork(i) = 0
         work(i)  = 0d0
 10   continue
c-----------------------------------------------------------------------
c     call of the subroutine GAM
c-----------------------------------------------------------------------
      call GAM(neqn,fevalr,t0,y,tf,h,
     &                  rtol,atol,itol,
     &                  jevalr ,ijac,mljac,mujac,
     &                  solout,iout,
     &                  work,lwork,iwork,liwork,rpar,ipar,idid)
      if (idid.ne.1) then
         write(6,*) 'ERROR: tabim returned idid =', idid
      endif
c-----------------------------------------------------------------------
c     print final solution
c-----------------------------------------------------------------------
      write(6,11)rtol,atol,h0,tf
   11 format(/,' we solved the problem with',//,
     +       ' relative error tolerance = ',d10.4,',',/,
     +       ' absolute error tolerance = ',d10.4,/,
     +       '     and initial stepsize = ',d10.4,//,
     +       ' the numerical solution at t = ',d22.14,' :',/)
      do 20 i=1,neqn
      write(6,21)i,y(i)
   20 continue
   21 format(' y(',i3,') = ',d22.14)
c-----------------------------------------------------------------------
c     print error with respect to reference solution
c-----------------------------------------------------------------------
      call solut(neqn,yref)
      write(6,29)
   29 format(/)
      do 30 i=1,neqn
   30 write(6,31)i,abs(y(i)-yref(i))
   31 format(' absolute error in component ',i3,' = ',d22.14)
c-----------------------------------------------------------------------
c     print statistics
c-----------------------------------------------------------------------
      NSTEPS = 0
      DO I=12,23
       NSTEPS = NSTEPS + iwork(I)
      END DO
      NACCEPT = IWORK(12)+IWORK(13)+IWORK(14)+IWORK(15)
      write(6,41) NSTEPS,NACCEPT,iwork(10),iwork(11),iwork(24)
   41 format(  ' # steps  =         ',i8,/,
     +         ' # accept =         ',i8,/,
     +         ' # f-eval =         ',i8,/,
     +         ' # Jac-eval =       ',i8,/,
     +         ' # LU-decomp =      ',i8/)
      end
c-----------------------------------------------------------------------
c     since in GAM  the format of the subroutines containing the
c     derivative function and its jacobian  differs from the format in
c     the testset, we transform them
c-----------------------------------------------------------------------
      subroutine fevalr(neqn,t,y,dy,rpar,ipar)
      double precision t,y,dy,rpar
      integer neqn,ipar
      dimension y(neqn),dy(neqn)
      call feval(neqn,t,y,dy)
      return
      end
      subroutine jevalr(neqn,t,y,jac,ldim,rpar,ipar)
      double precision t,y,jac,rpar
      integer neqn,ldim,ipar
      dimension y(neqn),jac(ldim,neqn)
      call jeval(neqn,t,y,jac,ldim)
      return
      end
      subroutine SOLOUT(R,TP,YP,F1,NT1,DBLK,ORD,RPAR,IPAR,IRTRN)
      implicit none
      integer R, DBLK, ORD, IPAR(*), IRTRN, NT1
      double precision TP(*),YP(R,*),RPAR(*),F1(R,*)
      RETURN
      END