1 2! KGEN-generated Fortran source file 3! 4! Filename : mcica_random_numbers.f90 5! Generated at: 2015-07-07 00:48:25 6! KGEN version: 0.4.13 7 8 9 10 MODULE mersennetwister 11 USE kgen_utils_mod, ONLY : kgen_dp, check_t, kgen_init_check, kgen_print_check 12 ! ------------------------------------------------------------- 13 USE shr_kind_mod, ONLY: r8 => shr_kind_r8 14 ! use parkind, only : jpim, jprb 15 IMPLICIT NONE 16 PRIVATE 17 ! Algorithm parameters 18 ! ------- 19 ! Period parameters 20 INTEGER, parameter :: blocksize = 624 21 INTEGER, parameter :: lmask = 2147483647 22 INTEGER, parameter :: umask = (-lmask) - 1 23 INTEGER, parameter :: m = 397 24 INTEGER, parameter :: matrix_a = -1727483681 25 ! constant vector a (0x9908b0dfUL) 26 ! least significant r bits (0x7fffffffUL) 27 ! most significant w-r bits (0x80000000UL) 28 ! Tempering parameters 29 INTEGER, parameter :: tmaskb= -1658038656 30 INTEGER, parameter :: tmaskc= -272236544 ! (0x9d2c5680UL) 31 ! (0xefc60000UL) 32 ! ------- 33 ! The type containing the state variable 34 TYPE randomnumbersequence 35 INTEGER :: currentelement ! = blockSize 36 INTEGER, dimension(0:blocksize -1) :: state ! = 0 37 END TYPE randomnumbersequence 38 39 INTERFACE new_randomnumbersequence 40 MODULE PROCEDURE initialize_scalar, initialize_vector 41 END INTERFACE new_randomnumbersequence 42 PUBLIC randomnumbersequence 43 PUBLIC new_randomnumbersequence, getrandomreal, getrandomint 44 ! ------------------------------------------------------------- 45 46 ! read interface 47 PUBLIC kgen_read 48 INTERFACE kgen_read 49 MODULE PROCEDURE kgen_read_randomnumbersequence 50 END INTERFACE kgen_read 51 52 PUBLIC kgen_verify 53 INTERFACE kgen_verify 54 MODULE PROCEDURE kgen_verify_randomnumbersequence 55 END INTERFACE kgen_verify 56 57 CONTAINS 58 59 ! write subroutines 60 ! No module extern variables 61 SUBROUTINE kgen_read_randomnumbersequence(var, kgen_unit, printvar) 62 INTEGER, INTENT(IN) :: kgen_unit 63 CHARACTER(*), INTENT(IN), OPTIONAL :: printvar 64 TYPE(randomnumbersequence), INTENT(out) :: var 65 READ(UNIT=kgen_unit) var%currentelement 66 IF ( PRESENT(printvar) ) THEN 67 print *, "** " // printvar // "%currentelement **", var%currentelement 68 END IF 69 READ(UNIT=kgen_unit) var%state 70 IF ( PRESENT(printvar) ) THEN 71 print *, "** " // printvar // "%state **", var%state 72 END IF 73 END SUBROUTINE 74 SUBROUTINE kgen_verify_randomnumbersequence(varname, check_status, var, ref_var) 75 CHARACTER(*), INTENT(IN) :: varname 76 TYPE(check_t), INTENT(INOUT) :: check_status 77 TYPE(check_t) :: dtype_check_status 78 TYPE(randomnumbersequence), INTENT(IN) :: var, ref_var 79 80 check_status%numTotal = check_status%numTotal + 1 81 CALL kgen_init_check(dtype_check_status) 82 CALL kgen_verify_integer("currentelement", dtype_check_status, var%currentelement, ref_var%currentelement) 83 CALL kgen_verify_integer_4_dim1("state", dtype_check_status, var%state, ref_var%state) 84 IF ( dtype_check_status%numTotal == dtype_check_status%numIdentical ) THEN 85 check_status%numIdentical = check_status%numIdentical + 1 86 ELSE IF ( dtype_check_status%numFatal > 0 ) THEN 87 check_status%numFatal = check_status%numFatal + 1 88 ELSE IF ( dtype_check_status%numWarning > 0 ) THEN 89 check_status%numWarning = check_status%numWarning + 1 90 END IF 91 END SUBROUTINE 92 SUBROUTINE kgen_verify_integer( varname, check_status, var, ref_var) 93 character(*), intent(in) :: varname 94 type(check_t), intent(inout) :: check_status 95 integer, intent(in) :: var, ref_var 96 check_status%numTotal = check_status%numTotal + 1 97 IF ( var == ref_var ) THEN 98 check_status%numIdentical = check_status%numIdentical + 1 99 if(check_status%verboseLevel > 1) then 100 WRITE(*,*) 101 WRITE(*,*) trim(adjustl(varname)), " is IDENTICAL( ", var, " )." 102 endif 103 ELSE 104 if(check_status%verboseLevel > 0) then 105 WRITE(*,*) 106 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 107 if(check_status%verboseLevel > 2) then 108 WRITE(*,*) "KERNEL: ", var 109 WRITE(*,*) "REF. : ", ref_var 110 end if 111 end if 112 check_status%numFatal = check_status%numFatal + 1 113 END IF 114 END SUBROUTINE kgen_verify_integer 115 116 SUBROUTINE kgen_verify_integer_4_dim1( varname, check_status, var, ref_var) 117 character(*), intent(in) :: varname 118 type(check_t), intent(inout) :: check_status 119 integer, intent(in), DIMENSION(:) :: var, ref_var 120 check_status%numTotal = check_status%numTotal + 1 121 IF ( ALL( var == ref_var ) ) THEN 122 123 check_status%numIdentical = check_status%numIdentical + 1 124 if(check_status%verboseLevel > 1) then 125 WRITE(*,*) 126 WRITE(*,*) "All elements of ", trim(adjustl(varname)), " are IDENTICAL." 127 !WRITE(*,*) "KERNEL: ", var 128 !WRITE(*,*) "REF. : ", ref_var 129 IF ( ALL( var == 0 ) ) THEN 130 if(check_status%verboseLevel > 2) then 131 WRITE(*,*) "All values are zero." 132 end if 133 END IF 134 end if 135 ELSE 136 if(check_status%verboseLevel > 0) then 137 WRITE(*,*) 138 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 139 WRITE(*,*) count( var /= ref_var), " of ", size( var ), " elements are different." 140 end if 141 142 check_status%numFatal = check_status%numFatal+1 143 END IF 144 END SUBROUTINE kgen_verify_integer_4_dim1 145 146 ! ------------------------------------------------------------- 147 ! Private functions 148 ! --------------------------- 149 150 FUNCTION mixbits(u, v) 151 INTEGER, intent( in) :: u 152 INTEGER, intent( in) :: v 153 INTEGER :: mixbits 154 mixbits = ior(iand(u, UMASK), iand(v, LMASK)) 155 END FUNCTION mixbits 156 ! --------------------------- 157 158 FUNCTION twist(u, v) 159 INTEGER, intent( in) :: u 160 INTEGER, intent( in) :: v 161 INTEGER :: twist 162 ! Local variable 163 INTEGER, parameter, dimension(0:1) :: t_matrix = (/ 0, matrix_a /) 164 twist = ieor(ishft(mixbits(u, v), -1), t_matrix(iand(v, 1))) 165 twist = ieor(ishft(mixbits(u, v), -1), t_matrix(iand(v, 1))) 166 END FUNCTION twist 167 ! --------------------------- 168 169 SUBROUTINE nextstate(twister) 170 TYPE(randomnumbersequence), intent(inout) :: twister 171 ! Local variables 172 INTEGER :: k 173 do k = 0, blockSize - M - 1 174 twister%state(k) = ieor(twister%state(k + M), & 175 twist(twister%state(k), twister%state(k + 1))) 176 end do 177 do k = blockSize - M, blockSize - 2 178 twister%state(k) = ieor(twister%state(k + M - blockSize), & 179 twist(twister%state(k), twister%state(k + 1))) 180 end do 181 twister%state(blockSize - 1) = ieor(twister%state(M - 1), & 182 twist(twister%state(blockSize - 1), twister%state(0))) 183 twister%currentElement = 0 184 END SUBROUTINE nextstate 185 ! --------------------------- 186 187 elemental FUNCTION temper(y) 188 INTEGER, intent(in) :: y 189 INTEGER :: temper 190 INTEGER :: x 191 ! Tempering 192 x = ieor(y, ishft(y, -11)) 193 x = ieor(x, iand(ishft(x, 7), TMASKB)) 194 x = ieor(x, iand(ishft(x, 15), TMASKC)) 195 temper = ieor(x, ishft(x, -18)) 196 END FUNCTION temper 197 ! ------------------------------------------------------------- 198 ! Public (but hidden) functions 199 ! -------------------- 200 201 FUNCTION initialize_scalar(seed) RESULT ( twister ) 202 INTEGER, intent(in ) :: seed 203 TYPE(randomnumbersequence) :: twister 204 INTEGER :: i 205 ! See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. In the previous versions, 206 ! MSBs of the seed affect only MSBs of the array state[]. 207 ! 2002/01/09 modified by Makoto Matsumoto 208 twister%state(0) = iand(seed, -1) 209 do i = 1, blockSize - 1 ! ubound(twister%state) ! ubound(twister%state) 210 twister%state(i) = 1812433253 * ieor(twister%state(i-1), & 211 ishft(twister%state(i-1), -30)) + i 212 twister%state(i) = iand(twister%state(i), -1) ! for >32 bit machines ! for >32 bit machines 213 end do 214 twister%currentElement = blockSize 215 END FUNCTION initialize_scalar 216 ! ------------------------------------------------------------- 217 218 FUNCTION initialize_vector(seed) RESULT ( twister ) 219 INTEGER, dimension(0:), intent(in) :: seed 220 TYPE(randomnumbersequence) :: twister 221 INTEGER :: nwraps 222 INTEGER :: nfirstloop 223 INTEGER :: k 224 INTEGER :: i 225 INTEGER :: j 226 nWraps = 0 227 twister = initialize_scalar(19650218) 228 nFirstLoop = max(blockSize, size(seed)) 229 do k = 1, nFirstLoop 230 i = mod(k + nWraps, blockSize) 231 j = mod(k - 1, size(seed)) 232 if(i == 0) then 233 twister%state(i) = twister%state(blockSize - 1) 234 twister%state(1) = ieor(twister%state(1), & 235 ieor(twister%state(1-1), & 236 ishft(twister%state(1-1), -30)) * 1664525) + & 237 seed(j) + j ! Non-linear 238 ! Non-linear 239 twister%state(i) = iand(twister%state(i), -1) ! for >32 bit machines ! for >32 bit machines 240 nWraps = nWraps + 1 241 else 242 twister%state(i) = ieor(twister%state(i), & 243 ieor(twister%state(i-1), & 244 ishft(twister%state(i-1), -30)) * 1664525) + & 245 seed(j) + j ! Non-linear 246 ! Non-linear 247 twister%state(i) = iand(twister%state(i), -1) ! for >32 bit machines ! for >32 bit machines 248 end if 249 end do 250 ! 251 ! Walk through the state array, beginning where we left off in the block above 252 ! 253 do i = mod(nFirstLoop, blockSize) + nWraps + 1, blockSize - 1 254 twister%state(i) = ieor(twister%state(i), & 255 ieor(twister%state(i-1), & 256 ishft(twister%state(i-1), -30)) * 1566083941) - i ! Non-linear 257 ! Non-linear 258 twister%state(i) = iand(twister%state(i), -1) ! for >32 bit machines ! for >32 bit machines 259 end do 260 twister%state(0) = twister%state(blockSize - 1) 261 do i = 1, mod(nFirstLoop, blockSize) + nWraps 262 twister%state(i) = ieor(twister%state(i), & 263 ieor(twister%state(i-1), & 264 ishft(twister%state(i-1), -30)) * 1566083941) - i ! Non-linear 265 ! Non-linear 266 twister%state(i) = iand(twister%state(i), -1) ! for >32 bit machines ! for >32 bit machines 267 end do 268 twister%state(0) = UMASK 269 twister%currentElement = blockSize 270 END FUNCTION initialize_vector 271 ! ------------------------------------------------------------- 272 ! Public functions 273 ! -------------------- 274 275 FUNCTION getrandomint(twister) 276 TYPE(randomnumbersequence), intent(inout) :: twister 277 INTEGER :: getrandomint 278 ! Generate a random integer on the interval [0,0xffffffff] 279 ! Equivalent to genrand_int32 in the C code. 280 ! Fortran doesn't have a type that's unsigned like C does, 281 ! so this is integers in the range -2**31 - 2**31 282 ! All functions for getting random numbers call this one, 283 ! then manipulate the result 284 if(twister%currentElement >= blockSize) call nextState(twister) 285 getRandomInt = temper(twister%state(twister%currentElement)) 286 twister%currentElement = twister%currentElement + 1 287 END FUNCTION getrandomint 288 ! -------------------- 289 290 ! -------------------- 291 !! mji - modified Jan 2007, double converted to rrtmg real kind type 292 293 FUNCTION getrandomreal(twister) 294 TYPE(randomnumbersequence), intent(inout) :: twister 295 ! double precision :: getRandomReal 296 REAL(KIND=r8) :: getrandomreal 297 ! Generate a random number on [0,1] 298 ! Equivalent to genrand_real1 in the C code 299 ! The result is stored as double precision but has 32 bit resolution 300 INTEGER :: localint 301 localInt = getRandomInt(twister) 302 if(localInt < 0) then 303 ! getRandomReal = dble(localInt + 2.0d0**32)/(2.0d0**32 - 1.0d0) 304 getRandomReal = (localInt + 2.0**32_r8)/(2.0**32_r8 - 1.0_r8) 305 else 306 ! getRandomReal = dble(localInt )/(2.0d0**32 - 1.0d0) 307 getRandomReal = (localInt )/(2.0**32_r8 - 1.0_r8) 308 end if 309 END FUNCTION getrandomreal 310 ! -------------------- 311 312 ! -------------------- 313 END MODULE mersennetwister 314 315 MODULE mcica_random_numbers 316 USE kgen_utils_mod, ONLY : kgen_dp, check_t, kgen_init_check, kgen_print_check 317 ! Generic module to wrap random number generators. 318 ! The module defines a type that identifies the particular stream of random 319 ! numbers, and has procedures for initializing it and getting real numbers 320 ! in the range 0 to 1. 321 ! This version uses the Mersenne Twister to generate random numbers on [0, 1]. 322 ! 323 ! The random number engine. 324 !! mji 325 !! use time_manager_mod, only: time_type, get_date 326 USE shr_kind_mod, ONLY: r8 => shr_kind_r8 327 ! use parkind, only : jpim, jprb 328 IMPLICIT NONE 329 PRIVATE 330 331 332 !! mji 333 !! initializeRandomNumberStream, getRandomNumbers, & 334 !! constructSeed 335 CONTAINS 336 337 ! write subroutines 338 ! No subroutines 339 ! No module extern variables 340 ! --------------------------------------------------------- 341 ! Initialization 342 ! --------------------------------------------------------- 343 344 ! --------------------------------------------------------- 345 346 ! --------------------------------------------------------- 347 ! Procedures for drawing random numbers 348 ! --------------------------------------------------------- 349 350 ! --------------------------------------------------------- 351 352 ! --------------------------------------------------------- 353 354 ! mji 355 ! ! --------------------------------------------------------- 356 ! ! Constructing a unique seed from grid cell index and model date/time 357 ! ! Once we have the GFDL stuff we'll add the year, month, day, hour, minute 358 ! ! --------------------------------------------------------- 359 ! function constructSeed(i, j, time) result(seed) 360 ! integer, intent( in) :: i, j 361 ! type(time_type), intent( in) :: time 362 ! integer, dimension(8) :: seed 363 ! 364 ! ! Local variables 365 ! integer :: year, month, day, hour, minute, second 366 ! 367 ! 368 ! call get_date(time, year, month, day, hour, minute, second) 369 ! seed = (/ i, j, year, month, day, hour, minute, second /) 370 ! end function constructSeed 371 END MODULE mcica_random_numbers 372