dqagpe Subroutine

    subroutine dqagpe(f, a, b, Npts2, Points, Epsabs, Epsrel, Limit, Result, &
                      Abserr, Neval, Ier, Alist, Blist, Rlist, Elist, Pts, &
                      Iord, Level, Ndin, Last)
        implicit none

        procedure(func) :: f
        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) :: Alist(Limit) !! vector of dimension at least `limit`, the first
                                              !! `last` elements of which are the left end points
                                              !! of the subintervals in the partition of the given
                                              !! integration range (a,b)
        real(wp), intent(out)  :: Blist(Limit) !! vector of dimension at least `limit`, the first
                                               !! `last` elements of which are the right end points
                                               !! of the subintervals in the partition of the given
                                               !! integration range (a,b)
        real(wp), intent(out)  :: Elist(Limit) !! vector of dimension at least `limit`, the first
                                               !! `last` elements of which are the moduli of the
                                               !! absolute error estimates on the subintervals
        real(wp), intent(out)  :: Rlist(Limit) !! vector of dimension at least `limit`, the first
                                               !! `last` elements of which are the integral
                                               !! approximations on the subintervals
        real(wp), intent(in) :: Epsabs !! absolute accuracy requested
        real(wp), intent(in) :: Epsrel !! relative accuracy requested
                                       !! if  `epsabs<=0`
                                       !! and `epsrel<max(50*rel.mach.acc.,0.5e-28)`,
                                       !! the routine will end with ier = 6.
        real(wp), intent(in) :: Points(Npts2) !! vector of dimension npts2, the first (npts2-2)
                                              !! elements of which are the user provided break
                                              !! points. if these points do not constitute an
                                              !! ascending sequence there will be an automatic
                                              !! sorting.
        real(wp), intent(out) :: Pts(Npts2) !! vector of dimension at least npts2, containing the
                                            !! integration limits and the break points of the
                                            !! interval in ascending sequence.
        integer, intent(out) :: Ier !! * ier = 0 normal and reliable termination of the
                                    !!   routine. it is assumed that the requested
                                    !!   accuracy has been achieved.
                                    !! * ier>0 abnormal termination of the routine.
                                    !!   the estimates for integral and error are
                                    !!   less reliable. it is assumed that the
                                    !!   requested accuracy has not been achieved.
                                    !!
                                    !! error messages:
                                    !!
                                    !! * ier = 1 maximum number of subdivisions allowed
                                    !!   has been achieved. one can allow more
                                    !!   subdivisions by increasing the value of
                                    !!   limit (and taking the according dimension
                                    !!   adjustments into account). however, if
                                    !!   this yields no improvement it is advised
                                    !!   to analyze the integrand in order to
                                    !!   determine the integration difficulties. if
                                    !!   the position of a local difficulty can be
                                    !!   determined (i.e. singularity,
                                    !!   discontinuity within the interval), it
                                    !!   should be supplied to the routine as an
                                    !!   element of the vector points. if necessary
                                    !!   an appropriate special-purpose integrator
                                    !!   must be used, which is designed for
                                    !!   handling the type of difficulty involved.
                                    !! * ier = 2 the occurrence of roundoff error is
                                    !!   detected, which prevents the requested
                                    !!   tolerance from being achieved.
                                    !!   the error may be under-estimated.
                                    !! * ier = 3 extremely bad integrand behaviour occurs
                                    !!   at some points of the integration
                                    !!   interval.
                                    !! * ier = 4 the algorithm does not converge.
                                    !!   roundoff error is detected in the
                                    !!   extrapolation table. it is presumed that
                                    !!   the requested tolerance cannot be
                                    !!   achieved, and that the returned result is
                                    !!   the best which can be obtained.
                                    !! * ier = 5 the integral is probably divergent, or
                                    !!   slowly convergent. it must be noted that
                                    !!   divergence can occur with any other value
                                    !!   of ier>0.
                                    !! * ier = 6 the input is invalid because
                                    !!   `npts2<2` or
                                    !!   break points are specified outside
                                    !!   the integration range or
                                    !!   `(epsabs<=0 and epsrel<max(50*rel.mach.acc.,0.5e-28))`
                                    !!   or `limit<npts2`.
                                    !!   `result`, `abserr`, `neval`, `last`, `rlist(1)`,
                                    !!   and elist(1) are set to zero. alist(1) and
                                    !!   blist(1) are set to `a` and `b` respectively.
        integer, intent(out) :: Iord(Limit) !! vector of dimension at least `limit`, the first `k`
                                            !! elements of which are pointers to the
                                            !! error estimates over the subintervals,
                                            !! such that `elist(iord(1)), ..., elist(iord(k))`
                                            !! form a decreasing sequence, with `k = last`
                                            !! if `last<=(limit/2+2)`, and `k = limit+1-last`
                                            !! otherwise
        integer, intent(out) :: Last !! number of subintervals actually produced in the
                                     !! subdivisions process
        integer, intent(in) :: Limit !! gives an upper bound on the number of subintervals
                                     !! in the partition of `(a,b)`, `limit>=npts2`
                                     !! if `limit<npts2`, the routine will end with
                                     !! ier = 6.
        integer, intent(in) :: Npts2 !! number equal to two more than the number of
                                     !! user-supplied break points within the integration
                                     !! range, `npts2>=2`.
                                     !! if `npts2<2`, the routine will end with ier = 6.
        integer, intent(out) :: Ndin(Npts2) !! vector of dimension at least npts2, after first
                                            !! integration over the intervals `(pts(i)),pts(i+1)`,
                                            !! `i = 0,1, ..., npts2-2`, the error estimates over
                                            !! some of the intervals may have been increased
                                            !! artificially, in order to put their subdivision
                                            !! forward. if this happens for the subinterval
                                            !! numbered `k`, `ndin(k)` is put to 1, otherwise
                                            !! `ndin(k) = 0`.
        integer, intent(out) :: Neval !! number of integrand evaluations
        integer, intent(out) :: Level(Limit) !! vector of dimension at least `limit`, containing the
                                             !! subdivision levels of the subinterval, i.e. if
                                             !! `(aa,bb)` is a subinterval of `(p1,p2)` where `p1` as
                                             !! well as `p2` is a user-provided break point or
                                             !! integration limit, then `(aa,bb)` has level `l` if
                                             !! `abs(bb-aa) = abs(p2-p1)*2**(-l)`.

        real(wp) :: a, abseps, b, correc, defabs, &
                    dres, ertest, resa, reseps, Result, &
                    res3la(3), sign, temp, resabs
        integer :: i, id, ierro, ind1, ind2, ip1, iroff1, &
                   iroff2, iroff3, j, jlow, jupbnd, k, ksgn, ktmin, &
                   levcur, levmax, nint, nintp1, npts, nrmax
        real(wp) :: area1, a1, b1, defab1, error1 !! variable for the left subinterval
        real(wp) :: area2, a2, b2, defab2, error2 !! variable for the right subinterval
        real(wp) :: area12 !! `area1 + area2`
        real(wp) :: erro12 !! `error1 + error2`
        real(wp) :: rlist2(limexp + 2) !! array of dimension at least `limexp+2`
                                       !! containing the part of the epsilon table which
                                       !! is still needed for further computations.
        real(wp) :: erlast !! error on the interval currently subdivided
                           !! (before that subdivision has taken place)
        real(wp) :: errsum !! sum of the errors over the subintervals
        real(wp) :: errbnd !! requested accuracy `max(epsabs,epsrel*abs(result))`
        real(wp) :: area !! sum of the integrals over the subintervals
        real(wp) :: erlarg !! sum of the errors over the intervals larger
                           !! than the smallest interval considered up to now
        real(wp) :: errmax !! `elist(maxerr)`
        logical :: extrap !! logical variable denoting that the routine
                          !! is attempting to perform extrapolation. i.e.
                          !! before subdividing the smallest interval we
                          !! try to decrease the value of `erlarg`.
        logical :: noext !! logical variable denoting that extrapolation is
                         !! no longer allowed (true-value)
        integer :: maxerr !! pointer to the interval with largest error estimate
        integer :: nres !! number of calls to the extrapolation routine
        integer :: numrl2 !! number of elements in `rlist2`. if an appropriate
                          !! approximation to the compounded integral has
                          !! been obtained, it is put in `rlist2(numrl2)` after
                          !! `numrl2` has been increased by one.

        ! test on validity of parameters

        Ier = 0
        Neval = 0
        Last = 0
        Result = 0.0_wp
        Abserr = 0.0_wp
        Alist(1) = a
        Blist(1) = b
        Rlist(1) = 0.0_wp
        Elist(1) = 0.0_wp
        Iord(1) = 0
        Level(1) = 0
        npts = Npts2 - 2
        if (Npts2 < 2 .or. Limit <= npts .or. &
            (Epsabs <= 0.0_wp .and. Epsrel < max(50.0_wp*epmach, 0.5e-28_wp))) &
            Ier = 6
        if (Ier == 6) return

        ! if any break points are provided, sort them into an
        ! ascending sequence.

        sign = 1.0_wp
        if (a > b) sign = -1.0_wp
        Pts(1) = min(a, b)
        if (npts /= 0) then
            do i = 1, npts
                Pts(i + 1) = Points(i)
            end do
        end if
        Pts(npts + 2) = max(a, b)
        nint = npts + 1
        a1 = Pts(1)
        if (npts /= 0) then
            nintp1 = nint + 1
            do i = 1, nint
                ip1 = i + 1
                do j = ip1, nintp1
                    if (Pts(i) > Pts(j)) then
                        temp = Pts(i)
                        Pts(i) = Pts(j)
                        Pts(j) = temp
                    end if
                end do
            end do
            if (Pts(1) /= min(a, b) .or. Pts(nintp1) /= max(a, b)) Ier = 6
            if (Ier == 6) return
        end if

        main : block

            ! compute first integral and error approximations.

            resabs = 0.0_wp
            do i = 1, nint
                b1 = Pts(i + 1)
                call dqk21(f, a1, b1, area1, error1, defabs, resa)
                Abserr = Abserr + error1
                Result = Result + area1
                Ndin(i) = 0
                if (error1 == resa .and. error1 /= 0.0_wp) Ndin(i) = 1
                resabs = resabs + defabs
                Level(i) = 0
                Elist(i) = error1
                Alist(i) = a1
                Blist(i) = b1
                Rlist(i) = area1
                Iord(i) = i
                a1 = b1
            end do
            errsum = 0.0_wp
            do i = 1, nint
                if (Ndin(i) == 1) Elist(i) = Abserr
                errsum = errsum + Elist(i)
            end do

            ! test on accuracy.

            Last = nint
            Neval = 21*nint
            dres = abs(Result)
            errbnd = max(Epsabs, Epsrel*dres)
            if (Abserr <= 100.0_wp*epmach*resabs .and. Abserr > errbnd) Ier = 2
            if (nint /= 1) then
            do i = 1, npts
                jlow = i + 1
                ind1 = Iord(i)
                do j = jlow, nint
                    ind2 = Iord(j)
                    if (Elist(ind1) <= Elist(ind2)) then
                        ind1 = ind2
                        k = j
                    end if
                end do
                if (ind1 /= Iord(i)) then
                    Iord(k) = Iord(i)
                    Iord(i) = ind1
                end if
            end do
            if (Limit < Npts2) Ier = 1
            end if
            if (Ier /= 0 .or. Abserr <= errbnd) exit main

            ! initialization

            rlist2(1) = Result
            maxerr = Iord(1)
            errmax = Elist(maxerr)
            area = Result
            nrmax = 1
            nres = 0
            numrl2 = 1
            ktmin = 0
            extrap = .false.
            noext = .false.
            erlarg = errsum
            ertest = errbnd
            levmax = 1
            iroff1 = 0
            iroff2 = 0
            iroff3 = 0
            ierro = 0
            Abserr = oflow
            ksgn = -1
            if (dres >= (1.0_wp - 50.0_wp*epmach)*resabs) ksgn = 1

            ! main do-loop

            loop: do Last = Npts2, Limit

                ! bisect the subinterval with the nrmax-th largest error
                ! estimate.

                levcur = Level(maxerr) + 1
                a1 = Alist(maxerr)
                b1 = 0.5_wp*(Alist(maxerr) + Blist(maxerr))
                a2 = b1
                b2 = Blist(maxerr)
                erlast = errmax
                call dqk21(f, a1, b1, area1, error1, resa, defab1)
                call dqk21(f, a2, b2, area2, error2, resa, defab2)

                ! improve previous approximations to integral
                ! and error and test for accuracy.

                Neval = Neval + 42
                area12 = area1 + area2
                erro12 = error1 + error2
                errsum = errsum + erro12 - errmax
                area = area + area12 - Rlist(maxerr)
                if (defab1 /= error1 .and. defab2 /= error2) then
                    if (abs(Rlist(maxerr) - area12) <= 0.1e-4_wp*abs(area12) .and. &
                        erro12 >= 0.99_wp*errmax) then
                        if (extrap) iroff2 = iroff2 + 1
                        if (.not. extrap) iroff1 = iroff1 + 1
                    end if
                    if (Last > 10 .and. erro12 > errmax) iroff3 = iroff3 + 1
                end if
                Level(maxerr) = levcur
                Level(Last) = levcur
                Rlist(maxerr) = area1
                Rlist(Last) = area2
                errbnd = max(Epsabs, Epsrel*abs(area))

                ! test for roundoff error and eventually set error flag.

                if (iroff1 + iroff2 >= 10 .or. iroff3 >= 20) Ier = 2
                if (iroff2 >= 5) ierro = 3

                ! set error flag in the case that the number of
                ! subintervals equals limit.

                if (Last == Limit) Ier = 1

                ! set error flag in the case of bad integrand behaviour
                ! at a point of the integration range

                if (max(abs(a1), abs(b2)) <= (1.0_wp + 100.0_wp*epmach) &
                    *(abs(a2) + 1000.0_wp*uflow)) Ier = 4

                ! append the newly-created intervals to the list.

                if (error2 > error1) then
                    Alist(maxerr) = a2
                    Alist(Last) = a1
                    Blist(Last) = b1
                    Rlist(maxerr) = area2
                    Rlist(Last) = area1
                    Elist(maxerr) = error2
                    Elist(Last) = error1
                else
                    Alist(Last) = a2
                    Blist(maxerr) = b1
                    Blist(Last) = b2
                    Elist(maxerr) = error1
                    Elist(Last) = error2
                end if

                ! call subroutine dqpsrt to maintain the descending ordering
                ! in the list of error estimates and select the subinterval
                ! with nrmax-th largest error estimate (to be bisected next).

                call dqpsrt(Limit, Last, maxerr, errmax, Elist, Iord, nrmax)
                ! ***jump out of do-loop
                if (errsum <= errbnd) then
                    ! compute global integral sum.
                    Result = sum(Rlist(1:Last))
                    Abserr = errsum
                    exit main
                end if
                ! ***jump out of do-loop
                if (Ier /= 0) exit loop
                if (.not. (noext)) then
                    erlarg = erlarg - erlast
                    if (levcur + 1 <= levmax) erlarg = erlarg + erro12
                    if (.not. (extrap)) then
                        ! test whether the interval to be bisected next is the
                        ! smallest interval.
                        if (Level(maxerr) + 1 <= levmax) cycle loop
                        extrap = .true.
                        nrmax = 2
                    end if
                    if (ierro /= 3 .and. erlarg > ertest) then
                        ! the smallest interval has the largest error.
                        ! before bisecting decrease the sum of the errors over
                        ! the larger intervals (erlarg) and perform extrapolation.
                        id = nrmax
                        jupbnd = Last
                        if (Last > (2 + Limit/2)) jupbnd = Limit + 3 - Last
                        do k = id, jupbnd
                            maxerr = Iord(nrmax)
                            errmax = Elist(maxerr)
                            ! ***jump out of do-loop
                            if (Level(maxerr) + 1 <= levmax) cycle loop
                            nrmax = nrmax + 1
                        end do
                    end if

                    ! perform extrapolation.

                    numrl2 = numrl2 + 1
                    rlist2(numrl2) = area
                    if (numrl2 > 2) then
                        call dqelg(numrl2, rlist2, reseps, abseps, res3la, nres)
                        ktmin = ktmin + 1
                        if (ktmin > 5 .and. Abserr < 0.1e-02_wp*errsum) Ier = 5
                        if (abseps < Abserr) then
                            ktmin = 0
                            Abserr = abseps
                            Result = reseps
                            correc = erlarg
                            ertest = max(Epsabs, Epsrel*abs(reseps))
                            ! ***jump out of do-loop
                            if (Abserr < ertest) exit loop
                        end if
                        ! prepare bisection of the smallest interval.
                        if (numrl2 == 1) noext = .true.
                        if (Ier >= 5) exit loop
                    end if
                    maxerr = Iord(1)
                    errmax = Elist(maxerr)
                    nrmax = 1
                    extrap = .false.
                    levmax = levmax + 1
                    erlarg = errsum
                end if

            end do loop

            ! set the final result.

            if (Abserr /= oflow) then
                if ((Ier + ierro) /= 0) then
                    if (ierro == 3) Abserr = Abserr + correc
                    if (Ier == 0) Ier = 3
                    if (Result == 0.0_wp .or. area == 0.0_wp) then
                        if (Abserr > errsum) then
                            ! compute global integral sum.
                            Result = sum(Rlist(1:Last))
                            Abserr = errsum
                            exit main
                        end if
                        if (area == 0.0_wp) exit main
                    elseif (Abserr/abs(Result) > errsum/abs(area)) then
                        ! compute global integral sum.
                        Result = sum(Rlist(1:Last))
                        Abserr = errsum
                        exit main
                    end if
                end if

                ! test on divergence.

                if (ksgn /= (-1) .or. max(abs(Result), abs(area)) &
                    > resabs*0.01_wp) then
                    if (0.01_wp > (Result/area) .or. (Result/area) > 100.0_wp .or. &
                        errsum > abs(area)) Ier = 6
                end if

            else
                ! compute global integral sum.
                Result = sum(Rlist(1:Last))
                Abserr = errsum
            end if

        end block main

        if (Ier > 2) Ier = Ier - 1
        Result = Result*sign

    end subroutine dqagpe