dsimpson Subroutine

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

    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 !! 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) :: bma,is,tol,fa,fm,fb
    real(wp),dimension(5) :: yy
    integer :: k !! number of calls to the recursive function

    real(wp),parameter :: eps = epsilon(1.0_wp)
    real(wp),dimension(5),parameter :: c = [.9501_wp, .2311_wp, .6068_wp, .4860_wp, .8913_wp]
    integer,parameter :: kmax = 10000 !! maximum number of calls to the recursive function (probably should be an input)

    k = 0
    ierr = 1
    bma = b-a
    tol = max(eps, error_tol)

    fa    = f(a)
    fm    = f((a+b)/2.0_wp)
    fb    = f(b)
    yy(1) = f(a+c(1)*bma )
    yy(2) = f(a+c(2)*bma )
    yy(3) = f(a+c(3)*bma )
    yy(4) = f(a+c(4)*bma )
    yy(5) = f(a+c(5)*bma )

    is = bma/8.0_wp * (fa+fm+fb+sum(yy))
    if (is==0.0_wp) is = bma
    is = is*tol/eps

    call adaptive_simpson_step(a,b,fa,fm,fb,is,ans)

    contains

    recursive subroutine adaptive_simpson_step (a,b,fa,fm,fb,is,ans)
        !!  Recursive function used by adaptive_simpson.
        !!  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)   :: fm
        real(wp),intent(in)   :: fb
        real(wp),intent(in)   :: is
        real(wp),intent(out)  :: ans

        real(wp) :: m,h,fml,fmr,i1,i2,q1,q2

        k = k + 1
        if (k>kmax) then
            ierr = 2
            ans = 0.0_wp
            return
        end if
        m = (a + b)/2.0_wp
        h = (b - a)/4.0_wp
        fml = f(a + h)
        fmr = f(b - h)
        i1 = h/1.5_wp * (fa + 4.0_wp*fm + fb)
        i2 = h/3.0_wp * (fa + 4.0_wp*(fml + fmr) + 2.0_wp*fm + fb)
        i1 = (16.0_wp*i2 - i1)/15.0_wp

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

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

        else

            if (ierr==1) call adaptive_simpson_step (a,m,fa,fml,fm,is,q1)
            if (ierr==1) call adaptive_simpson_step (m,b,fm,fmr,fb,is,q2)

            if (ierr==1) then
                ans = q1 + q2
            else
                ans = i1
            end if

        end if

        end subroutine adaptive_simpson_step
    !**************************************************************

    end subroutine dsimpson