model/src/port_rand.F

C $Header: /u/gcmpack/MITgcm/model/src/port_rand.F,v 1.6 2005/03/01 16:51:27 jmc Exp $
C $Name:  $

#undef _USE_INTEGERS

C--  File port_rand.F: Portable random number generator functions
C--   Contents
C--   o PORT_RAND
C--   o PORT_RANARR
C--   o PORT_RAND_NORM

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: PORT_RAND
C     !INTERFACE:
      real*8 FUNCTION PORT_RAND(seed)

C     !DESCRIPTION:
C     Portable random number generator
C      seed >=0 :: initialise using this seed ; and return 0
C      seed < 0 :: if first call then initialise using the default seed (=mseed)
C                  and always return a random number

C     !USES:
      IMPLICIT NONE

C     !INPUT PARAMETERS:
#ifdef _USE_INTEGERS
      INTEGER seed
#else
      real*8  seed
#endif
CEOP

C     !LOCAL VARIABLES:
      INTEGER nff,idum
      PARAMETER(nff=55)
      PARAMETER(idum=-2)
      real*8 fac
#ifdef _USE_INTEGERS
      INTEGER mbig,mseed,mZ
      PARAMETER (mbig=1000000000,mz=0,fac=1.d0/mbig)
      INTEGER mj,mk,ma(nff)
      DATA mseed/161803398/
#else
      real*8 mbig,mseed,mz
      PARAMETER (mbig=4000000.,mz=0.,fac=1.d0/mbig)
      real*8 mj,mk,ma(nff)
      DATA mseed/1618033./
#endif
      LOGICAL firstCall
      INTEGER i,ii,inext,inextp,k
      DATA firstCall /.true./
      SAVE firstCall,inext,inextp,ma

C-    Initialise the random number generator
      IF (firstCall .OR. seed.GE.mz) THEN
        IF (seed.GE.mz) mseed = seed
        firstCall=.false.
        mj=mseed-iabs(idum)
        mj=mod(mj,mbig)
        ma(nff)=mj
        mk=1
        DO i=1,nff-1
          ii=mod(21*i,nff)
          ma(ii)=mk
          mk=mj-mk
          IF (mk.LT.mz) mk=mk+mbig
          mj=ma(ii)
        ENDDO
        DO k=1,4
          DO i=1,nff
            ma(i)=ma(i)-ma(1+mod(i+30,nff))
            IF (ma(i).LT.mz) ma(i)=ma(i)+mbig
          ENDDO
        ENDDO
        inext=0
        inextp=31
      ENDIF

C-    Compute a random number (only if seed < 0)
      IF (seed.GE.mz) THEN
        port_rand=0.d0
      ELSE
        inext=mod(inext,nff)+1
        inextp=mod(inextp,nff)+1
        mj=ma(inext)-ma(inextp)
        IF (mj.LT.mz) mj=mj+mbig
        ma(inext)=mj
        port_rand=mj*fac
      ENDIF

      RETURN
      END

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP
C     !ROUTINE: PORT_RANARR
C     !INTERFACE:
      SUBROUTINE PORT_RANARR(n,arr)

C     !DESCRIPTION:
C     Portable random number array generator

C     !USES:
      IMPLICIT NONE

C     !INPUT PARAMETERS:
      INTEGER n
      real*8 arr(n)
CEOP

C     !LOCAL VARIABLES:
      INTEGER i
      real*8 port_rand
#ifdef _USE_INTEGERS
      INTEGER seed
      seed=-1
#else
      real*8  seed
      seed=-1.d0
#endif
c     seed=1618033.0d0
      DO i=1,n
       arr(i)=port_rand(seed)
      ENDDO

      RETURN
      END

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: PORT_RAND_NORM
C     !INTERFACE:
      real*8 FUNCTION PORT_RAND_NORM()

C     !DESCRIPTION:
C     This function generates a normally distributed random number with
C     the so called polar algorithm. The algorithm actually generates 2
C     numbers, but only 1 is returned for maximum compatibility with old
C     code.  The most obvious way to improve this function would be to
C     make sure that the second number is not wasted.

C     Changed: 2004.09.06 antti.westerlund@fimr.fi

C     !USES:
      IMPLICIT NONE
CEOP

C     !LOCAL VARIABLES:
      real*8 port_rand
      real*8 x1, x2, xs, t

#ifdef _USE_INTEGERS
      INTEGER seed
      seed=-1
#else
      real*8  seed
      seed=-1.d0
#endif
c     seed=1618033.0d0

C     first generate 2 equally distributed random numbers (-1,1)
      DO WHILE (1 .EQ. 1)
         x1=2.0*port_rand(seed)-1.0
         x2=2.0*port_rand(seed)-1.0
         xs=x1**2+x2**2
         IF (xs .LT. 1.0 .AND. xs .NE. 0.0) THEN
            GOTO 100
         ENDIF
      ENDDO
 100  CONTINUE

      t = SQRT(-2.0*LOG(xs)/xs)
      port_rand_norm = t*x1
C
C     also t*x2 would be a gaussian random number and could be returned

      RETURN
      END
1	jmc	1.7	C $Header: /u/gcmpack/MITgcm/model/src/port_rand.F,v 1.6 2005/03/01 16:51:27 jmc Exp $
2	cnh	1.3	C $Name: $
3	adcroft	1.2
4	jmc	1.7	#undef _USE_INTEGERS
5	adcroft	1.2
6	jmc	1.7	C-- File port_rand.F: Portable random number generator functions
7			C-- Contents
8			C-- o PORT_RAND
9			C-- o PORT_RANARR
10			C-- o PORT_RAND_NORM
11	adcroft	1.2
12	edhill	1.5	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
13	cnh	1.3	CBOP
14	jmc	1.7	C !ROUTINE: PORT_RAND
15	cnh	1.3	C !INTERFACE:
16	jmc	1.7	real*8 FUNCTION PORT_RAND(seed)
17	cnh	1.3
18	edhill	1.5	C !DESCRIPTION:
19	cnh	1.3	C Portable random number generator
20	jmc	1.6	C seed >=0 :: initialise using this seed ; and return 0
21			C seed < 0 :: if first call then initialise using the default seed (=mseed)
22			C and always return a random number
23	jmc	1.7
24	edhill	1.5	C !USES:
25	jmc	1.7	IMPLICIT NONE
26
27	jmc	1.6	C !INPUT PARAMETERS:
28			#ifdef _USE_INTEGERS
29	jmc	1.7	INTEGER seed
30	jmc	1.6	#else
31			real*8 seed
32			#endif
33	edhill	1.5	CEOP
34	cnh	1.3
35			C !LOCAL VARIABLES:
36	jmc	1.7	INTEGER nff,idum
37			PARAMETER(nff=55)
38			PARAMETER(idum=-2)
39	adcroft	1.2	real*8 fac
40			#ifdef _USE_INTEGERS
41	jmc	1.7	INTEGER mbig,mseed,mZ
42			PARAMETER (mbig=1000000000,mz=0,fac=1.d0/mbig)
43			INTEGER mj,mk,ma(nff)
44			DATA mseed/161803398/
45	adcroft	1.2	#else
46	jmc	1.6	real*8 mbig,mseed,mz
47	jmc	1.7	PARAMETER (mbig=4000000.,mz=0.,fac=1.d0/mbig)
48	jmc	1.6	real*8 mj,mk,ma(nff)
49	jmc	1.7	DATA mseed/1618033./
50	adcroft	1.2	#endif
51	jmc	1.7	LOGICAL firstCall
52			INTEGER i,ii,inext,inextp,k
53			DATA firstCall /.true./
54			SAVE firstCall,inext,inextp,ma
55	edhill	1.5
56	jmc	1.6	C- Initialise the random number generator
57	jmc	1.7	IF (firstCall .OR. seed.GE.mz) THEN
58			IF (seed.GE.mz) mseed = seed
59	adcroft	1.2	firstCall=.false.
60			mj=mseed-iabs(idum)
61			mj=mod(mj,mbig)
62			ma(nff)=mj
63			mk=1
64	jmc	1.7	DO i=1,nff-1
65	adcroft	1.2	ii=mod(21*i,nff)
66			ma(ii)=mk
67			mk=mj-mk
68	jmc	1.7	IF (mk.LT.mz) mk=mk+mbig
69	adcroft	1.2	mj=ma(ii)
70	jmc	1.7	ENDDO
71			DO k=1,4
72			DO i=1,nff
73	adcroft	1.2	ma(i)=ma(i)-ma(1+mod(i+30,nff))
74	jmc	1.7	IF (ma(i).LT.mz) ma(i)=ma(i)+mbig
75			ENDDO
76			ENDDO
77	adcroft	1.2	inext=0
78			inextp=31
79	jmc	1.7	ENDIF
80	jmc	1.6
81			C- Compute a random number (only if seed < 0)
82	jmc	1.7	IF (seed.GE.mz) THEN
83	jmc	1.6	port_rand=0.d0
84	jmc	1.7	ELSE
85	jmc	1.6	inext=mod(inext,nff)+1
86			inextp=mod(inextp,nff)+1
87			mj=ma(inext)-ma(inextp)
88	jmc	1.7	IF (mj.LT.mz) mj=mj+mbig
89	jmc	1.6	ma(inext)=mj
90			port_rand=mj*fac
91	jmc	1.7	ENDIF
92	jmc	1.6
93	jmc	1.7	RETURN
94			END
95	adcroft	1.2
96	edhill	1.5	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
97
98	jmc	1.7	CBOP
99			C !ROUTINE: PORT_RANARR
100			C !INTERFACE:
101			SUBROUTINE PORT_RANARR(n,arr)
102
103			C !DESCRIPTION:
104			C Portable random number array generator
105
106			C !USES:
107			IMPLICIT NONE
108
109			C !INPUT PARAMETERS:
110			INTEGER n
111			real*8 arr(n)
112			CEOP
113
114			C !LOCAL VARIABLES:
115			INTEGER i
116	jmc	1.6	real*8 port_rand
117			#ifdef _USE_INTEGERS
118	jmc	1.7	INTEGER seed
119	jmc	1.6	seed=-1
120			#else
121			real*8 seed
122			seed=-1.d0
123			#endif
124			c seed=1618033.0d0
125	jmc	1.7	DO i=1,n
126	molod	1.4	arr(i)=port_rand(seed)
127	jmc	1.7	ENDDO
128	adcroft	1.2
129	jmc	1.7	RETURN
130			END
131	edhill	1.5
132			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
133			CBOP
134	jmc	1.7	C !ROUTINE: PORT_RAND_NORM
135	edhill	1.5	C !INTERFACE:
136	jmc	1.7	real*8 FUNCTION PORT_RAND_NORM()
137	edhill	1.5
138			C !DESCRIPTION:
139			C This function generates a normally distributed random number with
140			C the so called polar algorithm. The algorithm actually generates 2
141			C numbers, but only 1 is returned for maximum compatibility with old
142			C code. The most obvious way to improve this function would be to
143			C make sure that the second number is not wasted.
144
145			C Changed: 2004.09.06 antti.westerlund@fimr.fi
146
147			C !USES:
148	jmc	1.7	IMPLICIT NONE
149	edhill	1.5	CEOP
150
151			C !LOCAL VARIABLES:
152	jmc	1.6	real*8 port_rand
153	jmc	1.7	real*8 x1, x2, xs, t
154	edhill	1.5
155	jmc	1.6	#ifdef _USE_INTEGERS
156	jmc	1.7	INTEGER seed
157	jmc	1.6	seed=-1
158			#else
159			real*8 seed
160			seed=-1.d0
161			#endif
162			c seed=1618033.0d0
163
164	jmc	1.7	C first generate 2 equally distributed random numbers (-1,1)
165			DO WHILE (1 .EQ. 1)
166	edhill	1.5	x1=2.0*port_rand(seed)-1.0
167			x2=2.0*port_rand(seed)-1.0
168			xs=x12+x22
169	jmc	1.7	IF (xs .LT. 1.0 .AND. xs .NE. 0.0) THEN
170			GOTO 100
171			ENDIF
172			ENDDO
173			100 CONTINUE
174	edhill	1.5
175	jmc	1.7	t = SQRT(-2.0*LOG(xs)/xs)
176	edhill	1.5	port_rand_norm = t*x1
177			C
178			C also t*x2 would be a gaussian random number and could be returned
179
180	jmc	1.7	RETURN
181			END