-
Notifications
You must be signed in to change notification settings - Fork 1
/
mod_rnorm.f90
194 lines (160 loc) · 4.35 KB
/
mod_rnorm.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
module mod_rnorm
!$$$ program documentation block
! . . .
! Program name: mod_randnorm
! programmer: Da,Cheng org: umd date: 2016-Jan-07
!
! Purpose:
! Generate random number from a normal distribution with mean of 0, and
! standard deviation of 1.
! This module is directly copied from the the NOAA GSI codes, and is named
! of rnorm.f90.
!
! !!!! Note !!!!!
! Credit should be given to NCEP GSI codes
!
! Revision history:
! 2015-MON-DD Da,Cheng -Creator
!
! File dependencies
!
! Attributes:
! language: fortran 90
! machine :
!
!$$$ end documentation block
use mod_type, only: r_single => rsp, i_kind => isp
private
public :: rnorm, rnorm1d, set_random_seed
contains
subroutine rnorm1d( nn, r1d )
implicit none
integer,intent(in) :: nn
real(r_single),intent(inout) :: r1d(nn)
integer :: i
do i =1 , nn
r1d(i) = rnorm()
enddo
endsubroutine
real(r_single) function rnorm() result( fn_val )
! Generate a random normal deviate using the polar method.
! Reference: Marsaglia,G. & Bray,T.A. 'A convenient method for generating
! normal variables', Siam Rev., vol.6, 260-264, 1964.
implicit none
! Local variables
real(r_single) :: u, sum
real(r_single), save :: v, sln
logicaL, save :: second = .false.
real(r_single), parameter :: one = 1.0_r_single, vsmall = TINY( one )
if (second) then
! If second, use the second random number generated on last call
second = .false.
fn_val = v*sln
else
! First call; generate a pair of random normals
second = .true.
do
call random_number( u )
call random_number( v )
u = scale( u, 1 ) - one
v = scale( v, 1 ) - one
sum = u*u + v*v + vsmall ! vsmall added to prevent LOG(zero) / zero
IF(sum < one) exit
end do
sln = sqrt(- scale( log(sum), 1 ) / sum)
fn_val = u*sln
end if
return
end function rnorm
subroutine iran(l,n,ran_int)
! generate an array of N random integers between 0 and L-1.
real(r_single) :: rnd
integer(i_kind) :: L,i,N
integer(i_kind) :: ran_int(N)
do i = 1,n
call random_number(rnd)
ran_int(i) = nint(float(l-1)*rnd)
end do
end subroutine iran
subroutine set_random_seed ( iseed , myrank)
!
!*******************************************************************************
!
!! SET_RANDOM_SEED initializes the FORTRAN 90 random number generator.
!
!
! Discussion:
!
! If ISEED is nonzero, then that value is used to construct a seed.
!
! If ISEED is zero, then the seed is determined by calling the date
! and time routine. Thus, if the code is run at different times,
! different seed values will be set.
!
! Parameters:
!
! Input, integer ISEED, is nonzero for a user seed, or 0 if the
! seed should be determined by this routine.
!
implicit none
!
integer(i_kind), intent(in), optional :: myrank
integer(i_kind) date_time(8)
logical, parameter :: debug = .false.
integer(i_kind) i,j,k
integer(i_kind) iseed
integer(i_kind), allocatable :: seed(:)
!
! Initialize the random seed routine.
!
call random_seed
!
! Request the size of a typical seed.
! (On the DEC ALPHA, K is returned as 2.)
!
call random_seed ( size = k )
if ( debug ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SET_RANDOM_SEED:'
write ( *, '(a,i6)' ) ' Random seed size is K = ', k
end if
!
! Set up space for a seed vector.
!
allocate ( seed(k) )
if ( iseed /= 0 ) then
seed(1:k) = iseed
else
!
! Make up a "random" value based on date and time information.
!
call date_and_time ( values = date_time )
do i = 1, k
seed(i) = 0
do j = 1, 8
if (present(myrank)) then
seed(i) = seed(i) + ( j + i ) * date_time(j) + myrank * 100
else
seed(i) = seed(i) + ( j + i ) * date_time(j)
endif
seed(i) = ishftc ( seed(i), 4 * ( j - 1 ) )
end do
end do
end if
if ( debug ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SET_RANDOM_SEED:'
write ( *, '(a)' ) ' The random seed vector:'
write ( *, '(a)' ) ' '
do i = 1, k
write ( *, '(i12)' ) seed(i)
end do
end if
!
! Send this random value back to the RANDOM_SEED routine, to be
! used as the seed of the random number generator.
!
call random_seed ( put = seed(1:k) )
deallocate ( seed )
end subroutine set_random_seed
end module