      SUBROUTINE DQAG(F,A,B,EPSABS,EPSREL,KEY,RESULT,ABSERR,NEVAL,IER,
     *    LIMIT,LENW,LAST,IWORK,WORK)
C***BEGIN PROLOGUE  DQAG
C***DATE WRITTEN   800101   (YYMMDD)
C***REVISION DATE  830518   (YYMMDD)
C***CATEGORY NO.  H2A1A1
C***KEYWORDS  AUTOMATIC INTEGRATOR, GENERAL-PURPOSE,
C             INTEGRAND EXAMINATOR, GLOBALLY ADAPTIVE,
C             GAUSS-KRONROD
C***AUTHOR  PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV - K.U.LEUVEN
C           DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
C***PURPOSE  THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
C            DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B),
C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
C            ABS(I-RESULT)LE.MAX(EPSABS,EPSREL*ABS(I)).
C***DESCRIPTION
C
C        COMPUTATION OF A DEFINITE INTEGRAL
C        STANDARD FORTRAN SUBROUTINE
C        DOUBLE PRECISION VERSION
C
C            F      - DOUBLE PRECISION
C                     FUNCTION SUBPROGAM DEFINING THE INTEGRAND
C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C            A      - DOUBLE PRECISION
C                     LOWER LIMIT OF INTEGRATION
C
C            B      - DOUBLE PRECISION
C                     UPPER LIMIT OF INTEGRATION
C
C            EPSABS - DOUBLE PRECISION
C                     ABSOLUTE ACCORACY REQUESTED
C            EPSREL - DOUBLE PRECISION
C                     RELATIVE ACCURACY REQUESTED
C                     IF  EPSABS.LE.0
C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
C                     THE ROUTINE WILL END WITH IER = 6.
C
C            KEY    - INTEGER
C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
C                     A GAUSS-KRONROD PAIR IS USED WITH
C                       7 - 15 POINTS IF KEY.LT.2,
C                      10 - 21 POINTS IF KEY = 2,
C                      15 - 31 POINTS IF KEY = 3,
C                      20 - 41 POINTS IF KEY = 4,
C                      25 - 51 POINTS IF KEY = 5,
C                      30 - 61 POINTS IF KEY.GT.5.
C
C         ON RETURN
C            RESULT - DOUBLE PRECISION
C                     APPROXIMATION TO THE INTEGRAL
C
C            ABSERR - DOUBLE PRECISION
C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
C
C            NEVAL  - INTEGER
C                     NUMBER OF INTEGRAND EVALUATIONS
C
C            IER    - INTEGER
C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
C                             ACCURACY HAS BEEN ACHIEVED.
C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
C                             THE ESTIMATES FOR RESULT AND ERROR ARE
C                             LESS RELIABLE. IT IS ASSUMED THAT THE
C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
C                      ERROR MESSAGES
C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE
C                             SUBDIVISIONS BY INCREASING THE VALUE OF
C                             LIMIT (AND TAKING THE ACCORDING DIMENSION
C                             ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF
C                             THIS YIELD NO IMPROVEMENT IT IS ADVISED
C                             TO ANALYZE THE INTEGRAND IN ORDER TO
C                             DETERMINE THE INTEGRATION DIFFICULATIES.
C                             IF THE POSITION OF A LOCAL DIFFICULTY CAN
C                             BE DETERMINED (I.E.SINGULARITY,
C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
C                             INTERVAL AT THIS POINT AND CALLING THE
C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
C                             SHOULD BE USED WHICH IS DESIGNED FOR
C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS
C                             DETECTED, WHICH PREVENTS THE REQUESTED
C                             TOLERANCE FROM BEING ACHIEVED.
C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS
C                             AT SOME POINTS OF THE INTEGRATION
C                             INTERVAL.
C                         = 6 THE INPUT IS INVALID, BECAUSE
C                             (EPSABS.LE.0 AND
C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
C                             OR LIMIT.LT.1 OR LENW.LT.LIMIT*4.
C                             RESULT, ABSERR, NEVAL, LAST ARE SET
C                             TO ZERO.
C                             EXCEPT WHEN LENW IS INVALID, IWORK(1),
C                             WORK(LIMIT*2+1) AND WORK(LIMIT*3+1) ARE
C                             SET TO ZERO, WORK(1) IS SET TO A AND
C                             WORK(LIMIT+1) TO B.
C
C         DIMENSIONING PARAMETERS
C            LIMIT - INTEGER
C                    DIMENSIONING PARAMETER FOR IWORK
C                    LIMIT DETERMINES THE MAXIMUM NUMBER OF SUBINTERVALS
C                    IN THE PARTITION OF THE GIVEN INTEGRATION INTERVAL
C                    (A,B), LIMIT.GE.1.
C                    IF LIMIT.LT.1, THE ROUTINE WILL END WITH IER = 6.
C
C            LENW  - INTEGER
C                    DIMENSIONING PARAMETER FOR WORK
C                    LENW MUST BE AT LEAST LIMIT*4.
C                    IF LENW.LT.LIMIT*4, THE ROUTINE WILL END WITH
C                    IER = 6.
C
C            LAST  - INTEGER
C                    ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS
C                    PRODUCED IN THE SUBDIVIOSION PROCESS, WHICH
C                    DETERMINES THE NUMBER OF SIGNIFICANT ELEMENTS
C                    ACTUALLY IN THE WORK ARRAYS.
C
C         WORK ARRAYS
C            IWORK - INTEGER
C                    VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K
C                    ELEMENTS OF WHICH CONTAIN POINTERS TO THE ERROR
C                    ESTIMATES OVER THE SUBINTERVALS, SUCH THAT
C                    WORK(LIMIT*3+IWORK(1)),... , WORK(LIMIT*3+IWORK(K))
C                    FORM A DECREASING SEQUENCE WITH K = LAST IF
C                    LAST.LE.(LIMIT/2+2), AND K = LIMIT+1-LAST OTHERWISE
C
C            WORK  - DOUBLE PRECISION
C                    VECTOR OF DIMENSION AT LEAST LENW
C                    ON RETURN
C                    WORK(1), ..., WORK(LAST) CONTAIN THE LEFT END
C                    POINTS OF THE SUBINTERVALS IN THE PARTITION OF
C                     (A,B),
C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN THE
C                     RIGHT END POINTS,
C                    WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) CONTAIN
C                     THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS,
C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST) CONTAIN
C                     THE ERROR ESTIMATES.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DQAGE,XERROR
C***END PROLOGUE  DQAG
      DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK
      INTEGER IER,IWORK,KEY,LAST,LENW,LIMIT,LVL,L1,L2,L3,NEVAL
C
      DIMENSION IWORK(LIMIT),WORK(LENW)
C
      EXTERNAL F
C
C         CHECK VALIDITY OF LENW.
C
C***FIRST EXECUTABLE STATEMENT  DQAG
      IER = 6
      NEVAL = 0
      LAST = 0
      RESULT = 0.0D+00
      ABSERR = 0.0D+00
      IF(LIMIT.LT.1.OR.LENW.LT.LIMIT*4) GO TO 10
C
C         PREPARE CALL FOR DQAGE.
C
      L1 = LIMIT+1
      L2 = LIMIT+L1
      L3 = LIMIT+L2
C
      CALL DQAGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR,NEVAL,
     *  IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST)
C
C         CALL ERROR HANDLER IF NECESSARY.
C
      LVL = 0
10    IF(IER.EQ.6) LVL = 1
      IF(IER.NE.0) CALL XERROR(26HABNORMAL RETURN FROM DQAG ,26,IER,LVL)
      RETURN
      END
