PROGRAM Compute_PI
   IMPLICIT NONE

   INTEGER N, i, num_threads, OMP_GET_NUM_THREADS, OMP_GET_THREAD_NUM
   DOUBLE PRECISION w, x, sum, pi, mypi, 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 = 50000000		!! Number of intervals
   w = 1.0d0/N  	!! width of each interval

   sum = 0.0d0

!$OMP    PARALLEL PRIVATE(i, num_threads, x, mypi)

   num_threads = OMP_GET_NUM_THREADS()
   mypi = 0.0d0;

   DO i = OMP_GET_THREAD_NUM(), N-1, num_threads
     x = w * (i + 0.5d0)
     mypi = mypi + w*f(x)
   END DO

!$OMP CRITICAL
   pi = pi + mypi
!$OMP END CRITICAL

!$OMP    END PARALLEL 

   PRINT *, "Pi = ", pi

   END PROGRAM