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 branches/2015/dev_r5177_CNRS4_stopar/NEMOGCM/NEMO/OPA_SRC – NEMO

source: branches/2015/dev_r5177_CNRS4_stopar/NEMOGCM/NEMO/OPA_SRC/storng.F90 @ 5304

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

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