double precision function urand(iseed)
      integer iseed
      double precision dl
      dl=dmod(16807.0d0*dble(iseed),2147483647.0d0)
      iseed=idint(dl)
      urand=dl*4.6566128752458d-10
      return
      end