PROGRAM Compute_PI
   IMPLICIT NONE

   INTEGER N, i, num_threads
   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(x, mypi)

   mypi = 0.0d0;

!$OMP    DO
   DO i = 0, N-1                !! Parallel Loop
     x = w * (i + 0.5d0)
     mypi = mypi + w*f(x)
   END DO
!$OMP    END DO

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

!$OMP    END PARALLEL 

   PRINT *, "Pi = ", pi

   END PROGRAM