dqk15w Subroutine

    subroutine dqk15w(f, w, p1, p2, p3, p4, Kp, a, b, Result, Abserr, Resabs, &
                      Resasc)
        implicit none

        procedure(func) :: f !! function subprogram defining the integrand function `f(x)`.
        procedure(weight_func) :: w !! function subprogram defining the integrand weight function `w(x)`.
        real(wp), intent(in) :: p1 !! parameter in the weight function
        real(wp), intent(in) :: p2 !! parameter in the weight function
        real(wp), intent(in) :: p3 !! parameter in the weight function
        real(wp), intent(in) :: p4 !! parameter in the weight function
        integer, intent(in) :: Kp !! key for indicating the type of weight function
        real(wp), intent(in) :: a !! lower limit of integration
        real(wp), intent(in) :: b !! upper limit of integration
        real(wp), intent(out) :: Result !! approximation to the integral i
                                        !! `result` is computed by applying the 15-point
                                        !! kronrod rule (resk) obtained by optimal addition
                                        !! of abscissae to the 7-point gauss rule (resg).
        real(wp), intent(out) :: Abserr !! estimate of the modulus of the absolute error,
                                        !! which should equal or exceed `abs(i-result)`
        real(wp), intent(out) :: Resabs !! approximation to the integral of `abs(f)`
        real(wp), intent(out) :: Resasc !! approximation to the integral of `abs(f-i/(b-a))`

        real(wp) :: absc1, absc2, dhlgth, fc, fsum, fv1(7), fv2(7)
        integer :: j, jtw, jtwm1
        real(wp) :: centr !! mid point of the interval
        real(wp) :: hlgth !! half-length of the interval
        real(wp) :: absc !! abscissa
        real(wp) :: fval1 !! function value
        real(wp) :: fval2 !! function value
        real(wp) :: resg !! result of the 7-point gauss formula
        real(wp) :: resk !! result of the 15-point kronrod formula
        real(wp) :: reskh !! approximation to the mean value of f*w over `(a,b)`, i.e. to `i/(b-a)`

        ! the abscissae and weights are given for the interval (-1,1).
        ! because of symmetry only the positive abscissae and their
        ! corresponding weights are given.

        real(wp), dimension(8), parameter :: xgk = [ &
                                             9.91455371120812639206854697526328516642e-1_wp, &
                                             9.49107912342758524526189684047851262401e-1_wp, &
                                             8.64864423359769072789712788640926201211e-1_wp, &
                                             7.41531185599394439863864773280788407074e-1_wp, &
                                             5.86087235467691130294144838258729598437e-1_wp, &
                                             4.05845151377397166906606412076961463347e-1_wp, &
                                             2.07784955007898467600689403773244913480e-1_wp, &
                                             0.00000000000000000000000000000000000000e0_wp] !! abscissae of the 15-point kronrod rule:
                                                                                            !!
                                                                                            !! * xgk(2), xgk(4), ... abscissae of the 7-point
                                                                                            !!   gauss rule
                                                                                            !! * xgk(1), xgk(3), ... abscissae which are optimally
                                                                                            !!   added to the 7-point gauss rule

        real(wp), dimension(8), parameter :: wgk = [ &
                                             2.29353220105292249637320080589695919936e-2_wp, &
                                             6.30920926299785532907006631892042866651e-2_wp, &
                                             1.04790010322250183839876322541518017444e-1_wp, &
                                             1.40653259715525918745189590510237920400e-1_wp, &
                                             1.69004726639267902826583426598550284106e-1_wp, &
                                             1.90350578064785409913256402421013682826e-1_wp, &
                                             2.04432940075298892414161999234649084717e-1_wp, &
                                             2.09482141084727828012999174891714263698e-1_wp] !! weights of the 15-point kronrod rule

        real(wp), dimension(4), parameter :: wg = [ &
                                             1.29484966168869693270611432679082018329e-1_wp, &
                                             2.79705391489276667901467771423779582487e-1_wp, &
                                             3.81830050505118944950369775488975133878e-1_wp, &
                                             4.17959183673469387755102040816326530612e-1_wp] !! weights of the 7-point gauss rule

        centr = 0.5_wp*(a + b)
        hlgth = 0.5_wp*(b - a)
        dhlgth = abs(hlgth)

        ! compute the 15-point kronrod approximation to the
        ! integral, and estimate the error.

        fc = f(centr)*w(centr, p1, p2, p3, p4, Kp)
        resg = wg(4)*fc
        resk = wgk(8)*fc
        Resabs = abs(resk)
        do j = 1, 3
            jtw = j*2
            absc = hlgth*xgk(jtw)
            absc1 = centr - absc
            absc2 = centr + absc
            fval1 = f(absc1)*w(absc1, p1, p2, p3, p4, Kp)
            fval2 = f(absc2)*w(absc2, p1, p2, p3, p4, Kp)
            fv1(jtw) = fval1
            fv2(jtw) = fval2
            fsum = fval1 + fval2
            resg = resg + wg(j)*fsum
            resk = resk + wgk(jtw)*fsum
            Resabs = Resabs + wgk(jtw)*(abs(fval1) + abs(fval2))
        end do
        do j = 1, 4
            jtwm1 = j*2 - 1
            absc = hlgth*xgk(jtwm1)
            absc1 = centr - absc
            absc2 = centr + absc
            fval1 = f(absc1)*w(absc1, p1, p2, p3, p4, Kp)
            fval2 = f(absc2)*w(absc2, p1, p2, p3, p4, Kp)
            fv1(jtwm1) = fval1
            fv2(jtwm1) = fval2
            fsum = fval1 + fval2
            resk = resk + wgk(jtwm1)*fsum
            Resabs = Resabs + wgk(jtwm1)*(abs(fval1) + abs(fval2))
        end do
        reskh = resk*0.5_wp
        Resasc = wgk(8)*abs(fc - reskh)
        do j = 1, 7
            Resasc = Resasc + wgk(j) &
                     *(abs(fv1(j) - reskh) + abs(fv2(j) - reskh))
        end do
        Result = resk*hlgth
        Resabs = Resabs*dhlgth
        Resasc = Resasc*dhlgth
        Abserr = abs((resk - resg)*hlgth)
        if (Resasc /= 0.0_wp .and. Abserr /= 0.0_wp) &
            Abserr = Resasc*min(1.0_wp, (200.0_wp*Abserr/Resasc)**1.5_wp)
        if (Resabs > uflow/(50.0_wp*epmach)) &
            Abserr = max((epmach*50.0_wp)*Resabs, Abserr)

    end subroutine dqk15w