source: NEMO/trunk/src/OCE/STO/storng.F90 @ 12649

Last change on this file since 12649 was 12649, checked in by smasson, 8 months ago

trunk: clean and unique 3rd dimension in iom_rstput, see #2432

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