dqagie Subroutine

    subroutine dqagie(f, Bound, Inf, Epsabs, Epsrel, Limit, Result, Abserr, &
                      Neval, Ier, Alist, Blist, Rlist, Elist, Iord, Last)
        implicit none

        procedure(func) :: f !! function subprogram defining the integrand function `f(x)`.
        integer, intent(in) :: Limit !! gives an upper bound on the number of subintervals
                                     !! in the partition of `(a,b)`, `limit>=1`
        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 transformed integration range (0,1).
        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 transformed integration range (0,1).
        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(out) :: Result !! approximation to the integral
        real(wp), intent(in) :: Bound !! finite bound of integration range
                                      !! (has no meaning if interval is doubly-infinite)
        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 result 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 (e.g. singularity,
                                    !!   discontinuity within the interval) one
                                    !!   will probably gain from splitting up the
                                    !!   interval at this point and calling the
                                    !!   integrator on the subranges. if possible,
                                    !!   an appropriate special-purpose integrator
                                    !!   should 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 assumed 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.
                                    !! * ier = 6 the input is invalid, because
                                    !!   `(epsabs<=0 and epsrel<max(50*rel.mach.acc.,0.5e-28)`,
                                    !!   `result`, `abserr`, `neval`, `last`, `rlist(1)`,
                                    !!   `elist(1)` and `iord(1)` are set to zero.
                                    !!   `alist(1)` and `blist(1)` are set to 0
                                    !!   and 1 respectively.
        integer, intent(in) :: Inf !! indicating the kind of integration range involved
                                   !! * inf = 1  corresponds to `(bound,+infinity)`
                                   !! * inf = -1 corresponds to `(-infinity,bound)`
                                   !! * inf = 2  corresponds to `(-infinity,+infinity)`
        integer, intent(out) :: Iord(Limit) !! vector of dimension `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 subdivision process
        integer, intent(out) :: Neval !! number of integrand evaluations

        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) :: errmax !! `elist(maxerr)`
        real(wp) :: erlast !! error on the interval currently subdivided
                           !! (before that subdivision has taken place)
        real(wp) :: area !! sum of the integrals over the subintervals
        real(wp) :: errsum !! sum of the errors over the subintervals
        real(wp) :: errbnd !! requested accuracy `max(epsabs,epsrel*abs(result))`
        real(wp) :: small !! length of the smallest interval considered up
                          !! to now, multiplied by 1.5
        real(wp) :: erlarg !! sum of the errors over the intervals larger
                           !! than the smallest interval considered up to now
        integer :: maxerr !! pointer to the interval with largest error estimate
        integer :: nres !! number of calls to the extrapolation routine
        integer :: numrl2 !! number of elements currently 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.
        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)
        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) :: abseps, boun, correc, defabs, dres, &
                    ertest, resabs, reseps, res3la(3)
        integer :: id, ierro, iroff1, iroff2, iroff3, &
                   jupbnd, k, ksgn, ktmin, nrmax

        ! test on validity of parameters

        Ier = 0
        Neval = 0
        Last = 0
        Result = 0.0_wp
        Abserr = 0.0_wp
        Alist(1) = 0.0_wp
        Blist(1) = 1.0_wp
        Rlist(1) = 0.0_wp
        Elist(1) = 0.0_wp
        Iord(1) = 0
        if (Epsabs <= 0.0_wp .and. Epsrel < max(50.0_wp*epmach, 0.5e-28_wp)) &
            Ier = 6
        if (Ier == 6) return

        main : block

            ! first approximation to the integral

            ! determine the interval to be mapped onto (0,1).
            ! if inf = 2 the integral is computed as i = i1+i2, where
            ! i1 = integral of f over (-infinity,0),
            ! i2 = integral of f over (0,+infinity).

            boun = Bound
            if (Inf == 2) boun = 0.0_wp
            call dqk15i(f, boun, Inf, 0.0_wp, 1.0_wp, Result, Abserr, defabs, &
                        resabs)

            ! test on accuracy

            Last = 1
            Rlist(1) = Result
            Elist(1) = Abserr
            Iord(1) = 1
            dres = abs(Result)
            errbnd = max(Epsabs, Epsrel*dres)
            if (Abserr <= 100.0_wp*epmach*defabs .and. Abserr > errbnd) Ier = 2
            if (Limit == 1) Ier = 1
            if (Ier /= 0 .or. (Abserr <= errbnd .and. Abserr /= resabs) .or. &
                Abserr == 0.0_wp) exit main

            ! initialization

            rlist2(1) = Result
            errmax = Abserr
            maxerr = 1
            area = Result
            errsum = Abserr
            Abserr = oflow
            nrmax = 1
            nres = 0
            ktmin = 0
            numrl2 = 2
            extrap = .false.
            noext = .false.
            ierro = 0
            iroff1 = 0
            iroff2 = 0
            iroff3 = 0
            ksgn = -1
            if (dres >= (1.0_wp - 50.0_wp*epmach)*defabs) ksgn = 1

            ! main do-loop

            loop: do Last = 2, Limit

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

                a1 = Alist(maxerr)
                b1 = 0.5_wp*(Alist(maxerr) + Blist(maxerr))
                a2 = b1
                b2 = Blist(maxerr)
                erlast = errmax
                call dqk15i(f, boun, Inf, a1, b1, area1, error1, resabs, defab1)
                call dqk15i(f, boun, Inf, a2, b2, area2, error2, resabs, defab2)

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

                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
                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 some points 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)
                if (errsum <= errbnd) then
                    ! compute global integral sum.
                    Result = sum(Rlist(1:Last))
                    Abserr = errsum
                    exit main
                end if
                if (Ier /= 0) exit
                if (Last == 2) then
                    small = 0.375_wp
                    erlarg = errsum
                    ertest = errbnd
                    rlist2(2) = area
                elseif (.not. (noext)) then
                    erlarg = erlarg - erlast
                    if (abs(b1 - a1) > small) erlarg = erlarg + erro12
                    if (.not. (extrap)) then
                        ! test whether the interval to be bisected next is the
                        ! smallest interval.
                        if (abs(Blist(maxerr) - Alist(maxerr)) > small) 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)
                            if (abs(Blist(maxerr) - Alist(maxerr)) > small) cycle loop
                            nrmax = nrmax + 1
                        end do
                    end if

                    ! perform extrapolation.

                    numrl2 = numrl2 + 1
                    rlist2(numrl2) = area
                    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))
                        if (Abserr <= ertest) exit
                    end if

                    ! prepare bisection of the smallest interval.

                    if (numrl2 == 1) noext = .true.
                    if (Ier == 5) exit
                    maxerr = Iord(1)
                    errmax = Elist(maxerr)
                    nrmax = 1
                    extrap = .false.
                    small = small*0.5_wp
                    erlarg = errsum
                end if

            end do loop

            ! set final result and error estimate.

            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)) > defabs*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

        Neval = 30*Last - 15
        if (Inf == 2) Neval = 2*Neval
        if (Ier > 2) Ier = Ier - 1

    end subroutine dqagie