source: LMDZ6/branches/Ocean_skin/libf/phylmd/ecrad/random_numbers_mix.F90 @ 5407

Last change on this file since 5407 was 4368, checked in by lguez, 2 years ago

Sync latest trunk changes to Ocean_skin
File size: 17.4 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=Z"3FFFFFFF"
242  INTEGER(KIND=JPIM) :: IVAR
243  DATA IVAR /Z"3FFFFFFF"/
244  TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM
245  REAL(KIND=JPRB), DIMENSION(:),     INTENT(  OUT) :: PX
246  INTEGER(KIND=JPIM)                :: JJ, JK, IN, IFILLED
247 
248  ! This test is a little dirty but Fortran 90 doesn't allow for the initialization
249  !   of components of derived types.
250  REAL(KIND=JPRB) :: ZHOOK_HANDLE
251! DR_HOOK removed to reduce overhead
252! IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:UNIFORM_DISTRIBUTION',0,ZHOOK_HANDLE)
253  IF(YD_STREAM%INITTEST /= INITVALUE) &
254    & CALL ABOR1 ('uniform_distribution called before initialize_random_numbers')
255 
256  !--------------------------------------------------------------------------------
257  ! Copy numbers that were generated during the last call, but not used.
258  !--------------------------------------------------------------------------------
259 
260  IN=SIZE(PX)
261  IFILLED=0
262 
263  DO JJ=YD_STREAM%IUSED+1,MIN(JPQ,IN+YD_STREAM%IUSED)
264    PX(JJ-YD_STREAM%IUSED) = YD_STREAM%IX(JJ)*YD_STREAM%ZRM
265    IFILLED=IFILLED+1
266  ENDDO
267 
268  YD_STREAM%IUSED=YD_STREAM%IUSED+IFILLED
269 
270  IF (IFILLED==IN)  THEN
271!   IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:UNIFORM_DISTRIBUTION',1,ZHOOK_HANDLE)
272! DR_HOOK removed to reduce overhead
273    RETURN
274  ENDIF
275 
276  !--------------------------------------------------------------------------------
277  ! Generate batches of jpq numbers until px has been filled
278  !--------------------------------------------------------------------------------
279 
280  DO WHILE (IFILLED<IN)
281 
282  !--------------------------------------------------------------------------------
283  ! Generate jpq numbers in vectorizable loops. The first loop is length jpp. The
284  ! remaining jpq-jpp elements are calculated in loops of length shorter than jpp.
285  !--------------------------------------------------------------------------------
286 
287  !OCL NOVREC
288    DO JJ=1,JPP
289!     yd_stream%ix(jj) = yd_stream%ix(jj) + yd_stream%ix(jj-jpp+jpq)
290!     if (yd_stream%ix(jj)>=jpm) yd_stream%ix(jj) = yd_stream%ix(jj)-jpm
291      YD_STREAM%IX(JJ) = IAND(IVAR,YD_STREAM%IX(JJ) + YD_STREAM%IX(JJ-JPP+JPQ))
292    ENDDO
293 
294    DO JK=1,JPNUMSPLIT
295  !OCL NOVREC
296      DO JJ=1+JPP+(JK-1)*JPLENSPLIT,MIN(JPQ,JPP+JK*JPLENSPLIT)
297!       yd_stream%ix(jj) = yd_stream%ix(jj) + yd_stream%ix(jj-jpp)
298!       if (yd_stream%ix(jj)>=jpm) yd_stream%ix(jj) = yd_stream%ix(jj)-jpm
299        YD_STREAM%IX(JJ) = IAND(IVAR,YD_STREAM%IX(JJ) + YD_STREAM%IX(JJ-JPP))
300      ENDDO
301    ENDDO
302 
303    YD_STREAM%IUSED = MIN(JPQ,IN-IFILLED)
304    PX(IFILLED+1:IFILLED+YD_STREAM%IUSED) = YD_STREAM%IX(1:YD_STREAM%IUSED)*YD_STREAM%ZRM
305    IFILLED = IFILLED+YD_STREAM%IUSED
306  ENDDO
307 
308!IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:UNIFORM_DISTRIBUTION',1,ZHOOK_HANDLE)
309! DR_HOOK removed to reduce overhead
310END SUBROUTINE UNIFORM_DISTRIBUTION
311!-------------------------------------------------------------------------------
312SUBROUTINE GAUSSIAN_DISTRIBUTION (PX, YD_STREAM)
313  TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM
314  REAL(KIND=JPRB),                   INTENT(  OUT) :: PX(:)
315  !--------------------------------------------------------------------------------
316  ! Generate normally-distributed random numbers using the Box-Muller method.
317  !
318  ! NB: this routine does not use buffering. This means that the following calls:
319  !     call gaussian_distribution (zx(1:k))
320  !     call gaussian_distribution (zx(k+1:n))
321  ! will produce different numbers for elements k+1 onwards than the single call:
322  !     call gaussian_distribution (zx(1:n))
323  !--------------------------------------------------------------------------------
324 
325  INTEGER(KIND=JPIM) :: ILEN, J
326  REAL(KIND=JPRB) :: ZFAC, ZTWOPI
327  REAL(KIND=JPRB) :: ZX(SIZE(PX)+1)
328 
329  !--------------------------------------------------------------------------------
330  ! Generate uniform random points in the range [0,1)
331  !--------------------------------------------------------------------------------
332
333    REAL(KIND=JPRB) :: ZHOOK_HANDLE
334    IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:GAUSSIAN_DISTRIBUTION',0,ZHOOK_HANDLE)
335    CALL UNIFORM_DISTRIBUTION (ZX, YD_STREAM)
336
337  !--------------------------------------------------------------------------------
338  ! Generate gaussian deviates using Box-Muller method
339  !--------------------------------------------------------------------------------
340 
341  ZTWOPI = 8.0_JPRB*ATAN(1.0_JPRB)
342  ILEN=SIZE(PX)
343 
344  DO J=1,ILEN-1,2
345    ZFAC = SQRT(-2.0_JPRB*LOG(1.0_JPRB-ZX(J)))
346    PX(J  ) = ZFAC*COS(ZTWOPI*ZX(J+1))
347    PX(J+1) = ZFAC*SIN(ZTWOPI*ZX(J+1))
348  ENDDO
349 
350  !--------------------------------------------------------------------------------
351  ! Generate the last point if ilen is odd
352  !--------------------------------------------------------------------------------
353 
354  IF (MOD(ILEN,2) /= 0) THEN
355    ZFAC = SQRT(-2.0_JPRB*LOG(1.0_JPRB-ZX(ILEN)))
356    PX(ILEN) = ZFAC*COS(ZTWOPI*ZX(ILEN+1))
357  ENDIF
358 
359IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:GAUSSIAN_DISTRIBUTION',1,ZHOOK_HANDLE)
360END SUBROUTINE GAUSSIAN_DISTRIBUTION
361!-------------------------------------------------------------------------------
362!!$SUBROUTINE RANDOM_NUMBER_RESTARTFILE( CDFNAME, CDACTION,YD_STREAM )
363!!$!--------------------------------------------------------------------------------
364!!$!
365!!$! read (cdaction='r') or write (cdaction='w') restart file
366!!$! for random number generator
367!!$!
368!!$!--------------------------------------------------------------------------------
369!!$CHARACTER (LEN=*),   INTENT(IN) :: CDFNAME
370!!$CHARACTER (LEN=1  ), INTENT(IN) :: CDACTION
371!!$TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM
372!!$ 
373!!$INTEGER(KIND=JPIM) :: IUNIT, IRET, IBYTES_IN_JPIM
374!!$
375!!$REAL(KIND=JPRB) :: ZHOOK_HANDLE
376!!$IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:RANDOM_NUMBER_RESTARTFILE',0,ZHOOK_HANDLE)
377!!$IBYTES_IN_JPIM= CEILING(REAL(BIT_SIZE(YD_STREAM%IUSED))/8.0_JPRB - TINY(1.0_JPRB))
378!!$
379!!$IF (IBYTES_IN_JPIM /= 4) THEN
380!!$  CALL ABOR1('random_number_restartfile: number of bytes for JPIM is not 4 ')       
381!!$ENDIF
382!!$
383!!$CALL PBOPEN(IUNIT, CDFNAME, CDACTION, IRET)
384!!$IF (IRET /= 0) THEN
385!!$  CALL ABOR1('random_number_restartfile: PBOPEN FAILED opening '//CDFNAME)   
386!!$ENDIF
387!!$
388!!$
389!!$IF (CDACTION=='r' .OR. CDACTION=='R') THEN
390!!$  CALL PBREAD(IUNIT, YD_STREAM%IX,    IBYTES_IN_JPIM*JPQ, IRET)
391!!$  IF (IRET < 0) THEN
392!!$    CALL ABOR1('random_number_restartfile: PBREAD could not read ix from '//CDFNAME)   
393!!$  ENDIF
394!!$  CALL PBREAD(IUNIT, YD_STREAM%IUSED, IBYTES_IN_JPIM    , IRET)
395!!$  IF (IRET < 0) THEN
396!!$    CALL ABOR1('random_number_restartfile: PBREAD could not read iused from '//CDFNAME)   
397!!$  ENDIF
398!!$
399!!$!  l_initialized = .TRUE.
400!!$  YD_STREAM%INITTEST = INITVALUE
401!!$  YD_STREAM%ZRM=1.0_JPRB/REAL(JPM,JPRB)
402!!$ELSEIF(CDACTION=='w' .OR. CDACTION=='W') THEN
403!!$  CALL PBWRITE(IUNIT, YD_STREAM%IX, IBYTES_IN_JPIM*JPQ, IRET)
404!!$  IF (IRET < 0) THEN
405!!$    CALL ABOR1('random_number_restartfile: PBWRITE could not write ix on '//CDFNAME)   
406!!$  ENDIF
407!!$  CALL PBWRITE(IUNIT, YD_STREAM%IUSED, IBYTES_IN_JPIM , IRET)
408!!$  IF (IRET < 0) THEN
409!!$    CALL ABOR1('random_number_restartfile: PBWRITE could not write iused on '//CDFNAME)   
410!!$  ENDIF
411!!$
412!!$ELSE
413!!$  CALL ABOR1 ('random_number_restartfile: cdaction = '//CDACTION//' is undefined.')
414!!$ENDIF
415!!$
416!!$CALL PBCLOSE(IUNIT, IRET)
417!!$IF (IRET /= 0) THEN
418!!$  CALL ABOR1('random_number_restartfile: PBCLOSE FAILED closing '//CDFNAME)   
419!!$ENDIF
420!!$
421!!$IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:RANDOM_NUMBER_RESTARTFILE',1,ZHOOK_HANDLE)
422!!$END SUBROUTINE RANDOM_NUMBER_RESTARTFILE
423
424
425SUBROUTINE WR_RANGEN_STATE( KUNIT, YD_STREAM )
426!--------------------------------------------------------------------------------
427! write state of random number generator to unit kunit
428!--------------------------------------------------------------------------------
429INTEGER(KIND=JPIM), INTENT(IN) :: KUNIT
430TYPE(RANDOMNUMBERSTREAM), INTENT(IN) :: YD_STREAM
431
432REAL(KIND=JPRB) :: ZHOOK_HANDLE
433IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:WR_RANGEN_STATE',0,ZHOOK_HANDLE)
434WRITE( KUNIT, * ) 'module random_numbers_mix, generator state is'
435WRITE( KUNIT, '(8I10)') YD_STREAM%IX
436WRITE( KUNIT, '(I10)')  YD_STREAM%IUSED
437
438IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:WR_RANGEN_STATE',1,ZHOOK_HANDLE)
439END SUBROUTINE WR_RANGEN_STATE
440
441END MODULE RANDOM_NUMBERS_MIX
Note: See TracBrowser for help on using the repository browser.