dqnc79 Subroutine

    subroutine dqnc79(fun,a,b,err,ans,ierr,k)

    implicit none

    procedure(func) :: fun !! 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 (may be less than `A`)
    real(wp),intent(in) :: err !! a requested error tolerance.  Normally, pick a value
                               !! `0 < ERR < 1.0e-8`.
    real(wp),intent(out) :: ans !! computed value of the integral.  Hopefully, `ANS` is
                                !! accurate to within `ERR *` integral of `ABS(FUN(X))`.
    integer,intent(out) :: ierr !! a status code:
                                !!
                                !!  * Normal codes
                                !!    * **1** `ANS` most likely meets requested error tolerance.
                                !!    * **-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.
    integer,intent(out) :: k !! the number of function evaluations actually used to do
                             !! the integration.  A value of `K > 1000` indicates a
                             !! difficult problem; other programs may be more efficient.
                             !! `DQNC79` will gracefully give up if `K` exceeds 5000.

    real(wp),parameter :: w1 = 41.0_wp/140.0_wp
    real(wp),parameter :: w2 = 216.0_wp/140.0_wp
    real(wp),parameter :: w3 = 27.0_wp/140.0_wp
    real(wp),parameter :: w4 = 272.0_wp/140.0_wp
    real(wp),parameter :: sq2 = sqrt(2.0_wp)
    real(wp),parameter :: ln2 = log(2.0_wp)
    integer,parameter :: nbits = int(d1mach(5)*i1mach14/0.30102000_wp)  !! is 0.30102000 supposed to be log10(2.0_wp) ???
    integer,parameter :: nlmx = min(99, (nbits*4)/5)
    integer,parameter :: nlmn = 2
    integer,parameter :: kml = 7
    integer,parameter :: kmx = 5000      !! JW : is this the max function evals? should be an input
    integer,parameter :: array_size = 99 !! JW : what is this magic number 99 array size ??
                                         !! does it depend on the number of function evals ?
                                         !! (see comment in revision history)

    real(wp) :: ae,area,bank,blocal,c,ce,ee,ef,eps,q13,q7,q7l,test,tol,vr
    integer :: i,l,lmn,lmx,nib
    real(wp),dimension(13) :: f
    real(wp),dimension(array_size) :: aa,f1,f2,f3,f4,f5,f6,f7,hh,q7r,vl
    integer,dimension(array_size) :: lr

    ans = 0.0_wp
    ierr = 1
    if ( a==b ) return ! JW : this was an error return in the original code

    ce = 0.0_wp
    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 ) then
                    ierr = -1
                    call xerror('a and b are too nearly equal to allow normal integration. '&
                                //'ans is set to zero and ierr to -1.',-1,-1)
                    return
                end if
                nib = 0.5_wp - log(c)/ln2
                lmx = min(nlmx,nbits-nib-4)
                if ( lmx<2 ) then
                    call xerror('a and b are too nearly equal to allow normal integration. '&
                                //'ans is set to zero and ierr to -1.',-1,-1)
                    return
                end if
                lmn = min(lmn,lmx)
            endif
        endif
    endif
    tol = max(abs(err),2.0_wp**(5-nbits))
    if ( err==0.0_wp ) tol = sqrt(epmach)
    eps = tol
    hh(1) = (b-a)/12.0_wp
    aa(1) = a
    lr(1) = 1
    do i = 1 , 11 , 2
        f(i) = fun(a+(i-1)*hh(1))
    enddo
    blocal = b
    f(13) = fun(blocal)
    k = 7
    l = 1
    area = 0.0_wp
    q7 = 0.0_wp
    ef = 256.0_wp/255.0_wp
    bank = 0.0_wp

    loop : do

        ! compute refined estimates, estimate the error, etc.
        do i = 2 , 12 , 2
            f(i) = fun(aa(l)+(i-1)*hh(l))
        enddo
        k = k + 6

        ! compute left and right half estimates
        q7l = hh(l)*((w1*(f(1)+f(7))+w2*(f(2)+f(6)))+(w3*(f(3)+f(5))+w4*f(4)))
        q7r(l) = hh(l)*((w1*(f(7)+f(13))+w2*(f(8)+f(12)))+(w3*(f(9)+f(11))+w4*f(10)))

        ! update estimate of integral of absolute value
        area = area + (abs(q7l)+abs(q7r(l))-abs(q7))

        ! do not bother to test convergence before minimum refinement level
        if ( l>=lmn ) then

            ! estimate the error in new value for whole interval, q13
            q13 = q7l + q7r(l)
            ee = abs(q7-q13)*ef

            ! compute nominal allowed error
            ae = eps*area

            ! borrow from bank account, but not too much
            test = min(ae+0.8_wp*bank,10.0_wp*ae)

            ! don't ask for excessive accuracy
            test = max(test,tol*abs(q13),0.00003_wp*tol*area)   ! jw : should change ?

            ! now, did this interval pass or not?
            if ( ee<=test ) then
                ! on good intervals accumulate the theoretical estimate
                ce = ce + (q7-q13)/255.0_wp
            else
                ! consider the left half of next deeper level
                if ( k>kmx ) lmx = min(kml,lmx)
                if ( l<lmx ) then
                    call f200()
                    cycle loop
                end if
                ! have hit maximum refinement level -- penalize the cumulative error
                ce = ce + (q7-q13)
            endif

            ! update the bank account.  don't go into debt.
            bank = bank + (ae-ee)
            if ( bank<0.0_wp ) bank = 0.0_wp

            ! did we just finish a left half or a right half?
            if ( lr(l)<=0 ) then
                ! proceed to right half at this level
                vl(l) = q13
                call f300()
                cycle loop
            else
                ! left and right halves are done, so go back up a level
                vr = q13
                do
                    if ( l<=1 ) then
                        !   exit
                        ans = vr
                        if ( abs(ce)>2.0_wp*tol*area ) then
                            ierr = 2
                            call xerror('ans is probably insufficiently accurate.',2,1)
                        endif
                        return
                    else
                        if ( l<=17 ) ef = ef*sq2
                        eps = eps*2.0_wp
                        l = l - 1
                        if ( lr(l)<=0 ) then
                            vl(l) = vl(l+1) + vr
                            call f300()
                            cycle loop
                        else
                            vr = vl(l+1) + vr
                        endif
                    endif
                end do
            endif
        endif

        call f200()

    end do loop

    contains

        subroutine f200()
            l = l + 1
            eps = eps*0.5_wp
            if ( l<=17 ) ef = ef/sq2
            hh(l) = hh(l-1)*0.5_wp
            lr(l) = -1
            aa(l) = aa(l-1)
            q7 = q7l
            f1(l) = f(7)
            f2(l) = f(8)
            f3(l) = f(9)
            f4(l) = f(10)
            f5(l) = f(11)
            f6(l) = f(12)
            f7(l) = f(13)
            f(13) = f(7)
            f(11) = f(6)
            f(9) = f(5)
            f(7) = f(4)
            f(5) = f(3)
            f(3) = f(2)
        end subroutine f200

        subroutine f300()
            q7 = q7r(l-1)
            lr(l) = 1
            aa(l) = aa(l) + 12.0_wp*hh(l)
            f(1) = f1(l)
            f(3) = f2(l)
            f(5) = f3(l)
            f(7) = f4(l)
            f(9) = f5(l)
            f(11) = f6(l)
            f(13) = f7(l)
        end subroutine f300

    end subroutine dqnc79