For those interested in this old thread, I just posted a code in this thread Spacing of the results of Fortran intrinsic random_number - #33 by RonShepard that uses this approach to return values from the full set of normal floating point numbers in [0.0,1.0). In the last version of the posted code, the random number codes are pure subroutines, so I think that code is thread-safe, but that feature is being discussed in that other thread.