!**********************************************************************
!   piArea.f90 - compute pi by counting random points in circle and
!                  circumscribing square.
!                The circle has centre in the origin and radius 1
!
!
!   Each process:
!    1) receives the number of rectangles used in the approximation
!    2) generates and count points in rectangle and circle
!    3) Synchronizes for a global summation.
!   Process 0 prints the result.
!
!  Variables:
!
!    pi          the calculated result
!    n           number of random points 
!    x1, x2      extremes of each rectangle
!    SqPoints    points in each rectangle
!    CiPoints    points in circle      
!
!****************************************************************************
program main
 implicit none

 include "mpif.h"

 real(8), parameter :: PI25DT=4.0D0*ATAN(1.0D0)
! (PI25DT = 3.141592653589793238462643d0)

 real(8) ::  mypi, pi, h
 real(8) ::  sTime, eTime, elapsed
 real(8) :: x, y, x1, x2, y1, y2
 integer :: n, myid, numprocs, i, rc, ierr, ls
 integer :: SqPoints, CiPoints, my_SqPoints, my_CiPoints
 character(MPI_MAX_PROCESSOR_NAME) :: server
 

 call MPI_INIT( ierr )
 call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
 call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr )
 call MPI_Get_processor_name(server,ls,ierr)
 print *, 'Process ', myid, ' of ', numprocs, ' is alive on ',server(1:ls)

 do
    call MPI_BARRIER(MPI_COMM_WORLD,ierr)
    if ( myid .eq. 0 ) then
       write(6,"('Enter the number of points: (0 quits)') ")
       read(5,"(i10)") n
       sTime = MPI_Wtime()
    endif

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

! check for quit signal
    if ( n .le. 0 ) exit

! calculate the side size
    h = 2.0d0/numprocs
    x1 = -1.0 + myid * h
    x2 = x1 + h
    y1 = -1.0d0
    y2 = 1.0d0

    my_SqPoints  = 0
    my_CiPoints  = 0
    do i = 1, n, numprocs
       call random_number(x); x = x1 + x * (x2 - x1)
       call random_number(y); y = y1 + y * (y2 - y1)
       my_SqPoints = my_SqPoints + 1
       if ( SQRT( x*x + y*y ) <= 1.0D0 ) my_CiPoints = my_CiPoints + 1
    enddo

! collect all the partial sums
    call MPI_REDUCE(my_SqPoints,SqPoints,1,MPI_INTEGER,MPI_SUM,0, &
                    MPI_COMM_WORLD,ierr)
    call MPI_REDUCE(my_CiPoints,CiPoints,1,MPI_INTEGER,MPI_SUM,0, &
                    MPI_COMM_WORLD,ierr)

! process 0 prints the answer.
    if (myid .eq. 0) then
        pi = 4.0D0 * ( DBLE(CiPoints) / DBLE(SqPoints) )
        write(6, "(' pi is approximately: ', F18.16, &
        &       '  Error is: ', F18.16)") pi, abs(pi - PI25DT)
       eTime = MPI_Wtime()
        write(6, "(' Wall clock time = ',F6.2/)") (eTime-sTime)
    endif

 enddo

 call MPI_FINALIZE(rc)
 stop
end