! udrnrt.f90  uniformly distributed random numbers
double precision function udrnrt()
  double precision :: modulus
  double precision :: multiplier
  double precision, save :: state
  double precision :: t
  integer :: i

  data modulus/2147483647.0d0/ ! 2**31-1 
  data multiplier/16807.0d0/   ! 7**5
  data state/30000000.0d0/     ! arbitrary initial value

  t = state * multiplier
  i = t / modulus
  state = t - (i * modulus)
  if(state .lt. 0.0d0) then              
    state = state + modulus
  end if
  udrnrt = state / modulus
end function udrnrt