1 ! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 2 ! C 3 ! Module name: Random Number Generation Utilities C 4 ! Purpose: Removed from interpolation mod and added built-in random 5 ! number routines instead of Pope's 6 ! C 7 ! C 8 ! Author: Sreekanth Pannala and Rahul Garg Date: 23-Oct-08 C 9 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 10 MODULE randomno 11 12 USE constant 13 14 IMPLICIT NONE 15 16 PRIVATE 17 PUBLIC :: uni_rno, nor_rno 18 19 20 CONTAINS 21 22 23 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv 24 SUBROUTINE UNI_RNO(Y) 25 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 26 IMPLICIT NONE 27 28 !----------------------------------------------- 29 ! Local variables 30 !----------------------------------------------- 31 double precision, intent(out), dimension(:) :: y 32 double precision rmean, variance, sigma 33 integer i, nsize 34 !----------------------------------------------- 35 36 nsize = size(y(:)) 37 38 call init_random_seed 39 call random_number(y) 40 41 rmean = sum(y(:))/nsize 42 43 ! write(*,*) 'Generating Uniform Random Variables for size', nsize 44 ! write(*,*) 'mean', rmean 45 46 variance = 0.0 47 do i = 1, nsize 48 ! write(20,*) i, y(i) 49 variance = variance + (y(i)-rmean)**2 50 end do 51 52 close(20) 53 54 variance = variance/nsize 55 sigma = sqrt(variance) 56 57 ! write(*,*) 'sigma', sigma 58 59 RETURN 60 END SUBROUTINE UNI_RNO 61 62 63 64 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv 65 SUBROUTINE NOR_RNO(Y, mean, sigma) 66 67 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 68 IMPLICIT NONE 69 70 !----------------------------------------------- 71 ! Local variables 72 !----------------------------------------------- 73 double precision, intent(out), dimension(:) :: y 74 double precision mean, sigma 75 76 double precision lmean, lvariance, lsigma 77 double precision x(2), w 78 integer i, nsize, n 79 ! no. of times this routine has been called 80 integer, save :: COUNTER = 0 81 ! so all components are written 82 integer fileunit 83 !----------------------------------------------- 84 COUNTER = COUNTER + 1 85 fileunit = 20 + COUNTER 86 87 nsize = size(y(:)) 88 89 call init_random_seed 90 91 do i = 1, ceiling(real(nsize/2.0)) 92 do n = 1,100000 93 call random_number(x) 94 x = 2.0 * x - 1.0 95 w = x(1)**2 + x(2)**2 96 if(w.lt.1.0) exit 97 end do 98 99 w = sqrt( (-2.0 * log( w ) ) / w ) 100 y(2*i-1) = x(1) * w * sigma + mean 101 if(2*i.lt.nsize) y(2*i) = x(2) * w * sigma + mean 102 end do 103 104 lmean = sum(y(:))/nsize 105 106 !write(*,'(7X,A)') 'Generating Normal Random Variables' 107 !write(*,'(7X,A,2X,ES15.5)') 'specified mean =', mean 108 !write(*,'(7X,A,2X,ES15.5)') 'computed mean =', lmean 109 110 !write(fileunit,'(A)') 'FROM NOR_RNO' 111 ! specific to the call from init_particles_jn 112 !write(fileunit,'(A,I5,A)') 'FOR DIRECTION = ', & 113 ! COUNTER, ' where (1=X,2=Y,3=Z)' 114 !write(fileunit,'(5X,A,5X,A)') 'particle no.', 'velocity component' 115 116 lvariance = 0.0 117 do i = 1, nsize 118 !write(fileunit,'(I10,5X,ES15.5)') i, y(i) 119 lvariance = lvariance + (y(i)-lmean)**2 120 end do 121 122 !close(fileunit) 123 124 lvariance = lvariance/nsize 125 lsigma = sqrt(lvariance) 126 127 !write(*,'(7X,A,2X,ES15.5)') 'specified sigma =', sigma 128 !write(*,'(7X,A,2X,ES15.5)') 'computed sigma =', lsigma 129 130 RETURN 131 END SUBROUTINE NOR_RNO 132 133 134 135 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv 136 SUBROUTINE init_random_seed 137 138 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 139 IMPLICIT NONE 140 141 !----------------------------------------------- 142 ! Local variables 143 !----------------------------------------------- 144 INTEGER :: isize,idate(8) 145 INTEGER,ALLOCATABLE :: iseed(:) 146 !----------------------------------------------- 147 148 CALL DATE_AND_TIME(VALUES=idate) 149 CALL RANDOM_SEED(SIZE=isize) 150 ALLOCATE( iseed(isize) ) 151 CALL RANDOM_SEED(GET=iseed) 152 iseed = iseed * (idate(8)-500) ! idate(8) contains millisecond 153 CALL RANDOM_SEED(PUT=iseed) 154 155 DEALLOCATE( iseed ) 156 157 END SUBROUTINE init_random_seed 158 159 160 END MODULE randomno 161 162