      program compute_pi
      parameter (n = 1000000)
      integer i
      double precision w, x, sum, pi, f, y
      f(y) = 4.d0/(1.d0+y*y)
      w = 1.0d0/n
      sum = 0.0d0;
!$omp parallel do private(x) shared(w)
!$omp& reduction(+:sum)
      do i=1,n
         x = w*(i-0.5d0)
         sum = sum + f(x)
      end do
      pi = w*sum
      print *, 'pi = ', pi
      end
