program main
      use MPI
      integer n, myid, numprocs, i, ierr
      double precision  PI25DT
      parameter (PI25DT = 3.141592653589793238462643d0)
      double precision  mypi, pi, h, sum, x, f, a
      double precision startwtime, endwtime
      integer namelen
      character*(MPI_MAX_PROCESSOR_NAME) name

c                                 function to integrate
      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 )
      call MPI_GET_PROCESSOR_NAME(name, namelen, ierr)
      print "('MPI Task ', I3,' on ',A)", myid, trim(name)
      call MPI_BARRIER(MPI_COMM_WORLD, ierr)
 10   if ( myid .eq. 0 ) then
         write(6,98)
 98      format('Enter the number of intervals: (0 quits)')
         read(5,99) n
 99      format(i10)
         startwtime=MPI_WTIME()
      endif
      call MPI_BCAST( n, 1, MPI_INTEGER, 0,
     +                MPI_COMM_WORLD, ierr)
c                             check for quit signal
      if ( n .le. 0 ) goto 30
c                    calculate the interval size
      h = 1.0d0/n
      sum  = 0.0d0
      do 20 i = myid+1, n, numprocs
        x   = h * (dble(i) - 0.5d0)
        sum = sum + f(x)
 20   continue
      mypi = h * sum
c                   collect all the partial sums
      call MPI_REDUCE( mypi, pi, 1, MPI_DOUBLE_PRECISION,
     +                MPI_SUM, 0, MPI_COMM_WORLD,ierr)
c                              node 0 prints the answer
      if (myid .eq. 0) then
         endwtime = MPI_WTIME()
         write(6, 96) pi, abs(pi - PI25DT)
 96      format('pi is approximately: ', F18.16,
     +          ', Error is: ', F18.16)
         write(6, 97) endwtime-startwtime
 97      format('wall clock time = ', F18.16)
      endif
      goto 10
 30   call MPI_FINALIZE(ierr)
      end