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

.