subroutine si(ni,y,ysc,t,h,eps,hdone,ierr)
C     
C  This routine does Semi-Implicit integration of an ODE whose r.h.s.
C  is given by an external function func1. Uses Gnedin & Gnedin method
C  from Astrophysical Journal 1998, 509, 11. The equation(s) to be solved
C  have to be in the following form:
C
C   dy/dt = -a*y + s
C
C  Parameters:
C  
C  ni      [in]:     number of equations.
C  y(ni)   [in/out]: current value of the vector of solutions.
C  ysc(ni) [in]:     error is scaled to ysc: abs(y) < eps*ysc. Has to
C                    be positive.
C  t       [in/out]: current values of the independent variable.
C  h       [in/out]: on the input, the trial size of the time-step. On the
C                    output, the new trial time-step for the next step.
C  eps     [in]:     the requested error relative to ysc.
C  hdone   [out]:    the time-step actually taken.
C  ierr    [out]:    error index: ierr=0 means success; ierr=-1 means
C                    ni is too large, the internal parameter NE needs to
C                    be increased; ierr>0 means a failure, in this case
C                    ierr gives the number of the offending equation.
C
C  This routine needs another routine with the following signature:
C
C      subroutine func1(ni,y,f,a,t)
C
C  where:
C  
C    ni    [in]:  number of equations.
C    y(ni) [in]:  current value of the vector of solutions.
C    f(ni) [out]: value of the r.h.s.
C    a(ni) [in]:  value of the incomplete Jacobian a
C    t     [in]:  current value of the independent variable.
C
C *******************************************************
C
C
      parameter(gam = 0.788675134594812882251)
      parameter(a21 = 1.0)
      parameter(a31 = 0.56698729810778067662)
      parameter(a32 = 0.25)
      parameter(c21 = -1.26794919243112270647)
      parameter(c31a = -1.5) 
      parameter(c32a = -1.1830127018922193234)
      parameter(b1a = 1.3779915320718537844)
      parameter(b2a = 0.9553418012614795489)
      parameter(b3a = 2.0/3.0)
C
      parameter(b3b = 0.21132486540518711775)
C
      parameter(c31b = -0.71891108675446473633 -
     .     0.52072594216369017578/b3b)
      parameter(c32b = -0.31698729810778067662 - 
     .     0.5773502691896257645/b3b)
      parameter(b1b = 1.654700538379251529 - 0.4150635094610966169*b3b)
      parameter(b2b = 1.077350269189625765 - 0.1830127018922193234*b3b)
C
C  itmax is the allowed maximum number of attempts to decrease the time-step
C  and redo the calculation is the error tolerance is not satisfied.
C
      parameter(itmax = 99)
C
C  NE is the maximum number of equation accepted.
C
      parameter(NE=1000)
C
C  Input arrays
C
      dimension y(*), ysc(*)
C
C  Internal storage
C
      dimension y0(NE), f0(NE), a0(NE), aa(NE), f(NE)
      dimension g1(NE), g2(NE), g3a(NE), g3b(NE)
C
C  First: set ierr to zero, check that the internal limits for the number
C         of equations is sufficient, and store the initial values of y
C         and t. Reset the iteration counter it.
C
      ierr = 0

      if(ni .gt. NE) then
         ierr = -1
         return
      endif

      do j=1,ni
         y0(j) = y(j)
      enddo

      t0 = t
      it = 0
C
C  Calculate the r.h.s. and incomplete Jacobian at the initial values 
C  of y and t. Needs to be done only once per call.
C
      call func1(ni,y0,f0,a0,t)
C
C  Start the iteration loop. If the initial h is not sufficient, we come
C  back here, if it is, we do not come back.
C
 10   it = it + 1
C
C  Set g1
C
      do j=1,ni
         g1(j) = h*f0(j)/(1.0+gam*h*a0(j))
      enddo
C
C  Set g2
C
      t = t0 + a21*h
      do j=1,ni
         y(j) = y0(j) + a21*g1(j)
      enddo
      call func1(ni,y,f,aa,t)
      do j=1,ni
         g2(j) = (h*f(j)+c21*g1(j))/(1.0+gam*h*a0(j))
      enddo
C
C  Set g3(a,b)
C
      t = t0 + (a31+a32)*h
      do j=1,ni
         y(j) = y0(j) + a31*g1(j) + a32*g2(j)
      enddo
      call func1(ni,y,f,aa,t)
      do j=1,ni
         g3a(j) = (h*f(j)+c31a*g1(j)+c32a*g2(j))/(1.0+gam*h*a0(j))
         g3b(j) = (h*f(j)+c31b*g1(j)+c32b*g2(j))/(1.0+gam*h*a0(j))
      enddo
C
C  Find the new y and at the same time find the maximum error errmax and 
C  the equation number jmax where it is reached.
C
      errmax = -1.
      jmax = 0
      do j=1,ni
         y(j) = y0(j) + b1a*g1(j) + b2a*g2(j) + b3a*g3a(j)
         errabs = (b1a-b1b)*g1(j) + (b2a-b2b)*g2(j) +
     .        b3a*g3a(j) - b3b*g3b(j)
         err = abs(errabs)/ysc(j)
         if(err .gt. errmax) then
            errmax = err
            jmax = j
         endif
      enddo
C
C  Check the maximum error and the tolerance.
C
      if(errmax .gt. eps) then
C
C  Error is too large. Reduce the time-step, reset t, and start again
C  unless we have reached the iteration limit. In that case report the 
C  failure.
C
         t = t0
         h = 0.9*h*(eps/errmax)**0.5
         if(it .lt. itmax) goto 10
         ierr = jmax
         hdone = 0.0
         return
      else
C
C  Error is acceptable. Update the independent variable, compute the 
C  estimate for the next time-step, and exit.
C
         t = t0 + h
         hdone = h
         h = 0.9*h*(eps/max(0.04*eps,errmax))**0.3333
      endif
C
C  We are done!
C
      return
      end