75
76
77
78
79
80
81
82
83
84
85
86
87
88
92
93
94
95 real(r8), dimension(:), intent(out) :: harvest
96
97
98
99 logical, save :: gaus_stored = .false.
100
101 logical, dimension(SIZE(harvest)) :: mask
102
103 integer(i8b), save :: last_allocated = 0
104 integer(i8b) :: i, ii, m, mc, n, ng, nn
105
106 real(r8), dimension(SIZE(harvest)) :: rsq, v1, v2, v3
107 real(r8), allocatable, dimension(:), save :: g
108
109
110
111
112
113
114
115 n=SIZE(harvest)
116 IF (n.ne.last_allocated) THEN
117 IF (last_allocated.ne.0) DEALLOCATE (g)
118 ALLOCATE ( g(n) )
119 last_allocated=n
120 gaus_stored=.false.
121 END IF
122
123
124
125
126 IF (gaus_stored) THEN
127 harvest=g
128 gaus_stored=.false.
129 ELSE
130 ng=1
131 DO
132 IF (ng.gt.n) EXIT
135 v1(ng:n)=2.0_r8*v1(ng:n)-1.0_r8
136 v2(ng:n)=2.0_r8*v2(ng:n)-1.0_r8
137
138
139
140
141
142 rsq(ng:n)=v1(ng:n)**2+v2(ng:n)**2
143 mask(ng:n)=((rsq(ng:n).gt.0.0_r8).and.(rsq(ng:n).lt.1.0_r8))
144 mc=count(mask(ng:n))
145
146 ii=0
147 DO i=ng,n
148 IF (mask(i)) THEN
149 ii=ii+1
150 v3(ii)=v1(i)
151 END IF
152 END DO
154 ii=ng-1
155 DO i=ng,n
156 IF (mask(i)) THEN
157 ii=ii+1
158 v2(ii)=v2(i)
159 rsq(ii)=rsq(i)
160 END IF
161 END DO
162 ng=ng+nn
163 END DO
164
165
166
167
168 rsq=sqrt(-2.0_r8*log(rsq)/rsq)
169 harvest=v1*rsq
170 g=v2*rsq
171 gaus_stored=.true.
172 END IF
173
174 RETURN