module gauss
contains
  function f(x)
    implicit none
    real(kind(1d0)) :: f
    real(kind(1d0)), intent(in) :: x
    f=exp(-x*x)
  end function
end module

program int

  use mpi

! For earlier version of NAG,
! use nag_f77_d_chapter
! For current version of NAG
  use nag_library, only: d01ahf
  use gauss

  implicit none
  integer :: nproc, rank, ierr, ifail, npts
  real (kind(1d0)) :: xmin,xmax,tol,total,res,relerr

  call mpi_init(ierr)

  call mpi_comm_size(mpi_comm_world, nproc, ierr)
  call mpi_comm_rank(mpi_comm_world, rank, ierr)

  xmin=rank*(10d0/nproc)
  xmax=(rank+1)*(10d0/nproc)
  tol=1d-4
  relerr=0
  ifail=0

  total=D01AHF(xmin, xmax, tol, npts, relerr, f, 10000, ifail)

  call MPI_Reduce(total,res,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD)

  if (rank.eq.0) then
    write(*,*)"Integral is ",res
  endif

  call mpi_finalize(ierr)
end program
