program adaptive_quadrature

implicit none
integer :: fun_eval
real :: a, b, integral, tol, f, f_a, f_b, adapt_trap
external f

write(*,*) 'Input endpoints a,b & error tolerance tol:'
read*, a, b, tol

f_a=f(a)
f_b=f(b)
write(*,*) '   a               b'
integral = adapt_trap(f,a,b,f_a,f_b,tol,fun_eval)
fun_eval = fun_eval+2

write(*,*) 'Integral = ', integral
write(*,*) 'Number of function evaluations = ', fun_eval

stop
end program adaptive_quadrature
!====================================================
function f(x) result(f_result)
! Integrand.
implicit none
real, intent (in) :: x
real :: f_result

!f_result=4./(1.+x**2)
if (x==0.0) then
   f_result = 0.0
else
   f_result=sqrt(x)*log(x)
endif

return
end function f
!=====================================================================
recursive function adapt_trap(f,a,b,f_a,f_b,tol,fun_eval) result(trap)

implicit none
real, intent (in):: a, b, tol, f_a, f_b
real :: f, h, mid, f_mid
real :: T2, T3, error_T3
real :: trap, trap_left, trap_right
integer :: fun_eval, fun_eval_left, fun_eval_right
external f

write(*,*) a, b
h=b-a
mid=(a+b)/2.0
f_mid = f(mid)
T2 = h*( f_a + f_b )/2.0
T3 = h*( f_a + 2.0*f_mid + f_b )/4.0

error_T3 = (T3-T2)/3.0
if (abs(error_T3) < tol) then
   trap = T3 + error_T3
   fun_eval = 1
else
   trap_left = adapt_trap(f,a,mid,f_a,f_mid,tol/2.0,fun_eval_left)
   trap_right = adapt_trap(f,mid,b,f_mid,f_b,tol/2.0,fun_eval_right)
   trap = trap_left + trap_right
   fun_eval = fun_eval_left + fun_eval_right
endif

return
end function adapt_trap