program mpi_square
use mpi_f08
implicit none

integer, parameter :: count=100000
integer :: nproc,rank,i,hit=0
integer, allocatable :: seed(:)
real(kind(1d0)) :: x(2),result,t

call mpi_init()
call mpi_comm_size(mpi_comm_world, nproc)
call mpi_comm_rank(mpi_comm_world, rank)

! If one needs to set the seeds explicitly and differently on each rank
! uncomment this section
!
! call random_seed(size=i)
! allocate(seed(i))
! i=rank+1
! call random_seed(put=seed)

do i=1,count
  call random_number(x)
  if ((x(1)**2+x(2)**2)<1d0) hit=hit+1
enddo

t=hit/real(count,kind(1d0))
write(*,101) rank, t

call mpi_reduce(t,result,1,mpi_double_precision,mpi_sum,0,mpi_comm_world)

if (rank.eq.0) write(*,102) result/nproc


101 format ('On rank ', i3,' ratio: ',f13.9)
102 format('Final ratio: ',f13.9)

call mpi_finalize()
end program
