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

  use kinds
  implicit none

  include 'mpif.h'

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

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

  call mpi_init(ierr)
  call mpi_comm_size(mpi_comm_world, ntasks, ierr)
  call mpi_comm_rank(mpi_comm_world, task_id, ierr)

  if (task_id == 0) then
     print *, 'Enter number of intervals: '
     read *, n
  endif
      
  call mpi_bcast(n, 1, mpi_integer, 0, mpi_comm_world, ierr)
  
  ! calculate the integral
  h = 1.0d0 / n
  localsum  = 0.0d0
  
  do i = task_id + 1, n, ntasks
     x = h * (dble(i) - 0.5d0)
     localsum = localsum + f(x)
  enddo

  call mpi_reduce(localsum, sum, 1, mpi_double_precision, &
       mpi_sum, 0, mpi_comm_world, ierr)
  
  ! print the result
  if (task_id .eq. 0) then
     pi = h * sum
     print *, 'pi computes to: ', pi
     print *, 'and is:         ', pi25dt
  endif
  
  call mpi_finalize(ierr)
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