46
47
48
55 USE format_mod , ONLY : lfield
56
57
58
59#include "implicit_f.inc"
60
61
62
63#include "com04_c.inc"
64#include "scr17_c.inc"
65#include "units_c.inc"
66#include "param_c.inc"
67#include "sphcom.inc"
68
69
70
71 TYPE (UNIT_TYPE_),INTENT(IN) ::UNITAB
73 . rnoise(nperturb,numelc+numeltg+numels+numsph),
74 . qp_rperturb(nperturb,4)
75 INTEGER IPART(LIPART1,*),IPARTC(*),
76 . IPARTG(*),IPM(NPROPMI,*),
77 . PERTURB(NPERTURB),OFFS,
78 . IDPERTURB(NPERTURB),INDEX(NUMELC+NUMELTG+NUMELS+NUMSPH),
79 . INDEX_ITYP(NUMELC+NUMELTG+NUMELS+NUMSPH),NPART_SHELL,
80 . QP_IPERTURB(NPERTURB,6)
81 TYPE(SUBMODEL_DATA) LSUBMODEL(*)
82
83 TYPE (GROUP_) , DIMENSION(NGRPART) :: IGRPART
84
85
86
87 INTEGER I,J,K,NUMA,I_METHOD,MAX_PART,
88 . CPT_PART,NB_RANDOM,I_SEED,DISTRIB(50),
89 . II,NB_INTERV,IGRPRT,N,IOK,SEED,SEED_RANDOM,
90 . ITYP,I_PERTURB_VAR,SIZEY,EMPTY
91 CHARACTER(LEN=NCHARTITLE) :: TITR
92 INTEGER, DIMENSION(:,:), ALLOCATABLE :: TAB_PART
93 INTEGER, DIMENSION(:), ALLOCATABLE :: A_SEED
94 INTEGER, DIMENSION(1:8) :: DT_SEED
96 . mean,sd,mean_input,sd_input,max_distrib,temp,min_value,
97 . max_value,minval,maxval,bid
98 my_real,
DIMENSION(:),
ALLOCATABLE :: array
99 CHARACTER*100 CHAR(100)
100 CHARACTER*100 CHAR1(100)
101 CHARACTER*100 CHAR2
102 CHARACTER MESS*40
103 CHARACTER(LEN=NCHARFIELD)::CHVAR
104 LOGICAL IS_AVAILABLE
105
106
107
108 DATA mess/'PERTURBATION DEFINITION '/
109
110
111
112
113
114 max_part = 0
115 ityp = 0
116 bid = 0
117 is_available = .false.
119
120
121
122 DO i=1+offs,npart_shell+offs
123
124
125 titr = ''
127 . option_id = idperturb(i),
128 . option_titr = titr)
129
130
131 ityp = 1
132
133 i_perturb_var = 0
134 cpt_part = 0
135 igrprt = 0
136 iok = 0
137 chvar = ''
138
139
140 CALL hm_get_intv(
'grpart_ID' ,igrprt ,is_available,lsubmodel)
142 IF (chvar(1:5) == 'thick' .OR. chvar(1:5) == 'THICK') i_perturb_var = 1
143
144 empty = 1
145 DO k = 1,lfield
146 IF(chvar(k:k) /= ' ') empty = 0
147 ENDDO
148
149
150 IF (i_perturb_var /= 1 .AND. empty == 0)
CALL ancmsg(msgid=1194,
151 . msgtype=msgerror,
152 . anmode=aninfo,
153 . i1=idperturb(i),
154 . c1=titr,
155 . c2=chvar)
156
157
158 IF (igrprt /= 0) THEN
159 DO n=1,ngrpart
160 IF (igrpart(n)%ID == igrprt) THEN
161 igrprt=n
162 iok = 1
163 ityp = 1
164 EXIT
165 END IF
166 END DO
167 ENDIF
168
169
170 perturb(i) = ityp
171
172
173 IF (iok == 0) THEN
175 . msgtype=msgerror,
176 . anmode=aninfo,
177 . i1=idperturb(i),
178 . c1=titr,
179 . i2=igrprt,
180 . c2='GROUP OF PART')
181 ELSEIF(iok == 1)THEN
182 cpt_part = igrpart(igrprt)%NENTITY
183 ENDIF
184 max_part =
max(max_part,cpt_part)
185 ENDDO
186
187 ALLOCATE(tab_part(nperturb,max_part))
188
189
190
191
193 DO i = 1+offs,npart_shell+offs
194
195
196 index(1:(numelc+numeltg+numels+numsph)) = 0
197 index_ityp(1:numelc+numeltg+numels+numsph) = 0
198
199
200 titr = ''
202 . option_id = idperturb(i),
203 . option_titr = titr)
204
205
206 ityp = 1
207
208
209 CALL hm_get_floatv(
'F_Mean' ,mean ,is_available, lsubmodel, unitab)
210 CALL hm_get_floatv(
'Deviation' ,sd ,is_available, lsubmodel, unitab)
211 CALL hm_get_floatv(
'Min_cut' ,minval ,is_available, lsubmodel, unitab)
212 CALL hm_get_floatv(
'Max_cut' ,maxval ,is_available, lsubmodel, unitab)
214 CALL hm_get_intv(
'Idistri' ,i_method,is_available, lsubmodel)
215
216
217 IF(i_method == 0) i_method = 2
218 IF(minval == zero .AND. maxval == zero) THEN
219 IF(i_method == 1) THEN
220 ELSEIF(i_method == 2)THEN
221 minval = -ep30
222 maxval = ep30
223 ENDIF
224 ENDIF
225 sd_input = sd
226 mean_input = mean
227
228
229 qp_iperturb(i,1) = idperturb(i)
230 qp_iperturb(i,2) = ityp
231 qp_iperturb(i,3) =
seed
232 qp_iperturb(i,4) = i_method
233 qp_rperturb(i,1) = mean
234 qp_rperturb(i,2) = sd
235 qp_rperturb(i,3) = minval
236 qp_rperturb(i,4) = maxval
237
238
239 cpt_part = 0
240 igrprt = 0
241 iok = 0
242
243
244 CALL hm_get_intv(
'grpart_ID',igrprt,is_available,lsubmodel)
245 qp_iperturb(i,5) = igrprt
247 IF (chvar(1:5) == 'thick' .OR. chvar(1:5) == 'THICK') qp_iperturb(i,6) = 1
248
249
250 IF (igrprt /= 0) THEN
251 DO n=1,ngrpart
252 IF (igrpart(n)%ID == igrprt) THEN
253 igrprt = n
254 iok = 1
255 ityp = 1
256 EXIT
257 END IF
258 END DO
259 ENDIF
260
261
262 IF (iok == 1) THEN
263 cpt_part = 0
264 DO j=1,igrpart(igrprt)%NENTITY
265 cpt_part = cpt_part + 1
266 numa = igrpart(igrprt)%ENTITY(j)
267 tab_part(i,cpt_part) = numa
268 END DO
269 ENDIF
270
271
272 IF(i_method == 2) THEN
273 WRITE (iout,1000)
274 . idperturb(i),
'GAUSSIAN',mean_input,sd_input,
seed
275 WRITE (iout,'(10I10)') ipart(4,tab_part(i,1:cpt_part))
276 WRITE(iout,*) ' '
277 WRITE(iout,*) ' '
278 ELSEIF(I_METHOD == 1) THEN
279 WRITE (IOUT,1100)
280 . IDPERTURB(I),'random',SEED
281 WRITE (IOUT,'(10i10)') IPART(4,TAB_PART(I,1:CPT_PART))
282 WRITE(IOUT,*) ' '
283 WRITE(IOUT,*) ' '
284 ENDIF
285
286
287 nb_random = 0
288 DO ii=1,numelc
289 DO k=1,cpt_part
290 IF (ipartc(ii) == tab_part(i,k)) THEN
291 nb_random = nb_random + 1
292 index(nb_random) = ii
293 index_ityp(nb_random) = 3
294 ENDIF
295 ENDDO
296 ENDDO
297 DO ii=1,numeltg
298 DO k=1,cpt_part
299 IF(ipartg(ii) == tab_part(i,k)) THEN
300 nb_random = nb_random + 1
301 index(nb_random) = ii
302 index_ityp(nb_random) = 7
303 ENDIF
304 ENDDO
305 ENDDO
306
307
309 CALL random_seed(size=i_seed)
310 ALLOCATE(a_seed(1:i_seed))
311 CALL random_seed(get=a_seed)
312 CALL date_and_time(values=dt_seed)
313 a_seed=dt_seed(8)*dt_seed(7)*dt_seed(6)
314 seed=dt_seed(8)*dt_seed(7)*dt_seed(6)
315 CALL random_seed(put=a_seed)
316 seed_random = 1
317 DEALLOCATE(a_seed)
318 ELSE
319 CALL random_seed(size=i_seed)
320 ALLOCATE(a_seed(1:i_seed))
322 CALL random_seed(put=a_seed)
323 seed_random = 0
324 DEALLOCATE(a_seed)
325 ENDIF
326
327
328 char=''
329 char1=''
330 char2=''
331 distrib(1:50) = 0
332 ALLOCATE(array(nb_random+2))
333 CALL random_number(array)
334
335
336 max_value = -ep30
337 min_value = ep30
338 IF ( i_method == 2) THEN
339 DO ii = 1, nb_random+1, 2
340 temp = sd * sqrt(-2.0*log(array(ii))) * cos(2*pi*array(ii+1)) +
341 . mean
342 array(ii+1) =
343 . sd * sqrt(-2.0*log(array(ii))) * sin(2*pi*array(ii+1)) + mean
344 array(ii) = temp
345 END DO
346 DO ii = 1, nb_random
347 array(ii) =
max(
min(maxval,array(ii)),minval)
348 max_value =
max(array(ii),max_value)
349 min_value =
min(array(ii),min_value)
350 END DO
351 ELSEIF(i_method == 1)THEN
352 DO ii = 1, nb_random
353 array(ii) = array(ii)*(maxval-minval)+minval
354 max_value =
max(array(ii),max_value)
355 min_value =
min(array(ii),min_value)
356 END DO
357 ENDIF
358
359
360 DO ii = 1, nb_random
361 IF (index_ityp(ii) == 3) THEN
362 rnoise(i,index(ii)) = array(ii)
363 ELSEIF (index_ityp(ii) == 7) THEN
364 rnoise(i,index(ii)+numelc) = array(ii)
365 ENDIF
366 ENDDO
367
368
369 mean = sum(array)/nb_random
370 sd = sqrt(sum((array - mean)**2)/nb_random)
371
372
373 IF(i_method == 2) THEN
374 max_distrib = one /(sd*sqrt(two * pi))
375 ELSEIF(i_method == 1) THEN
376 max_distrib = one /(max_value-min_value)
377 ENDIF
378 WRITE (iout,3000)
379 WRITE(iout,*) ' '
380 nb_interv = 50
381 sizey = 20
382 IF (minval /= -ep30 .AND. maxval /= ep30)THEN
383 min_value = minval
384 max_value = maxval
385 ENDIF
386 CALL plot_distrib( array,nb_random, nb_interv,sizey,min_value,
387 . max_value,max_distrib,'#')
388 IF(i_method == 2) THEN
389 WRITE (iout,2000) mean,sd
390 ELSEIF(i_method == 1) THEN
391 WRITE (iout,2050) mean
392 ENDIF
393 IF(seed_random == 1)
WRITE (iout,2100)
seed
394 WRITE(iout,*) ' '
395 WRITE(iout,*) ' '
396 IF (ALLOCATED(array)) DEALLOCATE(array)
397 ENDDO
398 DEALLOCATE(tab_part)
399
400 1000 FORMAT(/' PERTURBATION ID',i10/
401 + ' ---------------'/
402 + ' TYPE . . . . . . . . . . . . . . .',a/
403 + ' INPUT MEAN VALUE . . . . . . . . .',1pg20.13/
404 + ' INPUT STANDARD DEVIATION . . . . .',1pg20.13/
405 + ' INPUT SEED VALUE . . . . . . . . .',i10/
406 + ' SHELL THICKNESSES, PARTS:')
407 1100 FORMAT(/' PERTURBATION ID',i10/
408 + ' ---------------'/
409 + ' TYPE . . . . . . . . . . . . . . .',a/
410 + ' INPUT SEED VALUE . . . . . . . . .',i10/
411 + ' SHELL THICKNESSES, PARTS:')
412
413 2000 FORMAT(/
414 + ' GENERATED MEAN VALUE . . . . . . .',1pg20.13/
415 + ' GENERATED STANDARD DEVIATION . . .',1pg20.13)
416 2050 FORMAT(/
417 + ' GENERATED MEAN VALUE . . . . . . .',1pg20.13)
418 2100 FORMAT(/
419 + ' GENERATED SEED VALUE . . . . . . .',i10/)
420
421 3000 FORMAT(/
422 + ' DISTRIBUTION OF SCALE FACTORS APPLIED TO THICKNESSES OF SHELLS')
423
424
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)