      SUBROUTINE cnjgrd (n, x, stplim, f, h, ih, g, ifail, x0, g0, s)

      IMPLICIT UNDEFINED (A-Z)
      INTEGER n, ih, ifail
      REAL    h (ih, n), x (n), x0 (n), g0 (n), s (n), g (n)

C     Minimisation of a function by conjugate gradient (?) method of
C     R. Fletcher and M. Powell, Comput. J., vol 6 (1963) 163.

C     The function to be minimised should be defined in the 
C     subroutine 'func'.  The derivatives of the function are
C     also required to be defined in 'func'.  This subroutine
C     calls three auxilliary routines 'linmin', 'linzro' and 'linfunc', 
C     as well as 'func'.

C     Parameters to this routine are:
C
C          n         The number of variables
C          x (n)     On entry, initial estimate of location of
C                    minimum.  On exit holds location of the minimum
C                    discovered
C          stplim    On entry, maximum step length in any direction
C                    in single step.  This MUST be set! (making
C                    'stplim' small helps avoid overflows, etc.,
C                    in 'func').
C          f         On exit, contains the value of the function at
C                    the minimum discovered.
C          h (ih, n), 
C          ih >= n   On exit contains estimated value of inverse of
C                    matrix of second derivatives at the minimum.
C          g (n)     On exit contains gradient of function at
C                    "minimum" (should be rather small!).
C          ifail     On exit:
C                      ifail = 0 => Converged ok
C                      ifail = 1,...,4 => Less satisfactory convergence,
C                                         in increasing order of 
C                                         dodginess.
C                      ifail = 5 => n <= 0, ih < n or stplim <= 0.
C                                   Routine aborted.
C           x0 (n),
C           g0 (n),
C           s  (n)    Vectors used for workspace.
C
C     ***********************************************************
C     * 30/Jan/86.  Possible improvements, not yet implemented: *
C     * It ought to be possible to unify the two different      *
C     * linear search routines to get a more efficient search.  *
C     * Problem is that, by design, they do different things.   *
C     * 'linmin' is guaranteed to reduce the function, whereas  *
C     * 'linzero' is guaranteed to isolate a zero of the deriv- *
C     * ative.  In general these are incompatible goals due to  *
C     * numerical inaccuracies.  Also the treatment of the case *
C     * where denominators vanish in updatine H needs thinking  *
C     * about more carefully.                                   *
C     ***********************************************************
C
C     ***********************************************************
C     * 4/Jun/90.  This program has been resurrected from a     *
C     * listing.  I've taken the opportunity to improve some    *
C     * of the "lexicographic" conventions (whitespace and      *
C     * indentation, keywords in upper case, etc).  The only    *
C     * significant change (mainly to trap typing errors) is    *
C     * that the subroutines now have "implicit undefined" set, *
C     * so many variables have had to have declarations added.  *
C     * In retrospect the programming style looks a bit odd,    *
C     * but I don't really want to risk rewriting it.           *
C     ***********************************************************

        REAL    lambd0, linstl, lineps
        LOGICAL contin, sttnry, rptmin, test1, test2, test3

        LOGICAL debug
        COMMON /debug/ debug
        DATA    debug /.FALSE./

C       Defined constants:
C
C           eps      Required maximum estimated distance from 
C                    minimum on exit.
C           maxevn   Maximum number of calls to 'func' (per variable)
C           rnderr   A (generous) upper bound on the relative rounding
C                    error in real arithmetic

        REAL    eps, rnderr
        INTEGER maxevn
        DATA eps /1E-5/
        DATA maxevn /500/
        DATA rnderr /1E-4/

        contin = .TRUE.

        IF (n .LE. 0) THEN
          WRITE (6, *) '***conjgrd: n less than or equal to 0.'
          ifail  = 5
          contin = .FALSE.
        ELSEIF (ih .LT. n) THEN
          WRITE (6, *) '***conjgrd: ih is smaller than n.'
          ifail  = 5
          contin = .FALSE.
        ELSEIF (stplim .LE. 0.0) THEN
          WRITE (6, *) '***conjgrd: stplim less than or equal to 0.'
          ifail  = 5
          contin = .FALSE.
        ENDIF

        IF (contin) THEN
          ncount = n * maxevn

C         Initialise the matrix 'H' to unit matrix.

          DO 20 i = 1, n
            DO 10 j = 1, n
              h (i, j) = 0.0
 10         CONTINUE
            h (i, i) = 1.0
 20       CONTINUE

          hnorm  = n * 1.0
          lstrst = 0
          icount = 0
          CALL func (n, x, f, g)
          icount = icount + 1

          sttnry = .TRUE.
          DO 30 i = 1, n
            sttnry = sttnry .AND. (g (i) .EQ. 0.0)
 30       CONTINUE

          IF (stnnry) THEN
            WRITE (6, 1000)
 1000       FORMAT (' ***cnjgrd: Gradient vanished less than n',
     1              'iterations after resetting H.')
            ifail  = 1
            contin = .FALSE.
          ENDIF

        ENDIF

C       Start of main iteration... "WHILE (contin)"

 40     IF (.NOT. contin) GOTO 150

          IF (debug) THEN
            WRITE (6, *) 'h:'
            DO 45 i = 1, n
              WRITE (6, *) (h (i, j), j = 1, n)
 45         CONTINUE
          ENDIF

          DO 50 i = 1, n
            x0 (i) = x (i)
            g0 (i) = g (i)
 50       CONTINUE

          f0 = f

C         Calculate 's' and <g0 | s>.
C
C         <g0 | s> >= 0 implies H has lost its positivity, due to
C         either rounding errors (for example, caused by too large
C         an initial estimate of 'H' being suddenly greatly reduced
C         by subtraction), or else the approximate nature of the
C         minimisation in 'linfun'.  If this condition occurs,
C         'H' is reset to a constant matrix.

          DO 60 i = 1, n
            s (i) = 0.0
            DO 55 j = 1, n
              s (i) = s (i) - h (i, j) * g0 (j)
 55         CONTINUE
 60       CONTINUE

          g0s = 0.0
          DO 65 i = 1, n
            g0s = g0s + g0 (i) * s (i)
 65       CONTINUE

C         Calculate a norm for current 'H'

          hnorm0 = hnorm

          hnorm  = 0.0
          DO 75 i = 1, n
            DO 70 j = 1, n
              hnorm = hnorm + ABS (h (i, j))
 70         CONTINUE
 75       CONTINUE

          IF ((hnorm .LE. rnderr * hnorm0) .OR. (g0s .GE. 0.0)) THEN

            lstrst = 0

            IF (debug) THEN
              WRITE (6, 1005)
 1005         FORMAT (' ***cnjgrd: H reset to constant matrix.')
            ENDIF

C           'H' is reset to a hopefully intelligently chosen multiple
C           of the unit matrix.  This helps to avoid reoccurence of 
C           rounding errors that spoil positivity.

            g0s   = 0.0
            hnorm = rnderr * hnorm0

            DO 85 i = 1, n
              DO 80 j = 1, n
                h (i, j) = 0.0
 80           CONTINUE
              h (i, i) = hnorm / n
              s (i)    = - h (i, i) * g0 (i)
              g0s      = g0s + g0 (i) * s (i)
 85         CONTINUE

            IF (debug) THEN
              WRITE (6, *) 'h:'
              DO 90 i = 1, n
                WRITE (6, *) (h (i, j), j = 1, n)
 90           CONTINUE
            ENDIF

          ENDIF

          abssmx = 0.0
          DO 95 i = 1, n
            abssmx = MAX (ABS (s (i)), abssmx)
 95       CONTINUE

C         Is 's' vanishes routine exits.  If 'n' or more improvement
C         iterations have been made to the current 'H', this is 
C         accepted as proper convergence.  (There is a potential
C         problem if these iterations were not consecutive due
C         to vanishing denominators.  This isn't sorted yet).

          sttnry = abssmx .EQ. 0.0

          IF (sttnry .AND. (lstrst .LT. n)) THEN
            WRITE (6, 1010)
 1010       FORMAT (' ***cnjgrd: Predicted distance to minimum ',
     1              'vanished less than' /
     2              'n iterations after resetting H.')
            ifail  = 2
            contin = .FALSE.
          ELSEIF (sttnry) THEN
            ifail  = 0
            contin = .FALSE.
          ENDIF

C         Make an estimate of the minimum of the function along the
C         s-direction from the present position.  The cubic
C         interpolation method suggested by Fletcher and Powell is
C         used (in routine 'linmin').  Initial choice of 'alpha' is
C         fairly arbitrary.

          IF (contin) THEN
            linstl = stplim / abssmx
            alpha0 = MIN (1.0, linstl)

            IF (debug) THEN
              WRITE (6, *) 'icount =', icount
              WRITE (6, *) 'x0     =', x0
              WRITE (6, *) 'f      =', f
              WRITE (6, *) 's      =', s
              WRITE (6, *) 'alpha0 =', alpha0
            ENDIF

            CALL linfun (alpha0, f, g1s, n, x0, s, x, g)
            icount = icount + 1

            IF (abssmx .GT. eps) THEN
              lambd0 = 0.0
              alpha  = alpha0
              f1     = f0
              df1    = g0s
              f2     = f
              df2    = g1s
              CALL linmin (lambd0, f1, df1, alpha, f2, df2, linstl,
     1                     jfail, icount, ncount, n, x0, s, x, g)
              rptmin = ((ABS (alpha) * abssmx) .LE. eps) .OR.
     1                 ((ABS (alpha - lambd0) * abssmx) .LE. eps) .OR.
     2                 (jfail .GT. 0)
            ENDIF

          ENDIF

C         In the following case, the smallness of alpha * s is
C         confirmed or refuted by a routine which simply
C         searches for zeroes of the derivative along the line.
C         This has been found necessary because inconsistent
C         rounding errors in evaluating the function and its
C         derivative can otherwise lead to spuriously small
C         values of alpha being predicted by cubic interpolation.

          IF (contin .AND. ((abssmx .LE. eps) .OR. rptmin)) THEN
            lambd0 = 0.0
            alpha  = alpha0
            f1     = f0
            df1    = g0s
            f2     = f
            df2    = g1s
            lineps = eps / abssmx
            CALL linzro (lambd0, f1, df1, alpha, f2, df2, linstl,
     1                   lineps, jfail, icount, ncount, n,
     2                   x0, s, x, g)
            IF (jfail .NE. 0) THEN
              WRITE (6, *) '***cnjgrd: failure of linzro.'
              ifail  = 3
              contin = .FALSE.
            ENDIF
          ENDIF

          IF (contin) THEN
            f  = f2
            df = df2

C           Test for convergence of algorithm.

            sttnry = .TRUE.
            DO 100 i = 1, n
              sttnry = sttnry .AND. (g (i) .EQ. 0.0)
 100          CONTINUE

            IF (sttnry .AND. (lstrst .LT. n)) THEN
              WRITE (6, 1000)
              ifail  = 1
              contin = .FALSE.
            ELSEIF (sttnry) THEN
              ifail  = 0
              contin = .FALSE
            ENDIF
          ENDIF

          IF (contin) THEN
            test1 = abssmx .LT. eps
            test2 = (ABS (alpha) * abssmx) .LT. eps
            test3 = lstrst .GE. n

            IF (test1 .AND. test2 .AND. test3) THEN
              ifail  = 0
              contin = .FALSE.
            ELSEIF (icount .GE. ncount) THEN
              WRITE (6, *)
   1           '***cnjgrd: Too many function evaluations required.'
              ifail  = 4
              contin = .FALSE.
            ENDIF
          ENDIF

          IF (contin) THEN

C           Improvement of the matrix 'H'.
C           New usage of vectors (in notation of F & P) is:
C               g0  -->  y
C               x0  -->  H * y

            sy = 0.0
            DO 110 i = 1, n
              g0 (i) = g (i) - g0 (i)
              sy     = sy + (s (i) * g0 (i))
 110          CONTINUE

            yhy = 0.0
            DO 120 i = 1, n
              x0 (i) = 0.0
              DO 115 j = 1, n
                x0 (i) = x0 (i) + (h (i, j) * g0 (j))
 115            CONTINUE
              yhy = yhy + (g0 (i) * x0 (i))
 120          CONTINUE

C           If either of the denominators vanishes, H is not changed.

            IF ((sy .NE. 0.0) .AND. (yhy .NE. 0.0)) THEN

              lstrst = lstrst + 1
              asyinv = alpha / sy
              yhyinv = 1.0 / yhy

C             Update 'H'.  (NB. Order in which arithmetic is performed
C             here is important to maintain strict symmetry of 'H',
C             due to non-associativity of rounded arithmetic.
C             Hence brackets around s * s, x0 * x0).

              DO 140 i = 1, n
                DO 135 j = 1, n
                  h (i, j) = h (i, j) + asyinv * (s (i) * s (j)) -
     1                       yhyinv * (x0 (i) * x0 (j))
 135              CONTINUE
 140            CONTINUE

            ELSEIF (debug) THEN
              IF (sy .EQ. 0.0) THEN
                WRITE (6, *) '***cnjgrd: <s(i) | y(i)> vanished.'
              ELSEIF (yhy .EQ. 0.0) THEN
                WRITE (6, *) '***cnjgrd: <y(i) | H(i) | y(i)> vanished.'
              ENDIF
            ENDIF

          ENDIF

          GOTO 40
 150      CONTINUE

        RETURN
      END

      SUBROUTINE linmin (x1, f1, df1, x2, f2, df2, stplim, ifail,
     1                   icount, ncount, n, wk1, wk2, wk3, wk4)

      IMPLICIT UNDEFINED
      REAL    x1, x2, f1, f2, stplim
      INTEGER ifail, icount, ncount
      INTEGER n
      REAL    wk1 (n), wk2 (n), wk3 (n), wk4 (n)

C     Estimate location of minimum of function along a line.
C     A routine for evaluating the function and its derivative
C     must be provided (linfun).  The method is based on cubic
C     interpolation and starts from the given points 'x1', 'x2'.
C
C     Notes:
C         i) On entry 'f1', 'df1', 'f2' and 'df2' must be set, and
C            it is required that SIGN (df1) = - SIGN (x2 - x1), (ie,
C            the function slopes down into the [x1, x2] interval from
C            the 'x1' side).  The slope at 'x2' is unconstrained.
C        ii) On exit 'x2' is the best estimate of the location of the
C            minimum.  'x1' is a "second best estimate", and in general
C            is changed from entry value.  'f1', 'df1', 'f2', 'df2'
C            have appropriate values.
C       iii) Except for cubic or lower polynomials the location
C            of the minimum is only approximate.  The routine
C            can adaptively change the interpolation interval,
C            but exits as soon as cubic interpolation (or
C            extrapolation) predicts a minimum at a point
C            where the evaluated function is indeed smaller than
C            at the endpoints of the interval.
C        iv) On exit, the condition in i) on the signs of 'df1' and 
C            x2 - x1 still holds, so this routine may be called
C            iteratively to locate the minimum accurately.
C
C     Other parameters:
C         stplim     On entry, maximum stepsize to be taken in
C                    a single step
C         ifail      On exit:
C                      ifail = 0 => Ended OK.
C                      ifail > 0 => Didn't end OK.
C         icount     Counter of number of 'linfun' calls so far.
C         ncount     Maximum number of 'linfun' calls.
C         n,
C         wk1,
C         ...,
C         wk4,       Parameters and workspace passed to
C                    routine 'linfun'.

        REAL num, den
        LOGICAL contin, okextr

        LOGICAL debug
        COMMON /debug/debug

C       Defined constants:
C           expfac     Determines factor by which interval may be 
C                      expanded in a single search step.

        DATA expfac /5.0/

        IF (debug) THEN
          WRITE (6, *) 'entering linmin'
        ENDIF

        contin = .TRUE.

C       WHILE (contin)...

  10    IF (.NOT. contin) GOTO 20
          delint = x2 - x1

          IF (delint .EQ. 0.0) THEN
            WRITE (6, *) '***linmin: Interpolation interval is zero'
            ifail  = 1
            contin = .FALSE.
          ENDIF

          IF (contin) THEN
            sgnint = SIGN (1.0, delint)

            IF (.NOT. ((sgnint * df1) .LE. 0.0)) THEN
              WRITE (6, *) '***linmin: bad sign in linear minimisation'
              ifail  = 2
              contin = .FALSE.
            ENDIF

          ENDIF

          IF (contin .AND. (f2 .LT. f1) .AND. 
     1                    ((sgnint * df2) .LT. 0.0)) THEN

C           In this case the function appears to be monotonically
C           decreasing in the interval.  Cubic extrapolation is
C           attempted, but if this goes outside a range determined
C           by constant 'expfac', the interval is just moved
C           onto the 'x2' side and expanded by a factor 'expfac'.

            outlim = x2 + sgnint * MIN (expfac * ABS (delint), stplim)

C           (Note delint <> 0 was checked earlier in routine).

            z  = (3 * (f1 - f2) / delint) + df1 + df2
            w2 = (z * z) - (df1 * df2)
            okextr = w2 .GE. 0.0

            IF (okextr) THEN
              w   = sgnint * SQRT (ABS (w2))
              num = df2 + w - z
              den = df2 - df1 + (2 * w)
              okextr = den .NE. 0.0
            ENDIF

            IF (okextr) THEN
              sgnden = SIGN (1.0, den)

C             (Following test fails if extrapolation is to point outside
C             range between 'x2' and point delint * expfac to 'x2' side
C             of present interval).

              okextr = ((sgnden * num) .GT. - (sgnden * expfac * den))
     1                 .AND. ((sgnden * num) .LE. 0.0)
            ENDIF

            IF (okextr) THEN
              num = num / den
              x   = x1 + (delint * (1.0 - num))
            ELSE
              x   = outlim
            ENDIF

            IF (debug) THEN
              WRITE (6, *) 'x1 =', x1, 'x2 =', x2
              WRITE (6, *) 'f1 =', f1, 'f2 =', f2
              WRITE (6, *) 'df1 =', df1, 'df2 =', df2
              WRITE (6, *) 'x =', x
            ENDIF

            x1  = x2
            f1  = f2
            df1 = df2
            x2  = x
            CALL linfun (x2, f2, df2, n, wk1, wk2, wk3, wk4)
            icount = icount + 1

            IF (okextr .AND. (f2 .LE. f1)) THEN
              ifail  = 0
              contin = .FALSE.
            ENDIF

C           In any other case there must be a local minimum inside
C           interval (assuming continuity, differentiability).
C           Cubic interpolation is used to predict location, 'x', of
C           minimum.  If function at 'x' is less than (or equal to) 
C           function at either end of interval, 'x' is accepted as
C           estimate of true minimum.  Otherwise interval is changed
C           to that bounded by 'x' and the smaller of 'x1', 'x2'.
C
C           (The first conditional here tests for an ambiguous
C           exceptional case which would lead to a "divide by zero"
C           if cubic interpolation was applied).

          ELSEIF (contin .AND. (df1 .EQ. 0.0) .AND. (df2 .EQ. 0.0)
     1            .AND. (f2 .EQ. f1)) THEN
            ifail  = 0
            contin = .FALSE.

          ELSEIF (contin) THEN

            IF (f2 .LT. f1) THEN
              temp = x1
              x1   = x2
              x2   = temp
              temp = f1
              f1   = f2
              f2   = temp
              temp  = df1
              df1   = df2
              df2   = temp
              delint = - delint
              sgnint = - sgnint
            ENDIF

C           (Note delint <> 0 was checked earlier in routine).

            z = (3 * (f1 - f2) / delint) + df1 + df2

C           (In absence of rounding errors, argument of following square
C           root is +ve definite without ABS!).

            w2 = (z * z) - (df1 * df2)

            IF (w2 .LT. 0.0) THEN
              WRITE (6, *) '***linmin: Warning... w squared was less',
     1                     'than zero.'
            ENDIF

            w  = sgnint * SQRT (ABS (w2))
            temp = df2 + w - z
            den  = df2 - df1 + (2 * w)

C           (In absence of rounding errors, following error conditions
C           should be impossible!)

            IF (den .EQ. 0.0) THEN
              WRITE (6, *) '***linmin: Divide by zero in interpolation.'
              ifail  = 3
              contin = .FALSE.
            ENDIF

            IF (contin) THEN
              temp = temp / den
              x    = x1 + (delint * (1.0 - temp))

              IF (debug) THEN
                WRITE (6, *) 'x1 =', x1, 'x2 =', x2
                WRITE (6, *) 'f1 =', f1, 'f2 =', f2
                WRITE (6, *) 'df1 =', df1, 'df2 =', df2
                WRITE (6, *) 'x =', x
              ENDIF

              IF (((sgnint * (x - x1)) .LT. 0.0) .OR.
     1            ((sgnint * (x2 - x)) .LT. 0.0)) THEN
                WRITE (6, *) '***linmin: interpolation went outside',
     1                       'interval!'
                ifail  = 4
                contin = .FALSE.
              ENDIF

            ENDIF

            IF (contin) THEN
              x2 = x
              CALL linfun (x2, f2, df2, n, wk1, wk2, wk3, wk4)
              icount = icount + 1

              IF (f2 .LT. f1) THEN
                ifail  = 0
                contin = .FALSE.
              ENDIF

            ENDIF

          ENDIF

          IF (contin .AND. (icount .GE. ncount)) THEN
            WRITE (6, *) '***linmin: Failed to isolate minimum on line.'
            ifail  = 5
            contin = .FALSE.
          ENDIF

          GOTO 10
 20     CONTINUE

        IF (debug) THEN
          WRITE (6, *) 'exiting linmin'
        ENDIF

        RETURN
      END

      SUBROUTINE linzro (x1, f1, df1, x2, f2, df2, stplim, eps, ifail,
     1                   icount, ncount, n, wk1, wk2, wk3, wk4)
      
      IMPLICIT UNDEFINED
      REAL    x1, f1, df1, x2, df2, stplim, eps
      INTEGER ifail, icount, ncount
      INTEGER n
      REAL    wk1 (n), wk2 (n), wk3 (n), wk4 (n)

C     Locate minimum of a function along a line.
C     A routine for evaluating the function and its derivative
C     must be provided ('linfun').  The method is based on linear
C     interpolation, and starts from the given points 'x1', 'x2'.
C
C     Parameters:
C         x1, x2   On entry, initial intervale from which search 
C                  starts.  On exit 'x2' is a best estimate of
C                  minimum, 'x1' is a second best estimate.
C         f1, df1,
C         f2, df2  Corresponding values of function and derivative.
C                  On entry df1 > 0 is required.  'f1' and 'f2' are
C                  just carried along for the ride in this routine.
C         stplim   On entry, maximum stepsize to be taken in
C                  a single step.
C         eps      Routine terminates when ABS (x2 - x1) < eps.
C         ifail    On exit:
C                    ifail = 0 => Ended OK.
C                    ifail > 0 => Didn't end OK.
C         icount   Counter for number of 'linfun' calls.
C         ncount   Maximum number of 'linfun' calls.
C         n,
C         wk1,
C         ...,
C         wk4,       Parameters and workspace passed to
C                    routine 'linfun'.

        LOGICAL contin, gotmin, okextr, bisect

        LOGICAL debug
        COMMON /debug/ debug

        REAL outlim, inlim

C       Defined constants:
C           expfac     Determines factor by which interval may be 
C                      expanded in a single search step.
C           nintc      Maximum number of interval changes allowed
C                      in search.

        DATA expfac /5.0/
        DATA nintch /40/

        IF (debug) THEN
          WRITE (6, *) 'entering linzro'
        ENDIF

        ifail  = 0
        contin = .TRUE.

        IF ((x2 - x1) .EQ. 0.0) THEN
          WRITE (6, *) '***linzro: Interpolation interval is zero.'
          ifail  = 2
          contin = .FALSE.
        ELSEIF (df1 .GT. 0.0)
          WRITE (6, *) '***linzro: bad sign in linear minimisation.'
          ifail  = 2
          contin = .FALSE.
        ENDIF

        IF (contin .AND. (df1 .EQ. 0.0)) THEN
          temp = x1
          x1   = x2
          x2   = temp
          temp = f1
          f1   = f2
          f2   = temp
          temp  = df1
          df1   = df2
          df2   = temp
          contin = .FALSE.
        ENDIF

C       First isolate minimum by finding range over which
C       the derivative changes sign.

        IF (contin) THEN
          gotmin = df2 .GT. 0.0
        ENDIF

C       WHILE (contin .AND. (.NOT. gotmin))...

 10     IF ((.NOT. contin) .OR. gotmin) GOTO 20

C         The function is monotonically decreasing in the interval.
C         Linear extrapolation is attempted, but if this goes
C         outside a range determined by constant 'expfac', the 
C         interval is just moved onto the 'x2' side and expanded
C         by a factor 'expfac'.

          outlim = x2 + MIN (expfac * (x2 - x1), stplim)
          inlim  = x2 + MIN (eps, stplim)
          deldf  = df2 - df1
          okextr = deldf .NE. 0.0

          IF (okextr) THEN
            x = ((x1 * df2) - (x2 * df1)) / deldf
            okextr = (x .GT. x2) .AND. (x .LE. outlim)
          ENDIF

          IF (okextr .AND. (x .LT. inlim)) THEN
            x = inlim
          ELSEIF (.NOT. okextr) THEN
            x = outlim
          ENDIF

          IF (debug) THEN
            WRITE (6, *) 'x1 =', x1, 'x2 =', x2
            WRITE (6, *) 'f1 =', f1, 'f2 =', f2
            WRITE (6, *) 'df1 =', df1, 'df2 =', df2
            WRITE (6, *) 'x =', x
          ENDIF

          x1  = x2
          f1  = f2
          df1 = df2
          x2  = x
          CALL linfun (x2, f2, df2, n, wk1, wk2, wk3, wk4)
          icount = icount + 1
          gotmin = df2 .GE. 0.0

          IF (df2 .EQ. 0.0) THEN
            contin = .FALSE.
          ELSEIF ((.NOT. gotmin) .AND. (lintch .EQ. nintch)) THEN
            WRITE (6, *) '***linzro: Too many interval changes.'
            ifail  = 3
            contin = .FALSE.
          ENDIF

          GOTO 10
 20     CONTINUE

C       Now refine minimum by successive linear interpolation,
C       alternated with interval bisection.

        IF (contin) THEN
          gotmin = ABS (x2 - x1) .LT. eps
          sgn1   = -1.0
          sgn2   = 1.0
          bisect = .FALSE.
        ENDIF

C       WHILE (contin .AND. (.NOT. gotmin))...

 30     IF ((.NOT. contin) .OR. gotmin) GOTO 40

          IF (bisect) THEN
            x      = 0.5 * (x1 + x2)
            bisect = .FALSE.
          ELSE
            deldf  = df2 - df1
            x      = ((x1 * df2) - (x2 * df1)) / deldf
            bisect = .TRUE.
          ENDIF

          IF (debug) THEN
            WRITE (6, *) 'x1 =', x1, 'x2 =', x2
            WRITE (6, *) 'f1 =', f1, 'f2 =', f2
            WRITE (6, *) 'df1 =', df1, 'df2 =', df2
            WRITE (6, *) 'x =', x
          ENDIF

          CALL linfun (x, f, df, n, wk1, wk2, wk3, wk4)
          icount = icount + 1
          gotmin = df .EQ. 0.0
          sgn    = SIGN (1.0, df)

          IF (gotmin .OR. (sgn .eq. sgn1)) THEN
            x1   = x2
            f1   = f2
            df1  = df2
            sgn1 = sgn2
          ENDIF

          x2   = x
          f2   = f
          df2  = df
          sgn2 = sgn
          gotmin = gotmin .OR. (ABS (x2 - x1) .LT. eps)

          IF ((.NOT. gotmin) .AND. (icount .GE. ncount)) THEN
            WRITE (6, *) '***linzro: Too many interval changes.'
            ifail  = 3
            contin = .FALSE.
          ENDIF

          GOTO 30
 40     CONTINUE

        IF (debug) THEN
          WRITE (6, *) 'exiting linzro'
        ENDIF

        RETURN
      END

      SUBROUTINE linfun (lambda, f, df, n, x0, s, x, g)

      IMPLICIT UNDEFINED
      INTEGER n
      REAL    x0 (n), s (n), x (n), g (n)
      REAL    lambda

C     Evaluate function along line x = x0 + (lambda * s), where 'x0' is
C     initial point, 's' is some vector, and 'lambda' is parameter.
C     Also derivative, 'df', of 'f' with respect to parameter is
C     returned.  'g' returns the full gradient.

        DO 10 i = 1, n
          x (i) = x0 (i) + (lambda * s (i))
 10     CONTINUE

        CALL func (n, x, f, g)

        df = 0
        DO 20 i = 1, n
          df = df + (s (i) * g (i))
 20     CONTINUE

        RETURN
      END

      SUBROUTINE func (n, v, p, g)

      IMPLICIT UNDEFINED
      INTEGER n
      REAL    v (n), g (n)
      REAL    p

C       Chi-squared and derivatives.

        COMMON /data/ x (40), y (40), s (40), npnts

        REAL    df (10)
        LOGICAL debug
        COMMON /debug/ debug

        DATA ldf /10/

        IF (n .GT. ldf) THEN
          WRITE (6, *) '***func: Array needs redimensioning.'
          STOP
        ENDIF

        IF (debug) THEN
          WRITE (6, *) 'Entering func.'
          WRITE (6, *) 'v =', v
        ENDIF

        p = 0.0
        DO 10 j = 1, n
          g (j) = 0.0
 10     CONTINUE

        DO 30 i = 1, npts
          CALL model (x (i), v, n, thy, df)
          IF (s (i) .EQ. 0.0) THEN
            WRITE (6, *) '***func: Standard deviation = 0.  i =', i
            STOP
          ENDIF
          dev = (thy - y (i)) / s (i)
          p   = p + dev**2
          DO 20 j = 1, n
            g (j) = g (j) + ((2 * dev * df (j)) / s (i))
 20       CONTINUE
 30     CONTINUE
        
        IF (debug) THEN
          WRITE (6, *) 'p =', p
          WRITE (6, *) 'g:', g
          WRITE (6, *) 'Exiting func'
        ENDIF

        RETURN
      END

      SUBROUTINE model (t, v, n, y, dy)

      IMPLICIT UNDEFINED
      INTEGER n
      REAL    v (n), dy (n)
      REAL    t, y

C     Curve to be fitted, and derivatives.

        b     = v (1)
        tcrit = v (2)
  
        IF (t .GT. tcrit) THEN
          x = t - tcrit
          y = x / b**2
          dy (1) = - (2 * x) / b**3
          dy (2) = - 1.0 / b**2
        ELSE
          y = 0.0
          dy (1) = 0.0
          dy (2) = 0.0
        ENDIF
  
        RETURN
      END
