c   example program for epsode -- diurnal chemistry problem
c     taken from
c        hindmarsh and byrne, applications of episode: an experimental
c           package for the integration of systems of ordinary differ-
c           ential equations, from numerical methods for differential
c           systems, academic press 1976.
c
      implicit real*8 (a-h,o-z)
      dimension y0(2)
      common/epcom9/ hused, nqused, nstep, nfe, nje
      common/pcom/ q1, q2, c3, c4, freq, y3
      external f, jacob
c
c     call errset (208, 256, -1, 1, 0, 0)
      n = 2
      t0 = 0.0d0
      h0 = 1.0d-8
      y0(1) = 1.0d6
      y0(2) = 1.0d12
      eps = 1.0d-6
      ierror = 2
      mf = 21
      q1 = 1.63d-16
      q2 = 4.66d-16
      c3 = 22.62d0
      c4 = 7.601d0
      freq = 3.141592653589793d0 / 43200.0d0
      y3 = 3.7d16
      write(6,1000) mf, eps
      nout = 26
      tout = 7200.0d0
      index = 1
      do 60 iout = 1, nout
         call epsode(f, jacob, n, t0, h0, y0, tout, eps,
     *               ierror, mf, index)
         write(6,1001) t0, nstep, nfe, nje, nqused, h0, y0(1), y0(2)
         if (index .ne. 0) stop
         if (iout .lt. 6) tout = tout + 7200.0d0
         if (iout .eq. 6) tout = tout + 21600.0d0
         if (iout .gt. 6) tout = tout + 43200.0d0
   60 continue
      stop
 1000 format(' diurnal chemistry problem'/
     *       '0mf =',i3,'   eps =',1pd10.1/
     *       '-   t',8x,'nstep',3x,'nfe',3x,'nje',2x,'nq',4x,'h',10x,
     *       'y1',10x,'y2'//)
 1001 format(1x,1pd11.3,3i6,i4,d10.2,2d12.4)
      end
      subroutine f(n, t, y, ydot)
      implicit real*8 (a-h,o-z)
      dimension y(2), ydot(2)
      common/epcom2/ ymax(20)
      common/pcom/ q1, q2, c3, c4, freq, y3
c
      call diurn (t, q3, q4)
      qq1 = q1 * y(1) * y3
      qq2 = q2 * y(1) * y(2)
      qq3 = q3 * y3
      qq4 = q4 * y(2)
      ydot(1) = -qq1 - qq2 + 2.0d0 * qq3 + qq4
      ydot(2) = qq1 - qq2 - qq4
      ymax(1) = dmax1(dabs(y(1)), 1.0d-20)
      ymax(2) = dmax1(dabs(y(2)), 1.0d-20)
      return
      end
      subroutine jacob(n, t, y, pd, n0)
      implicit real*8 (a-h,o-z)
      dimension y(2), pd(n0,2)
      common/pcom/ q1, q2, c3, c4, freq, y3
c
      call diurn (t, q3, q4)
      pd(1,1) = -q1 * y3 - q2 * y(2)
      pd(1,2) = -q2 * y(1) + q4
      pd(2,1) = q1 * y3 - q2 * y(2)
      pd(2,2) = -q2 * y(1) - q4
      return
      end
      subroutine diurn (t, q3, q4)
      implicit real*8 (a-h,o-z)
      common/pcom/ q1, q2, c3, c4, freq, y3
c
      s = dsin(freq * t)
      if (s .le. 0.0d0) go to 20
      q3 = dexp(-c3 / s)
      q4 = dexp(-c4 / s)
      return
   20 q3 = 0.0d0
      q4 = 0.0d0
      return
      end