*******************************************************************
      PROGRAM XSURAND2D
*
*     This program is written to test
*     the SURAND random number generator
*     in 2D (See Chapter 11:Monte Carlo).
*
      integer i,j
      double precision ranu
      double precision random, random_new
      open(unit=10,file='output_surand2.dat',status='new')
      random=ranu()
      do 20 j=1,10
      do 10 i=1,500000000
      random_new=ranu()
      if ((random_new .le. 0.001) .and. (random .le. 0.001)) then
         write(10,*), random, random_new
      endif
      random=random_new
 10   continue
 20   continue
*
*     Now view the 2D plot in 'output_surand2.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
*******************************************************************