Program useNRrandom

IMPLICIT NONE

INTERFACE 
   FUNCTION ranNR(idum)
     REAL ranNR
     INTEGER, INTENT(inout) :: idum
   END FUNCTION ranNR
END INTERFACE


INTEGER, PARAMETER :: K4B=selected_int_kind(9)
INTEGER(K4B) :: idum
INTEGER i, f90seedsize
INTEGER, ALLOCATABLE :: F90seed(:)
REAL rnd

idum=-293829
!  this is our seed for the NR routine 

print *, "random numbers from NR routine"
do i =1 ,10
   print *, ranNR(idum)
enddo

CALL srand(3456)
print *, "random numbers from built-in rand"
do i =1 ,10
   print *, rand()
enddo


print *, " "
print *, " Fancy seed for F90 random number generator "
CALL RANDOM_SEED(size=f90seedsize)
allocate(F90seed(f90seedsize))
CALL RANDOM_SEED(get=F90seed)
            write (*, *) F90seed
print *, "random numbers from updated built-in generator"
do i =1 ,10
   CALL RANDOM_NUMBER(rnd)
   print *, rnd
enddo



end Program useNRrandom





FUNCTION ranNR(idum)
IMPLICIT NONE
INTEGER, PARAMETER :: K4B=selected_int_kind(9) 
INTEGER(K4B), INTENT(INOUT) :: idum
REAL :: ranNR
!Minimal random number generator of Park and Miller combined with a 
!Marsaglia shift sequence. Returns a uniform random deviate between 
!0.0 and 1.0 (exclusive of the endpoint values). This fully portable, 
!scalar generator has the “traditional” (not Fortran 90) calling sequence 
!with a random deviate as the returned function value: call with idum 
!a negative integer to initialize; thereafter, do not alter idum except 
!to reinitialize. The period of this generator is about 3.1 x 10^18.

INTEGER(K4B), PARAMETER :: IA=16807,IM=2147483647,IQ=127773,IR=2836
REAL, SAVE :: am
INTEGER(K4B), SAVE :: ix=-1,iy=-1,k
if (idum <= 0 .or. iy < 0) then
   am=nearest(1.0,-1.0)/IM 
   iy=ior(ieor(888889999,abs(idum)),1) 
   ix=ieor(777755555,abs(idum)) 
   idum=abs(idum)+1
end if
ix=ieor(ix,ishft(ix,13)) 
ix=ieor(ix,ishft(ix,-17)) 
ix=ieor(ix,ishft(ix,5))
k=iy/IQ
iy=IA*(iy-k*IQ)-IR*k
if (iy < 0) iy=iy+IM 
ranNR=am*ior(iand(IM,ieor(ix,iy)),1) 
END FUNCTION ranNR