      SUBROUTINE DQAGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,
     *  NEVAL,IER,LIMIT,LENW,LAST,IWORK,WORK)
C ......................................................................
C
C 1.     DQAGS
C        COMPUTATION OF A DEFINITE INTEGRAL
C           STANDARD FORTRAN SUBROUTINE
C           DOUBLE PRECISION VERSION
C
C 2.     PURPOSE
C           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
C 3.     CALLING SEQUENCE
C           CALL DQAGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,
C                      NEVAL,IER,LIMIT,LENW,LAST,IWORK,WORK)
C
C        PARAMETERS
C         ON ENTRY
C            F      - DOUBLE PRECISION
C                     FUNCTION SUBPROGRAM 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 ACCURACY 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         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 INTEGRAL AND ERROR ARE
C                             LESS RELIABLE. IT IS ASSUMED THAT THE
C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
C
C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB-
C                             DIVISIONS BY INCREASING THE VALUE OF
C                             LIMIT (AND TAKING THE ACCORDING
C                             DIMENSION ADJUSTMENTS INTO ACCOUNT).
C                             HOWEVER, IF THIS YIELDS NO IMPROVEMENT
C                             IT IS ADVISED TO ANALYZE THE INTEGRAND
C                             IN ORDER TO DETERMINE THE INTEGRATION
C                             DIFFICULTIES. IF THE POSITION OF A
C                             LOCAL DIFFICULTY CAN BE DETERMINED (E.G.
C                             SINGULARITY, DISCONTINUITY WITHIN THE
C                             INTERVAL) ONE WILL PROBABLY GAIN FROM
C                             SPLITTING UP THE INTERVAL AT THIS POINT
C                             AND CALLING THE INTEGRATOR ON THE SUB-
C                             RANGES. IF POSSIBLE, AN APPROPRIATE
C                             SPECIAL-PURPOSE INTEGRATOR SHOULD BE USED,
C                             WHICH IS DESIGNED FOR HANDLING THE TYPE
C                             OF DIFFICULTY INVOLVED.
C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETEC-
C                             TED, WHICH PREVENTS THE REQUESTED
C                             TOLERANCE FROM BEING ACHIEVED.
C                             THE ERROR MAY BE UNDER-ESTIMATED.
C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR
C                             OCCURS AT SOME  POINTS OF THE INTEGRATION
C                             INTERVAL.
C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
C                             IT IS PRESUMED THAT THE REQUESTED
C                             TOLERANCE CANNOT BE ACHIEVED, AND THAT THE
C                             RETURNED RESULT IS THE BEST WHICH CAN BE
C                             OBTAINED.
C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
C                             SLOWLY CONVERGENT. IT MUST BE NOTED
C                             THAT DIVERGENCE CAN OCCUR WITH ANY OTHER
C                             VALUE OF IER.
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
C                             SET TO ZERO. EXEPT WHEN LIMIT OR
C                             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
C                    OF SUBINTERVALS IN THE PARTITION OF THE GIVEN
C                    INTEGRATION INTERVAL (A,B), LIMIT.GE.1.
C                    IF LIMIT.LT.1, THE ROUTINE WILL END WITH
C                    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
C                    WITH IER = 6.
C
C            LAST  - INTEGER
C                    ON RETURN, LAST EQUALS THE NUMBER OF
C                    SUBINTERVALS PRODUCED IN THE SUBDIVISION
C                    PROCESS, WHICH DETERMINES THE NUMBER OF
C                    SIGNIFICANT ELEMENTS ACTUALLY IN THE WORK
C                    ARRAYS.
C
C         WORK ARRAYS
C            IWORK - INTEGER
C                    VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K
C                    ELEMENTS OF WHICH CONTAIN POINTERS
C                    TO THE ERROR ESTIMATES OVER THE SUBINTERVALS
C                    SUCH THAT WORK(LIMIT*3+IWORK(1)),... ,
C                    WORK(LIMIT*3+IWORK(K)) FORM A DECREASING
C                    SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2), AND
C                    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
C                     END POINTS OF THE SUBINTERVALS IN THE
C                     PARTITION OF (A,B),
C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN
C                     THE RIGHT END POINTS,
C                    WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST)
C                     CONTAIN THE INTEGRAL APPROXIMATIONS OVER
C                     THE SUBINTERVALS,
C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
C                     CONTAIN THE ERROR ESTIMATES.
C
C
C 4.      SUBROUTINES OR FUNCTIONS NEEDED
C              - DQAGSE
C              - DQK21
C              - DQPSRT
C              - DQELG
C              - F (USER-PROVIDED FUNCTION)
C              - DPMPAR
C              - FORTRAN DABS, DMAX1, DMIN1
C
C ......................................................................
C
      DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK
      INTEGER IER,IWORK,LENW,LIMIT,L1,L2,L3,NEVAL
C
      DIMENSION IWORK(LIMIT),WORK(LENW)
C
      EXTERNAL F
C
C         CHECK VALIDITY OF LIMIT AND LENW.
C
      IER = 6
      NEVAL = 0
      LAST = 0
      RESULT = 0.0D+00
      ABSERR = 0.0D+00
      IF(LIMIT.LT.1.OR.LENW.LT.LIMIT*4) RETURN
C
C         PREPARE CALL FOR DQAGSE
C
      L1 = LIMIT+1
      L2 = LIMIT+L1
      L3 = LIMIT+L2
C
      CALL DQAGSE(F,A,B,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,NEVAL,
     *  IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST)
      RETURN
      END
