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

Last change on this file since 12377 was 12377, checked in by acc, 6 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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