!**********************************************************************
!   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

!! MPI Fortran include file

 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

 
!! MPI initialization: define numprocs and myid
 myid = 0; numprocs = 1;

 print *, 'Process ', myid, ' of ', numprocs, ' is alive '

 do
    if ( myid .eq. 0 ) then
       write(6,"('Enter the number of points: (0 quits)') ")
       read(5,"(i10)") n
    endif

    !! Send number of points to all processes

! 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
	
    !! Sum up SqPoints from all processes: result in SqPoints
    !! Sum up CiPoints from all processes: result in CiPoints
    SqPoints = my_SqPoints;
    CiPoints = my_CiPoints;

! 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)
    endif

 enddo
    !! MPI finalization    
 stop
end