!  GeNpois.f90
!
!  Poisson distributed random number genarator,
!  adapted from npois written by PG.
!
!  Detlef Irmscher, Thomas Ullrich, Uni Heidelberg, September 1993
!  Last update: 16.3.95 tu
!
integer function GeNpois (amean)
   use f90_kind
   implicit none
   real(kind=double), intent(in)   :: amean
   real(kind=double)               :: rand, p, psum
   integer                         :: nmax, n

   call random_number(rand)
   p    = exp(-amean)
   psum = p
   nmax = max(amean*10., 1._double)

   do n = 1, nmax
      if (psum > rand) then
          GeNpois = n - 1
          return
      endif
      p    = p*amean/n
      psum = psum + p
   enddo

   GeNpois = nmax
end function GeNpois