dqk41 Subroutine

    subroutine dqk41(f, a, b, Result, Abserr, Resabs, Resasc)
        implicit none

        procedure(func) :: f !! function subprogram defining the integrand function `f(x)`.
        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 41-point
                                        !! gauss-kronrod rule (resk) obtained by optimal
                                        !! addition of abscissae to the 20-point gauss
                                        !! rule (resg).
        real(wp), intent(out) :: Abserr !! estimate of the modulus of the absolute error,
                                        !! which should not exceed `abs(i-result)`
        real(wp), intent(out) :: Resabs !! approximation to the integral j
        real(wp), intent(out) :: Resasc !! approximation to the integral of abs(f-i/(b-a))
                                        !! over `(a,b)`

        real(wp) :: dhlgth, fc, fsum, fv1(20), fv2(20)
        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 20-point gauss formula
        real(wp) :: resk !! result of the 41-point kronrod formula
        real(wp) :: reskh !! approximation to mean value of `f` 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(10), parameter :: wg = [ &
                                              1.76140071391521183118619623518528163621e-2_wp, &
                                              4.06014298003869413310399522749321098791e-2_wp, &
                                              6.26720483341090635695065351870416063516e-2_wp, &
                                              8.32767415767047487247581432220462061002e-2_wp, &
                                              1.01930119817240435036750135480349876167e-1_wp, &
                                              1.18194531961518417312377377711382287005e-1_wp, &
                                              1.31688638449176626898494499748163134916e-1_wp, &
                                              1.42096109318382051329298325067164933035e-1_wp, &
                                              1.49172986472603746787828737001969436693e-1_wp, &
                                              1.52753387130725850698084331955097593492e-1_wp] !! weights of the 20-point gauss rule

        real(wp), dimension(21), parameter :: xgk = [ &
                                              9.98859031588277663838315576545863010000e-1_wp, &
                                              9.93128599185094924786122388471320278223e-1_wp, &
                                              9.81507877450250259193342994720216944567e-1_wp, &
                                              9.63971927277913791267666131197277221912e-1_wp, &
                                              9.40822633831754753519982722212443380274e-1_wp, &
                                              9.12234428251325905867752441203298113049e-1_wp, &
                                              8.78276811252281976077442995113078466711e-1_wp, &
                                              8.39116971822218823394529061701520685330e-1_wp, &
                                              7.95041428837551198350638833272787942959e-1_wp, &
                                              7.46331906460150792614305070355641590311e-1_wp, &
                                              6.93237656334751384805490711845931533386e-1_wp, &
                                              6.36053680726515025452836696226285936743e-1_wp, &
                                              5.75140446819710315342946036586425132814e-1_wp, &
                                              5.10867001950827098004364050955250998425e-1_wp, &
                                              4.43593175238725103199992213492640107840e-1_wp, &
                                              3.73706088715419560672548177024927237396e-1_wp, &
                                              3.01627868114913004320555356858592260615e-1_wp, &
                                              2.27785851141645078080496195368574624743e-1_wp, &
                                              1.52605465240922675505220241022677527912e-1_wp, &
                                              7.65265211334973337546404093988382110048e-2_wp, &
                                              0.00000000000000000000000000000000000000e0_wp] !! abscissae of the 41-point gauss-kronrod rule:
                                                                                             !!
                                                                                             !! * xgk(2), xgk(4), ...  abscissae of the 20-point
                                                                                             !!   gauss rule
                                                                                             !! * xgk(1), xgk(3), ...  abscissae which are optimally
                                                                                             !!   added to the 20-point gauss rule

        real(wp), dimension(21), parameter :: wgk = [ &
                                              3.07358371852053150121829324603098748803e-3_wp, &
                                              8.60026985564294219866178795010234725213e-3_wp, &
                                              1.46261692569712529837879603088683561639e-2_wp, &
                                              2.03883734612665235980102314327547051228e-2_wp, &
                                              2.58821336049511588345050670961531429995e-2_wp, &
                                              3.12873067770327989585431193238007378878e-2_wp, &
                                              3.66001697582007980305572407072110084875e-2_wp, &
                                              4.16688733279736862637883059368947380440e-2_wp, &
                                              4.64348218674976747202318809261075168421e-2_wp, &
                                              5.09445739237286919327076700503449486648e-2_wp, &
                                              5.51951053482859947448323724197773291948e-2_wp, &
                                              5.91114008806395723749672206485942171364e-2_wp, &
                                              6.26532375547811680258701221742549805858e-2_wp, &
                                              6.58345971336184221115635569693979431472e-2_wp, &
                                              6.86486729285216193456234118853678017155e-2_wp, &
                                              7.10544235534440683057903617232101674129e-2_wp, &
                                              7.30306903327866674951894176589131127606e-2_wp, &
                                              7.45828754004991889865814183624875286161e-2_wp, &
                                              7.57044976845566746595427753766165582634e-2_wp, &
                                              7.63778676720807367055028350380610018008e-2_wp, &
                                              7.66007119179996564450499015301017408279e-2_wp] !! weights of the 41-point gauss-kronrod rule

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

        ! compute the 41-point gauss-kronrod approximation to
        ! the integral, and estimate the absolute error.

        resg = 0.0_wp
        fc = f(centr)
        resk = wgk(21)*fc
        Resabs = abs(resk)
        do j = 1, 10
            jtw = j*2
            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)*(abs(fval1) + abs(fval2))
        end do
        do j = 1, 10
            jtwm1 = j*2 - 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)*(abs(fval1) + abs(fval2))
        end do
        reskh = resk*0.5_wp
        Resasc = wgk(21)*abs(fc - reskh)
        do j = 1, 20
            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 dqk41