      SUBROUTINE DQK21(F,A,B,RESULT,ABSERR,RESABS,RESASC,ISIG)
C     ..................................................................
C
C 1.        DQK21
C           INTEGRATION RULES
C              STANDARD FORTRAN SUBROUTINE
C              DOUBLE PRECISION VERSION
C
C 2.        PURPOSE
C              TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR
C                             ESTIMATE
C                         J = INTEGRAL OF ABS(F) OVER (A,B)
C
C 3.        CALLING SEQUENCE
C              CALL DQK21(F,A,B,RESULT,ABSERR,RESABS,RESASC,ISIG)
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
C                       TO BE DECLARED E X T E R N A L IN THE
C                       CALLING PROGRAM.
C
C              A      - DOUBLE PRECISION
C                       LOWER LIMIT OF INTEGRATION
C
C              B      - DOUBLE PRECISION
C                       UPPER LIMIT OF INTEGRATION
C
C            ON RETURN
C              RESULT - DOUBLE PRECISION
C                       APPROXIMATION TO THE INTEGRAL I
C                       RESULT IS COMPUTED BY APPLYING THE 21-POINT
C                       KRONROD RULE (RESK) OBTAINED BY OPTIMAL
C                       ADDITION OF ABSCISSAE TO THE 10-POINT GAUSS
C                       RULE (RESG).
C
C              ABSERR - DOUBLE PRECISION
C                       ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
C                       WHICH SHOULD NOT EXCEED ABS(I-RESULT)
C
C              RESABS - DOUBLE PRECISION
C                       APPROXIMATION TO THE INTEGRAL J
C
C              RESASC - DOUBLE PRECISION
C                       APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A))
C                       OVER (A,B)
C
C              ISIG   - INTEGER
C                       ISIG=0  THE INTEGRAL WAS APPROXIMATED.
C                       ISIG=5  THE INTERVAL (A,B) IS TOO SHORT.
C                               THE INTEGRAL CANNOT BE COMPUTED.
C
C 4.        SUBROUTINES OR FUNCTIONS NEEDED
C                 - F (USER-PROVIDED FUNCTION)
C                 - DPMPAR
C                 - FORTRAN DABS, DMAX1, DMIN1
C
C     ..................................................................
C
      DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DHLGTH,DPMPAR,EPMACH,F,FC,
     *  FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESG,RESK,RESKH,RESULT,
     *  UFLOW,WG,WGK,XGK,RESASC
      INTEGER ISIG,J,JTW,JTWM1
      EXTERNAL F
C
      DIMENSION FV1(10),FV2(10),WG(5),WGK(11),XGK(11)
C
C           THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1).
C           BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
C           CORRESPONDING WEIGHTS ARE GIVEN.
C
C           XGK    - ABSCISSAE OF THE 21-POINT KRONROD RULE
C                    XGK(2), XGK(4), ...  ABSCISSAE OF THE 10-POINT
C                    GAUSS RULE
C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
C                    ADDED TO THE 10-POINT GAUSS RULE
C
C           WGK    - WEIGHTS OF THE 21-POINT KRONROD RULE
C
C           WG     - WEIGHTS OF THE 10-POINT GAUSS RULE
C
      DATA XGK(1),XGK(2),XGK(3),XGK(4),XGK(5),XGK(6),XGK(7),
     *  XGK(8),XGK(9),XGK(10),XGK(11)/
     *     0.995657163025808080735527280689D+00,
     *     0.973906528517171720077964012084D+00,
     *     0.930157491355708226001207180060D+00,
     *     0.865063366688984510732096688423D+00,
     *     0.780817726586416897063717578345D+00,
     *     0.679409568299024406234327365115D+00,
     *     0.562757134668604683339000099273D+00,
     *     0.433395394129247190799265943166D+00,
     *     0.294392862701460198131126603104D+00,
     *     0.148874338981631210884826001130D+00,
     *     0.000000000000000000000000000000D+00/
C
      DATA WGK(1),WGK(2),WGK(3),WGK(4),WGK(5),WGK(6),WGK(7),
     *  WGK(8),WGK(9),WGK(10),WGK(11)/
     *     0.116946388673718742780643960622D-01,
     *     0.325581623079647274788189724594D-01,
     *     0.547558965743519960313813002446D-01,
     *     0.750396748109199527670431409162D-01,
     *     0.931254545836976055350654650834D-01,
     *     0.109387158802297641899210590326D+00,
     *     0.123491976262065851077958109831D+00,
     *     0.134709217311473325928054001772D+00,
     *     0.142775938577060080797094273139D+00,
     *     0.147739104901338491374841515972D+00,
     *     0.149445554002916905664936468390D+00/
C
      DATA WG(1),WG(2),WG(3),WG(4),WG(5)/
     *     0.666713443086881375935688098933D-01,
     *     0.149451349150580593145776339658D+00,
     *     0.219086362515982043995534934228D+00,
     *     0.269266719309996355091226921569D+00,
     *     0.295524224714752870173892994651D+00/
C
C
C           LIST OF MAJOR VARIABLES
C           -----------------------
C
C           CENTR  - MID POINT OF THE INTERVAL
C           HLGTH  - HALF-LENGTH OF THE INTERVAL
C           ABSC   - ABSCISSA
C           FVAL*  - FUNCTION VALUE
C           RESG   - RESULT OF THE 10-POINT GAUSS FORMULA
C           RESK   - RESULT OF THE 21-POINT KRONROD FORMULA
C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B),
C                    I.E. TO I/(B-A)
C
C
C           MACHINE DEPENDENT CONSTANTS
C           ---------------------------
C
C           EPMACH IS THE LARGEST RELATIVE SPACING.
C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
C
      EPMACH = DPMPAR(1)
      UFLOW = DPMPAR(2)
C
      CENTR = 0.5D+00*(A+B)
      HLGTH = 0.5D+00*(B-A)
      DHLGTH = DABS(HLGTH)
C
C             CHECK IF THE INTERVAL (A,B) IS TOO SHORT
C
      ISIG = 5
      ABSC = DABS(CENTR) + DHLGTH*0.14D+00
      IF (ABSC .EQ. DABS(CENTR)) RETURN
C
C           COMPUTE THE 21-POINT KRONROD APPROXIMATION TO
C           THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
C
      ISIG = 0
      RESG = 0.0D+00
      FC = F(CENTR)
      RESK = WGK(11)*FC
      RESABS = DABS(RESK)
      DO 10 J=1,5
        JTW = 2*J
        ABSC = HLGTH*XGK(JTW)
        FVAL1 = F(CENTR-ABSC)
        FVAL2 = F(CENTR+ABSC)
        FV1(JTW) = FVAL1
        FV2(JTW) = FVAL2
        FSUM = FVAL1+FVAL2
        RESG = RESG+WG(J)*FSUM
        RESK = RESK+WGK(JTW)*FSUM
        RESABS = RESABS+WGK(JTW)*(DABS(FVAL1)+DABS(FVAL2))
   10 CONTINUE
      DO 15 J = 1,5
        JTWM1 = 2*J-1
        ABSC = HLGTH*XGK(JTWM1)
        FVAL1 = F(CENTR-ABSC)
        FVAL2 = F(CENTR+ABSC)
        FV1(JTWM1) = FVAL1
        FV2(JTWM1) = FVAL2
        FSUM = FVAL1+FVAL2
        RESK = RESK+WGK(JTWM1)*FSUM
        RESABS = RESABS+WGK(JTWM1)*(DABS(FVAL1)+DABS(FVAL2))
   15 CONTINUE
      RESKH = RESK*0.5D+00
      RESASC = WGK(11)*DABS(FC-RESKH)
      DO 20 J=1,10
        RESASC = RESASC+WGK(J)*(DABS(FV1(J)-RESKH)+DABS(FV2(J)-RESKH))
   20 CONTINUE
      RESULT = RESK*HLGTH
      RESABS = RESABS*DHLGTH
      RESASC = RESASC*DHLGTH
      ABSERR = DABS((RESK-RESG)*HLGTH)
      IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00)
     *  ABSERR = RESASC*DMIN1(0.1D+01,
     *  (0.2D+03*ABSERR/RESASC)**1.5D+00)
      IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = DMAX1
     *  ((EPMACH*0.5D+02)*RESABS,ABSERR)
      RETURN
      END
