c     microsized version of an ``s-process" network. For teaching
c     purposes only.  One hundred isotopes are included
c     and the initial abundance of the first, nominally fe56 is taken
c     equal to one. Abundances are calculated for given exposures
c     tau. Though the ``y" notation is preserved, the abundances may 
c     equally well be thought of as ``n"s. The exposure may either be 
c     thought of as tau, the integral of neutron density, velocity, 
c     and time.

      dimension y(100), ynew(100), dy(100), signg(100)

    5 do 10 i=1,100
   10 y(i)=0.

c     for now take a very simple representation of the cross section
c     with fake `magic numbers' at 30 and 60 and odd-even effects

      do 11  i=2,100,2
   11 signg(i)=1.
      do 12  i=1,99,2
   12 signg(i)=2.

      signg(100)=0.
      signg(30)=.1
      signg(60)=.1
      y(1)=1.

      write (6,15)
   15 format (' input tau:')
      read (5,*) tau
      if (tau.lt.0.) go to 101
      dt=1.e-6*tau
      t=0.
      ifin=0

   20 if (t.gt.tau) go to 99
      ynew(1)=y(1)/(1.+signg(1)*dt) 
      dy(1)=y(1)-ynew(1)
      y(1)=ynew(1)
      do 30 i=2,100
      ynew(i)=(y(i)/dt+y(i-1)*signg(i-1))/(1./dt+signg(i))
      dy(i)=ynew(i)-y(i)
   30 y(i)=ynew(i)
      if (ifin.ne.0.) go to 999
      call dtnuc(dt,y,dy)
      t=t+dt
      go to 20
   99 ifin=1
      t=t-dt
      dt=tau-t
      go to 20

  999 t=t+dt
      write (6,40) t,tau
   40 format (2f10.2)
      do 50 i=1,100
      write (6,60) i,y(i)
   50 continue
   60 format (i5,1pe11.2)
      go to 5
  101 call exit(1)

      end   
cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

      subroutine dtnuc (dt,y,dy)

      dimension y(1),dy(1)
c
c
c
c     delchi is maximum fractional change allowed in number density
c     chimin is smallest y that can affect the time step
c     fdtn is maximum enlargement factor for new time step
c
      delchi=0.15
      chimin=.001
      fdtn=2.

      tau = fdtn*dt
      do 100 j = 1,100
      if(y(j)-chimin) 100,100,10
   10 taug = abs(y(j)/dy(j))*delchi*dt
      if(taug) 100,100,20
   20 if(tau-taug) 100,100,30
   30 tau = taug
  100 continue
      dt = tau

      return
      end