program pi_mpi
include 'mpif.h'
integer          :: i,n,ierr,myrank,numprocs
double precision :: x,sum,pi,h,mypi

call MPI_Init(ierr)
call MPI_Comm_rank(MPI_COMM_WORLD, myrank, ierr)
call MPI_Comm_size(MPI_COMM_WORLD, numprocs, ierr)

if ( myrank == 0 ) then
   open(1, file="pi.dat")
   read(1,*) n
end if

call MPI_Bcast(n, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)

h = 1.0d0 / n
sum = 0.0d0
do i = myrank+1, n, numprocs
   x = (i - 0.5d0)*h
   sum = sum + (4.d0/(1.d0 + x*x))
end do
mypi = h * sum

call MPI_Reduce(mypi, pi, 1, MPI_DOUBLE_PRECISION, &
                MPI_SUM, 0, MPI_COMM_WORLD, ierr)

if ( myrank == 0 ) then
  write(*, fmt="(A, F16.12)") "Value of pi is ", pi
endif

call MPI_Finalize(ierr)
end program