**************************************************************
* This program performs a simple histogram analysis of ten million calls
* to RAN1.  If RAN1 is correctly yielding a uniformly distributed value
* between 0.0 and slightly less than 1.0, then the counts in each bin
* should be approximately 100,000, plus or minus a few tenths of a
* percent.  For the last bin on the last line of output, the count
* should be zero.  If you see anything wrong with your output, please 
* let me know and I'll try to help figure out why this doesn't work on 
* your system.
*
* Once you have confirmed that RAN1 appears to be working correctly on
* your system, then you can comment out or delete PROGRAM TEST and copy
* RAN1 directly into the source code for your Monte Carlo program.

      program test

      integer icount(101)
      data icount/101*0/
      real x, ran1
      external ran1
      integer idum

* initialize with negative number
      idum = -1 
      x = ran1(idum)

* tabulate values returned by many calls
      do i=1,10000000
         x = ran1(idum)
         ibin = (x*100) + 1
         icount(ibin) = icount(ibin) + 1
      enddo
      write(*,'(10i7)')icount
      end

***************************************************************
*    "Minimal" random number generator of Park and Miller with
*    Bays-Durham shuffle and added safeguards.  Returns a uniform random
*    deviate between 0.0 and 1.0 (exclusive of the endpoint values).
*    Call with IDUM a negative integer to initialize; thereafter, do not
*    alter IDUM between successive deviates in a sequence.  RNMX should
*    approximate the largest floating value that is less than 1.

      function ran1(idum)
      integer idum,ia,im,iq,ir,ntab,ndiv
      real ran1,am,eps,rnmx
      parameter (ia=16807,im=2147483647,am=1./im,iq=127773,ir=2836,
     &     ntab=32,ndiv=1+(im-1)/ntab,eps=1.2e-7,rnmx=1.-eps)
      integer j,k,iv(ntab),iy
      save iv,iy
      data iv /ntab*0/, iy /0/
      if (idum.le.0.or.iy.eq.0) then
         idum=max(-idum,1)
         do 11 j=ntab+8,1,-1
            k=idum/iq
            idum=ia*(idum-k*iq)-ir*k
            if (idum.lt.0) idum=idum+im
            if (j.le.ntab) iv(j)=idum
 11      continue
         iy=iv(1)
      endif
      k=idum/iq
      idum=ia*(idum-k*iq)-ir*k
      if (idum.lt.0) idum=idum+im
      j=1+iy/ndiv
      iy=iv(j)
      iv(j)=idum
      ran1=min(am*iy,rnmx)
      return
      end
*
*     After Press et al., Numerical Recipes for Fortran