program calculate_pi
! This program calculates pi with a Monte Carlo method i.e.
! integrate f(x) = 4/(1+x**2) on the range [0,1] by approximating 
! with a sum of f(x), with N random x values on the interval, 
! divided by N. 
  implicit none
  integer, parameter :: wp = SELECTED_REAL_KIND(15,307)
  integer, parameter :: n = 1000000
  real (kind=wp), dimension (n) :: x
  real (kind=wp) :: pi

!Use intrinsic subroutine "RANDOM_NUMBER" and instrinsic function "SUM"
  call RANDOM_NUMBER(x) ! x vector of psuedorandom numbers in [0,1[

  pi = 4.0_wp*SUM(1.0_wp/(1.0_wp+x**2))/REAL( n,wp )

! Print the result and compare it with the "exact" value (4arctan(1)) 
! using the intrinsic function ATAN
  print *, pi, pi-4.0_wp*ATAN(1.0_wp)

end program calculate_pi