program gamd

      integer md
      parameter (md=400)
      integer lwork, liwork
      parameter (lwork=2*md*md+40*md+18,liwork=md+24)

      integer neqn,ndisc,mljac,mujac,mlmas,mumas,ind(md),
     +        iwork(liwork),ipar(2),idid
      double precision y(md),dy(md),t(0:100),
     +                 h0,rtol(md),atol(md),
     +                 work(lwork),rpar(1)
      logical numjac,nummas,consis,tolvec

      integer itol,ijac,imas,iout
      double precision h

      double precision yref(md)
      character fullnm*40, problm*8, type*3
      character driver*8, solver*8
      parameter (driver = 'gamd', solver='GAM')
      external odef,odejac,odemas
      external solout

      double precision solmax
      real gettim, timtim, cputim
      external gettim
      double precision scd, mescd

      integer i,j,icount(10:24)

      do 10 i=1,18
         iwork(i) = 0
         work(i)  = 0d0
   10 continue
c
c     NMAX , the maximal number of steps
c
      iwork(2) = 1000000
      iout = 0

c-----------------------------------------------------------------------
c     check the problem definition interface date
c-----------------------------------------------------------------------

      call chkdat(driver,solver,20030829)

c-----------------------------------------------------------------------
c     get the problem dependent parameters
c-----------------------------------------------------------------------

      call prob(fullnm,problm,type,
     +          neqn,ndisc,t,
     +          numjac,mljac,mujac,
     +          nummas,mlmas,mumas,
     +          ind)
      if (type.eq.'IDE') then
         print *, 'GAMD: ERROR: ',
     +            'GAM can not solve IDEs'
         stop
      elseif (type.eq.'ODE') then
         imas = 0
      elseif (type.eq.'DAE') then
         print *, 'GAMD: ERROR: ',
     +            'GAM can not solve DAEs'
         stop
      else
         print *, 'GAMD: ERROR: ',
     +            'unknown Test Set problem type', type
         stop
      endif
      if (numjac) then
         ijac = 0
      else
         ijac = 1
      endif

c-----------------------------------------------------------------------
c     get the initial values
c-----------------------------------------------------------------------

      call init(neqn,t(0),y,dy,consis)

c-----------------------------------------------------------------------
c     read the tolerances and initial stepsize
c-----------------------------------------------------------------------

      call getinp(driver,problm,solver,fullnm,
     +            tolvec,rtol,atol,h0,solmax)

      if (tolvec) then
         itol = 1
      else
         itol = 0
      endif

c-----------------------------------------------------------------------
c     call of the subroutine GAM
c-----------------------------------------------------------------------

      do 40 j=10,24
         icount(j) = 0
   40 continue

      timtim = gettim()
      timtim = gettim() - timtim

      cputim = gettim()

      do 60 i=0,ndisc

         h = h0

         if (type.eq.'ODE') then
            call gam(neqn,odef,t(i),y,t(i+1),h,
     +               rtol,atol,itol,
     +               odejac ,ijac,mljac,mujac,
     +               solout,iout,
     +               work,lwork,iwork,liwork,rpar,ipar,idid)
         endif

         if (idid.ne.1) then
            print *, 'GAMD: ERROR: ',
     +               'GAM returned IDID = ', idid
            stop
         endif

         do 50 j=10,24
            icount(j) = icount(j) + iwork(j)
   50    continue

   60 continue

      cputim = gettim() - cputim - timtim

      do 70 j=10,24
         iwork(j) = icount(j)
   70 continue

c-----------------------------------------------------------------------
c     print numerical solution in endpoint and integration statistics
c-----------------------------------------------------------------------

      call solut(neqn,t(ndisc+1),yref)
c     
 call getscd(mescd,scd,neqn,yref,y,problm,tolvec,atol,rtol)
      call report(
     +   driver,problm,solver,
     +   rtol,atol,h0,solmax,
     +   iwork,cputim,scd,mescd
     +)

      end

      subroutine solout(nr,xold,x,y,cont,lrc,n,rpar,ipar,irtrn)
      integer nr,lrc,n,ipar(*),irtrn
      double precision xold,x,y(n),cont(lrc),rpar(*)
      return
      end

C=======================================================================
C     `Test Set for IVP Solvers' ODE wrappers for GAM
C=======================================================================
c
c     since in GAM the format of the subroutines for the
c     function f and the Jacobian differ from the format
c     in the testset, we transform them
c
c-----------------------------------------------------------------------
      subroutine odef(n,x,y,f,ierr,rpar,ipar)
      integer n,ipar(*)
      double precision x,y(n),f(n),rpar(*)
      integer ierr
   
      call feval(n,x,y,y,f,ierr,rpar,ipar)
    
      return
      end

      subroutine odejac(n,x,y,dfy,ldfy,rpar,ipar)
      integer ldfy,n,ipar(*)
      double precision x,y(n),dfy(ldfy,n),rpar(*)
      integer ierr
      double precision dy
      ierr = 0
      call jeval(ldfy,n,x,y,dy,dfy,ierr,rpar,ipar)
      if (ierr.ne.0) then
         print *, 'GAMD: ERROR: ',
     +            'GAM can not handle JEVAL IERR'
         stop
      endif
      return
      end

      subroutine odemas(n,am,lmas,rpar,ipar)
      integer n,lmas,ipar(*)
      double precision am(lmas,n),rpar(*)
      return
      end