PROGRAM Compute_PI
   INTEGER n, i
   DOUBLE PRECISION w, x, sum, pi, f, a

! Integrate(1/SQRT(1.d0 - x*x)) = arcsin(x)
! arcsin(1) = Pi/2
!
! Function to integrate

   f(a) = 2.d0 / SQRT(1.d0 - a*a)

   n = 100000000	!! Number of intervals
   w = 1.0d0/n  	!! width of each interval

   sum = 0.0d0

!$OMP    PARALLEL DO   PRIVATE(x), SHARED(w), REDUCTION(+ : sum)

   DO i = 1, n
     x = w * (i - 0.5d0)
     sum = sum + f(x)
   END DO

!$OMP    END PARALLEL DO

   pi = w * sum
   PRINT *, "Pi = ", pi

   END PROGRAM