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/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/STO – NEMO

source: branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/STO/storng.F90

Last change on this file was 7470, checked in by jcastill, 8 years ago

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