source: LMDZ6/branches/contrails/libf/phylmd/ecrad.v1.5.1/random_numbers_mix.F90 @ 5452

Last change on this file since 5452 was 4489, checked in by lguez, 21 months ago

Merge LMDZ_ECRad branch back into trunk!
File size: 17.3 KB
Line 
1!**** *RANDOM_NUMBERS_MIX*  - Portable Random Number Generator
2
3! (C) Copyright 2002- ECMWF.
4!
5! This software is licensed under the terms of the Apache Licence Version 2.0
6! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
7!
8! In applying this licence, ECMWF does not waive the privileges and immunities
9! granted to it by virtue of its status as an intergovernmental organisation
10! nor does it submit to any jurisdiction.
11
12!     Purpose.
13!     --------
14!           Generate machine-independent pseudo-random numbers
15
16!**   Interface.
17!     ----------
18!        CALL initialize_random_numbers (kseed, yd_stream)
19!        CALL uniform_distribution      (px   , yd_stream)
20!        CALL gaussian_distribution     (px   , yd_stream)
21!        CALL random_number_restartfile (fname, action)
22!        CALL wr_rangen_state           (nunit)
23
24!        Explicit arguments :
25!        --------------------
26!        kseed  (input)    : integer seed in the range [0,HUGE(kseed)]
27!        yd_stream (optional) : the state of the random number generator
28!        px     (output)   : array to receive random numbers in the range
29
30!        In the case of uniform_distribution, px has values in the range [0.0,1.0)
31
32!        Implicit arguments :
33!        --------------------
34!        None
35
36!     Method.
37!     -------
38!        Based loosly on ZUFALL (Petersen, 1994).
39
40!        The main difference between this generator and ZUFALL is that integer arithmetic
41!        is used. This ensures portability to vector machines that implement different
42!        real arithmetic. In particular, vector machines often implement non-IEEE
43!        arithmetic for their vector pipes. This routine will give identical results for
44!        any integer type with at least 32 bits.
45
46!        The generator is a lagged-Fibonacci generator: x(i) = x(i-p) + x(i-q) mod 2**m.
47!        Lagged-Fibonacci generators have very long repeat periods: (2**q -1) * 2**(m-1)
48!        (i.e about 2.85E191 for q=607, m=30). They pass most tests for randomness.
49
50!        p and q must be chosen carefully. Values from the following table are OK.
51!        Larger values give better random numbers, but smaller values are more
52!        cache-friendly.
53
54!          q         p
55!        9689      4187
56!        4423      2098
57!        2281      1029
58!        1279       418
59!         607       273
60!         521       168
61!         250       103
62!         127        63
63!          97        33
64!          55        24
65!          43        22
66!          31        13
67!          24        10
68
69!        The initial q values of x are set using the binary shirt register method of
70!        Burns and Pryor 1999.
71
72!        Mascagni et al (1995) show how to choose different sets of initial values that
73!        are guaranteed to be drawn from different maximal-length cycles. This requires
74!        the initial values of x(1)...x(q) to be in "canonical form". Specifically,
75!        x(1) must be zero and all but a particular one or two values of x must be
76!        even. For q=607 and p=273, only one element (jpq-jps) must be odd.
77
78!     Externals.
79!     ----------
80!        None
81
82!     Reference.
83!     ----------
84!        Burns P.J. and Pryor D.V. 1999,
85!                             Surface Radiative Transport at Large Scale via Monte Carlo.
86!                             Annual Review of Heat Transfer, Vol 9.
87!
88!        Petersen W.P., 1994, Lagged Fibonacci Series Random Number Generator
89!                             for the NEC SX-3. International Journal of High Speed Computing
90!                             Vol. 6, No. 3, pp387-398.
91!
92!        Mascagni M., Cuccaro S.A., Pryor D.V., Robinson M.L., 1995,
93!                             A Fast, High Quality and Reproducible Parallel Lagged-Fibonacci
94!                             Pseudorandom Number Generator. Journal of Computational Physics
95!                             Vol 119. pp211-219.
96
97!     Author.
98!     -------
99!        Mike Fisher *ECMWF*
100
101!     Modifications.
102!     --------------
103!        Original : 2002-09-25
104!        Made parallel friendly: 2003-08-11 Robert Pincus
105!        M Leutbecher: 2004-05-10 restart capability
106!        M Fisher:     2005-03-30 replaced LCG initialization with shift register
107!     ------------------------------------------------------------------
108
109#ifdef RS6K
110@PROCESS HOT(NOVECTOR) NOSTRICT
111#endif
112MODULE RANDOM_NUMBERS_MIX
113USE YOMHOOK,  ONLY : LHOOK, DR_HOOK
114USE PARKIND1, ONLY : JPIM, JPRB
115
116IMPLICIT NONE
117
118SAVE
119
120PRIVATE
121PUBLIC RANDOMNUMBERSTREAM,WR_RANGEN_STATE, &
122     & INITIALIZE_RANDOM_NUMBERS, UNIFORM_DISTRIBUTION, GAUSSIAN_DISTRIBUTION ! ,&
123!     & RANDOM_NUMBER_RESTARTFILE,
124
125INTEGER(KIND=JPIM), PARAMETER      :: JPP=273, JPQ=607, JPS=105
126INTEGER(KIND=JPIM), PARAMETER      :: JPMM=30
127INTEGER(KIND=JPIM), PARAMETER      :: JPM=2**JPMM
128INTEGER(KIND=JPIM), PARAMETER      :: JPNUMSPLIT=(JPQ-2)/(JPP-1)
129INTEGER(KIND=JPIM), PARAMETER      :: JPLENSPLIT=(JPQ-JPP+JPNUMSPLIT-1)/JPNUMSPLIT
130INTEGER(KIND=JPIM), PARAMETER      :: INITVALUE = 12345678
131
132TYPE RANDOMNUMBERSTREAM
133  PRIVATE
134  INTEGER(KIND=JPIM)                 :: IUSED
135  INTEGER(KIND=JPIM)                 :: INITTEST ! Should initialize to zero, but can't in F90
136  INTEGER(KIND=JPIM), DIMENSION(JPQ) :: IX
137  REAL(KIND=JPRB)                    :: ZRM
138END TYPE RANDOMNUMBERSTREAM
139
140CONTAINS
141!-------------------------------------------------------------------------------
142SUBROUTINE INITIALIZE_RANDOM_NUMBERS (KSEED, YD_STREAM)
143  !-------------------------------------------------------------------------------
144  ! Initialize fibgen
145  !-------------------------------------------------------------------------------
146  INTEGER(KIND=JPIM),                INTENT(IN   ) :: KSEED
147  TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM
148 
149  INTEGER, PARAMETER :: JPMASK=123459876
150  INTEGER(KIND=JPIM), PARAMETER     :: JPWARMUP_SHFT=64, JPWARMUP_LFG=999
151  INTEGER(KIND=JPIM)                :: IDUM,JK,JJ,JBIT
152  REAL(KIND=JPRB), DIMENSION(JPWARMUP_LFG)   :: ZWARMUP
153
154  !-------------------------------------------------------------------------------
155  ! Initialize the buffer using a binary shift register (Burns and Pryor, 1999).
156  ! The Galois representation is used for the shift register as it is more
157  ! efficient than the Fibonacci representation. The magic numbers 31 and 87
158  ! define the shift register primitive polynomial=(32,7,5,3,2,1,0).
159  !
160  ! To ensure that different seeds produce distinct initial buffer states in
161  ! canonical form, bits 0...jpmm-2 of the initial seed (after XORing with jpmask
162  ! and spinning up using the linear congruential generator) are used to construct
163  ! x(2), and the remaining bits are used to construct x(jpq).
164  !-------------------------------------------------------------------------------
165 
166  REAL(KIND=JPRB) :: ZHOOK_HANDLE
167  IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:INITIALIZE_RANDOM_NUMBERS',0,ZHOOK_HANDLE)
168  IDUM = ABS(IEOR(KSEED,JPMASK))
169  IF (IDUM==0) IDUM=JPMASK
170
171  DO JJ=1,JPWARMUP_SHFT
172    IF (BTEST(IDUM,31)) THEN
173      IDUM=IBSET(ISHFT(IEOR(IDUM,87),1),0)
174    ELSE
175      IDUM=IBCLR(ISHFT(IDUM,1),0)
176    ENDIF
177  ENDDO
178
179  YD_STREAM%IX(1:JPQ-1)= 0
180  YD_STREAM%IX(2)      = ISHFT(IBITS(IDUM,0,JPMM-1),1)
181  YD_STREAM%IX(JPQ)    = IBITS(IDUM,JPMM-1,BIT_SIZE(IDUM)+1-JPMM)
182
183  DO JBIT=1,JPMM-1
184    DO JJ=3,JPQ-1
185      IF (BTEST(IDUM,31)) THEN
186        IDUM=IBSET(ISHFT(IEOR(IDUM,87),1),0)
187        YD_STREAM%IX(JJ)=IBSET(YD_STREAM%IX(JJ),JBIT)
188      ELSE
189        IDUM=IBCLR(ISHFT(IDUM,1),0)
190      ENDIF
191    ENDDO
192  ENDDO
193
194  YD_STREAM%IX(JPQ-JPS) = IBSET(YD_STREAM%IX(JPQ-JPS),0)
195 
196  !-------------------------------------------------------------------------------
197  ! Initialize some constants
198  !-------------------------------------------------------------------------------
199 
200  YD_STREAM%IUSED=JPQ
201  YD_STREAM%ZRM=1.0_JPRB/REAL(JPM,JPRB)
202 
203  !-------------------------------------------------------------------------------
204  ! Check the calculation of jpnumsplit and jplensplit.
205  !-------------------------------------------------------------------------------
206 
207  IF (JPP+JPNUMSPLIT*JPLENSPLIT < JPQ) THEN
208    CALL ABOR1 ('initialize_random_numbers: upper limit of last loop < jpq')
209  ENDIF
210 
211  IF (JPLENSPLIT >=JPP) THEN
212    CALL ABOR1 ('initialize_random_numbers: loop length > jpp')
213  ENDIF
214 
215  IF (JPNUMSPLIT>1) THEN
216    IF ((JPQ-JPP+JPNUMSPLIT-2)/(JPNUMSPLIT-1) < JPP) THEN
217      CALL ABOR1 ('initialize_random_numbers: jpnumsplit is bigger than necessary')
218    ENDIF
219  ENDIF
220
221  !-------------------------------------------------------------------------------
222  ! Set initTest to show that the stream is initialized.
223  !-------------------------------------------------------------------------------
224
225  YD_STREAM%INITTEST = INITVALUE
226 
227  !-------------------------------------------------------------------------------
228  ! Warm up the generator.
229  !-------------------------------------------------------------------------------
230
231  CALL UNIFORM_DISTRIBUTION (ZWARMUP, YD_STREAM)
232
233IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:INITIALIZE_RANDOM_NUMBERS',1,ZHOOK_HANDLE)
234END SUBROUTINE INITIALIZE_RANDOM_NUMBERS
235
236!@PROCESS HOT NOSTRICT
237SUBROUTINE UNIFORM_DISTRIBUTION (PX,YD_STREAM)
238  !--------------------------------------------------------------------------------
239  ! Generate uniformly distributed random numbers in the range 0.0<= px < 1.0
240  !--------------------------------------------------------------------------------
241  INTEGER(KIND=JPIM), PARAMETER :: IVAR = INT(Z"3FFFFFFF",JPIM)
242  TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM
243  REAL(KIND=JPRB), DIMENSION(:),     INTENT(  OUT) :: PX
244
245  INTEGER(KIND=JPIM)                :: JJ, JK, IN, IFILLED
246 
247  ! This test is a little dirty but Fortran 90 doesn't allow for the initialization
248  !   of components of derived types.
249  REAL(KIND=JPRB) :: ZHOOK_HANDLE
250! DR_HOOK removed to reduce overhead
251! IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:UNIFORM_DISTRIBUTION',0,ZHOOK_HANDLE)
252  IF(YD_STREAM%INITTEST /= INITVALUE) &
253    & CALL ABOR1 ('uniform_distribution called before initialize_random_numbers')
254
255  !--------------------------------------------------------------------------------
256  ! Copy numbers that were generated during the last call, but not used.
257  !--------------------------------------------------------------------------------
258 
259  IN=SIZE(PX)
260  IFILLED=0
261 
262  DO JJ=YD_STREAM%IUSED+1,MIN(JPQ,IN+YD_STREAM%IUSED)
263    PX(JJ-YD_STREAM%IUSED) = YD_STREAM%IX(JJ)*YD_STREAM%ZRM
264    IFILLED=IFILLED+1
265  ENDDO
266 
267  YD_STREAM%IUSED=YD_STREAM%IUSED+IFILLED
268 
269  IF (IFILLED==IN)  THEN
270!   IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:UNIFORM_DISTRIBUTION',1,ZHOOK_HANDLE)
271! DR_HOOK removed to reduce overhead
272    RETURN
273  ENDIF
274 
275  !--------------------------------------------------------------------------------
276  ! Generate batches of jpq numbers until px has been filled
277  !--------------------------------------------------------------------------------
278 
279  DO WHILE (IFILLED<IN)
280 
281  !--------------------------------------------------------------------------------
282  ! Generate jpq numbers in vectorizable loops. The first loop is length jpp. The
283  ! remaining jpq-jpp elements are calculated in loops of length shorter than jpp.
284  !--------------------------------------------------------------------------------
285 
286  !OCL NOVREC
287    DO JJ=1,JPP
288!     yd_stream%ix(jj) = yd_stream%ix(jj) + yd_stream%ix(jj-jpp+jpq)
289!     if (yd_stream%ix(jj)>=jpm) yd_stream%ix(jj) = yd_stream%ix(jj)-jpm
290      YD_STREAM%IX(JJ) = IAND(IVAR,YD_STREAM%IX(JJ) + YD_STREAM%IX(JJ-JPP+JPQ))
291    ENDDO
292 
293    DO JK=1,JPNUMSPLIT
294  !OCL NOVREC
295      DO JJ=1+JPP+(JK-1)*JPLENSPLIT,MIN(JPQ,JPP+JK*JPLENSPLIT)
296!       yd_stream%ix(jj) = yd_stream%ix(jj) + yd_stream%ix(jj-jpp)
297!       if (yd_stream%ix(jj)>=jpm) yd_stream%ix(jj) = yd_stream%ix(jj)-jpm
298        YD_STREAM%IX(JJ) = IAND(IVAR,YD_STREAM%IX(JJ) + YD_STREAM%IX(JJ-JPP))
299      ENDDO
300    ENDDO
301 
302    YD_STREAM%IUSED = MIN(JPQ,IN-IFILLED)
303    PX(IFILLED+1:IFILLED+YD_STREAM%IUSED) = YD_STREAM%IX(1:YD_STREAM%IUSED)*YD_STREAM%ZRM
304    IFILLED = IFILLED+YD_STREAM%IUSED
305  ENDDO
306 
307!IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:UNIFORM_DISTRIBUTION',1,ZHOOK_HANDLE)
308! DR_HOOK removed to reduce overhead
309END SUBROUTINE UNIFORM_DISTRIBUTION
310!-------------------------------------------------------------------------------
311SUBROUTINE GAUSSIAN_DISTRIBUTION (PX, YD_STREAM)
312  TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM
313  REAL(KIND=JPRB),                   INTENT(  OUT) :: PX(:)
314  !--------------------------------------------------------------------------------
315  ! Generate normally-distributed random numbers using the Box-Muller method.
316  !
317  ! NB: this routine does not use buffering. This means that the following calls:
318  !     call gaussian_distribution (zx(1:k))
319  !     call gaussian_distribution (zx(k+1:n))
320  ! will produce different numbers for elements k+1 onwards than the single call:
321  !     call gaussian_distribution (zx(1:n))
322  !--------------------------------------------------------------------------------
323 
324  INTEGER(KIND=JPIM) :: ILEN, J
325  REAL(KIND=JPRB) :: ZFAC, ZTWOPI
326  REAL(KIND=JPRB) :: ZX(SIZE(PX)+1)
327 
328  !--------------------------------------------------------------------------------
329  ! Generate uniform random points in the range [0,1)
330  !--------------------------------------------------------------------------------
331
332    REAL(KIND=JPRB) :: ZHOOK_HANDLE
333    IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:GAUSSIAN_DISTRIBUTION',0,ZHOOK_HANDLE)
334    CALL UNIFORM_DISTRIBUTION (ZX, YD_STREAM)
335
336  !--------------------------------------------------------------------------------
337  ! Generate gaussian deviates using Box-Muller method
338  !--------------------------------------------------------------------------------
339 
340  ZTWOPI = 8.0_JPRB*ATAN(1.0_JPRB)
341  ILEN=SIZE(PX)
342 
343  DO J=1,ILEN-1,2
344    ZFAC = SQRT(-2.0_JPRB*LOG(1.0_JPRB-ZX(J)))
345    PX(J  ) = ZFAC*COS(ZTWOPI*ZX(J+1))
346    PX(J+1) = ZFAC*SIN(ZTWOPI*ZX(J+1))
347  ENDDO
348 
349  !--------------------------------------------------------------------------------
350  ! Generate the last point if ilen is odd
351  !--------------------------------------------------------------------------------
352 
353  IF (MOD(ILEN,2) /= 0) THEN
354    ZFAC = SQRT(-2.0_JPRB*LOG(1.0_JPRB-ZX(ILEN)))
355    PX(ILEN) = ZFAC*COS(ZTWOPI*ZX(ILEN+1))
356  ENDIF
357 
358IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:GAUSSIAN_DISTRIBUTION',1,ZHOOK_HANDLE)
359END SUBROUTINE GAUSSIAN_DISTRIBUTION
360!-------------------------------------------------------------------------------
361!!$SUBROUTINE RANDOM_NUMBER_RESTARTFILE( CDFNAME, CDACTION,YD_STREAM )
362!!$!--------------------------------------------------------------------------------
363!!$!
364!!$! read (cdaction='r') or write (cdaction='w') restart file
365!!$! for random number generator
366!!$!
367!!$!--------------------------------------------------------------------------------
368!!$CHARACTER (LEN=*),   INTENT(IN) :: CDFNAME
369!!$CHARACTER (LEN=1  ), INTENT(IN) :: CDACTION
370!!$TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM
371!!$ 
372!!$INTEGER(KIND=JPIM) :: IUNIT, IRET, IBYTES_IN_JPIM
373!!$
374!!$REAL(KIND=JPRB) :: ZHOOK_HANDLE
375!!$IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:RANDOM_NUMBER_RESTARTFILE',0,ZHOOK_HANDLE)
376!!$IBYTES_IN_JPIM= CEILING(REAL(BIT_SIZE(YD_STREAM%IUSED))/8.0_JPRB - TINY(1.0_JPRB))
377!!$
378!!$IF (IBYTES_IN_JPIM /= 4) THEN
379!!$  CALL ABOR1('random_number_restartfile: number of bytes for JPIM is not 4 ')       
380!!$ENDIF
381!!$
382!!$CALL PBOPEN(IUNIT, CDFNAME, CDACTION, IRET)
383!!$IF (IRET /= 0) THEN
384!!$  CALL ABOR1('random_number_restartfile: PBOPEN FAILED opening '//CDFNAME)   
385!!$ENDIF
386!!$
387!!$
388!!$IF (CDACTION=='r' .OR. CDACTION=='R') THEN
389!!$  CALL PBREAD(IUNIT, YD_STREAM%IX,    IBYTES_IN_JPIM*JPQ, IRET)
390!!$  IF (IRET < 0) THEN
391!!$    CALL ABOR1('random_number_restartfile: PBREAD could not read ix from '//CDFNAME)   
392!!$  ENDIF
393!!$  CALL PBREAD(IUNIT, YD_STREAM%IUSED, IBYTES_IN_JPIM    , IRET)
394!!$  IF (IRET < 0) THEN
395!!$    CALL ABOR1('random_number_restartfile: PBREAD could not read iused from '//CDFNAME)   
396!!$  ENDIF
397!!$
398!!$!  l_initialized = .TRUE.
399!!$  YD_STREAM%INITTEST = INITVALUE
400!!$  YD_STREAM%ZRM=1.0_JPRB/REAL(JPM,JPRB)
401!!$ELSEIF(CDACTION=='w' .OR. CDACTION=='W') THEN
402!!$  CALL PBWRITE(IUNIT, YD_STREAM%IX, IBYTES_IN_JPIM*JPQ, IRET)
403!!$  IF (IRET < 0) THEN
404!!$    CALL ABOR1('random_number_restartfile: PBWRITE could not write ix on '//CDFNAME)   
405!!$  ENDIF
406!!$  CALL PBWRITE(IUNIT, YD_STREAM%IUSED, IBYTES_IN_JPIM , IRET)
407!!$  IF (IRET < 0) THEN
408!!$    CALL ABOR1('random_number_restartfile: PBWRITE could not write iused on '//CDFNAME)   
409!!$  ENDIF
410!!$
411!!$ELSE
412!!$  CALL ABOR1 ('random_number_restartfile: cdaction = '//CDACTION//' is undefined.')
413!!$ENDIF
414!!$
415!!$CALL PBCLOSE(IUNIT, IRET)
416!!$IF (IRET /= 0) THEN
417!!$  CALL ABOR1('random_number_restartfile: PBCLOSE FAILED closing '//CDFNAME)   
418!!$ENDIF
419!!$
420!!$IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:RANDOM_NUMBER_RESTARTFILE',1,ZHOOK_HANDLE)
421!!$END SUBROUTINE RANDOM_NUMBER_RESTARTFILE
422
423
424SUBROUTINE WR_RANGEN_STATE( KUNIT, YD_STREAM )
425!--------------------------------------------------------------------------------
426! write state of random number generator to unit kunit
427!--------------------------------------------------------------------------------
428INTEGER(KIND=JPIM), INTENT(IN) :: KUNIT
429TYPE(RANDOMNUMBERSTREAM), INTENT(IN) :: YD_STREAM
430
431REAL(KIND=JPRB) :: ZHOOK_HANDLE
432IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:WR_RANGEN_STATE',0,ZHOOK_HANDLE)
433WRITE( KUNIT, * ) 'module random_numbers_mix, generator state is'
434WRITE( KUNIT, '(8I10)') YD_STREAM%IX
435WRITE( KUNIT, '(I10)')  YD_STREAM%IUSED
436
437IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:WR_RANGEN_STATE',1,ZHOOK_HANDLE)
438END SUBROUTINE WR_RANGEN_STATE
439
440END MODULE RANDOM_NUMBERS_MIX
Note: See TracBrowser for help on using the repository browser.