************************************************************************
*   FUNCTION:  Find a zero of a single nonlinear equation.             *
************************************************************************
      SUBROUTINE ZEROS(F,B,C,ABSERR,RELERR,RESIDL,FLAG)
C
C  ZERO computes a root of the nonlinear equation F(X) = 0 when F(X)
C  is a continuous real function of a single real variable X.  The
C  method used is the secant rule.
C
      DOUBLE PRECISION F,B,C,ABSERR,RELERR,RESIDL
      INTEGER FLAG
      EXTERNAL F
C
C  Input parameters:
C     F      = function program defining F(X).  This subprogram
C              must have the form:
C                 DOUBLE PRECISION FUNCTION F(X)
C                 DOUBLE PRECISION X
C                 F = ...
C                 RETURN
C                 END
C              The function name F must appear in an EXTERNAL
C              statement in the calling program.
C     ABSERR,RELERR = absolute and relative error tolerances.
C              The stopping criterion is:
C                ABS(B-C) .LE. 2.0*MAX(ABSERR,ABS(B)*RELERR).
C  Output parameters:
C     B,C    = see FLAG returns.
C     RESIDL = value of final residual F(B).
C     FLAG   = 0 for normal return; F(B)*F(C) < 0 and the
C                stopping criterion is met (or F(B)=0).  B is
C                always selected so that ABS(F(B)) .LE. ABS(F(C)).
C            = 1 if too many function evaluations were made; in
C                this version 500 are allowed.
C            = 2 if ABS(F(B)) is more than 100 times the larger
C                of the residuals for the input B and C.  B and C
C                probably approximate a pole rather than a root.
C            =-1 if ABSERR is less than zero, or RELERR too small.
C            =-2 if F(B)*F(C) is positive on input.
C
C  Local variables:
      DOUBLE PRECISION RES
      SAVE             RES
      INTEGER          COUNT,MAXF,NUMF
      PARAMETER       (MAXF = 500)
      DOUBLE PRECISION A,ACMB,CMB,FA,FB,FC,P,Q,TOL,U,WIDTH,
     &  C0,HALF,ONE,EIGHT,TEN,HUNDRD
      DATA C0/0.D0/,HALF/0.5D0/,ONE/1.D0/,EIGHT/8.D0/,TEN/10.D0/,
     &     HUNDRD/100.D0/
************************************************************************
*  Machine dependent constant:                                         *
*                                                                      *
*  Set U to the appropriate unit roundoff level for your machine;      *
*  here we assume IEEE double precision.                               *
      DATA U/2.22D-16/
************************************************************************
C
C     Test the input tolerances.
C
      IF ((RELERR .LT. TEN*U) .OR. (ABSERR .LT. C0)) THEN
         FLAG = -1
         RETURN
      ENDIF
C
C     Initialization.
C
      COUNT = 0
      WIDTH = ABS(B-C)
      A = C
      FA = F(A)
      NUMF = 1
      IF (ABS(FA) .EQ. C0) THEN
         FLAG = 0
         B = A
         RESIDL = C0
         RETURN
      ENDIF
      FB = F(B)
      NUMF = 2
      IF (ABS(FB) .EQ. C0) THEN
         FLAG = 0
         RESIDL = C0
         RETURN
      ENDIF
      IF (SIGN(ONE,FA) .EQ. SIGN(ONE,FB)) THEN
         FLAG = -2
c        RETURN
      ENDIF
      RES = MAX(ABS(FA),ABS(FB))
      FC = FA
   20 IF (ABS(FC) .LT. ABS(FB)) THEN
C
C        Interchange B and C so that ABS(F(B)) .LE. ABS(F(C)).
C
         A = B
         FA = FB
         B = C
         FB = FC
         C = A
         FC = FA
      ENDIF
      CMB = HALF*(C-B)
      ACMB = ABS(CMB)
      TOL = MAX(ABSERR,ABS(B)*RELERR)
C
C     Test the stopping criterion and function count.
C
      IF (ACMB .LE. TOL) THEN
         RESIDL = FB
         IF(ABS(RESIDL) .GT. HUNDRD*RES) THEN
            FLAG = 2
         ELSE
            FLAG = 0
         ENDIF
         RETURN
      ENDIF
      IF (NUMF .GE. MAXF) THEN
         FLAG = 1
         RETURN
      ENDIF
C
C     Calculate new iterate implicitly as B+P/Q where we arrange
C     P .GE. 0.  The implicit form is used to prevent overflow.
C
      P = (B-A)*FB
      Q = FA-FB
      IF (P .LT. C0) THEN
         P = -P
         Q = -Q
      ENDIF
C
C     Update A; check if reduction in the size of bracketing
C     interval is being reduced at a reasonable rate.  If not,
C     bisect until it is.
C
      A = B
      FA = FB
      COUNT = COUNT+1
      BISECT = .FALSE.
      IF (COUNT .GE. 4) THEN
         IF (EIGHT*ACMB .GE. WIDTH) THEN
            BISECT = .TRUE.
         ELSE
            COUNT = 0
            WIDTH = ACMB
         ENDIF
      ENDIF
C
C     Test for too small a change.
C
      IF (P .LE. ABS(Q)*TOL) THEN
C
C        If the change is too small, increment by TOL.
C
         B = B+SIGN(TOL,CMB)
      ELSE
C
C        Root ought to be between B and (C+B)/2.
C
         IF (P .LT. CMB*Q) THEN
C
C           Use secant rule.
C
            B = B+P/Q
         ELSE
C
C           Use bisection.
C
            BISECT = .TRUE.
         ENDIF
      ENDIF
      IF (BISECT) THEN
         B = C - CMB
      ENDIF
C
C     The computation for new iterate B has been completed.
C
      FB = F(B)
      NUMF = NUMF+1
      IF (ABS(FB) .EQ. C0) THEN
         FLAG = 0
         C = B
         RESIDL = FB
         RETURN
      ENDIF
      IF (SIGN(ONE,FB) .EQ. SIGN(ONE,FC)) THEN
         C = A
         FC = FA
      ENDIF
      GOTO 20
      END
C=======================================================================
