FUNCTION RANDOM(ISEED)
IMPLICIT NONE
    INTEGER(4), PARAMETER:: K=16087, L=0, NB=31
    REAL(8), PARAMETER:: F=2.0D0**(-NB)
    INTEGER(4), INTENT(INOUT) :: ISEED
    REAL(8) :: RANDOM
    ISEED=ISEED*K+L
    ISEED=IBCLR(ISEED,NB)
    RANDOM=F*ISEED
END FUNCTION RANDOM