45
46
47
54
55
56
57#include "implicit_f.inc"
58
59
60
61#include "com04_c.inc"
62#include "scr17_c.inc"
63#include "units_c.inc"
64#include "param_c.inc"
65#include "sphcom.inc"
66
67
68
69 TYPE (UNIT_TYPE_),INTENT(IN) ::UNITAB
71 . rnoise(nperturb,numelc+numeltg+numels+numsph),
72 . qp_rperturb(nperturb,4)
73 INTEGER IPART(LIPART1,*),
74 . IPM(NPROPMI,*),OFFS,
75 . IPARTS(*),PERTURB(NPERTURB),
76 . IDPERTURB(NPERTURB),INDEX(NUMELC+NUMELTG+NUMELS+NUMSPH),
77 . INDEX_ITYP(NUMELC+NUMELTG+NUMELS+NUMSPH),NPART_SOLID,
78 . QP_IPERTURB(NPERTURB,6)
79 TYPE(SUBMODEL_DATA) LSUBMODEL(*)
80
81 TYPE (GROUP_) , DIMENSION(NGRPART) :: IGRPART
82
83
84
85 INTEGER I,J,K,NUMA,I_METHOD,MAX_PART,
86 . CPT_PART,NB_RANDOM,I_SEED,DISTRIB(50),
87 . II,NB_INTERV,N,IOK,SEED,SEED_RANDOM,
88 . ITYP,I_PERTURB_VAR,IGRPRTS,SIZEY
89 CHARACTER(LEN=NCHARTITLE) :: TITR
90 INTEGER, DIMENSION(:,:), ALLOCATABLE :: TAB_PART
91 INTEGER, DIMENSION(:), ALLOCATABLE :: A_SEED
92 INTEGER, DIMENSION(1:8) :: DT_SEED
94 . mean,sd,mean_input,sd_input,max_distrib,temp,min_value,
95 . max_value,minval,maxval,bid
96 my_real,
DIMENSION(:),
ALLOCATABLE :: array
97 CHARACTER*100 CHAR(100)
98 CHARACTER*100 CHAR1(100)
99 CHARACTER*100 CHAR2
100 CHARACTER MESS*40
101 CHARACTER(LEN=NCHARFIELD)::CHVAR
102 LOGICAL IS_AVAILABLE
103
104
105
106 DATA mess/'PERTURBATION DEFINITION '/
107
108
109
110
111
112 max_part = 0
113 ityp = 0
114 bid = 0
115 is_available = .false.
117
118
119
120 DO i=1+offs,npart_solid+offs
121
122
123 titr = ''
125 . option_id = idperturb(i),
126 . option_titr = titr)
127
128
129 ityp = 3
130
131 i_perturb_var = 0
132 cpt_part = 0
133 igrprts = 0
134 iok = 0
135
136
137 CALL hm_get_intv(
'grpart_ID',igrprts,is_available,lsubmodel)
139 IF (chvar(1:4) == 'dens' .OR. chvar(1:4) == 'DENS') i_perturb_var = 1
140
141
142 IF (i_perturb_var /= 1)
CALL ancmsg(msgid=1194,
143 . msgtype=msgerror,
144 . anmode=aninfo,
145 . i1=idperturb(i),
146 . c1=titr,
147 . c2=chvar)
148
149
150 IF (igrprts /= 0)THEN
151 DO n=1,ngrpart
152 IF (igrpart(n)%ID == igrprts) THEN
153 igrprts = n
154 iok = 1
155 ityp = 3
156 EXIT
157 END IF
158 END DO
159 ENDIF
160
161
162 perturb(i) = ityp
163
164
165 IF (iok == 0) THEN
167 . msgtype=msgerror,
168 . anmode=aninfo,
169 . i1=idperturb(i),
170 . c1=titr,
171 . i2=igrprts,
172 . c2='GROUP OF PART')
173 ELSEIF (iok == 1)THEN
174 cpt_part = igrpart(igrprts)%NENTITY
175 ENDIF
176 max_part =
max(max_part,cpt_part)
177 ENDDO
178
179 ALLOCATE(tab_part(nperturb,max_part))
180
181
182
183
185 DO i=1+offs,npart_solid+offs
186
187
188 index(1:(numelc+numeltg+numels+numsph)) = 0
189 index_ityp(1:numelc+numeltg+numels+numsph) = 0
190
191
192 titr = ''
194 . option_id = idperturb(i),
195 . option_titr = titr)
196
197
198 ityp = 3
199
200
201 CALL hm_get_floatv(
'F_Mean' ,mean ,is_available, lsubmodel, unitab)
202 CALL hm_get_floatv(
'Deviation' ,sd ,is_available, lsubmodel, unitab)
203 CALL hm_get_floatv(
'Min_cut' ,minval ,is_available, lsubmodel, unitab)
204 CALL hm_get_floatv(
'Max_cut' ,maxval ,is_available, lsubmodel, unitab)
206 CALL hm_get_intv(
'Idistri' ,i_method,is_available, lsubmodel)
207
208
209 IF(i_method == 0) i_method = 2
210 IF(minval == zero .AND. maxval == zero) THEN
211 IF(i_method == 1) THEN
212 ELSEIF(i_method == 2)THEN
213 minval = -ep30
214 maxval = ep30
215 ENDIF
216 ENDIF
217 sd_input = sd
218 mean_input = mean
219
220
221 qp_iperturb(i,1) = idperturb(i)
222 qp_iperturb(i,2) = ityp
223 qp_iperturb(i,3) =
seed
224 qp_iperturb(i,4) = i_method
225 qp_rperturb(i,1) = mean
226 qp_rperturb(i,2) = sd
227 qp_rperturb(i,3) = minval
228 qp_rperturb(i,4) = maxval
229
230
231 cpt_part = 0
232 igrprts = 0
233 iok = 0
234
235
236 CALL hm_get_intv(
'grpart_ID' ,igrprts ,is_available,lsubmodel)
237 qp_iperturb(i,5) = igrprts
239 IF (chvar(1:4) == 'dens' .OR. chvar(1:4) == 'DENS') qp_iperturb(i,6) = 1
240
241
242 IF (igrprts /= 0)THEN
243 DO n=1,ngrpart
244 IF (igrpart(n)%ID == igrprts) THEN
245 igrprts = n
246 iok = 1
247 ityp = 3
248 EXIT
249 END IF
250 END DO
251 ENDIF
252
253
254 IF (iok == 1) THEN
255 cpt_part = 0
256 DO j=1,igrpart(igrprts)%NENTITY
257 cpt_part = cpt_part + 1
258 numa = igrpart(igrprts)%ENTITY(j)
259 tab_part(i,cpt_part) = numa
260 END DO
261 ENDIF
262
263
264 IF(i_method == 2) THEN
265 WRITE (iout,6000)
266 . idperturb(i),
'GAUSSIAN',mean_input,sd_input,
seed
267 WRITE (iout,'(10I10)') ipart(4,tab_part(i,1:cpt_part))
268 WRITE(iout,*) ' '
269 WRITE(iout,*) ' '
270 ELSEIF(i_method == 1) THEN
271 WRITE (iout,6100)
272 . idperturb(i),
'RANDOM',
seed
273 WRITE (iout,'(10I10)') ipart(4,tab_part(i,1:cpt_part))
274 WRITE(iout,*) ' '
275 WRITE(iout,*) ' '
276 ENDIF
277
278
279 nb_random = 0
280 DO ii=1,numels
281 DO k=1,cpt_part
282 IF(iparts(ii) == tab_part(i,k)) THEN
283 nb_random = nb_random + 1
284 index(nb_random) = ii
285 index_ityp(nb_random) = 1
286 ENDIF
287 ENDDO
288 ENDDO
289
290
292 CALL random_seed(size=i_seed)
293 ALLOCATE(a_seed(1:i_seed))
294 CALL random_seed(get=a_seed)
295 CALL date_and_time(values=dt_seed)
296 a_seed=dt_seed(8)*dt_seed(7)*dt_seed(6)
297 seed=dt_seed(8)*dt_seed(7)*dt_seed(6)
298 CALL random_seed(put=a_seed)
299 seed_random = 1
300 DEALLOCATE(a_seed)
301 ELSE
302 CALL random_seed(size=i_seed)
303 ALLOCATE(a_seed(1:i_seed))
305 CALL random_seed(put=a_seed)
306 seed_random = 0
307 DEALLOCATE(a_seed)
308 ENDIF
309
310
311 char=''
312 char1=''
313 char2=''
314 distrib(1:50) = 0
315 ALLOCATE(array(nb_random+2))
316 CALL random_number(array)
317
318
319 max_value = -ep30
320 min_value = ep30
321 IF ( i_method == 2) THEN
322 DO ii = 1, nb_random+1, 2
323 temp = sd * sqrt(-2.0*log(array(ii))) * cos(2*pi*array(ii+1)) + mean
324 array(ii+1) = sd * sqrt(-2.0*log(array(ii))) * sin(2*pi*array(ii+1)) + mean
325 array(ii) = temp
326 END DO
327 DO ii = 1, nb_random
328 array(ii) =
max(
min(maxval,array(ii)),minval)
329 max_value =
max(array(ii),max_value)
330 min_value =
min(array(ii),min_value)
331 END DO
332 ELSEIF(i_method == 1)THEN
333 DO ii = 1, nb_random
334 array(ii) = array(ii)*(maxval-minval)+minval
335 max_value =
max(array(ii),max_value)
336 min_value =
min(array(ii),min_value)
337 END DO
338 ENDIF
339
340
341 DO ii = 1, nb_random
342 IF (index_ityp(ii) == 1) THEN
343 rnoise(i,index(ii)+numelc+numeltg) = array(ii)
344 ENDIF
345 ENDDO
346
347
348 mean = sum(array)/nb_random
349 sd = sqrt(sum((array - mean)**2)/nb_random)
350
351
352 IF(i_method == 2) THEN
353 max_distrib = one /(sd*sqrt(two * pi))
354 ELSEIF(i_method == 1) THEN
355 max_distrib = one /(max_value-min_value)
356 ENDIF
357 WRITE (iout,4500)
358 WRITE(iout,*) ' '
359 nb_interv = 50
360 sizey = 20
361 IF (minval /= -ep30 .AND. maxval /= ep30)THEN
362 min_value = minval
363 max_value = maxval
364 ENDIF
365 CALL plot_distrib( array,nb_random, nb_interv,sizey,min_value,
366 . max_value,max_distrib,'#')
367 IF(i_method == 2) THEN
368 WRITE (iout,2000) mean,sd
369 ELSEIF(i_method == 1) THEN
370 WRITE (iout,2050) mean
371 ENDIF
372 IF(seed_random == 1)
WRITE (iout,2100)
seed
373 WRITE(iout,*) ' '
374 WRITE(iout,*) ' '
375 IF (ALLOCATED(array)) DEALLOCATE(array)
376 ENDDO
377 DEALLOCATE(tab_part)
378
379 6000 FORMAT(/' PERTURBATION ID',i10/
380 + ' ---------------'/
381 + ' TYPE . . . . . . . . . . . . . . .',a/
382 + ' INPUT MEAN VALUE . . . . . . . . .',1pg20.13/
383 + ' INPUT STANDARD DEVIATION . . . . .',1pg20.13/
384 + ' INPUT SEED VALUE . . . . . . . . .',i10/
385 + ' SOLID DENSITIES, PARTS:')
386 6100 FORMAT(/' PERTURBATION ID',i10/
387 + ' ---------------'/
388 + ' TYPE . . . . . . . . . . . . . . .',a/
389 + ' INPUT SEED VALUE . . . . . . . . .',i10/
390 + ' SOLID DENSITIES, PARTS:')
391
392 2000 FORMAT(/
393 + ' GENERATED MEAN VALUE . . . . . . .',1pg20.13/
394 + ' GENERATED STANDARD DEVIATION . . .',1pg20.13)
395 2050 FORMAT(/
396 + ' GENERATED MEAN VALUE . . . . . . .',1pg20.13)
397 2100 FORMAT(/
398 + ' GENERATED SEED VALUE . . . . . . .',i10/)
399
400 4500 FORMAT(/
401 + ' DISTRIBUTION OF SCALE FACTORS APPLIED TO DENSITIES OF SOLIDS')
402
403
subroutine hm_get_floatv(name, rval, is_available, lsubmodel, unitab)
subroutine hm_get_intv(name, ival, is_available, lsubmodel)
subroutine hm_get_string(name, sval, size, is_available)
subroutine hm_option_start(entity_type)
integer, parameter nchartitle
integer, parameter ncharkey
integer, parameter ncharfield
subroutine plot_distrib(array, s_array, nb_interv, sizey, x_minvalue, x_maxvalue, y_maxvalue, ecrit)
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)