C
      REAL FUNCTION rnd(seed)
*     linear congruential random number generator with shuffling
*     based on ran1 in "Numerical Recipes" second edition
      IMPLICIT NONE
      INTEGER a,m,q,p,n,ndiv,j,k,seed
      REAL rm,rmax
*     2^31 - 1
*     m = a*q + p
      PARAMETER (a = 16807, m = 2147483647, rm = 1.0/m)
      PARAMETER (q = 127773, p = 2836, n = 32, ndiv = 1 + (m-1)/n)
      PARAMETER (rmax = 1.0 - 1.2e-7)
      INTEGER r(n),r0,r1
      SAVE r,r0,r1
      DATA r/n*0/
*     initialize table of random numbers
      IF(seed .NE. 0) THEN
        IF(seed .LT. 0) THEN
          r1 = 584287
        ELSE
          r1 = seed
        ENDIF
        DO j = n+8,1,-1
          k = r1/q
          r1 = a*(r1-k*q) - p*k
          IF (r1 .LT. 0) r1 = r1 + m
          IF (j .LE. n) r(j) = r1
        ENDDO
        r0 = r(1)
      END IF
*     beginning when not initializing
*     compute r1 = mod(a*r1,m) without overflows
      k = r1/q
      r1 = a*(r1 - k*q) - p*k
      IF (r1 .LT. 0) r1 = r1 + m
      j = 1 + r0/ndiv
      r0 = r(j)
      r(j) = r1
      rnd = min(rm*r0,rmax)
      END