dlobatto Subroutine

    recursive subroutine dlobatto (f, a, b, error_tol, ans, ierr)

    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 !! relative error tolerance
    real(wp),intent(out) :: ans !! computed value of integral
    integer,intent(out)  :: ierr !! status code:
                                 !!
                                 !! * 1 = success
                                 !! * 2 = requested accuracy may not be satisfied

    real(wp) :: m,h,s,erri1,erri2,is,tol,fa,fb,i1,i2,r
    real(wp),dimension(13) :: x,y
    integer :: i
    integer :: k !! number of calls to the recursive function

    integer,parameter :: kmax = 10000 !! maximum number of calls to the recursive function (probably should be an input)
    real(wp),parameter :: eps  = epsilon(1.0_wp)
    real(wp),parameter :: alpha = sqrt(2.0_wp/3.0_wp)
    real(wp),parameter :: beta  = 1.0_wp/sqrt(5.0_wp)
    real(wp),parameter :: x1   = .94288241569547971905635175843185720232_wp
    real(wp),parameter :: x2   = .64185334234578130578123554132903188354_wp
    real(wp),parameter :: x3   = .23638319966214988028222377349205292599_wp
    real(wp),dimension(7) :: c = [ .015827191973480183087169986733305510591_wp, &
                                   .094273840218850045531282505077108171960_wp, &
                                   .15507198733658539625363597980210298680_wp, &
                                   .18882157396018245442000533937297167125_wp, &
                                   .19977340522685852679206802206648840246_wp, &
                                   .22492646533333952701601768799639508076_wp, &
                                   .24261107190140773379964095790325635233_wp ]
    k = 0
    ierr = 1
    tol = max(eps, error_tol)
    m = (a+b)/2.0_wp
    h = (b-a)/2.0_wp

    x = [a, m-x1*h, m-alpha*h, m-x2*h, m-beta*h, m-x3*h, m, m+x3*h, m+beta*h, m+x2*h, m+alpha*h, m+x1*h, b]
    do i=1,13
        y(i) = f(x(i))
    end do

    fa=y(1)
    fb=y(13)
    i2=(h/6.0_wp)*(y(1)+y(13)+5.0_wp*(y(5)+y(9)))
    i1=(h/1470.0_wp)*(77.0_wp*(y(1)+y(13))+432.0_wp*(y(3)+y(11))+625.0_wp*(y(5)+y(9))+672.0_wp*y(7))

    is = h*(c(1)*(y(1)+y(13)) + &
            c(2)*(y(2)+y(12)) + &
            c(3)*(y(3)+y(11)) + &
            c(4)*(y(4)+y(10)) + &
            c(5)*(y(5)+y(9))  + &
            c(6)*(y(6)+y(8))  + &
            c(7)*y(7))

    s = sign(1.0_wp,is)
    if (s==0.0_wp) s = 1.0_wp
    erri1 = abs(i1-is)
    erri2 = abs(i2-is)
    r = 1.0_wp
    if (erri2/=0.0_wp) r=erri1/erri2
    if (r>0.0_wp .and. r<1.0_wp) tol=tol/r
    is=s*abs(is)*tol/eps
    if (is==0.0_wp) is=b-a

    call adaptive_lobatto_step(a,b,fa,fb,is,ans)

    contains

    recursive subroutine adaptive_lobatto_step(a,b,fa,fb,is,ans)

        !!  Recursive function used by adaptive_lobatto.
        !!  Tries to approximate the integral of f(x) from a to b
        !!  to an appropriate relative error.

        implicit none

        real(wp),intent(in)   :: a
        real(wp),intent(in)   :: b
        real(wp),intent(in)   :: fa
        real(wp),intent(in)   :: fb
        real(wp),intent(in)   :: is
        real(wp),intent(out)  :: ans

        real(wp) :: h,m,mll,ml,mr,mrr,fmll,fml,fm,fmr,fmrr,i2,i1
        real(wp),dimension(6) :: q

        k = k + 1
        if (k>kmax) then
            ierr = 2
            ans = 0.0_wp
            return
        end if
        h   = (b-a)/2.0_wp
        m   = (a+b)/2.0_wp
        mll = m-alpha*h
        ml  = m-beta*h
        mr  = m+beta*h
        mrr = m+alpha*h

        fmll = f(mll)
        fml  = f(ml)
        fm   = f(m)
        fmr  = f(mr)
        fmrr = f(mrr)

        i2 = (h/6.0_wp)*(fa+fb+5.0_wp*(fml+fmr))
        i1 = (h/1470.0_wp)*(77.0_wp*(fa+fb)+432.0_wp*(fmll+fmrr)+625.0_wp*(fml+fmr)+672.0_wp*fm)

        if ( (is+(i1-i2)==is) .or. (mll<=a) .or. (b<=mrr) ) then

            if (((m <= a) .or. (b<=m)) .and. (ierr==1)) ierr = 2
            ans = i1

        else

            if (ierr==1) call adaptive_lobatto_step(a,mll,fa,fmll,    is,q(1))
            if (ierr==1) call adaptive_lobatto_step(mll,ml,fmll,fml,  is,q(2))
            if (ierr==1) call adaptive_lobatto_step(ml,m,fml,fm,      is,q(3))
            if (ierr==1) call adaptive_lobatto_step(m,mr,fm,fmr,      is,q(4))
            if (ierr==1) call adaptive_lobatto_step(mr,mrr,fmr,fmrr,  is,q(5))
            if (ierr==1) call adaptive_lobatto_step(mrr,b,fmrr,fb,    is,q(6))

            if (ierr==1) then
                ans = sum(q)
            else
                ans = i1
            end if

        end if

    end subroutine adaptive_lobatto_step
!**************************************************************

    end subroutine dlobatto