*******************************************************************
      PROGRAM XSURAND3D
*
*     This program is written to test
*     the SURAND random number generator
*     in 3D (See Chapter 11:Monte Carlo).
*
      integer i,j
      double precision ranu
      double precision random_low, random_med, random_high
      open(unit=10,file='output_surand3.dat',status='new')
      random_low=ranu()
      random_med=ranu()
      do 20 j=1,10
      do 10 i=1,500000000
      random_high=ranu()
      if ((random_low .le. 0.01) .and. (random_med.le. 0.01) .and.
     > (random_high.le. 0.01)) then
         write(10,*), random_low, random_med, random_high
      endif
      random_low=random_med
      random_med=random_high
 10   continue
 20   continue
*
*     Now view the 3D plot in 'output_surand3.dat'
*     using Matlab.
*
      end

*******************************************************************
      double precision function ranu ()
*
* This function is given in Chapter 11: Monte Carlo.
* Based on good implementation of SURAND. See
* Park & Miller, reference given in Chapter 11.
* 
      integer a,m,q,r,seed
      double precision rm
      parameter (a=16807, m=2147483647, q=127773, r=2836, rm=1d0/m)
      common /random/ seed
      save /random/
      data seed /1/
      seed = a*mod(seed,q)-r*(seed/q)
      if (seed .le. 0) seed=seed+m
      ranu=seed*rm
      return
      end
*******************************************************************