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