      FUNCTION QSUBA(F,A,B,EPSIL,MCOUNT,RELERR,IND)
C   THIS FUNCTION ROUTINE PERFORMS AUTOMATIC INTEGRATION
C OVER A FINITE INTERVAL USING THE BASIC INTEGRATION
C ALGORITHM QUAD TOGETHER WITH,  IF NECESSARY AN ADAPTIVE
C SUBDIVISION PROCESS.  IT IS GENERALLY MORE EFFICIENT THAN
C THE NON-ADAPTIVE ALGORITHM QSUB BUT IS LIKELY TO BE LESS
C RELIABLE(SEE COMP.J.,14,189,1971).
C   THE CALL TAKES THE FORM
C           QSUBA(F,A,B,EPSIL,MCOUNT,RELERR,IND)
C AND CAUSES F(X) TO BE INTEGRATED OVER (A,B) WITH RELATIVE
C ERROR HOPEFULLY NOT EXCEEDING EPSIL.  SHOULD QUAD CONVERGE
C (ICHECK=0) THEN QSUBA WILL RETURN THE VALUE OBTAINED BY IT
C OTHERWISE SUBDIVISION WILL BE INVOKED AS A RESCUE
C OPERATION IN AN ADAPTIVE MANNER. THE ARGUMENT RELERR GIVES
C A CRUDE ESTIMATE OF THE ACTUAL RELATIVE ERROR OBTAINED.
C THE ARGUMENT MCOUNT IS THE MAXIMUM NUMBER OF POINTS AT WHICH
C THE INTEGRAND MAY BE EVALUATED.
C   THE SUBDIVISION STRATEGY IS AS FOLLOWS
C AT EACH STAGE OF THE PROCESS AN INTERVAL IS PRESENTED FOR
C SUBDIVISION (INITIALLY THIS WILL BE THE WHOLE INTERVAL
C (A,B)).  THE INTERVAL IS HALVED AND QUAD APPLIED TO EACH
C SUBINTERVAL.  SHOULD QUAD FAIL ON THE FIRST SUBINTERVAL
C THE SUBINTERVAL IS STACKED FOR FUTURE SUBDIVISION AND THE
C SECOND SUBINTERVAL IMMEDIATELY EXAMINED.  SHOULD QUAD FAIL
C ON THE SECOND SUBINTERVAL THE SUBINTERVAL IS
C IMMEDIATELY SUBDIVIDED AND THE WHOLE PROCESS REPEATED.
C EACH TIME A CONVERGED RESULT IS OBTAINED IT IS
C ACCUMULATED AS THE PARTIAL VALUE OF THE INTEGRAL. WHEN
C QUAD CONVERGES ON BOTH SUBINTERVALS THE INTERVAL LAST
C STACKED IS CHOSEN NEXT FOR SUBDIVISION AND THE PROCESS
C REPEATED.  A SUBINTERVAL IS NOT EXAMINED AGAIN ONCE A
C CONVERGED RESULT IS OBTAINED FOR IT SO THAT A SPURIOUS
C CONVERGENCE IS MORE LIKELY TO SLIP THROUGH THAN FOR THE
C NON-ADAPTIVE ALGORITHM QSUB.
C   THE VARIABLE IND IS SET BY QSUBA. IND MAY HAVE THE VALUES ...
C   IND=0   QSUBA IS SATISFIED THAT THE INTEGRAL HAS BEEN
C           COMPUTED TO THE DESIRED ACCURACY.
C   IND=1   THE INTEGRAL HAS BEEN COMPUTED, BUT QSUBA IS
C           NOT CERTAIN OF THE ACCURACY OF THE RESULT.
C   IND=2   THE INTEGRAND HAS BEEN EVALUATED AT MCOUNT
C           POINTS. IF MORE EVALUATIONS ARE NEEDED THEN
C           QSUBA TERMINATES.
C   IND=3   QSUBA CANNOT COMPUTE THE INTEGRAL TO THE DESIRED
C           ACCURACY. IND IS SET TO 3 WHENEVER THE STACK OF
C           INTERVALS BECOMES FULL (IT CURRENTLY CAN HOLD 50
C           INTERVALS). A RESULT IS OBTAINED BY CONTINUING
C           THE INTEGRATION IGNORING CONVERGENCE FAILURES
C           WHICH CANNOT BE ACCOMMODATED ON THE STACK.
C  THE RELIABILITY OF THE ALGORITHM WILL DECREASE FOR LARGE
C  VALUES OF EPSIL.  IT IS RECOMMENDED THAT EPSIL SHOULD
C  GENERALLY BE LESS THAN ABOUT 0.001.
      DIMENSION RESULT(8), STACK(100)
      EXTERNAL F
      DATA ISMAX/100/
      IND=0
      RELERR = 0.0
      CALL QUAD(A, B, RESULT, K, EPSIL, NPTS, ICHECK, F)
      QSUBA = RESULT(K)
      IF (QSUBA.NE.0.0) RELERR =
     * ABS((RESULT(K)-RESULT(K-1))/QSUBA)
C CHECK IF SUBDIVISION IS NEEDED
      IF (ICHECK.EQ.0) RETURN
C SUBDIVIDE
      RELERR = 0.0
      QSUBA = 0.0
      IS = 1
      SUB1 = A
      SUB3 = B
   10 IF (NPTS.GE.MCOUNT) GO TO 60
      SUB2 = (SUB1+SUB3)*0.5
      CALL QUAD(SUB1, SUB2, RESULT, K, EPSIL, NF, ICHECK, F)
      NPTS = NPTS + NF
      COMP = ABS(RESULT(K)-RESULT(K-1))
      SUM = QSUBA + RESULT(K)
      IF ((ICHECK.EQ.0).OR.(COMP.LE.ABS(EPSIL*SUM))) GO TO 30
      IF (IS.GE.ISMAX) GO TO 20
C STACK SUBINTERVAL (SUB1,SUB2) FOR FUTURE EXAMINATION
      STACK(IS) = SUB1
      IS = IS + 1
      STACK(IS) = SUB2
      IS = IS + 1
      GO TO 40
C UPDATE QSUBA AND CHECK IF ANY SIGNIFICANT DIGITS ARE LOST
   20 IND=3
   30 X=QSUBA
      QSUBA=SUM
      RELERR = RELERR + COMP
      IF ((IND.NE.0).OR.(X*RESULT(K).GE.0.)) GO TO 40
      X=.1*AMIN1(ABS(X),ABS(RESULT(K)))
      IF (ABS(QSUBA).LE.X) IND=1
C CALL QUAD FOR THE INTERVAL (SUB2,SUB3)
   40 IF (NPTS.GE.MCOUNT) GO TO 60
      CALL QUAD(SUB2, SUB3, RESULT, K, EPSIL, NF, ICHECK, F)
      NPTS = NPTS + NF
      COMP = ABS(RESULT(K)-RESULT(K-1))
      SUM = QSUBA + RESULT(K)
      IF ((ICHECK.EQ.0).OR.(COMP.LE.ABS(EPSIL*SUM))) GO TO 50
C SUBDIVIDE INTERVAL (SUB2,SUB3)
      SUB1 = SUB2
      GO TO 10
C UPDATE QSUBA AND CHECK IF ANY SIGNIFICANT DIGITS ARE LOST
   50 X=QSUBA
      QSUBA=SUM
      RELERR = RELERR + COMP
      IF ((IND.NE.0).OR.(X*RESULT(K).GE.0.)) GO TO 51
      X=.1*AMIN1(ABS(X),ABS(RESULT(K)))
      IF (ABS(QSUBA).LE.X) IND=1
   51 IF (IS.EQ.1) GO TO 61
C SUBDIVIDE THE DELINQUENT INTERVAL LAST STACKED
      IS = IS - 1
      SUB3 = STACK(IS)
      IS = IS - 1
      SUB1 = STACK(IS)
      GO TO 10
C SUBDIVISION RESULT
   60 IND=2
   61 IF (QSUBA.NE.0.) RELERR=RELERR/ABS(QSUBA)
      RETURN
      END
