[3880] | 1 | !**** *RANDOM_NUMBERS_MIX* - Portable Random Number Generator |
---|
| 2 | |
---|
| 3 | ! (C) Copyright 2002- ECMWF. |
---|
| 4 | ! |
---|
| 5 | ! This software is licensed under the terms of the Apache Licence Version 2.0 |
---|
| 6 | ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. |
---|
| 7 | ! |
---|
| 8 | ! In applying this licence, ECMWF does not waive the privileges and immunities |
---|
| 9 | ! granted to it by virtue of its status as an intergovernmental organisation |
---|
| 10 | ! nor does it submit to any jurisdiction. |
---|
| 11 | |
---|
| 12 | ! Purpose. |
---|
| 13 | ! -------- |
---|
| 14 | ! Generate machine-independent pseudo-random numbers |
---|
| 15 | |
---|
| 16 | !** Interface. |
---|
| 17 | ! ---------- |
---|
| 18 | ! CALL initialize_random_numbers (kseed, yd_stream) |
---|
| 19 | ! CALL uniform_distribution (px , yd_stream) |
---|
| 20 | ! CALL gaussian_distribution (px , yd_stream) |
---|
| 21 | ! CALL random_number_restartfile (fname, action) |
---|
| 22 | ! CALL wr_rangen_state (nunit) |
---|
| 23 | |
---|
| 24 | ! Explicit arguments : |
---|
| 25 | ! -------------------- |
---|
| 26 | ! kseed (input) : integer seed in the range [0,HUGE(kseed)] |
---|
| 27 | ! yd_stream (optional) : the state of the random number generator |
---|
| 28 | ! px (output) : array to receive random numbers in the range |
---|
| 29 | |
---|
| 30 | ! In the case of uniform_distribution, px has values in the range [0.0,1.0) |
---|
| 31 | |
---|
| 32 | ! Implicit arguments : |
---|
| 33 | ! -------------------- |
---|
| 34 | ! None |
---|
| 35 | |
---|
| 36 | ! Method. |
---|
| 37 | ! ------- |
---|
| 38 | ! Based loosly on ZUFALL (Petersen, 1994). |
---|
| 39 | |
---|
| 40 | ! The main difference between this generator and ZUFALL is that integer arithmetic |
---|
| 41 | ! is used. This ensures portability to vector machines that implement different |
---|
| 42 | ! real arithmetic. In particular, vector machines often implement non-IEEE |
---|
| 43 | ! arithmetic for their vector pipes. This routine will give identical results for |
---|
| 44 | ! any integer type with at least 32 bits. |
---|
| 45 | |
---|
| 46 | ! The generator is a lagged-Fibonacci generator: x(i) = x(i-p) + x(i-q) mod 2**m. |
---|
| 47 | ! Lagged-Fibonacci generators have very long repeat periods: (2**q -1) * 2**(m-1) |
---|
| 48 | ! (i.e about 2.85E191 for q=607, m=30). They pass most tests for randomness. |
---|
| 49 | |
---|
| 50 | ! p and q must be chosen carefully. Values from the following table are OK. |
---|
| 51 | ! Larger values give better random numbers, but smaller values are more |
---|
| 52 | ! cache-friendly. |
---|
| 53 | |
---|
| 54 | ! q p |
---|
| 55 | ! 9689 4187 |
---|
| 56 | ! 4423 2098 |
---|
| 57 | ! 2281 1029 |
---|
| 58 | ! 1279 418 |
---|
| 59 | ! 607 273 |
---|
| 60 | ! 521 168 |
---|
| 61 | ! 250 103 |
---|
| 62 | ! 127 63 |
---|
| 63 | ! 97 33 |
---|
| 64 | ! 55 24 |
---|
| 65 | ! 43 22 |
---|
| 66 | ! 31 13 |
---|
| 67 | ! 24 10 |
---|
| 68 | |
---|
| 69 | ! The initial q values of x are set using the binary shirt register method of |
---|
| 70 | ! Burns and Pryor 1999. |
---|
| 71 | |
---|
| 72 | ! Mascagni et al (1995) show how to choose different sets of initial values that |
---|
| 73 | ! are guaranteed to be drawn from different maximal-length cycles. This requires |
---|
| 74 | ! the initial values of x(1)...x(q) to be in "canonical form". Specifically, |
---|
| 75 | ! x(1) must be zero and all but a particular one or two values of x must be |
---|
| 76 | ! even. For q=607 and p=273, only one element (jpq-jps) must be odd. |
---|
| 77 | |
---|
| 78 | ! Externals. |
---|
| 79 | ! ---------- |
---|
| 80 | ! None |
---|
| 81 | |
---|
| 82 | ! Reference. |
---|
| 83 | ! ---------- |
---|
| 84 | ! Burns P.J. and Pryor D.V. 1999, |
---|
| 85 | ! Surface Radiative Transport at Large Scale via Monte Carlo. |
---|
| 86 | ! Annual Review of Heat Transfer, Vol 9. |
---|
| 87 | ! |
---|
| 88 | ! Petersen W.P., 1994, Lagged Fibonacci Series Random Number Generator |
---|
| 89 | ! for the NEC SX-3. International Journal of High Speed Computing |
---|
| 90 | ! Vol. 6, No. 3, pp387-398. |
---|
| 91 | ! |
---|
| 92 | ! Mascagni M., Cuccaro S.A., Pryor D.V., Robinson M.L., 1995, |
---|
| 93 | ! A Fast, High Quality and Reproducible Parallel Lagged-Fibonacci |
---|
| 94 | ! Pseudorandom Number Generator. Journal of Computational Physics |
---|
| 95 | ! Vol 119. pp211-219. |
---|
| 96 | |
---|
| 97 | ! Author. |
---|
| 98 | ! ------- |
---|
| 99 | ! Mike Fisher *ECMWF* |
---|
| 100 | |
---|
| 101 | ! Modifications. |
---|
| 102 | ! -------------- |
---|
| 103 | ! Original : 2002-09-25 |
---|
| 104 | ! Made parallel friendly: 2003-08-11 Robert Pincus |
---|
| 105 | ! M Leutbecher: 2004-05-10 restart capability |
---|
| 106 | ! M Fisher: 2005-03-30 replaced LCG initialization with shift register |
---|
| 107 | ! ------------------------------------------------------------------ |
---|
| 108 | |
---|
| 109 | #ifdef RS6K |
---|
| 110 | @PROCESS HOT(NOVECTOR) NOSTRICT |
---|
| 111 | #endif |
---|
| 112 | MODULE RANDOM_NUMBERS_MIX |
---|
| 113 | USE YOMHOOK, ONLY : LHOOK, DR_HOOK |
---|
| 114 | USE PARKIND1, ONLY : JPIM, JPRB |
---|
| 115 | |
---|
| 116 | IMPLICIT NONE |
---|
| 117 | |
---|
| 118 | SAVE |
---|
| 119 | |
---|
| 120 | PRIVATE |
---|
| 121 | PUBLIC RANDOMNUMBERSTREAM,WR_RANGEN_STATE, & |
---|
| 122 | & INITIALIZE_RANDOM_NUMBERS, UNIFORM_DISTRIBUTION, GAUSSIAN_DISTRIBUTION ! ,& |
---|
| 123 | ! & RANDOM_NUMBER_RESTARTFILE, |
---|
| 124 | |
---|
| 125 | INTEGER(KIND=JPIM), PARAMETER :: JPP=273, JPQ=607, JPS=105 |
---|
| 126 | INTEGER(KIND=JPIM), PARAMETER :: JPMM=30 |
---|
| 127 | INTEGER(KIND=JPIM), PARAMETER :: JPM=2**JPMM |
---|
| 128 | INTEGER(KIND=JPIM), PARAMETER :: JPNUMSPLIT=(JPQ-2)/(JPP-1) |
---|
| 129 | INTEGER(KIND=JPIM), PARAMETER :: JPLENSPLIT=(JPQ-JPP+JPNUMSPLIT-1)/JPNUMSPLIT |
---|
| 130 | INTEGER(KIND=JPIM), PARAMETER :: INITVALUE = 12345678 |
---|
| 131 | |
---|
| 132 | TYPE RANDOMNUMBERSTREAM |
---|
| 133 | PRIVATE |
---|
| 134 | INTEGER(KIND=JPIM) :: IUSED |
---|
| 135 | INTEGER(KIND=JPIM) :: INITTEST ! Should initialize to zero, but can't in F90 |
---|
| 136 | INTEGER(KIND=JPIM), DIMENSION(JPQ) :: IX |
---|
| 137 | REAL(KIND=JPRB) :: ZRM |
---|
| 138 | END TYPE RANDOMNUMBERSTREAM |
---|
| 139 | |
---|
| 140 | CONTAINS |
---|
| 141 | !------------------------------------------------------------------------------- |
---|
| 142 | SUBROUTINE INITIALIZE_RANDOM_NUMBERS (KSEED, YD_STREAM) |
---|
| 143 | !------------------------------------------------------------------------------- |
---|
| 144 | ! Initialize fibgen |
---|
| 145 | !------------------------------------------------------------------------------- |
---|
| 146 | INTEGER(KIND=JPIM), INTENT(IN ) :: KSEED |
---|
| 147 | TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM |
---|
| 148 | |
---|
| 149 | INTEGER, PARAMETER :: JPMASK=123459876 |
---|
| 150 | INTEGER(KIND=JPIM), PARAMETER :: JPWARMUP_SHFT=64, JPWARMUP_LFG=999 |
---|
| 151 | INTEGER(KIND=JPIM) :: IDUM,JK,JJ,JBIT |
---|
| 152 | REAL(KIND=JPRB), DIMENSION(JPWARMUP_LFG) :: ZWARMUP |
---|
| 153 | |
---|
| 154 | !------------------------------------------------------------------------------- |
---|
| 155 | ! Initialize the buffer using a binary shift register (Burns and Pryor, 1999). |
---|
| 156 | ! The Galois representation is used for the shift register as it is more |
---|
| 157 | ! efficient than the Fibonacci representation. The magic numbers 31 and 87 |
---|
| 158 | ! define the shift register primitive polynomial=(32,7,5,3,2,1,0). |
---|
| 159 | ! |
---|
| 160 | ! To ensure that different seeds produce distinct initial buffer states in |
---|
| 161 | ! canonical form, bits 0...jpmm-2 of the initial seed (after XORing with jpmask |
---|
| 162 | ! and spinning up using the linear congruential generator) are used to construct |
---|
| 163 | ! x(2), and the remaining bits are used to construct x(jpq). |
---|
| 164 | !------------------------------------------------------------------------------- |
---|
| 165 | |
---|
| 166 | REAL(KIND=JPRB) :: ZHOOK_HANDLE |
---|
| 167 | IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:INITIALIZE_RANDOM_NUMBERS',0,ZHOOK_HANDLE) |
---|
| 168 | IDUM = ABS(IEOR(KSEED,JPMASK)) |
---|
| 169 | IF (IDUM==0) IDUM=JPMASK |
---|
| 170 | |
---|
| 171 | DO JJ=1,JPWARMUP_SHFT |
---|
| 172 | IF (BTEST(IDUM,31)) THEN |
---|
| 173 | IDUM=IBSET(ISHFT(IEOR(IDUM,87),1),0) |
---|
| 174 | ELSE |
---|
| 175 | IDUM=IBCLR(ISHFT(IDUM,1),0) |
---|
| 176 | ENDIF |
---|
| 177 | ENDDO |
---|
| 178 | |
---|
| 179 | YD_STREAM%IX(1:JPQ-1)= 0 |
---|
| 180 | YD_STREAM%IX(2) = ISHFT(IBITS(IDUM,0,JPMM-1),1) |
---|
| 181 | YD_STREAM%IX(JPQ) = IBITS(IDUM,JPMM-1,BIT_SIZE(IDUM)+1-JPMM) |
---|
| 182 | |
---|
| 183 | DO JBIT=1,JPMM-1 |
---|
| 184 | DO JJ=3,JPQ-1 |
---|
| 185 | IF (BTEST(IDUM,31)) THEN |
---|
| 186 | IDUM=IBSET(ISHFT(IEOR(IDUM,87),1),0) |
---|
| 187 | YD_STREAM%IX(JJ)=IBSET(YD_STREAM%IX(JJ),JBIT) |
---|
| 188 | ELSE |
---|
| 189 | IDUM=IBCLR(ISHFT(IDUM,1),0) |
---|
| 190 | ENDIF |
---|
| 191 | ENDDO |
---|
| 192 | ENDDO |
---|
| 193 | |
---|
| 194 | YD_STREAM%IX(JPQ-JPS) = IBSET(YD_STREAM%IX(JPQ-JPS),0) |
---|
| 195 | |
---|
| 196 | !------------------------------------------------------------------------------- |
---|
| 197 | ! Initialize some constants |
---|
| 198 | !------------------------------------------------------------------------------- |
---|
| 199 | |
---|
| 200 | YD_STREAM%IUSED=JPQ |
---|
| 201 | YD_STREAM%ZRM=1.0_JPRB/REAL(JPM,JPRB) |
---|
| 202 | |
---|
| 203 | !------------------------------------------------------------------------------- |
---|
| 204 | ! Check the calculation of jpnumsplit and jplensplit. |
---|
| 205 | !------------------------------------------------------------------------------- |
---|
| 206 | |
---|
| 207 | IF (JPP+JPNUMSPLIT*JPLENSPLIT < JPQ) THEN |
---|
| 208 | CALL ABOR1 ('initialize_random_numbers: upper limit of last loop < jpq') |
---|
| 209 | ENDIF |
---|
| 210 | |
---|
| 211 | IF (JPLENSPLIT >=JPP) THEN |
---|
| 212 | CALL ABOR1 ('initialize_random_numbers: loop length > jpp') |
---|
| 213 | ENDIF |
---|
| 214 | |
---|
| 215 | IF (JPNUMSPLIT>1) THEN |
---|
| 216 | IF ((JPQ-JPP+JPNUMSPLIT-2)/(JPNUMSPLIT-1) < JPP) THEN |
---|
| 217 | CALL ABOR1 ('initialize_random_numbers: jpnumsplit is bigger than necessary') |
---|
| 218 | ENDIF |
---|
| 219 | ENDIF |
---|
| 220 | |
---|
| 221 | !------------------------------------------------------------------------------- |
---|
| 222 | ! Set initTest to show that the stream is initialized. |
---|
| 223 | !------------------------------------------------------------------------------- |
---|
| 224 | |
---|
| 225 | YD_STREAM%INITTEST = INITVALUE |
---|
| 226 | |
---|
| 227 | !------------------------------------------------------------------------------- |
---|
| 228 | ! Warm up the generator. |
---|
| 229 | !------------------------------------------------------------------------------- |
---|
| 230 | |
---|
| 231 | CALL UNIFORM_DISTRIBUTION (ZWARMUP, YD_STREAM) |
---|
| 232 | |
---|
| 233 | IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:INITIALIZE_RANDOM_NUMBERS',1,ZHOOK_HANDLE) |
---|
| 234 | END SUBROUTINE INITIALIZE_RANDOM_NUMBERS |
---|
| 235 | |
---|
| 236 | !@PROCESS HOT NOSTRICT |
---|
| 237 | SUBROUTINE UNIFORM_DISTRIBUTION (PX,YD_STREAM) |
---|
| 238 | !-------------------------------------------------------------------------------- |
---|
| 239 | ! Generate uniformly distributed random numbers in the range 0.0<= px < 1.0 |
---|
| 240 | !-------------------------------------------------------------------------------- |
---|
| 241 | INTEGER(KIND=JPIM), PARAMETER :: IVAR=Z"3FFFFFFF" |
---|
| 242 | TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM |
---|
| 243 | REAL(KIND=JPRB), DIMENSION(:), INTENT( OUT) :: PX |
---|
| 244 | |
---|
| 245 | INTEGER(KIND=JPIM) :: JJ, JK, IN, IFILLED |
---|
| 246 | |
---|
| 247 | ! This test is a little dirty but Fortran 90 doesn't allow for the initialization |
---|
| 248 | ! of components of derived types. |
---|
| 249 | REAL(KIND=JPRB) :: ZHOOK_HANDLE |
---|
| 250 | ! DR_HOOK removed to reduce overhead |
---|
| 251 | ! IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:UNIFORM_DISTRIBUTION',0,ZHOOK_HANDLE) |
---|
| 252 | IF(YD_STREAM%INITTEST /= INITVALUE) & |
---|
| 253 | & CALL ABOR1 ('uniform_distribution called before initialize_random_numbers') |
---|
| 254 | |
---|
| 255 | !-------------------------------------------------------------------------------- |
---|
| 256 | ! Copy numbers that were generated during the last call, but not used. |
---|
| 257 | !-------------------------------------------------------------------------------- |
---|
| 258 | |
---|
| 259 | IN=SIZE(PX) |
---|
| 260 | IFILLED=0 |
---|
| 261 | |
---|
| 262 | DO JJ=YD_STREAM%IUSED+1,MIN(JPQ,IN+YD_STREAM%IUSED) |
---|
| 263 | PX(JJ-YD_STREAM%IUSED) = YD_STREAM%IX(JJ)*YD_STREAM%ZRM |
---|
| 264 | IFILLED=IFILLED+1 |
---|
| 265 | ENDDO |
---|
| 266 | |
---|
| 267 | YD_STREAM%IUSED=YD_STREAM%IUSED+IFILLED |
---|
| 268 | |
---|
| 269 | IF (IFILLED==IN) THEN |
---|
| 270 | ! IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:UNIFORM_DISTRIBUTION',1,ZHOOK_HANDLE) |
---|
| 271 | ! DR_HOOK removed to reduce overhead |
---|
| 272 | RETURN |
---|
| 273 | ENDIF |
---|
| 274 | |
---|
| 275 | !-------------------------------------------------------------------------------- |
---|
| 276 | ! Generate batches of jpq numbers until px has been filled |
---|
| 277 | !-------------------------------------------------------------------------------- |
---|
| 278 | |
---|
| 279 | DO WHILE (IFILLED<IN) |
---|
| 280 | |
---|
| 281 | !-------------------------------------------------------------------------------- |
---|
| 282 | ! Generate jpq numbers in vectorizable loops. The first loop is length jpp. The |
---|
| 283 | ! remaining jpq-jpp elements are calculated in loops of length shorter than jpp. |
---|
| 284 | !-------------------------------------------------------------------------------- |
---|
| 285 | |
---|
| 286 | !OCL NOVREC |
---|
| 287 | DO JJ=1,JPP |
---|
| 288 | ! yd_stream%ix(jj) = yd_stream%ix(jj) + yd_stream%ix(jj-jpp+jpq) |
---|
| 289 | ! if (yd_stream%ix(jj)>=jpm) yd_stream%ix(jj) = yd_stream%ix(jj)-jpm |
---|
| 290 | YD_STREAM%IX(JJ) = IAND(IVAR,YD_STREAM%IX(JJ) + YD_STREAM%IX(JJ-JPP+JPQ)) |
---|
| 291 | ENDDO |
---|
| 292 | |
---|
| 293 | DO JK=1,JPNUMSPLIT |
---|
| 294 | !OCL NOVREC |
---|
| 295 | DO JJ=1+JPP+(JK-1)*JPLENSPLIT,MIN(JPQ,JPP+JK*JPLENSPLIT) |
---|
| 296 | ! yd_stream%ix(jj) = yd_stream%ix(jj) + yd_stream%ix(jj-jpp) |
---|
| 297 | ! if (yd_stream%ix(jj)>=jpm) yd_stream%ix(jj) = yd_stream%ix(jj)-jpm |
---|
| 298 | YD_STREAM%IX(JJ) = IAND(IVAR,YD_STREAM%IX(JJ) + YD_STREAM%IX(JJ-JPP)) |
---|
| 299 | ENDDO |
---|
| 300 | ENDDO |
---|
| 301 | |
---|
| 302 | YD_STREAM%IUSED = MIN(JPQ,IN-IFILLED) |
---|
| 303 | PX(IFILLED+1:IFILLED+YD_STREAM%IUSED) = YD_STREAM%IX(1:YD_STREAM%IUSED)*YD_STREAM%ZRM |
---|
| 304 | IFILLED = IFILLED+YD_STREAM%IUSED |
---|
| 305 | ENDDO |
---|
| 306 | |
---|
| 307 | !IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:UNIFORM_DISTRIBUTION',1,ZHOOK_HANDLE) |
---|
| 308 | ! DR_HOOK removed to reduce overhead |
---|
| 309 | END SUBROUTINE UNIFORM_DISTRIBUTION |
---|
| 310 | !------------------------------------------------------------------------------- |
---|
| 311 | SUBROUTINE GAUSSIAN_DISTRIBUTION (PX, YD_STREAM) |
---|
| 312 | TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM |
---|
| 313 | REAL(KIND=JPRB), INTENT( OUT) :: PX(:) |
---|
| 314 | !-------------------------------------------------------------------------------- |
---|
| 315 | ! Generate normally-distributed random numbers using the Box-Muller method. |
---|
| 316 | ! |
---|
| 317 | ! NB: this routine does not use buffering. This means that the following calls: |
---|
| 318 | ! call gaussian_distribution (zx(1:k)) |
---|
| 319 | ! call gaussian_distribution (zx(k+1:n)) |
---|
| 320 | ! will produce different numbers for elements k+1 onwards than the single call: |
---|
| 321 | ! call gaussian_distribution (zx(1:n)) |
---|
| 322 | !-------------------------------------------------------------------------------- |
---|
| 323 | |
---|
| 324 | INTEGER(KIND=JPIM) :: ILEN, J |
---|
| 325 | REAL(KIND=JPRB) :: ZFAC, ZTWOPI |
---|
| 326 | REAL(KIND=JPRB) :: ZX(SIZE(PX)+1) |
---|
| 327 | |
---|
| 328 | !-------------------------------------------------------------------------------- |
---|
| 329 | ! Generate uniform random points in the range [0,1) |
---|
| 330 | !-------------------------------------------------------------------------------- |
---|
| 331 | |
---|
| 332 | REAL(KIND=JPRB) :: ZHOOK_HANDLE |
---|
| 333 | IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:GAUSSIAN_DISTRIBUTION',0,ZHOOK_HANDLE) |
---|
| 334 | CALL UNIFORM_DISTRIBUTION (ZX, YD_STREAM) |
---|
| 335 | |
---|
| 336 | !-------------------------------------------------------------------------------- |
---|
| 337 | ! Generate gaussian deviates using Box-Muller method |
---|
| 338 | !-------------------------------------------------------------------------------- |
---|
| 339 | |
---|
| 340 | ZTWOPI = 8.0_JPRB*ATAN(1.0_JPRB) |
---|
| 341 | ILEN=SIZE(PX) |
---|
| 342 | |
---|
| 343 | DO J=1,ILEN-1,2 |
---|
| 344 | ZFAC = SQRT(-2.0_JPRB*LOG(1.0_JPRB-ZX(J))) |
---|
| 345 | PX(J ) = ZFAC*COS(ZTWOPI*ZX(J+1)) |
---|
| 346 | PX(J+1) = ZFAC*SIN(ZTWOPI*ZX(J+1)) |
---|
| 347 | ENDDO |
---|
| 348 | |
---|
| 349 | !-------------------------------------------------------------------------------- |
---|
| 350 | ! Generate the last point if ilen is odd |
---|
| 351 | !-------------------------------------------------------------------------------- |
---|
| 352 | |
---|
| 353 | IF (MOD(ILEN,2) /= 0) THEN |
---|
| 354 | ZFAC = SQRT(-2.0_JPRB*LOG(1.0_JPRB-ZX(ILEN))) |
---|
| 355 | PX(ILEN) = ZFAC*COS(ZTWOPI*ZX(ILEN+1)) |
---|
| 356 | ENDIF |
---|
| 357 | |
---|
| 358 | IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:GAUSSIAN_DISTRIBUTION',1,ZHOOK_HANDLE) |
---|
| 359 | END SUBROUTINE GAUSSIAN_DISTRIBUTION |
---|
| 360 | !------------------------------------------------------------------------------- |
---|
| 361 | !!$SUBROUTINE RANDOM_NUMBER_RESTARTFILE( CDFNAME, CDACTION,YD_STREAM ) |
---|
| 362 | !!$!-------------------------------------------------------------------------------- |
---|
| 363 | !!$! |
---|
| 364 | !!$! read (cdaction='r') or write (cdaction='w') restart file |
---|
| 365 | !!$! for random number generator |
---|
| 366 | !!$! |
---|
| 367 | !!$!-------------------------------------------------------------------------------- |
---|
| 368 | !!$CHARACTER (LEN=*), INTENT(IN) :: CDFNAME |
---|
| 369 | !!$CHARACTER (LEN=1 ), INTENT(IN) :: CDACTION |
---|
| 370 | !!$TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM |
---|
| 371 | !!$ |
---|
| 372 | !!$INTEGER(KIND=JPIM) :: IUNIT, IRET, IBYTES_IN_JPIM |
---|
| 373 | !!$ |
---|
| 374 | !!$REAL(KIND=JPRB) :: ZHOOK_HANDLE |
---|
| 375 | !!$IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:RANDOM_NUMBER_RESTARTFILE',0,ZHOOK_HANDLE) |
---|
| 376 | !!$IBYTES_IN_JPIM= CEILING(REAL(BIT_SIZE(YD_STREAM%IUSED))/8.0_JPRB - TINY(1.0_JPRB)) |
---|
| 377 | !!$ |
---|
| 378 | !!$IF (IBYTES_IN_JPIM /= 4) THEN |
---|
| 379 | !!$ CALL ABOR1('random_number_restartfile: number of bytes for JPIM is not 4 ') |
---|
| 380 | !!$ENDIF |
---|
| 381 | !!$ |
---|
| 382 | !!$CALL PBOPEN(IUNIT, CDFNAME, CDACTION, IRET) |
---|
| 383 | !!$IF (IRET /= 0) THEN |
---|
| 384 | !!$ CALL ABOR1('random_number_restartfile: PBOPEN FAILED opening '//CDFNAME) |
---|
| 385 | !!$ENDIF |
---|
| 386 | !!$ |
---|
| 387 | !!$ |
---|
| 388 | !!$IF (CDACTION=='r' .OR. CDACTION=='R') THEN |
---|
| 389 | !!$ CALL PBREAD(IUNIT, YD_STREAM%IX, IBYTES_IN_JPIM*JPQ, IRET) |
---|
| 390 | !!$ IF (IRET < 0) THEN |
---|
| 391 | !!$ CALL ABOR1('random_number_restartfile: PBREAD could not read ix from '//CDFNAME) |
---|
| 392 | !!$ ENDIF |
---|
| 393 | !!$ CALL PBREAD(IUNIT, YD_STREAM%IUSED, IBYTES_IN_JPIM , IRET) |
---|
| 394 | !!$ IF (IRET < 0) THEN |
---|
| 395 | !!$ CALL ABOR1('random_number_restartfile: PBREAD could not read iused from '//CDFNAME) |
---|
| 396 | !!$ ENDIF |
---|
| 397 | !!$ |
---|
| 398 | !!$! l_initialized = .TRUE. |
---|
| 399 | !!$ YD_STREAM%INITTEST = INITVALUE |
---|
| 400 | !!$ YD_STREAM%ZRM=1.0_JPRB/REAL(JPM,JPRB) |
---|
| 401 | !!$ELSEIF(CDACTION=='w' .OR. CDACTION=='W') THEN |
---|
| 402 | !!$ CALL PBWRITE(IUNIT, YD_STREAM%IX, IBYTES_IN_JPIM*JPQ, IRET) |
---|
| 403 | !!$ IF (IRET < 0) THEN |
---|
| 404 | !!$ CALL ABOR1('random_number_restartfile: PBWRITE could not write ix on '//CDFNAME) |
---|
| 405 | !!$ ENDIF |
---|
| 406 | !!$ CALL PBWRITE(IUNIT, YD_STREAM%IUSED, IBYTES_IN_JPIM , IRET) |
---|
| 407 | !!$ IF (IRET < 0) THEN |
---|
| 408 | !!$ CALL ABOR1('random_number_restartfile: PBWRITE could not write iused on '//CDFNAME) |
---|
| 409 | !!$ ENDIF |
---|
| 410 | !!$ |
---|
| 411 | !!$ELSE |
---|
| 412 | !!$ CALL ABOR1 ('random_number_restartfile: cdaction = '//CDACTION//' is undefined.') |
---|
| 413 | !!$ENDIF |
---|
| 414 | !!$ |
---|
| 415 | !!$CALL PBCLOSE(IUNIT, IRET) |
---|
| 416 | !!$IF (IRET /= 0) THEN |
---|
| 417 | !!$ CALL ABOR1('random_number_restartfile: PBCLOSE FAILED closing '//CDFNAME) |
---|
| 418 | !!$ENDIF |
---|
| 419 | !!$ |
---|
| 420 | !!$IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:RANDOM_NUMBER_RESTARTFILE',1,ZHOOK_HANDLE) |
---|
| 421 | !!$END SUBROUTINE RANDOM_NUMBER_RESTARTFILE |
---|
| 422 | |
---|
| 423 | |
---|
| 424 | SUBROUTINE WR_RANGEN_STATE( KUNIT, YD_STREAM ) |
---|
| 425 | !-------------------------------------------------------------------------------- |
---|
| 426 | ! write state of random number generator to unit kunit |
---|
| 427 | !-------------------------------------------------------------------------------- |
---|
| 428 | INTEGER(KIND=JPIM), INTENT(IN) :: KUNIT |
---|
| 429 | TYPE(RANDOMNUMBERSTREAM), INTENT(IN) :: YD_STREAM |
---|
| 430 | |
---|
| 431 | REAL(KIND=JPRB) :: ZHOOK_HANDLE |
---|
| 432 | IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:WR_RANGEN_STATE',0,ZHOOK_HANDLE) |
---|
| 433 | WRITE( KUNIT, * ) 'module random_numbers_mix, generator state is' |
---|
| 434 | WRITE( KUNIT, '(8I10)') YD_STREAM%IX |
---|
| 435 | WRITE( KUNIT, '(I10)') YD_STREAM%IUSED |
---|
| 436 | |
---|
| 437 | IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:WR_RANGEN_STATE',1,ZHOOK_HANDLE) |
---|
| 438 | END SUBROUTINE WR_RANGEN_STATE |
---|
| 439 | |
---|
| 440 | END MODULE RANDOM_NUMBERS_MIX |
---|