program compute_pi 
      integer n, i 
      double precision w, x, sum, pi, f, a 
c function to integrate
      f(a) = 4.d0 / (1.d0 + a*a) 
      print *, 'Enter number of intervals: ' 
      read (5,*) n 
c calculate the interval size
      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) 
      enddo 
      pi = w * sum 
      print *, 'computed pi = ', pi 
      print *, 'exact pi =    ', 
     x    '   3.14159_26535_89793_23846_26433_83279_50288_41971_...'
      end