module dpres
   integer, parameter :: dp = kind(1.d0)
end module dpres

subroutine fct(x,f)
   use dpres
   real(kind=dp) :: x,f
   f = 4.D0/(1.D0+x*x)
end subroutine fct

program calcul_pi
   
   use dpres
   implicit none
   include 'mpif.h'

   ! Variables globales
   integer :: N=1000000,i
   real(kind=dp) :: x0 = 0.d0,xN = 1.d0
   real(kind=dp) :: dx
   real(kind=dp) :: pi = 0.d0

   ! Variables mpi
   integer :: erreur, totalCPU, monCPU, NP, reste
   integer :: source,dest,flag
   real(kind=dp) :: x0p , int_partielle = 0.d0
   real(kind=dp) :: xg,xd,fg,fd
   integer :: statut(MPI_STATUS_SIZE)

   ! Demmarrage parallele
   call MPI_INIT(erreur)
   call MPI_COMM_SIZE(MPI_COMM_WORLD, totalCPU, erreur)
   call MPI_COMM_RANK(MPI_COMM_WORLD, monCPU, erreur)

   dx = (xN-x0)/dfloat(N)

   ! Decoupage parallele
   NP = N/totalCPU
   ! Reste des points a traiter
   reste = mod(N,totalCPU)

   x0p = x0+dfloat(NP*monCPU)*dx
   ! Rerepartition des points

   if (monCPU<reste) then
      NP = NP+1
      x0p = x0p+dfloat(monCPU)*dx
   else
      x0p = x0p+dfloat(reste)*dx
   endif

   xg = x0p
   xd = x0p+dx

   do i=0,NP-1
      call fct(xg,fg)
      call fct(xd,fd)
      int_partielle = int_partielle+0.5d0*(fg+fd)*dx
      xg = xg+dx
      xd = xd+dx
   enddo

   ! Affichage des resultats de chaque CPU
   print*, "ID =",monCPU,"integrale=",int_partielle

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! COMMUNICATION
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   if (monCPU>0) then
      dest = 0
      flag = 100+monCPU
      call MPI_SEND(int_partielle,1,MPI_DOUBLE_PRECISION,dest,flag,MPI_COMM_WORLD,erreur)
   else
      pi = int_partielle
      do i=1,totalCPU-1
         source = i
         flag = 100+i
         call MPI_RECV(int_partielle,1,MPI_DOUBLE_PRECISION,source,flag,MPI_COMM_WORLD,statut,erreur)
         pi = pi+int_partielle
      enddo
      print*, "Pi =",pi
   endif 


   call MPI_FINALIZE(erreur)


end program calcul_pi
