dgauss_generic Subroutine

    recursive subroutine dgauss_generic (me, lb, ub, error_tol, ans, ierr, err)

    implicit none

    class(integration_class_1d),intent(inout)  :: me
    real(wp),intent(in)   :: lb         !! lower bound of the integration
    real(wp),intent(in)   :: ub         !! upper bound of the integration
    real(wp),intent(in)   :: error_tol  !! is a requested pseudorelative error tolerance.  normally
                                        !! pick a value of abs(error_tol) so that
                                        !! dtol < abs(error_tol) <= 1.0e-3 where dtol is the larger
                                        !! of 1.0e-18 and the real(wp) unit roundoff d1mach(4).
                                        !! ans will normally have no more error than abs(error_tol)
                                        !! times the integral of the absolute value of fun(x).  usually,
                                        !! smaller values of error_tol yield more accuracy and require
                                        !! more function evaluations.
    real(wp),intent(out)  :: ans        !! computed value of integral
    integer,intent(out)   :: ierr       !! status code:
                                        !!
                                        !!  * normal codes:
                                        !!    * 1 : `ans` most likely meets requested error tolerance,
                                        !!      or `lb=ub`.
                                        !!    * -1 : `lb` and `ub` are too nearly equal to allow normal
                                        !!      integration. `ans` is set to zero.
                                        !!  * abnormal code:
                                        !!    * 2 : `ans` probably does not meet requested error tolerance.
    real(wp),intent(out)  :: err        !! an estimate of the absolute error in `ans`.
                                        !! the estimated error is solely for information to the user and
                                        !! should not be used as a correction to the computed integral.

    real(wp),parameter  :: sq2      = sqrt(two)
    real(wp),parameter  :: ln2      = log(two)
    integer,parameter   :: nlmn     = 1                   !! ??
    integer,parameter   :: kmx      = 5000                !! ??
    integer,parameter   :: kml      = 6                   !! ??
    real(wp),parameter  :: magic    = 0.30102000_wp       !! ??
    integer,parameter   :: iwork    = 60                  !! size of the work arrays. ?? Why 60 ??
    real(wp),parameter  :: bb       = radix(one)          !! machine constant
    real(wp),parameter  :: d1mach4  = bb**(1-digits(one)) !! machine constant
    real(wp),parameter  :: d1mach5  = log10(bb)           !! machine constant

    integer                   :: k,l,lmn,lmx,mxl,nbits,nib,nlmx
    real(wp)                  :: ae,anib,area,c,ee,ef,eps,est,gl,glr,tol
    real(wp),dimension(iwork) :: aa,hh,vl,gr
    integer,dimension(iwork)  :: lr

    ans = zero
    ierr = 1
    err = zero
    if (lb == ub) return
    aa = zero
    hh = zero
    vl = zero
    gr = zero
    lr = 0
    k = digits(one)
    anib = d1mach5*k/magic
    nbits = anib
    nlmx = min(60,(nbits*5)/8)         ! ... is this the same 60 as iwork???
    lmx = nlmx
    lmn = nlmn
    if (ub /= zero) then
        if (sign(one,ub)*lb > zero) then
            c = abs(one-lb/ub)
            if (c <= 0.1_wp) then
                if (c <= zero) return
                anib = one_half - log(c)/ln2
                nib = anib
                lmx = min(nlmx,nbits-nib-7)
                if (lmx < 1) then
                    ! lb and ub are too nearly equal to allow
                    ! normal integration [ans is set to zero]
                    ierr = -1
                    return
                end if
                lmn = min(lmn,lmx)
            end if
        end if
    end if
    tol = max(abs(error_tol),two**(5-nbits))/two
    if (error_tol == zero) tol = sqrt(d1mach4)
    eps = tol
    hh(1) = (ub-lb)/four
    aa(1) = lb
    lr(1) = 1
    l = 1
    est = me%g(aa(l)+two*hh(l),two*hh(l))
    k = 8
    area = abs(est)
    ef = one_half
    mxl = 0

    !compute refined estimates, estimate the error, etc.
    main : do

        gl = me%g(aa(l)+hh(l),hh(l))
        gr(l) = me%g(aa(l)+three*hh(l),hh(l))
        k = k + 16
        area = area + (abs(gl)+abs(gr(l))-abs(est))
        glr = gl + gr(l)
        ee = abs(est-glr)*ef
        ae = max(eps*area,tol*abs(glr))
        if (ee-ae > zero) then
            !consider the left half of this level
            if (k > kmx) lmx = kml
            if (l >= lmx) then
                mxl = 1
            else
                l = l + 1
                eps = eps*one_half
                ef = ef/sq2
                hh(l) = hh(l-1)*one_half
                lr(l) = -1
                aa(l) = aa(l-1)
                est = gl
                cycle main
            end if
        end if

        err = err + (est-glr)
        if (lr(l) > 0) then
            !return one level
            ans = glr
            do
                if (l <= 1) exit main ! finished
                l = l - 1
                eps = eps*two
                ef = ef*sq2
                if (lr(l) <= 0) then
                    vl(l) = vl(l+1) + ans
                    est = gr(l-1)
                    lr(l) = 1
                    aa(l) = aa(l) + four*hh(l)
                    cycle main
                end if
                ans = vl(l+1) + ans
            end do
        else
            !proceed to right half at this level
            vl(l) = glr
            est = gr(l-1)
            lr(l) = 1
            aa(l) = aa(l) + four*hh(l)
            cycle main
        end if

    end do main

    if ((mxl/=0) .and. (abs(err)>two*tol*area)) ierr = 2 ! ans is probably insufficiently accurate

    end subroutine dgauss_generic