program pidemo
  ! compute pi by integrating f(x) = 4/(1 + x*x)     
  ! example program for colloq
  ! scalar version

  use kinds
  implicit none

  real(kind=REAL8), parameter   :: PI25DT = 3.141592653589793238462643d0

  real(kind=REAL8)   :: h, x, sum, pi
  real(kind=REAL8)   :: f
  integer(kind=INT4) :: n, i


  print *, 'Enter number of intervals: '
  read *, n
      
  ! calculate the integral
  h = 1.0d0 / n
  sum  = 0.0d0

  !$OMP PARALLEL DO PRIVATE(x) REDUCTION(+: sum)
  do i = 1, n
     x = h * (dble(i) - 0.5d0)
     sum = sum + f(x)
  enddo
  !$OMP END PARALLEL DO

  pi = h * sum
  
  !     print the result
  print *, 'pi computes to: ', pi
  print *, 'and is:         ', PI25DT

end program pidemo


function f(x)
  use kinds
  implicit none

  real(kind=REAL8)   :: f, x
  f = 4.d0 / (1.d0 + x*x)
end function f