Integração - Soma de Riemann

program riemann
implicit none

integer :: i,N
real :: integ,a,b,dx,x

! zerando as variaveis
integ = 0.
x=0.
dx=0.

print *,"Digite os limites a,b e o numero N de intervalos:"
print *,"Use N grande, ex: 1000"
read *,a,b,N

dx = (b-a)/real(N)

!integ = 0.5*(f(a)+f(b))

do i = 1,N
    x = a + i*dx
    integ = integ + f(x)*dx
enddo

print *,'Integral = ', integ

stop '>>> Fim do Programa'

 contains
real function f(x)
real :: x
! f(x) = x^2
 f = x**2.
return
end function f

end program riemann