dgauss8 Subroutine

    subroutine dgauss8( f, a, b, error_tol, ans, ierr, err)

    implicit none

    procedure(func) :: f !! function subprogram defining the integrand function `f(x)`.
    real(wp),intent(in)   :: a          !! lower bound of the integration
    real(wp),intent(in)   :: b          !! 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 unit roundoff `d1mach(4)`.
                                        !! `ans` will normally have no more error than `abs(error_tol)`
                                        !! times the integral of the absolute value of `f(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 `a=b`.
                                        !!    * -1 : `a` and `b` 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.

    ! note: see also dqnc79 for some clues about the purpose of these numbers...
    real(wp),parameter  :: sq2    = sqrt(2.0_wp)
    real(wp),parameter  :: ln2    = log(2.0_wp)
    integer,parameter   :: kmx    = 5000
    integer,parameter   :: kml    = 6
    real(wp),parameter  :: magic  = 0.30102000_wp   !! is 0.30102000 supposed to be log10(2.0_wp) ???
    integer,parameter   :: iwork  = 60              !! size of the work arrays. ?? Why 60 ??
    integer,parameter   :: nbits  = int(d1mach(5)*i1mach14/magic)
    integer,parameter   :: nlmn   = 1
    integer,parameter   :: nlmx   = min(60,(nbits*5)/8)

    integer :: k !! number of function evaluations
    integer                   :: l,lmn,lmx,mxl,nib
    real(wp)                  :: ae,area,c,ee,ef,eps,est,gl,glr,tol
    real(wp),dimension(iwork) :: aa,hh,vl,gr
    integer,dimension(iwork)  :: lr

    ans = 0.0_wp
    ierr = 1
    err = 0.0_wp
    if (a == b) return

    aa = 0.0_wp
    hh = 0.0_wp
    vl = 0.0_wp
    gr = 0.0_wp
    lr = 0
    lmx = nlmx
    lmn = nlmn
    if (b /= 0.0_wp) then
        if (sign(1.0_wp,b)*a > 0.0_wp) then
            c = abs(1.0_wp-a/b)
            if (c <= 0.1_wp) then
                if (c <= 0.0_wp) return
                nib = int(0.5_wp - log(c)/ln2)
                lmx = min(nlmx,nbits-nib-7)
                if (lmx < 1) then
                    ! a and b 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
    if (error_tol == 0.0_wp) then
        tol = sqrt(epmach)
    else
        tol = max(abs(error_tol),2.0_wp**(5-nbits))/2.0_wp
    end if
    eps = tol
    hh(1) = (b-a)/4.0_wp
    aa(1) = a
    lr(1) = 1
    l = 1
    est = g(aa(l)+2.0_wp*hh(l),2.0_wp*hh(l))
    k = 8
    area = abs(est)
    ef = 0.5_wp
    mxl = 0

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

        gl = g(aa(l)+hh(l),hh(l))
        gr(l) = g(aa(l)+3.0_wp*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 > 0.0_wp) 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*0.5_wp
                ef = ef/sq2
                hh(l) = hh(l-1)*0.5_wp
                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*2.0_wp
                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) + 4.0_wp*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) + 4.0_wp*hh(l)
            cycle main
        end if

    end do main

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

    contains

    !************************************************************************************
    !>
    !  This is the 8-point formula from the original SLATEC routine
    !  [DGAUS8](http://www.netlib.org/slatec/src/dgaus8.f).
    !
    !@note replaced coefficients with high-precision ones from:
    !      http://processingjs.nihongoresources.com/bezierinfo/legendre-gauss-values.php

        function g(x, h)

        implicit none

        real(wp),intent(in) :: x
        real(wp),intent(in) :: h
        real(wp)            :: g

        !> abscissae:
        real(wp),parameter :: x1 = 0.18343464249564980493947614236018398066675781291297378231718847_wp
        real(wp),parameter :: x2 = 0.52553240991632898581773904918924634904196424312039285775085709_wp
        real(wp),parameter :: x3 = 0.79666647741362673959155393647583043683717173161596483207017029_wp
        real(wp),parameter :: x4 = 0.96028985649753623168356086856947299042823523430145203827163977_wp

        !> weights:
        real(wp),parameter :: w1 = 0.36268378337836198296515044927719561219414603989433054052482306_wp
        real(wp),parameter :: w2 = 0.31370664587788728733796220198660131326032899900273493769026394_wp
        real(wp),parameter :: w3 = 0.22238103445337447054435599442624088443013087005124956472590928_wp
        real(wp),parameter :: w4 = 0.10122853629037625915253135430996219011539409105168495705900369_wp

        g = h * ( w1*( f(x-x1*h) + f(x+x1*h) ) + &
                  w2*( f(x-x2*h) + f(x+x2*h) ) + &
                  w3*( f(x-x3*h) + f(x+x3*h) ) + &
                  w4*( f(x-x4*h) + f(x+x4*h) ) )

        end function g
    !************************************************************************************

    end subroutine dgauss8