program pi
  implicit none
  
  real(8) :: trap
  
  ! Serial version
  print*,'Pi=',trap(0.0d0,1.0d0,128)
  
end program pi

! Trapezoidal rule
real(8) function trap(a,b,n)
  real(8) :: a,b,x
  integer :: n,i
  real(8) :: h
  
  h=(b-a)/real(n,8)
  trap=0.5d0*(f(a)+f(b))
  
  do i=1,n-1
     x=a+real(i,8)*h
     trap=trap+f(x)
  end do
  
  trap=trap*h
  
contains
  
  ! Our function
  real(8) function f(x)
    real(8) :: x
    f=4.0d0/(1.0d0+x*x)
  end function f
  
end function trap