nrand Function

  real(8) function nrand(dseed_)
    integer,optional        :: dseed_
    integer                 :: dseed
    integer,parameter       :: IM1=2147483563, IM2=2147483399, IMM1=IM1-1, IA1=40014, &
         & IA2=40692, IQ1=53668, IQ2=52774, IR1=12211, IR2=3791,  &
         & NTAB=32, NDIV=1+IMM1/NTAB
    real(kind=8), parameter :: AM=1.0d0/IM1, EPS=1.2e-7, RNMX=1.-EPS
    integer                 :: dseed2, j, k, iv(NTAB), iy
    save iv, iy, dseed2
    data dseed2/123456789/, iv/NTAB*0/, iy/0/
    dseed = 123456 ;if(present(dseed_))dseed=dseed_
    if(dseed .le. 0) then
       dseed = max(-dseed,1)
       dseed2 = dseed
       do j=NTAB+8, 1, -1
          k = dseed/IQ1
          dseed = IA1*(dseed-k*IQ1)-k*IR1
          if(dseed .lt. 0) dseed = dseed+IM1
          if(j .le. NTAB) iv(j) = dseed
       enddo
       iy=iv(1)
    endif
    k = dseed/IQ1
    dseed = IA1*(dseed-k*IQ1)-k*IR1
    if(dseed .lt. 0) dseed = dseed+IM1
    k = dseed2/IQ2
    dseed2 = IA2*(dseed2-k*IQ2)-k*IR2
    if(dseed2 .lt. 0) dseed2 = dseed2+IM2
    j = 1+iy/NDIV
    iy = iv(j)-dseed2
    iv(j) = dseed
    if(iy .lt. 1) iy = iy+IMM1
    nrand = min(AM*iy,RNMX)
  end function nrand