module Random

contains

  double precision function drand(ix)

    implicit none

    !     This is the random number generator of Schrage:
    !
    !     L. Schrage, A more portable Fortran random number generator, ACM
    !     Transactions on Mathematical Software 5 (1979), 132-138.

    double precision ix

    double precision a,p,b15,b16,xhi,xalo,leftlo,fhi,k
    data a/16807.d0/,b15/32768.d0/,b16/65536.d0/,p/2147483647.d0/

    xhi= ix/b16
    xhi= xhi - dmod(xhi,1.d0)
    xalo= (ix-xhi*b16)*a
    leftlo= xalo/b16
    leftlo= leftlo - dmod(leftlo,1.d0)
    fhi= xhi*a + leftlo
    k= fhi/b15
    k= k - dmod(k,1.d0)
    ix= (((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k
    if (ix.lt.0.0d0) ix= ix + p
    drand= ix*4.656612875d-10

    return

  end function drand

end module Random