program hybpi
  use omp_lib
  implicit none
  include "mpif.h"  

  integer    maxcpu
  parameter (maxcpu = 10)

  real*8  mypi,pi,h,pil,x,f,a,t0,t1,t2,PI25DT
  integer istat(MPI_STATUS_SIZE),ierr
  integer myid,numprocs,j
  integer(kind=8), parameter:: n=80000000000
  integer(kind=8)  i
  integer(kind=8)  ibeg(0:maxcpu),iend(0:maxcpu)

  parameter  (PI25DT = 3.141592653589793238462643d0)

  f(a) = 4.d0 / (1.d0 + a*a)

  call MPI_INIT(ierr)
  call MPI_COMM_RANK(MPI_COMM_WORLD,myid,ierr)
  call MPI_COMM_SIZE(MPI_COMM_WORLD,numprocs,ierr)

  h = 1.0d0/dble(n)
  do j=0,numprocs-1
     ibeg(j)=int(j*n/numprocs,8)
     iend(j)=int((j+1)*n/numprocs-1,8)
  enddo

  write(*,*) 'Process ', myid, 'of ', numprocs, ' is alive'
  write(*,*) '              ib','                 ie' 
  write(*,*) ibeg(myid),iend(myid)
  call MPI_BARRIER(MPI_COMM_WORLD,ierr)

  t0  = MPI_WTIME()
  pil = 0.0d0
  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,x) SCHEDULE(STATIC) REDUCTION(+:pil)
  do i=ibeg(myid),iend(myid)
     x   = h * (dble(i) + 0.5d0)
     pil = pil + f(x)
  enddo
  !$OMP END PARALLEL DO
  mypi = h * pil
      
  call MPI_REDUCE(mypi,pi,1,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,ierr)

  t1  = MPI_WTIME()
  t2  = t1-t0
  if(myid.eq.0) write(*,*)'elapsed time: ',t2

  if (myid .eq. 0) then
     write(*,*) 'pi is approximately:', pi, 'Error is:',abs(pi - PI25DT)
  endif

  call MPI_FINALIZE(ierr)
end program hybpi
  
