New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
storng.F90 in trunk/NEMOGCM/NEMO/OPA_SRC – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/storng.F90 @ 5366

Last change on this file since 5366 was 5366, checked in by pabouttier, 9 years ago

Add AGRIF_DO_NOT_TREAT directve in storng to compile AGRIF configuration in SETTE; See Ticket #1527
File size: 13.1 KB
Line 
1MODULE storng
2!$AGRIF_DO_NOT_TREAT
3   !!======================================================================
4   !!                       ***  MODULE  storng  ***
5   !! Random number generator, used in NEMO stochastic parameterization
6   !!
7   !!=====================================================================
8   !! History :  3.3  ! 2011-10 (J.-M. Brankart)  Original code
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !! The module is based on (and includes) the
13   !! 64-bit KISS (Keep It Simple Stupid) random number generator
14   !! distributed by George Marsaglia :
15   !! http://groups.google.com/group/comp.lang.fortran/
16   !!        browse_thread/thread/a85bf5f2a97f5a55
17   !!----------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   kiss          : 64-bit KISS random number generator (period ~ 2^250)
21   !!   kiss_seed     : Define seeds for KISS random number generator
22   !!   kiss_state    : Get current state of KISS random number generator
23   !!   kiss_save     : Save current state of KISS (for future restart)
24   !!   kiss_load     : Load the saved state of KISS
25   !!   kiss_reset    : Reset the default seeds
26   !!   kiss_check    : Check the KISS pseudo-random sequence
27   !!   kiss_uniform  : Real random numbers with uniform distribution in [0,1]
28   !!   kiss_gaussian : Real random numbers with Gaussian distribution N(0,1)
29   !!   kiss_gamma    : Real random numbers with Gamma distribution Gamma(k,1)
30   !!   kiss_sample   : Select a random sample from a set of integers
31   !!
32   !!   ---CURRENTLY NOT USED IN NEMO :
33   !!   kiss_save, kiss_load, kiss_check, kiss_gamma, kiss_sample
34   !!----------------------------------------------------------------------
35   IMPLICIT NONE
36   PRIVATE
37
38   ! Public functions/subroutines
39   PUBLIC :: kiss, kiss_seed, kiss_state, kiss_save, kiss_load, kiss_reset, kiss_check
40   PUBLIC :: kiss_uniform, kiss_gaussian, kiss_gamma, kiss_sample
41
42   ! Default/initial seeds
43   INTEGER(KIND=8) :: x=1234567890987654321_8
44   INTEGER(KIND=8) :: y=362436362436362436_8
45   INTEGER(KIND=8) :: z=1066149217761810_8
46   INTEGER(KIND=8) :: w=123456123456123456_8
47
48   ! Parameters to generate real random variates
49   REAL(KIND=8), PARAMETER :: huge64=9223372036854775808.0_8  ! +1
50   REAL(KIND=8), PARAMETER :: zero=0.0_8, half=0.5_8, one=1.0_8, two=2.0_8
51
52   ! Variables to store 2 Gaussian random numbers with current index (ig)
53   INTEGER(KIND=8), SAVE :: ig=1
54   REAL(KIND=8), SAVE :: gran1, gran2
55
56   !!----------------------------------------------------------------------
57   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
58   !! $Id: dynhpg.F90 2528 2010-12-27 17:33:53Z rblod $
59   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   FUNCTION kiss()
64      !! --------------------------------------------------------------------
65      !!                  ***  FUNCTION kiss  ***
66      !!
67      !! ** Purpose :   64-bit KISS random number generator
68      !!
69      !! ** Method  :   combine several random number generators:
70      !!                (1) Xorshift (XSH), period 2^64-1,
71      !!                (2) Multiply-with-carry (MWC), period (2^121+2^63-1)
72      !!                (3) Congruential generator (CNG), period 2^64.
73      !!
74      !!                overall period:
75      !!                (2^250+2^192+2^64-2^186-2^129)/6
76      !!                            ~= 2^(247.42) or 10^(74.48)
77      !!
78      !!                set your own seeds with 'kiss_seed'
79      ! --------------------------------------------------------------------
80      IMPLICIT NONE
81      INTEGER(KIND=8) :: kiss, t
82
83      t = ISHFT(x,58) + w
84      IF (s(x).eq.s(t)) THEN
85         w = ISHFT(x,-6) + s(x)
86      ELSE
87         w = ISHFT(x,-6) + 1 - s(x+t)
88      ENDIF
89      x = t + x
90      y = m( m( m(y,13_8), -17_8 ), 43_8 )
91      z = 6906969069_8 * z + 1234567_8
92
93      kiss = x + y + z
94
95      CONTAINS
96
97         FUNCTION s(k)
98            INTEGER(KIND=8) :: s, k
99            s = ISHFT(k,-63)
100         END FUNCTION s
101
102         FUNCTION m(k, n)
103            INTEGER(KIND=8) :: m, k, n
104            m =  IEOR(k, ISHFT(k, n) )
105         END FUNCTION m
106
107   END FUNCTION kiss
108
109
110   SUBROUTINE kiss_seed(ix, iy, iz, iw)
111      !! --------------------------------------------------------------------
112      !!                  ***  ROUTINE kiss_seed  ***
113      !!
114      !! ** Purpose :   Define seeds for KISS random number generator
115      !!
116      !! --------------------------------------------------------------------
117      IMPLICIT NONE
118      INTEGER(KIND=8) :: ix, iy, iz, iw
119
120      x = ix
121      y = iy
122      z = iz
123      w = iw
124
125   END SUBROUTINE kiss_seed
126
127
128   SUBROUTINE kiss_state(ix, iy, iz, iw)
129      !! --------------------------------------------------------------------
130      !!                  ***  ROUTINE kiss_state  ***
131      !!
132      !! ** Purpose :   Get current state of KISS random number generator
133      !!
134      !! --------------------------------------------------------------------
135      IMPLICIT NONE
136      INTEGER(KIND=8) :: ix, iy, iz, iw
137
138      ix = x
139      iy = y
140      iz = z
141      iw = w
142
143   END SUBROUTINE kiss_state
144
145
146   SUBROUTINE kiss_reset()
147      !! --------------------------------------------------------------------
148      !!                  ***  ROUTINE kiss_reset  ***
149      !!
150      !! ** Purpose :   Reset the default seeds for KISS random number generator
151      !!
152      !! --------------------------------------------------------------------
153      IMPLICIT NONE
154
155      x=1234567890987654321_8
156      y=362436362436362436_8
157      z=1066149217761810_8
158      w=123456123456123456_8
159
160    END SUBROUTINE kiss_reset
161
162
163    SUBROUTINE kiss_check(check_type)
164      !! --------------------------------------------------------------------
165      !!                  ***  ROUTINE kiss_check  ***
166      !!
167      !! ** Purpose :   Check the KISS pseudo-random sequence
168      !!
169      !! ** Method  :   Check that it reproduces the correct sequence
170      !!                from the default seed
171      !!
172      !! --------------------------------------------------------------------
173      IMPLICIT NONE
174      INTEGER(KIND=8) :: iter, niter, correct, iran
175      CHARACTER(LEN=*) :: check_type
176      LOGICAL :: print_success
177
178      ! Save current state of KISS
179      CALL kiss_save()
180      ! Reset the default seed
181      CALL kiss_reset()
182
183      ! Select check type
184      SELECT CASE(check_type)
185      CASE('short')
186         niter = 5_8
187         correct = 542381058189297533_8
188         print_success = .FALSE.
189      CASE('long')
190         niter = 100000000_8
191         correct = 1666297717051644203_8 ! Check provided by G. Marsaglia
192         print_success = .TRUE.
193      CASE('default')
194      CASE DEFAULT
195         STOP 'Bad check type in kiss_check'
196      END SELECT
197
198      ! Run kiss for the required number of iterations (niter)
199      DO iter=1,niter
200         iran = kiss()
201      ENDDO
202
203      ! Check that last iterate is correct
204      IF (iran.NE.correct) THEN
205         STOP 'Check failed: KISS internal error !!'
206      ELSE
207         IF (print_success) PRINT *, 'Check successful: 100 million calls to KISS OK'
208      ENDIF
209
210      ! Reload the previous state of KISS
211      CALL kiss_load()
212
213   END SUBROUTINE kiss_check
214
215
216   SUBROUTINE kiss_save
217      !! --------------------------------------------------------------------
218      !!                  ***  ROUTINE kiss_save  ***
219      !!
220      !! ** Purpose :   Save current state of KISS random number generator
221      !!
222      !! --------------------------------------------------------------------
223      IMPLICIT NONE
224
225      OPEN(UNIT=30,FILE='.kiss_restart')
226      WRITE(30,*) x
227      WRITE(30,*) y
228      WRITE(30,*) z
229      WRITE(30,*) w
230      CLOSE(30)
231
232    END SUBROUTINE kiss_save
233
234
235    SUBROUTINE kiss_load
236      !! --------------------------------------------------------------------
237      !!                  ***  ROUTINE kiss_load  ***
238      !!
239      !! ** Purpose :   Load the saved state of KISS random number generator
240      !!
241      !! --------------------------------------------------------------------
242      IMPLICIT NONE
243      LOGICAL :: filexists
244
245      INQUIRE(FILE='.kiss_restart',EXIST=filexists)
246      IF (filexists) THEN
247         OPEN(UNIT=30,FILE='.kiss_restart')
248         READ(30,*) x
249         READ(30,*) y
250         READ(30,*) z
251         READ(30,*) w
252         CLOSE(30)
253      ENDIF
254
255   END SUBROUTINE kiss_load
256
257
258   SUBROUTINE kiss_uniform(uran)
259      !! --------------------------------------------------------------------
260      !!                  ***  ROUTINE kiss_uniform  ***
261      !!
262      !! ** Purpose :   Real random numbers with uniform distribution in [0,1]
263      !!
264      !! --------------------------------------------------------------------
265      IMPLICIT NONE
266      REAL(KIND=8) :: uran
267
268      uran = half * ( one + REAL(kiss(),8) / huge64 )
269
270   END SUBROUTINE kiss_uniform
271
272
273   SUBROUTINE kiss_gaussian(gran)
274      !! --------------------------------------------------------------------
275      !!                  ***  ROUTINE kiss_gaussian  ***
276      !!
277      !! ** Purpose :   Real random numbers with Gaussian distribution N(0,1)
278      !!
279      !! ** Method  :   Generate 2 new Gaussian draws (gran1 and gran2)
280      !!                from 2 uniform draws on [-1,1] (u1 and u2),
281      !!                using the Marsaglia polar method
282      !!                (see Devroye, Non-Uniform Random Variate Generation, p. 235-236)
283      !!
284      !! --------------------------------------------------------------------
285      IMPLICIT NONE
286      REAL(KIND=8) :: gran, u1, u2, rsq, fac
287
288      IF (ig.EQ.1) THEN
289         rsq = two
290         DO WHILE ( (rsq.GE.one).OR. (rsq.EQ.zero) )
291            u1 = REAL(kiss(),8) / huge64
292            u2 = REAL(kiss(),8) / huge64
293            rsq = u1*u1 + u2*u2
294         ENDDO
295         fac = SQRT(-two*LOG(rsq)/rsq)
296         gran1 = u1 * fac
297         gran2 = u2 * fac
298      ENDIF
299
300      ! Output one of the 2 draws
301      IF (ig.EQ.1) THEN
302         gran = gran1 ; ig = 2
303      ELSE
304         gran = gran2 ; ig = 1
305      ENDIF
306
307   END SUBROUTINE kiss_gaussian
308
309
310   SUBROUTINE kiss_gamma(gamr,k)
311      !! --------------------------------------------------------------------
312      !!                  ***  ROUTINE kiss_gamma  ***
313      !!
314      !! ** Purpose :   Real random numbers with Gamma distribution Gamma(k,1)
315      !!
316      !! --------------------------------------------------------------------
317      IMPLICIT NONE
318      REAL(KIND=8), PARAMETER :: p1 = 4.5_8
319      REAL(KIND=8), PARAMETER :: p2 = 2.50407739677627_8  ! 1+LOG(9/2)
320      REAL(KIND=8), PARAMETER :: p3 = 1.38629436111989_8  ! LOG(4)
321      REAL(KIND=8) :: gamr, k, u1, u2, b, c, d, xx, yy, zz, rr, ee
322      LOGICAL :: accepted
323
324      IF (k.GT.one) THEN
325         ! Cheng's rejection algorithm
326         ! (see Devroye, Non-Uniform Random Variate Generation, p. 413)
327         b = k - p3 ; d = SQRT(two*k-one) ; c = k + d
328
329         accepted=.FALSE.
330         DO WHILE (.NOT.accepted)
331            CALL kiss_uniform(u1)
332            yy = LOG(u1/(one-u1)) / d  ! Mistake in Devroye: "* k" instead of "/ d"
333            xx = k * EXP(yy)
334            rr = b + c * yy - xx
335            CALL kiss_uniform(u2)
336            zz = u1 * u1 * u2
337
338            accepted = rr .GE. (zz*p1-p2)
339            IF (.NOT.accepted) accepted =  rr .GE. LOG(zz)
340         ENDDO
341
342         gamr = xx
343
344      ELSEIF (k.LT.one) THEN
345        ! Rejection from the Weibull density
346        ! (see Devroye, Non-Uniform Random Variate Generation, p. 415)
347        c = one/k ; d = (one-k) * EXP( (k/(one-k)) * LOG(k) )
348
349        accepted=.FALSE.
350        DO WHILE (.NOT.accepted)
351           CALL kiss_uniform(u1)
352           zz = -LOG(u1)
353           xx = EXP( c * LOG(zz) )
354           CALL kiss_uniform(u2)
355           ee = -LOG(u2)
356
357           accepted = (zz+ee) .GE. (d+xx)  ! Mistake in Devroye: "LE" instead of "GE"
358        ENDDO
359
360        gamr = xx
361
362      ELSE
363         ! Exponential distribution
364         CALL kiss_uniform(u1)
365         gamr = -LOG(u1)
366
367      ENDIF
368
369   END SUBROUTINE kiss_gamma
370
371
372   SUBROUTINE kiss_sample(a,n,k)
373      !! --------------------------------------------------------------------
374      !!                  ***  ROUTINE kiss_sample  ***
375      !!
376      !! ** Purpose :   Select a random sample of size k from a set of n integers
377      !!
378      !! ** Method  :   The sample is output in the first k elements of a
379      !!                Set k equal to n to obtain a random permutation
380      !!                  of the whole set of integers
381      !!
382      !! --------------------------------------------------------------------
383      IMPLICIT NONE
384      INTEGER(KIND=8), DIMENSION(:) :: a
385      INTEGER(KIND=8) :: n, k, i, j, atmp
386      REAL(KIND=8) :: uran
387
388      ! Select the sample using the swapping method
389      ! (see Devroye, Non-Uniform Random Variate Generation, p. 612)
390      DO i=1,k
391         ! Randomly select the swapping element between i and n (inclusive)
392         CALL kiss_uniform(uran)
393         j = i - 1 + CEILING( REAL(n-i+1,8) * uran )
394         ! Swap elements i and j
395         atmp = a(i) ; a(i) = a(j) ; a(j) = atmp
396      ENDDO
397
398   END SUBROUTINE kiss_sample
399!$AGRIF_END_DO_NOT_TREAT
400END MODULE storng
Note: See TracBrowser for help on using the repository browser.