      DOUBLE PRECISION FUNCTION URAND (IY)
      IMPLICIT NONE
C
C UNIFORM RANDOM NUMBER GENERATOR BASED ON THEORY AND SUGGESTIONS GIVEN IN
C  D.E. KNUTH (1969), VOL 2.  THE INTEGER IY SHOULD BE INITIALIZED TO AN
C  ARBITRARY INTEGER PRIOR TO THE FIRST CALL TO URAND.  THE CALLING PROGRAM
C  SHOULD NOT ALTER THE VALUE OF IY BETWEEN SUBSEQUENT CALLS TO URAND.
C  VALUES OF URAND WILL BE RETURNED IN THE INTERVAL (0,1).
C
      INTEGER IY
C
      INTEGER IA, IC, ITWO, M2, M, MIC
      DOUBLE PRECISION HALFM, S
C
      DOUBLE PRECISION DATAN, DSQRT, DFLOAT
C
      DATA M2 / 0 /, ITWO / 2 /
C
C
      IF (M2.NE.0) GO TO 20
C
C ***      IF FIRST ENTRY, COMPUTE MACHINE INTEGER WORD LENGTH
      M = 1
10    CONTINUE
      M2 = M
      M = ITWO*M2
      IF (M.GT.M2) GO TO 10
      HALFM = M2
C
C ***      COMPUTE MULTIPLIER AND INCREMENT FOR LINEAR CONGRUENTIAL METHOD
      IA = 8*IDINT (HALFM*DATAN (1.D0)/8.D0)+5
      IC = 2*IDINT (HALFM*(0.5D0-DSQRT (3.D0)/6.D0))+1
      MIC = (M2-IC)+M2
C
C ***      S IS THE SCALE FACTOR FOR CONVERTING TO FLOATING POINT
      S = 0.5D0/HALFM
C
C ***      COMPUTE NEXT RANDOM NUMBER
20    CONTINUE
      IY = IY*IA
C
C ***      THE FOLLOWING STATEMENT IS FOR COMPUTERS WHICH DO NOT ALLOW
C           INTEGER OVERFLOW ON ADDITION
      IF (IY.GT.MIC) IY = (IY-M2)-M2
      IY = IY+IC
C
C ***      THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE THE WORD LENGTH
C           FOR ADDITION IS GREATER THAN FOR MULTIPLICATION
      IF (IY/2.GT.M2) IY = (IY-M2)-M2
C
C ***      THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE INTEGER OVERFLOW
C           AFFECTS THE SIGN BIT
      IF (IY.LT.0) IY = (IY+M2)+M2
      URAND = DFLOAT (IY)*S
C
      RETURN
      END
