MOM_random.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Provides gridded random number capability
6module mom_random
7
9use mom_time_manager, only : time_type, set_date, get_date
10
11use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit
12use iso_fortran_env, only : int32
13
14implicit none ; private
15
17public :: random_01
18public :: random_01_cb
19public :: random_norm
21public :: random_2d_01
22public :: random_2d_norm
23public :: random_unit_tests
24
25! Private period parameters for the Mersenne Twister
26integer, parameter :: &
27 blocksize = 624, & !< Size of the state vector
28 m = 397, & !< Pivot element in state vector
29 matrix_a = -1727483681, & !< constant vector a (0x9908b0dfUL)
30 umask = ibset(0, 31), & !< most significant w-r bits (0x80000000UL)
31 lmask = 2147483647 !< least significant r bits (0x7fffffffUL)
32
33! Private tempering parameters for the Mersenne Twister
34integer, parameter :: tmaskb= -1658038656, & !< (0x9d2c5680UL)
35 tmaskc= -272236544 !< (0xefc60000UL)
36
37!> A private type used by the Mersenne Twistor
39 integer :: currentelement !< Index into state vector
40 integer, dimension(0:blockSize -1) :: state !< State vector
42
43!> Container for pseudo-random number generators
44type, public :: prng ; private
45
46 !> Scalar random number generator for whole model
47 type(randomnumbersequence) :: stream0d
48
49 !> Random number generator for each cell on horizontal grid
50 type(randomnumbersequence), dimension(:,:), allocatable :: stream2d
51
52end type prng
53
54contains
55
56!> Returns a random number between 0 and 1
57real function random_01(CS)
58 type(prng), intent(inout) :: cs !< Container for pseudo-random number generators
59
60 random_01 = getrandomreal(cs%stream0d)
61
62end function random_01
63
64!> Returns a random number between 0 and 1
65!! See https://arxiv.org/abs/2004.06278. Not an exact reproduction of "squares" because Fortran
66!! doesn't have a uint64 type, and not all compilers provide integers with > 64 bits...
67real function random_01_cb(ctr, key)
68 use iso_fortran_env, only : int64
69 integer, intent(in) :: ctr !< ctr should be incremented each time you call the function
70 integer, intent(in) :: key !< key is like a seed: use a different key for each random stream
71 integer(kind=int64) :: x, y, z ! Follows "Squares" naming convention
72
73 x = (ctr + 1) * (key + 65536) ! 65536 added because keys below that don't work.
74 y = (ctr + 1) * (key + 65536)
75 z = y + (key + 65536)
76 x = x*x + y
77 x = ior(ishft(x,32),ishft(x,-32))
78 x = x*x + z
79 x = ior(ishft(x,32),ishft(x,-32))
80 x = x*x + y
81 x = ior(ishft(x,32),ishft(x,-32))
82 x = x*x + z
83 random_01_cb = .5*(1. + .5*real(int(ishft(x,-32)))/real(2**30))
84
85end function
86
87!> Returns an approximately normally distributed random number with mean 0 and variance 1
88real function random_norm(CS)
89 type(prng), intent(inout) :: cs !< Container for pseudo-random number generators
90 ! Local variables
91 integer :: i
92
93 random_norm = getrandomreal(cs%stream0d) - 0.5
94 do i = 1,11
95 random_norm = random_norm + ( getrandomreal(cs%stream0d) - 0.5 )
96 enddo
97
98end function random_norm
99
100!> Generates random numbers between 0 and 1 for each cell of the model grid
101subroutine random_2d_01(CS, HI, rand)
102 type(prng), intent(inout) :: cs !< Container for pseudo-random number generators
103 type(hor_index_type), intent(in) :: hi !< Horizontal index structure
104 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), intent(out) :: rand !< Random numbers between 0 and 1 [nondim]
105 ! Local variables
106 integer :: i,j
107
108 do j = hi%jsd,hi%jed
109 do i = hi%isd,hi%ied
110 rand(i,j) = getrandomreal( cs%stream2d(i,j) )
111 enddo
112 enddo
113
114end subroutine random_2d_01
115
116!> Returns an approximately normally distributed random number with mean 0 and variance 1
117!! for each cell of the model grid
118subroutine random_2d_norm(CS, HI, rand)
119 type(prng), intent(inout) :: cs !< Container for pseudo-random number generators
120 type(hor_index_type), intent(in) :: hi !< Horizontal index structure
121 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), intent(out) :: rand !< Random numbers between 0 and 1 [nondim]
122 ! Local variables
123 integer :: i,j,n
124
125 do j = hi%jsd,hi%jed
126 do i = hi%isd,hi%ied
127 rand(i,j) = getrandomreal( cs%stream2d(i,j) ) - 0.5
128 enddo
129 do n = 1,11
130 do i = hi%isd,hi%ied
131 rand(i,j) = rand(i,j) + ( getrandomreal( cs%stream2d(i,j) ) - 0.5 )
132 enddo
133 enddo
134 enddo
135
136end subroutine random_2d_norm
137
138!> Constructor for scalar PRNG. Can be used to reset the sequence.
139subroutine random_0d_constructor(CS, Time, seed)
140 type(prng), intent(inout) :: cs !< Container for pseudo-random number generators
141 type(time_type), intent(in) :: time !< Current model time
142 integer, intent(in) :: seed !< Seed for PRNG
143 ! Local variables
144 integer :: tseed
145
146 tseed = seed_from_time(time)
147 tseed = ieor(tseed, seed)
148 cs%stream0d = new_randomnumbersequence(tseed)
149
150end subroutine random_0d_constructor
151
152!> Constructor for gridded PRNG. Can be used to reset the sequence.
153subroutine random_2d_constructor(CS, HI, Time, seed)
154 type(prng), intent(inout) :: cs !< Container for pseudo-random number generators
155 type(hor_index_type), intent(in) :: hi !< Horizontal index structure
156 type(time_type), intent(in) :: time !< Current model time
157 integer, intent(in) :: seed !< Seed for PRNG
158 ! Local variables
159 integer :: i,j,sseed,tseed
160
161 if (.not. allocated(cs%stream2d)) allocate( cs%stream2d(hi%isd:hi%ied,hi%jsd:hi%jed) )
162
163 tseed = seed_from_time(time)
164
165 tseed = ieor(tseed*9007, seed)
166 do j = hi%jsd,hi%jed
167 do i = hi%isd,hi%ied
168 sseed = seed_from_index(hi, i, j)
169 sseed = ieor(tseed, sseed*7993)
170 cs%stream2d(i,j) = new_randomnumbersequence(sseed)
171 enddo
172 enddo
173
174end subroutine random_2d_constructor
175
176!> Return a seed derived as hash of values in Time
177integer function seed_from_time(Time)
178 type(time_type), intent(in) :: time !< Current model time
179 ! Local variables
180 integer :: yr,mo,dy,hr,mn,sc,s1,s2
181
182 call get_date(time,yr,mo,dy,hr,mn,sc)
183 s1 = sc + 61*(mn + 61*hr) + 379 ! Range 379 .. 89620
184 ! Fun fact: 2147483647 is the eighth Mersenne prime.
185 ! This is not the reason for using 2147483647 here. It is the
186 ! largest integer of kind=4.
187 s2 = modulo(dy + 32*(mo + 13*yr), 2147483647_4) ! Range 0 .. 2147483646
188 seed_from_time = ieor(s1*4111, s2)
189
190end function seed_from_time
191
192!> Create seed from position index
193integer function seed_from_index(HI, i, j)
194 type(hor_index_type), intent(in) :: hi !< Horizontal index structure
195 integer, intent(in) :: i !< i-index (of h-cell)
196 integer, intent(in) :: j !< j-index (of h-cell)
197 ! Local variables
198 integer :: ig, jg, ni, nj
199
200 ni = hi%niglobal
201 nj = hi%njglobal
202 ! Periodicity is assumed here but does not break non-periodic models
203 ig = mod(hi%idg_offset + i - 1 + ni, ni)+1
204 jg = max(hi%jdg_offset + j, 0)
205 if (jg>nj) then ! Tri-polar hard-coded until we put needed info in HI **TODO**
206 jg = 2*nj+1-jg
207 ig = ni+1-ig
208 endif
209 seed_from_index = ig + ni*(jg-1)
210
211end function seed_from_index
212
213!> Destructor for PRNG
214subroutine random_destruct(CS)
215 type(prng), pointer :: CS !< Container for pseudo-random number generators
216
217 if (allocated(cs%stream2d)) deallocate(cs%stream2d)
218 !deallocate(CS)
219end subroutine random_destruct
220
221!> Return an initialized twister using seed
222!!
223!! Code was based on initialize_scaler() from the FMS implementation of the Mersenne Twistor
224function new_randomnumbersequence(seed) result(twister)
225 integer, intent(in) :: seed !< Seed to initialize twister
226 type(randomnumbersequence) :: twister !< The Mersenne Twister container
227 ! Local variables
228 integer :: i
229
230 twister%state(0) = iand(seed, -1)
231 do i = 1, blocksize - 1 ! ubound(twister%state)
232 twister%state(i) = 1812433253 * ieor(twister%state(i-1), &
233 ishft(twister%state(i-1), -30)) + i
234 twister%state(i) = iand(twister%state(i), -1) ! for >32 bit machines
235 enddo
236 twister%currentElement = blocksize
237end function new_randomnumbersequence
238
239!> Return a random integer on interval [0,0xffffffff]
240!!
241!! Code was based on getRandomInt() from the FMS implementation of the Mersenne Twistor
242integer function getrandomint(twister)
243 type(randomnumbersequence), intent(inout) :: twister !< The Mersenne Twister container
244
245 if (twister%currentElement >= blocksize) call nextstate(twister)
246 getrandomint = temper(twister%state(twister%currentElement))
247 twister%currentElement = twister%currentElement + 1
248
249end function getrandomint
250
251!> Return a random real number on interval [0,1]
252!!
253!! Code was based on getRandomReal() from the FMS implementation of the Mersenne Twistor
254double precision function getrandomreal(twister)
255 type(randomnumbersequence), intent(inout) :: twister
256 ! Local variables
257 integer :: localint
258
259 localint = getrandomint(twister)
260 if (localint < 0) then
261 getrandomreal = dble(localint + 2.0d0**32)/(2.0d0**32 - 1.0d0)
262 else
263 getrandomreal = dble(localint )/(2.0d0**32 - 1.0d0)
264 endif
265end function getrandomreal
266
267!> Merge bits of u and v
268integer function mixbits(u, v)
269 integer, intent(in) :: u !< An integer
270 integer, intent(in) :: v !< An integer
271
272 mixbits = ior(iand(u, umask), iand(v, lmask))
273end function mixbits
274
275!> Twist bits of u and v
276integer function twist(u, v)
277 integer, intent(in) :: u !< An integer
278 integer, intent(in) :: v !< An integer
279 ! Local variable
280 integer, parameter, dimension(0:1) :: t_matrix = (/ 0, matrix_a /)
281
282 twist = ieor(ishft(mixbits(u, v), -1), t_matrix(iand(v, 1)))
283 twist = ieor(ishft(mixbits(u, v), -1), t_matrix(iand(v, 1)))
284end function twist
285
286!> Update internal state of twister to the next state in the sequence
287subroutine nextstate(twister)
288 type(randomnumbersequence), intent(inout) :: twister !< Container for the Mersenne Twister
289 ! Local variables
290 integer :: k
291
292 do k = 0, blocksize - m - 1
293 twister%state(k) = ieor(twister%state(k + m), &
294 twist(twister%state(k), twister%state(k + 1)))
295 enddo
296 do k = blocksize - m, blocksize - 2
297 twister%state(k) = ieor(twister%state(k + m - blocksize), &
298 twist(twister%state(k), twister%state(k + 1)))
299 enddo
300 twister%state(blocksize - 1) = ieor(twister%state(m - 1), &
301 twist(twister%state(blocksize - 1), twister%state(0)))
302 twister%currentElement = 0
303end subroutine nextstate
304
305!> Tempering of bits in y
306elemental integer function temper(y)
307 integer, intent(in) :: y !< An integer
308 ! Local variables
309 integer :: x
310
311 x = ieor(y, ishft(y, -11))
312 x = ieor(x, iand(ishft(x, 7), tmaskb))
313 x = ieor(x, iand(ishft(x, 15), tmaskc))
314 temper = ieor(x, ishft(x, -18))
315end function temper
316
317!> Runs some statistical tests on the PRNG
318logical function random_unit_tests(verbose)
319 logical :: verbose !< True if results should be written to stdout
320 ! Local variables
321 type(prng) :: test_rng ! Generator
322 type(time_type) :: time ! Model time
323 real :: r1, r2, r3 ! Some random numbers and re-used work variables [nondim]
324 real :: mean, var, ar1, std ! Some statistics [nondim]
325 integer :: stdunit ! For messages
326 integer, parameter :: n_samples = 800
327 integer :: i, j, ni, nj
328 ! Fake being on a decomposed domain
329 type(hor_index_type), pointer :: hi => null() !< Not the real HI
330 real, dimension(:,:), allocatable :: r2d ! Random numbers [nondim]
331
332 ! Fake a decomposed domain
333 ni = 6
334 nj = 9
335 allocate(hi)
336 hi%isd = 0
337 hi%ied = ni+1
338 hi%jsd = 0
339 hi%jed = nj+1
340 hi%niglobal = ni
341 hi%njglobal = nj
342 hi%idg_offset = 0
343 hi%jdg_offset = 0
344
345 random_unit_tests = .false.
346 stdunit = stdout
347 write(stdunit,'(1x,a)') '==== MOM_random: random_unit_tests ======================='
348
349 if (verbose) write(stdunit,'(1x,"random: ",a)') '-- Time-based seeds ---------------------'
350 ! Check time-based seed generation
351 time = set_date(1903, 11, 21, 13, 47, 29)
352 i = seed_from_time(time)
354 test_fn(verbose, i==212584341, 'time seed 1903/11/21 13:47:29', ivalue=i)
355 time = set_date(1903, 11, 22, 13, 47, 29)
356 i = seed_from_time(time)
358 test_fn(verbose, i==212584342, 'time seed 1903/11/22 13:47:29', ivalue=i)
359 time = set_date(1903, 11, 21, 13, 47, 30)
360 i = seed_from_time(time)
362 test_fn(verbose, i==212596634, 'time seed 1903/11/21 13:47:30', ivalue=i)
363
364 if (verbose) write(stdunit,'(1x,"random: ",a)') '-- PRNG tests ---------------------------'
365 ! Generate a random number, r1
366 call random_0d_constructor(test_rng, time, 1)
367 r1 = random_01(test_rng)
369 test_fn(verbose, abs(r1-4.75310122e-2)<1.e-9, 'first call', r1)
370
371 ! Check that we get a different number, r2, on a second call
372 r2 = random_01(test_rng)
374 test_fn(verbose, abs(r2-2.71289742e-1)<1.e-9, 'consecutive test', r2)
375
376 ! Check that we can reproduce r1 by resetting the seed
377 call random_0d_constructor(test_rng, time, 1)
378 r2 = random_01(test_rng)
380 test_fn(verbose, abs(r2-r1)==0., 'reproduce test', r2)
381
382 ! Check that we get a different number, r2, with a different seed but same date
383 call random_0d_constructor(test_rng, time, 2)
384 r2 = random_01(test_rng)
386 test_fn(verbose, abs(r2-7.15508473e-1)<1.e-9, 'different seed test', r2)
387
388 ! Check that we get a different number, r2, for a different date but same seed
389 time = set_date(1903, 11, 21, 13, 0, 29)
390 call random_0d_constructor(test_rng, time, 1)
391 r2 = random_01(test_rng)
393 test_fn(verbose, abs(r2-9.56667163e-1)<1.e-9, 'different date test', r2)
394
395 if (verbose) write(stdunit,'(1x,"random: ",a)') '-- index-based seeds --------------------'
396 ! Check index-based seed
397 i = seed_from_index(hi,1,1)
398 random_unit_tests = random_unit_tests .or. test_fn(verbose, i==1, 'seed from index (1,1)', ivalue=i)
399 j = seed_from_index(hi,ni+1,1)
400 random_unit_tests = random_unit_tests .or. test_fn(verbose, j==i, 'seed from index (n+1,1)', ivalue=j)
401 i = seed_from_index(hi,ni,1)
402 random_unit_tests = random_unit_tests .or. test_fn(verbose, i==6, 'seed from index (n,1)', ivalue=i)
403 j = seed_from_index(hi,0,1)
404 random_unit_tests = random_unit_tests .or. test_fn(verbose, j==i, 'seed from index (0,1)', ivalue=j)
405 i = seed_from_index(hi,1,nj)
406 random_unit_tests = random_unit_tests .or. test_fn(verbose, i==49, 'seed from index (1,n)', ivalue=i)
407 j = seed_from_index(hi,ni,nj+1)
408 random_unit_tests = random_unit_tests .or. test_fn(verbose, j==i, 'seed from index (n,n+1)', ivalue=j)
409 i = seed_from_index(hi,ni,nj)
410 random_unit_tests = random_unit_tests .or. test_fn(verbose, i==54, 'seed from index (n,n)', ivalue=i)
411 j = seed_from_index(hi,1,nj+1)
412 random_unit_tests = random_unit_tests .or. test_fn(verbose, j==i, 'seed from index (1,n+1)', ivalue=j)
413
414 if (.not.random_unit_tests) write(stdunit,'(1x,a)') 'Passed unit tests'
415 ! The rest of these are not unit tests but statistical tests and as such
416 ! could fail for different sample sizes but happen to pass here.
417
418 ! Check statistics of large samples for uniform generator
419 mean = 0. ; var = 0. ; ar1 = 0. ; r2 = 0.
420 do i = 1, n_samples
421 r1 = random_01(test_rng) - 0.5
422 mean = mean + r1
423 var = var + r1**2
424 ar1 = ar1 + r1*r2
425 r2 = r1 ! Keep copy of last value
426 enddo
427 mean = mean / real(n_samples) ! Expected mean is 0
428 var = var / real(n_samples) ! Expected variance is 1/12
429 ar1 = ar1 / real(n_samples-1) ! Autocovariance
430 std = sqrt(var) ! Expected std is sqrt(1/12)
431 r2 = mean*sqrt(real(12*n_samples)) ! Normalized error in mean
432 r3 = std*sqrt(12.) ! Normalized standard deviation
433 r1 = ( ar1 * sqrt(real(n_samples-1)) ) / var
434 if (verbose) then
435 write(stdunit,'(1x,"random: ",a)') '-- Uniform -0.5 .. 0.5 generator --------'
436 write(stdunit,'(1x,"random: ",a,f12.9)') 'mean =',mean,'std =',std,'AR1 =',ar1
437 write(stdunit,'(1x,"random: ",a,f12.9)') 'norm. mean =',r2, &
438 'norm. std =',r3,'norm. AR1 =',r1
439 endif
441 test_fn(verbose, abs(r2)<2., &
442 'n>>1, mean within 2 sigma [uniform]', r2)
444 test_fn(verbose, abs(r3-1.)<1./sqrt(real(n_samples)), &
445 'n>>1, std ~ 1/sqrt(12) [uniform]', r3-1.)
447 test_fn(verbose, abs(r1)<2., &
448 'n>>1, AR1 < std/sqrt(n) [uniform]', r1)
449
450 ! Check statistics of large samples for normal generator
451 mean = 0. ; var = 0. ; ar1 = 0. ; r2 = 0.
452 do i = 1, n_samples
453 r1 = random_norm(test_rng)
454 mean = mean + r1
455 var = var + r1**2
456 ar1 = ar1 + r1*r2
457 r2 = r1 ! Keep copy of last value for AR calculation
458 enddo
459 mean = mean / real(n_samples)
460 var = var / real(n_samples)
461 ar1 = ar1 / real(n_samples)
462 std = sqrt(var)
463 r3 = 1./sqrt(real(n_samples)) ! Standard error of mean
464 r2 = mean*sqrt(real(n_samples)) ! Normalized error in mean
465 r3 = std ! Normalized standard deviation
466 r1 = ( ar1 * sqrt(real(n_samples-1)) ) / var
467 if (verbose) then
468 write(stdunit,'(1x,"random: ",a)') '-- Normal distribution generator --------'
469 write(stdunit,'(1x,"random: ",a,f12.9)') 'mean =',mean,'std =',std,'AR1 =',ar1
470 write(stdunit,'(1x,"random: ",a,f12.9)') 'norm. error in mean =',r2, &
471 'norm. standard deviation =',r3,'norm. AR1 =',r1
472 endif
474 test_fn(verbose, abs(r2)<2., &
475 'n>>1, mean within 2 sigma [norm]', r2)
477 test_fn(verbose, abs(r3-1.)<1./sqrt(real(n_samples)), &
478 'n>>1, std ~ 1 [norm]', r3-1.)
480 test_fn(verbose, abs(r1)<2., &
481 'n>>1, AR1 < std/sqrt(n) [norm]', r1)
482
483 if (verbose) write(stdunit,'(1x,"random: ",a)') '-- 2d PRNG ------------------------------'
484 ! Check 2d random number generator 0..1
485 allocate( r2d(hi%isd:hi%ied,hi%jsd:hi%jed) )
486 call random_2d_constructor(test_rng, hi, time, 123)
487 r2d(:,:) = -999. ! Use -9. to detect unset values
488 call random_2d_01(test_rng, hi, r2d)
489 if (any(abs(r2d(:,:)+999.)<=0.)) random_unit_tests=.true.
490 r1 = minval(r2d)
491 r2 = maxval(r2d)
492 random_unit_tests = random_unit_tests .or. test_fn(verbose, r1>=0., '2d all set', r1)
493 random_unit_tests = random_unit_tests .or. test_fn(verbose, r2<=1., '2d all valid', r2)
494 mean = sum( r2d(1:ni,1:nj) - 0.5 )/real(ni*nj)
495 var = sum( (r2d(1:ni,1:nj) - 0.5 - mean)**2 )/real(ni*nj)
496 std = sqrt(var)
497 r3 = 1./sqrt(real(12*ni*nj)) ! Standard error of mean
498 r2 = mean*sqrt(real(12*ni*nj)) ! Normalized error in mean
499 r3 = std*sqrt(12.) ! Normalized standard deviation
500 if (verbose) then
501 write(stdunit,'(1x,"random: ",a)') '2D uniform 0..1 generator'
502 write(stdunit,'(1x,"random: ",a,f12.9)') 'mean =',mean,'std =',std
503 write(stdunit,'(1x,"random: ",a,f12.9)') 'norm. error in mean =',r2
504 write(stdunit,'(1x,"random: ",a,f12.9)') 'norm. standard deviation =',r3
505 endif
507 test_fn(verbose, abs(r2)<2., &
508 '2d, mean within 2 sigma [uniform]', r2)
510 test_fn(verbose, abs(r3-1.)<1./sqrt(real(ni*nj)), &
511 '2d, std ~ 1/sqrt(12) [uniform]', r3-1.)
512 if (verbose) then
513 write(stdunit,'(1x,"random:")')
514 write(stdunit,'(1x,"random:",8f8.5)') r2d
515 write(stdunit,'(1x,"random:")')
516 endif
517
518 ! Check 2d normal random number generator
519 call random_2d_norm(test_rng, hi, r2d)
520 mean = sum( r2d(1:ni,1:nj) )/real(ni*nj)
521 var = sum( r2d(1:ni,1:nj)**2 )/real(ni*nj)
522 std = sqrt(var)
523 r3 = 1./sqrt(real(ni*nj)) ! Standard error of mean
524 r2 = mean*sqrt(real(ni*nj)) ! Normalized error in mean
525 r3 = std ! Normalized standard deviation
526 if (verbose) then
527 write(stdunit,'(1x,"random: ",a)') '2D normal generator'
528 write(stdunit,'(1x,"random: ",a,f12.9)') 'mean =',mean,'std =',std
529 write(stdunit,'(1x,"random: ",a,f12.9)') 'norm. error in mean =',r2
530 write(stdunit,'(1x,"random: ",a,f12.9)') 'norm. standard deviation =',r3
531 endif
533 test_fn(verbose, abs(r2)<2., &
534 '2d, mean within 2 sigma [norm]', r2)
536 test_fn(verbose, abs(r3-1.)<1./sqrt(real(ni*nj)), &
537 '2d, std ~ 1/sqrt(12) [norm]', r3-1.)
538
539 ! Clean up
540 deallocate(r2d)
541 deallocate(hi)
542
543 if (.not.random_unit_tests) write(stdunit,'(1x,a)') 'Passed statistical tests'
544
545end function random_unit_tests
546
547!> Convenience function for reporting result of test
548logical function test_fn(verbose, good, label, rvalue, ivalue)
549 logical, intent(in) :: verbose !< Verbosity
550 logical, intent(in) :: good !< True if pass, false otherwise
551 character(len=*), intent(in) :: label !< Label for messages
552 real, intent(in) :: rvalue !< Result of calculation [nondim]
553 integer, intent(in) :: ivalue !< Result of calculation
554 optional :: rvalue, ivalue
555
556 if (present(ivalue)) then
557 if (.not. good) then
558 write(stdout,'(1x,a,i10,1x,a,a)') 'random: result =',ivalue,label,' <------- FAIL!'
559 write(stderr,'(1x,a,i10,1x,a,a)') 'random: result =',ivalue,label,' <------- FAIL!'
560 elseif (verbose) then
561 write(stdout,'(1x,a,i10,1x,a)') 'random: result =',ivalue,label
562 endif
563 else
564 if (.not. good) then
565 write(stdout,'(1x,a,1pe15.8,1x,a,a)') 'random: result =',rvalue,label,' <------- FAIL!'
566 write(stderr,'(1x,a,1pe15.8,1x,a,a)') 'random: result =',rvalue,label,' <------- FAIL!'
567 elseif (verbose) then
568 write(stdout,'(1x,a,1pe15.8,1x,a)') 'random: result =',rvalue,label
569 endif
570 endif
571 test_fn = .not. good
572
573end function test_fn
574
575end module mom_random
576
577!> \namespace mom_random
578!!
579!! Provides MOM6 implementation of the Mersenne Twistor, copied from the FMS implementation
580!! which was originally written by Robert Pincus (Robert.Pincus@colorado.edu).
581!! We once used the FMS implementation directly but since random numers do not need to be
582!! infrastructure specific, and because MOM6 should be infrastructure agnostic, we have copied
583!! the parts of MT that we used here.
584!!
585!! Example usage:
586!! \code
587!! type(PRNG) :: rng
588!! real :: rn
589!! call random_0d_constructor(rng, Time, seed) ! Call this each time-step
590!! rn = random_01(rng)
591!! rn = random_norm(rng)
592!!
593!! type(PRNG) :: rng
594!! real, dimension(:,:) :: rn2d
595!! call random_2d_constructor(rng, HI, Time, seed) ! Call this each time-step
596!! call random_2d_01(rng, HI, rn2d)
597!! call random_2d_norm(rng, HI, rn2d)
598!!
599!! Note: reproducibility across restarts is implemented by using time-derived
600!! seeds to pass to the Mersenne twister. It is therefore important that any
601!! PRNG type be re-initialized each time-step.
602!! \endcode