V3FIT
ran3.f
1  SUBROUTINE ran3(idum,rand)
2 c#######################################################################
3 c
4 c Returns a uniform random deviate between 0.0 and 1.0. Set idum to
5 c ANY negative value to initialize or reinitialize the sequence.
6 c This FUNCTION is taken from W.H. Press, "Numerical Recipes" p. 199.
7 c
8  USE stel_kinds
9  IMPLICIT NONE
10  INTEGER :: idum, IFf, inext, inextp, i, ii, k
11  REAL(rprec) :: rand, ma, mj, mk
12  SAVE
13  REAL(rprec), PARAMETER :: mbig=4000000, mseed=1618033,
14  1 mz=0, fac=1._dp/mbig
15 c PARAMETER (mbig=1000000000,mseed=161803398,mz=0,fac=1./mbig)
16 c
17 c According to Knuth, ANY large mbig, and ANY smaller (but still large)
18 c mseed can be substituted for the above values.
19  dimension ma(55)
20  data iff /0/
21  IF (idum.lt.0 .or. iff.eq.0) THEN
22  iff=1
23  mj=mseed-abs(idum)
24  mj=mod(mj,mbig)
25  ma(55)=mj
26  mk=1
27  DO 11 i=1,54
28  ii=mod(21*i,55)
29  ma(ii)=mk
30  mk=mj-mk
31  IF(mk.lt.mz) mk=mk+mbig
32  mj=ma(ii)
33  11 CONTINUE
34  DO 13 k=1,4
35  DO 12 i=1,55
36  ma(i)=ma(i)-ma(1+mod(i+30,55))
37  IF(ma(i).lt.mz) ma(i)=ma(i)+mbig
38  12 CONTINUE
39  13 CONTINUE
40  inext=0
41  inextp=31
42  idum=1
43  ENDif
44  inext=inext+1
45  IF(inext.eq.56) inext=1
46  inextp=inextp+1
47  IF(inextp.eq.56) inextp=1
48  mj=ma(inext)-ma(inextp)
49  IF(mj.lt.mz) mj=mj+mbig
50  ma(inext)=mj
51  rand=mj*fac
52 
53  END SUBROUTINE ran3