47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
67 USE pblast_mod
72
73
74
75#include "implicit_f.inc"
76
77
78
79#include "param_c.inc"
80#include "com04_c.inc"
81#include "units_c.inc"
82#include "tabsiz_c.inc"
83
84
85
86 INTEGER ITAB(NUMNOD), ITABM1(NUMNOD), NUMLOADP, ILOADP(SIZLOADP,NLOADP), LLOADP(SLLOADP)
87 my_real facloadp(lfacload,nloadp)
88 my_real,
INTENT(IN) :: x(3,numnod)
89 my_real,
INTENT(INOUT) :: bufsf(sbufsf)
90 my_real,
INTENT(IN) :: rtrans(ntransf,nrtrans)
91 TYPE (SURF_),TARGET, DIMENSION(NSURF) :: IGRSURF
92 TYPE (UNIT_TYPE_),INTENT(IN) ::UNITAB
93 TYPE(SUBMODEL_DATA), DIMENSION(NSUBMOD), INTENT(IN) :: LSUBMODEL
94 TYPE(pblast_),INTENT(INOUT) :: PBLAST
95
96
97
98 LOGICAL :: IS_AVAILABLE,lFOUND, IS_AVAILABLE_TSTOP, BOOL_COORD, BOOL_SKIP_CALC
99 LOGICAL :: BOOL_UNDERGROUND_CURRENT_LOAD, BOOL_UNDERGROUND_CURRENT_SEG
100 LOGICAL :: IS_ITA_SHIFT,IS_DECAY_TO_BE_COMPUTED
101
102 CHARACTER(LEN=NCHARLINE) :: MSGOUT1
103 CHARACTER(LEN=NCHARLINE) :: MSGOUT2
104 CHARACTER(LEN=NCHARTITLE)::TITR
105
106 my_real xdet,ydet,zdet,tdet,wtnt,pmin
107 my_real dx,dy,dz, z, w13, dnorm,lambda,cos_theta,t_a
108 my_real fac_m_bb,fac_l_bb,fac_t_bb,fac_p_bb,fac_i_bb,
norm, nn2, tmp(3)
109 my_real nx_seg,ny_seg,nz_seg, norm_,npt
110 my_real nx_surf,ny_surf,nz_surf, base_x,base_y,base_z
114 my_real lg,zg,tmp_var,z_struct,zx,zy,zz,projz(3),projdet(3)
115 my_real p_inci,p_refl,decay_inci,decay_refl,i_inci,i_refl
119 my_real,
ALLOCATABLE,
DIMENSION(:,:) :: output_user_params
120
121 INTEGER I,J,K,NUMSEG,K_BLAST
122 INTEGER IAD,ID,,internal_SURF_ID
123 INTEGER IABAC,IADPL,itmp,IS,SUB_ID
124 INTEGER ITA_SHIFT
125 INTEGER N1,N2,N3,N4,IERR1,NDT,ILD_PBLAST,IZ_UPDATE,IMODEL,NODE_ID
126 INTEGER SURF_ID_ground,internal_SURF_ID_GROUND
127 INTEGER curve_id1,curve_id2,SEG_UNDERGROUND
128 INTEGER, DIMENSION(:), POINTER :: INGR2USR
129 INTEGER :: ISHAPE
130 INTEGER :: NWARN
131 TYPE(FRIEDLANDER_PARAMS_) :: FRIEDLANDER_PARAMS
132
133
134
135 CHARACTER :: MESS*40
138 my_real :: cst_255_div_ln_z1_on_zn, log10_, fac_unit
139 DATA cst_180 /180.000000000000000000/
140 DATA pi_/3.141592653589793238462643d0/
141 DATA cst_255_div_ln_z1_on_zn/-38.147316611455952998/
142 DATA log10_ /2.30258509299405000000/
143 DATA z1_/0.50000000000000000/
144 DATA fac_unit/3.966977216838196139019/
145 DATA mess/'BLAST LOADING DEFINITION '/
146
147
148
149 INTEGER,EXTERNAL :: DICHOTOMIC_SEARCH_R_DESC,DICHOTOMIC_SEARCH_R_ASC,USR2SYS,NGR2USR
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
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217 base_x = -huge(base_x)
218 base_y = -huge(base_y)
219 base_z = -huge(base_z)
220
221
222 CALL pblast_init_tables(pblast%PBLAST_DATA)
223 tmp(1:3) = huge(tmp(1))
224
225 pblast%PBLAST_DT%TA_INF = ep20
226 numseg = 0
227 ild_pblast = 0
228 alpha_zc = zero
229 hc = zero
230 bool_skip_calc = .false.
231 is_ita_shift = .false.
232 nwarn = 0
233
234 ALLOCATE( output_user_params(6,pblast%NLOADP_B))
235 output_user_params(:,:) = zero
236
238
239 DO k = 1+nloadp_f, nloadp_f+pblast%NLOADP_B
240 k_blast = k-nloadp_f
242 iad = numloadp + 1
243 surf_id = 0
244 surf_id_ground = 0
245 internal_surf_id = 0
246 seg_underground = 0
247 bool_underground_current_load = .false.
248
249 xdet = zero
250 ydet = zero
251 zdet = zero
252 tdet = zero
253 tstop = zero
254 wtnt = zero
255 ta_pblast = ep20
256
257
258
259
260
261
262 CALL hm_get_intv(
'surf_ID', surf_id, is_available, lsubmodel)
263 CALL hm_get_intv(
'Exp_data', iabac, is_available, lsubmodel)
264 CALL hm_get_intv('i_tshift
', ITA_SHIFT, IS_AVAILABLE, LSUBMODEL)
265 CALL HM_GET_INTV('ndt', NDT, IS_AVAILABLE, LSUBMODEL)
266 CALL HM_GET_INTV('iz', IZ_UPDATE, IS_AVAILABLE, LSUBMODEL)
267 CALL HM_GET_INTV('imodel', IMODEL, IS_AVAILABLE, LSUBMODEL)
268 CALL HM_GET_INTV('node_id', NODE_ID, IS_AVAILABLE, LSUBMODEL)
269 !---line-2
270 CALL HM_GET_FLOATV('xdet', XDET, IS_AVAILABLE, LSUBMODEL, UNITAB)
271 CALL HM_GET_FLOATV('ydet', YDET, IS_AVAILABLE, LSUBMODEL, UNITAB)
272 CALL HM_GET_FLOATV('zdet', ZDET, IS_AVAILABLE, LSUBMODEL, UNITAB)
273 CALL HM_GET_FLOATV('tdet', TDET, IS_AVAILABLE, LSUBMODEL, UNITAB)
274 CALL HM_GET_FLOATV('wtnt', WTNT, IS_AVAILABLE, LSUBMODEL, UNITAB)
275 !---line-3
276 CALL HM_GET_FLOATV('pmin', PMIN, IS_AVAILABLE, LSUBMODEL, UNITAB)
277 CALL HM_GET_FLOATV('tstop',TSTOP,IS_AVAILABLE_TSTOP, LSUBMODEL, UNITAB)
278 !---line-4
279 CALL HM_GET_INTV('surf_ground_id', SURF_ID_GROUND, IS_AVAILABLE, LSUBMODEL)
280 CALL HM_GET_INTV('ishape', ISHAPE, IS_AVAILABLE, LSUBMODEL)
281
282 OUTPUT_USER_PARAMS(1,K_BLAST)=SURF_ID
283 OUTPUT_USER_PARAMS(2,K_BLAST)=SURF_ID_GROUND
284
285 BOOL_COORD=.FALSE.
286.OR..OR. IF(XDET/=ZERO YDET/=ZERO ZDET/=ZERO)THEN
287 BOOL_COORD=.TRUE.
288 ENDIF
289
290 IF(SUB_ID /= 0)CALL SUBROTPOINT(XDET,YDET,ZDET,RTRANS,SUB_ID,LSUBMODEL)
291
292 !units
293 FAC_M_bb = UNITAB%FAC_MASS*EP03
294 FAC_L_bb = UNITAB%FAC_LENGTH*EP02
295 FAC_T_bb = UNITAB%FAC_TIME*EP06
296 FAC_P_bb = FAC_M_bb/FAC_L_bb/FAC_T_bb/FAC_T_bb
297 FAC_I_bb = FAC_P_bb*FAC_T_bb !/FAC_M_bb**THIRD
298 W13 = (WTNT*FAC_M_bb)**THIRD
299
300 !--------------------------------------------------------------
301 ! D e f a u l t V a l u e s
302 !--------------------------------------------------------------
303
304 IF(TSTOP == ZERO)TSTOP=EP20
305 IF(NDT == 0) NDT=100
306
307 IF(PMIN < ZERO)THEN
308 MSGOUT1=''
309 MSGOUT1='pmin
PARAMETER is negative. please check
the load history
for negative pressure,
'
310 MSGOUT2=''
311 MSGOUT2='as
the pblast model is fitted only
for the positive phase
'
312 CALL ANCMSG(MSGID=1907, MSGTYPE=MSGWARNING, ANMODE=ANINFO,C1=TRIM(TITR), I1=ID, C2=MSGOUT1, C3=MSGOUT2)
313 ENDIF
314
315 IF(WTNT == ZERO)THEN
316 MSGOUT1=''
317 MSGOUT1='tnt mass must be defined.'
318 MSGOUT2=''
319 MSGOUT2='please set
VALUE for wtnt
PARAMETER'
320 CALL ANCMSG(MSGID=1907, MSGTYPE=MSGWARNING, ANMODE=ANINFO,C1=TRIM(TITR), I1=ID, C2=MSGOUT1, C3=MSGOUT2)
321 ENDIF
322
323 IF(Node_id > 0)THEN
324 NODE_ID=USR2SYS(NODE_ID,ITABM1,MESS,ID)
325 IF(NODE_ID > 0)THEN
326 XDET = X(1,NODE_ID)
327 YDET = X(2,NODE_ID)
328 ZDET = X(3,NODE_ID)
329 ENDIF
330 OUTPUT_USER_PARAMS(5,K_BLAST)=NODE_ID
331 IF(BOOL_COORD)THEN
332 MSGOUT1=''
333 MSGOUT1='node identifier is provided.'
334 MSGOUT2=''
335 MSGOUT2='coordinates(xdet,ydet,zdet) will be ignored'
336 CALL ANCMSG(MSGID=1907, MSGTYPE=MSGWARNING, ANMODE=ANINFO,C1=TRIM(TITR), I1=ID, C2=MSGOUT1, C3=MSGOUT2)
337 ENDIF
338 ENDIF
339
340 INGR2USR => IGRSURF(1:NSURF)%ID
341 internal_SURF_ID = NGR2USR(SURF_ID,INGR2USR,NSURF)
342 internal_SURF_ID_GROUND = NGR2USR(SURF_ID_GROUND,INGR2USR,NSURF)
343
344.AND. IF(IABAC >= 2 SURF_ID_GROUND /= 0)THEN
345 IF(internal_SURF_ID_GROUND == 0)THEN
346 !user surface not found in input file
347 !display error message : ground_id not found
348 MSGOUT1=''
349 MSGOUT1='surface identifier
for ground definition was not found
'
350 MSGOUT2=''
352 WRITE(MSGOUT2(9:19),FMT='(i10)') SURF_ID_GROUND
353 CALL ANCMSG(MSGID=2047,MSGTYPE=MSGERROR,ANMODE=ANINFO,C1=TRIM(TITR),I1=ID,C2=MSGOUT1,C3=MSGOUT2)
354 BOOL_SKIP_CALC = .TRUE.
355 ENDIF
356 ENDIF
357
358 IF(SURF_ID == 0)THEN
359 !user forgot to input surf_id
360 ! display error message : surf_id missing
361 MSGOUT1=''
362 MSGOUT1='surface identifier is not provided'
363 MSGOUT2=''
364 MSGOUT2='blast calculation will be skipped'
365 CALL ANCMSG(MSGID=1907,MSGTYPE=MSGWARNING,ANMODE=ANINFO,C1=TRIM(TITR),I1=ID,C2=MSGOUT1,C3=MSGOUT2)
366 BOOL_SKIP_CALC = .TRUE.
367 ELSE
368 IF(internal_SURF_ID == 0)THEN
369 !surface is not existing
370 ! display error message : surf_id not found
371 MSGOUT1=''
372 MSGOUT1='invalid surface identifier'
373 MSGOUT2=''
375 WRITE(MSGOUT2(9:19),FMT='(i10)') SURF_ID
376 CALL ANCMSG(MSGID=2047,MSGTYPE=MSGERROR,ANMODE=ANINFO,C1=TRIM(TITR),I1=ID,C2=MSGOUT1,C3=MSGOUT2)
377 BOOL_SKIP_CALC = .TRUE.
378 ENDIF
379 ENDIF
380
381 !ITA_SHIFT : BLAST ARRIVAL
382 !---1:no shift (default)
383 !---2:shift loading starts from the first cycle T=0 is shift to time on which first targeted segment is reached (within all pblast options)
384 IF(ITA_SHIFT==0) ITA_SHIFT = 2
385.OR. IF(ITA_SHIFT < 0 ITA_SHIFT >= 3)ITA_SHIFT = 1
386 IF(ITA_SHIFT == 2)IS_ITA_SHIFT=.TRUE.
387
388 ! IMODEL : flag for blast model
389 !---1:Friedlander
390 !---2:modified Friedlander (default)
391.OR. IF(IMODEL <= 0 IMODEL > 2)IMODEL=2
392 IS_DECAY_TO_BE_COMPUTED = .FALSE.
393 IF(IMODEL == 2)IS_DECAY_TO_BE_COMPUTED = .TRUE.
394
395 ! IZ : Z update during Engine
396 !---1 : do not update (default)
397 !---2 : update blast characteristing during Engine (because target may have move)
398 IF(IZ_UPDATE == 0)IZ_UPDATE=1
399 IF(IZ_UPDATE >=3)IZ_UPDATE=1
400
401 !IABAC flag : MODEL TYPE
402 ! 1: FREE AIR (default)
403 ! 2: SURFACE BURST
404 ! 3: AIR BURST
405 IF(IABAC <= 0)IABAC=1
406 IF(IABAC >= 4)IABAC=1
407
408 !ISHAPE
409 ! 1: spherical (with ground) - default
410 ! 2: spherical (without ground)
411 ! 3: cylindrical (with ground)
412.OR. IF(ISHAPE <= 0 ISHAPE >= 4) ISHAPE=1
413
414 !--------------------------------------------------------------
415 ! S u r f a c e S e g m e n t s
416 !--------------------------------------------------------------
417 IF(internal_SURF_ID /= 0)THEN
418 NUMSEG = IGRSURF(internal_SURF_ID)%NSEG
419 DO J=1,NUMSEG
420 LLOADP(IAD+4*(J-1)) = IGRSURF(internal_SURF_ID)%NODES(J,1)
421 LLOADP(IAD+4*(J-1)+1) = IGRSURF(internal_SURF_ID)%NODES(J,2)
422 LLOADP(IAD+4*(J-1)+2) = IGRSURF(internal_SURF_ID)%NODES(J,3)
423 LLOADP(IAD+4*(J-1)+3) = IGRSURF(internal_SURF_ID)%NODES(J,4)
424 IF(IGRSURF(internal_SURF_ID)%ELTYP(J)==7)LLOADP(IAD+4*(J-1)+3) = 0
425 ENDDO
426 ELSE
427 NUMSEG = 0
428 ENDIF
429
430 !--------------------------------------------------------------
431 ! Normal Vector for ground (IABAC==2 or 3)
432 ! Scaled Height of Charge (IABAC==3)
433 !--------------------------------------------------------------
434 ZC = ZERO
435 NORM = ZERO
436 IF(IABAC >= 2 )THEN !AIR BURST ONLY
437 IF(SURF_ID_GROUND == 0)THEN
438 Base_X=ZERO
439 Base_y=ZERO
440 Base_z=ZERO
441 Nx_SURF=ZERO
442 Ny_SURF=ZERO
443 Nz_SURF=ZERO
444 IF(IABAC == 3)THEN
445 Nz_SURF=-Zdet
446 NN2 = Nx_SURF*Nx_SURF+Ny_SURF*Ny_SURF+Nz_SURF*Nz_SURF
447 NORM = SQRT(NN2)
448 Nz_SURF=ONE
449 ELSE
450 Base_z=Zdet
451 Nz_SURF=ONE
452 NN2=ONE
453 NORM=ONE
454 MSGOUT1=''
455 MSGOUT1='missing ground_id identifier'
456 MSGOUT2=''
457 MSGOUT2='assuming ground with basis=(0,0,zdet) and normal=(0,0,1.)'
458 CALL ANCMSG(MSGID=1907,MSGTYPE=MSGWARNING,ANMODE=ANINFO,C1=TRIM(TITR),I1=ID,C2=MSGOUT1,C3=MSGOUT2)
459 ENDIF
460.AND. IF(ZDET == ZERO IABAC == 3)THEN
461 MSGOUT1=''
462 MSGOUT1='zdet is null. it is not possible to determine default surface
for ground definition
'
463 MSGOUT2=''
464 MSGOUT2='surface identifier must be input
for ground definition
'
465 CALL ANCMSG(MSGID=2047,MSGTYPE=MSGERROR,ANMODE=ANINFO,C1=TRIM(TITR),I1=ID,C2=MSGOUT1,C3=MSGOUT2)
466 BOOL_SKIP_CALC = .TRUE.
467 Nz_SURF=-ONE ! Starter must go on up to normal termination
468 NN2=ONE
469 NORM=ONE
470 ELSE
471 IF(IABAC == 3)THEN
472 MSGOUT1=''
473 MSGOUT1='missing ground_id identifier'
474 MSGOUT2=''
475 MSGOUT2='assuming ground with basis=(0,0,0) and normal=(0,0,zdet'
476 CALL ANCMSG(MSGID=1907, MSGTYPE=MSGWARNING, ANMODE=ANINFO,C1=TRIM(TITR), I1=ID, C2=MSGOUT1,C3=MSGOUT2)
477 ENDIF
478
479 !NORM IS STRICTLY POSITIVE : BECAUSE THERE IS A STARTER CHECK FROM READER, ERROR MESSAGE 891 OTHERWISE
480
481 !vector to find ground basis point from charge
482 Nx_SURF=Nx_SURF/NORM !lambda*NX_
483 Ny_SURF=Ny_SURF/NORM !lambda*NY_
484 Nz_SURF=Nz_SURF/NORM !lambda*NZ_
485
486 IF(IABAC == 3)THEN
487 !Determine Height
488 ! find Projection along line generated by n (ground vector) and over the gound plan
489 ! Proj=lambda.n
490 ! <z-Proj,n>=0
491 lambda=(Nx_SURF*XDET + Ny_SURF*YDET + Nz_SURF*ZDET - Nx_SURF*Base_X-Ny_SURF*Base_y-Nz_SURF*Base_z)/NN2
492 !Height is length Proj->Det. Storing Det->Proj into NN array
493 HC=lambda*NORM
494 Nx_SURF=HC*Nx_SURF
495 Ny_Surf=HC*Ny_SURF
496 Nz_Surf=HC*Nz_SURF
497 ENDIF
498
499 !Scaled Height of Charge
500 ZC = HC/W13
501 ZC = ZC * FAC_L_bb
502
503 ENDIF
504 ELSE
505 lfound = .FALSE.
506 DO IS=1,NSURF
507 IF (SURF_ID_GROUND == IGRSURF(IS)%ID)THEN
508 SELECT CASE(IGRSURF(IS)%TYPE)
509 CASE DEFAULT
510 !Surface type is not matching. Setting Default ground definition
511 IF (Zdet <= ZERO)THEN
512 MSGOUT1=''
513 MSGOUT1='surface
TYPE for ground definition is not matching
'
514 MSGOUT2=''
515 MSGOUT2='expected surface TYPE:/SURF/PLANE'
516 CALL ANCMSG(MSGID=2047,MSGTYPE=MSGERROR,ANMODE=ANINFO,C1=TRIM(TITR),I1=ID,C2=MSGOUT1,C3=MSGOUT2)
517 BOOL_SKIP_CALC = .TRUE.
518 ELSE
519 MSGOUT1=''
520 MSGOUT1='EXPECTED TYPE FOR GROUND SURFACE IS /SURF/PLANE.'
521 MSGOUT2=''
522 MSGOUT2='ASSUMING GROUND WITH BASIS=(0,0,0) and normal=(0,0,zdet)'
523 CALL ANCMSG(MSGID=1907, MSGTYPE=MSGWARNING, ANMODE=ANINFO,C1=TRIM(TITR), I1=ID, C2=MSGOUT1,C3=MSGOUT2)
524 ENDIF
525 IADPL = 0
526 lfound = .TRUE.
527 Base_X = ZERO
528 Base_y = ZERO
529 Base_z = ZERO
530 HC = Zdet
531 ZC = HC * FAC_L_bb/W13
532 Nx_SURF=ZERO
533 Ny_SURF=ZERO
534 Nz_SURF=-Zdet
535 IF(ZDET == ZERO)THEN
536 Nz_SURF = -ONE
537 ENDIF
538 EXIT
539 CASE(200)
540 !Surface type is matching:setting ground from surface definition
541 IADPL = IGRSURF(IS)%IAD_BUFR
542 lfound = .TRUE.
543 Base_X = BUFSF(IADPL+1)
544 Base_y = BUFSF(IADPL+2)
545 Base_z = BUFSF(IADPL+3)
546 Nx_SURF = BUFSF(IADPL+4)-Base_X
547 Ny_SURF = BUFSF(IADPL+5)-Base_Y
548 Nz_SURF = BUFSF(IADPL+6)-Base_Z
549 NN2 = Nx_SURF*Nx_SURF+Ny_SURF*Ny_SURF+Nz_SURF*Nz_SURF
550 NORM = SQRT(NN2)
551 !NORM IS STRICTLY POSITIVE : BECAUSE THERE IS A STARTER CHECK FROM READER, ERROR MESSAGE 891 OTHERWISE
552
553 IF(NN2 == ZERO)THEN
554 !Starter must go on up to normal termination.
555 !ERROR ID : 891 already printed (surface definition)
556 Nz_SURF=-ONE
557 NN2=ONE
558 NORM=ONE
559 ENDIF
560
561 !vector to find ground basis point from charge
562 Nx_SURF=Nx_SURF/NORM !NX_
563 Ny_SURF=Ny_SURF/NORM !NY_
564 Nz_SURF=Nz_SURF/NORM !NZ_
565 IF(IABAC == 3)THEN
566 !Determine Height
567 ! find Projection along line generated by n (ground vector) and over the gound plan
568 ! Proj=lambda.n
569 ! <z-Proj,n>=0
570 lambda=(Nx_SURF*XDET + Ny_SURF*YDET + Nz_SURF*ZDET - Nx_SURF*Base_X-Ny_SURF*Base_y-Nz_SURF*Base_z)
571 !Height is length Proj->Det. Storing Det->Proj into NN array
572 HC = lambda ! |N_SURF|=1.0
573 ! N_Surf becomes vector from detonation point to ground surface
574 Nx_SURF=-HC*Nx_SURF !lambda*NX_
575 Ny_SURF=-HC*Ny_SURF !lambda*NY_
576 Nz_SURF=-HC*Nz_SURF !lambda*NZ_
577 ENDIF
578 !Scaled Height of Charge
579 ZC = HC * FAC_L_bb/W13
580 EXIT
581 END SELECT
582 ENDIF
583 ENDDO
584 ! Provided ID is referring to an nonexistent surface
585.NOT. IF (lFOUND) THEN
586 Nx_Surf=Zero
587 Ny_Surf=Zero
588 Nz_Surf=-one
589 NN2=ONE
590 NORM=ONE
591 HC=ONE
592 ZC=ONE
593 ENDIF
594 ENDIF
595 ELSE
596 !initialized to zero but not used
597 Nx_SURF=ZERO
598 Ny_SURF=ZERO
599 Nz_SURF=ZERO
600 ZC = ZERO
601 Base_X = ZERO
602 Base_y = ZERO
603 Base_z = ZERO
604 ENDIF
605
606 !check that DEtonation Point belongs the ground
607 IF(IABAC == 2)THEN
608 tmp_VAR = Nx_SURF*(Xdet-Base_X)+Ny_SURF*(Ydet-Base_Y)+Nz_SURF*(Zdet-Base_Z)
609 tmp_VAR = ABS(tmp_VAR/NORM)
610 IF(tmp_VAR > EM06)THEN
611 lambda = (Base_X-Xdet)*Nx_SURF + (Base_Y-Ydet)*Ny_SURF + (Base_Z-Zdet)*Nz_SURF
612 Xdet = Xdet+lambda*Nx_SURF
613 Ydet = Ydet+lambda*Ny_SURF
614 Zdet = Zdet+lambda*Nz_SURF
615 MSGOUT1=''
616 MSGOUT1=' detonation center must be on
the ground
'
617 MSGOUT2=' projecting(xdet,ydet,zdet) to
the ground( , , )
'
618 WRITE(MSGOUT2(47:56),FMT='(e10.3)')Xdet
619 WRITE(MSGOUT2(58:67),FMT='(e10.3)')Ydet
620 WRITE(MSGOUT2(69:78),FMT='(e10.3)')Zdet
621 CALL ANCMSG(MSGID=1907,MSGTYPE=MSGWARNING,ANMODE=ANINFO,C1=TRIM(TITR),I1=ID,C2=MSGOUT1,C3=MSGOUT2)
622 ENDIF
623 ENDIF
624
625 IF(IABAC == 3)THEN
626 !Scaled Height of triple point
627 ! curve_id 1 2 3 4 5 6 7 8 9 10
628 ! 10 curves for ZC : {1.0; 1.5; 2.0; 2.5; 3.0; 3.5; 4.0; 5.0; 6.0; 7.0}
629 ! <--------curve_id=[2x-1]---------><----[x+3]----->
630 !
631 !get curve_id and interpolation factor. No matter unit system (curve_id1 and alpha_zc are dimensionless)
632
633 IF(HC<=ZERO)THEN
634 MSGOUT1=''
635 MSGOUT1=' detonation center must be above
the ground
'
636 MSGOUT2=''
637 CALL ANCMSG(MSGID=2047,MSGTYPE=MSGERROR,ANMODE=ANINFO,C1=TRIM(TITR),I1=ID,C2=MSGOUT1,C3=MSGOUT2)
638 BOOL_SKIP_CALC = .TRUE.
639 HC = ONE
640 ENDIF
641
642 tmp(1)=ZC
643 ZC=ZC/FAC_unit
644 IF(ZC < 4)THEN
645 itmp = INT( 2.*(ZC)-1. )
646 curve_id1 = max(1,itmp)
647 curve_id2 = curve_id1+1
648 IF(itmp < 1)THEN
649 !message out of bounds. curve 1 will be used. no extrapolation
650 alpha_zc = ZERO
651 ELSE
652 alpha_zc = (ZC - PBLAST%PBLAST_DATA%Curve_val_2_13(curve_id1)) /
653 . ( PBLAST%PBLAST_DATA%Curve_val_2_13(curve_id2) - PBLAST%PBLAST_DATA%Curve_val_2_13(curve_id1) )
654 ENDIF
655 ELSE
656 itmp = INT( ZC+3. )
657 curve_id1 = INT( min(10,itmp) )
658 curve_id2 = curve_id1+1
659 IF(curve_id1 == 10)THEN
660 !message out of bounds. curve 10 will be used. no extrapolation
661 alpha_zc = ZERO
662 ELSE
663 alpha_zc = (ZC - PBLAST%PBLAST_DATA%Curve_val_2_13(curve_id1)) /
664 . ( PBLAST%PBLAST_DATA%Curve_val_2_13(curve_id2) - PBLAST%PBLAST_DATA%Curve_val_2_13(curve_id1) )
665 ENDIF
666 ENDIF
667 alpha_zc = curve_id1+alpha_zc !integer part is curve id1, digits are standing for interpolation factor between curve_id1 & curve_id2
668 ENDIF
669 ZC = tmp(1)
670
671 !--------------------------------------------------------------
672 ! B u f f e r S t o r a g e (Transmitted to Engine)
673 !--------------------------------------------------------------
674 IF(internal_SURF_ID /= 0)THEN
675 ILOADP( 1,K) = 4*NUMSEG
676 ILOADP( 2,K) = internal_SURF_ID
677 ILOADP( 3,K) = ISHAPE
678 ILOADP( 4,K) = IAD
679 ILOADP( 5,K) = 3
680 ILOADP( 6,K) = IZ_UPDATE
681 ILOADP( 7,K) = IABAC
682 ILOADP( 8,K) = ID
683 ILOADP( 9,K) = ITA_SHIFT
684 ILOADP(10,K) = NDT
685 ILOADP(11,K) = IMODEL
686 !ILOADP(12,K) = curve_id1 ! getting interpolated SHTP from figure 2_13
687 IAD = IAD + 3*NUMSEG
688 FACLOADP( 1,K) = TDET
689 FACLOADP( 2,K) = XDET
690 FACLOADP( 3,K) = YDET
691 FACLOADP( 4,K) = ZDET
692 FACLOADP( 5,K) = WTNT
693 FACLOADP( 6,K) = PMIN
694 FACLOADP( 7,K) = ZERO !TA_PBLAST initialized below
695 FACLOADP( 8,K) = Nx_SURF
696 FACLOADP( 9,K) = Ny_SURF
697 FACLOADP(10,K) = Nz_SURF
698 HC = ABS(HC)
699 FACLOADP(11,K) = HC
700 FACLOADP(12,K) = alpha_zc ! inteprolation factor for curve if on figure 2_13 in [0,1[ + curve_id1 (no matter are units ft/lb**1/3 or cm/g**1/3)
701 FACLOADP(13,K) = TSTOP
702 ENDIF
703
704 !--------------------------------------------------------------
705 ! C o m p u t e W a v e P a r a m e t e r s
706 !--------------------------------------------------------------
707.AND..NOT. IF(internal_SURF_ID /= 0 BOOL_SKIP_CALC)THEN
708 NUMSEG = IGRSURF(internal_SURF_ID)%NSEG
709 IAD = ILOADP(4,K)
710 ILD_PBLAST = ILD_PBLAST + 1
711 ALLOCATE (PBLAST%PBLAST_TAB(ILD_PBLAST)%cos_theta(NUMSEG) ,STAT=IERR1); IF (IERR1 /= 0) GOTO 1000
712 ALLOCATE (PBLAST%PBLAST_TAB(ILD_PBLAST)%P_inci(NUMSEG) ,STAT=IERR1); IF (IERR1 /= 0) GOTO 1000
713 ALLOCATE (PBLAST%PBLAST_TAB(ILD_PBLAST)%P_refl(NUMSEG) ,STAT=IERR1); IF (IERR1 /= 0) GOTO 1000
714 ALLOCATE (PBLAST%PBLAST_TAB(ILD_PBLAST)%ta(NUMSEG) ,STAT=IERR1); IF (IERR1 /= 0) GOTO 1000
715 ALLOCATE (PBLAST%PBLAST_TAB(ILD_PBLAST)%t0(NUMSEG) ,STAT=IERR1); IF (IERR1 /= 0) GOTO 1000
716 ALLOCATE (PBLAST%PBLAST_TAB(ILD_PBLAST)%decay_inci(NUMSEG) ,STAT=IERR1); IF (IERR1 /= 0) GOTO 1000
717 ALLOCATE (PBLAST%PBLAST_TAB(ILD_PBLAST)%decay_refl(NUMSEG) ,STAT=IERR1); IF (IERR1 /= 0) GOTO 1000
718 ALLOCATE (PBLAST%PBLAST_TAB(ILD_PBLAST)%TAGMSG(NUMSEG) ,STAT=IERR1); IF (IERR1 /= 0) GOTO 1000
719 PBLAST%PBLAST_TAB(ILD_PBLAST)%SIZ=NUMSEG
720
721 ! Coordinates of detonation point projection = DETPOINT + HC*n !N_SURF is here -HC*n where n is normal unitary vector of the plane
722 IF(IABAC == 3)THEN
723 ProjDet(1)=Xdet+ Nx_SURF
724 ProjDet(2)=Ydet+ Ny_SURF
725 ProjDet(3)=Zdet+ Nz_SURF
726 ELSE
727 ProjDet(1) = Xdet
728 ProjDet(2) = Ydet
729 ProjDet(3) = Zdet
730 ENDIF
731
732 b_min = EP20
733 b_max = -EP20
734 DTMIN = EP20
735 DO I=1,NUMSEG
736 N1=IGRSURF(internal_SURF_ID)%NODES(I,1)
737 N2=IGRSURF(internal_SURF_ID)%NODES(I,2)
738 N3=IGRSURF(internal_SURF_ID)%NODES(I,3)
739 N4=IGRSURF(internal_SURF_ID)%NODES(I,4)
740.AND..AND..AND..AND. IF(N4 /= 0 N1 /= N2 N2 /= N3 N3 /= N4 N4 /= N1 )THEN
741 ! 4 NODE SEGMENT
742 NPt = FOUR
743 ! Segment Zentrum
744 Zx = X(1,N1)+X(1,N2)+X(1,N3)+X(1,N4)
745 Zy = X(2,N1)+X(2,N2)+X(2,N3)+X(2,N4)
746 Zz = X(3,N1)+X(3,N2)+X(3,N3)+X(3,N4)
747 Zx = Zx*FOURTH
748 Zy = Zy*FOURTH
749 Zz = Zz*FOURTH
750 ! Normal vectors
751 Nx_SEG = (X(2,N3)-X(2,N1))*(X(3,N4)-X(3,N2)) - (X(3,N3)-X(3,N1))*(X(2,N4)-X(2,N2))
752 Ny_SEG = (X(3,N3)-X(3,N1))*(X(1,N4)-X(1,N2)) - (X(1,N3)-X(1,N1))*(X(3,N4)-X(3,N2))
753 Nz_SEG = (X(1,N3)-X(1,N1))*(X(2,N4)-X(2,N2)) - (X(2,N3)-X(2,N1))*(X(1,N4)-X(1,N2))
754 NORM_= SQRT(Nx_SEG*Nx_SEG+Ny_SEG*Ny_SEG+Nz_SEG*Nz_SEG)
755 ELSE
756 ! 3 NODE SEGMENT
757 NPt = THREE
758 IF(N4==0)THEN
759 ! nothing to do
760 ELSEIF(N1 == N2)THEN
761 N2 = N3
762 N3 = N4
763 N4 = 0
764 ELSEIF(N2 == N3)THEN
765 N3 = N4
766 N4 = 0
767 ELSEIF(N3 == N4)THEN
768 N4 = 0
769 ELSEIF(N4 == N1)THEN
770 N4 = 0
771 ENDIF
772 ! Segment Zentrum
773 Zx = X(1,N1)+X(1,N2)+X(1,N3)
774 Zy = X(2,N1)+X(2,N2)+X(2,N3)
775 Zz = X(3,N1)+X(3,N2)+X(3,N3)
776 Zx = Zx*THIRD
777 Zy = Zy*THIRD
778 Zz = Zz*THIRD
779 Nx_SEG = (X(2,N3)-X(2,N1))*(X(3,N3)-X(3,N2)) - (X(3,N3)-X(3,N1))*(X(2,N3)-X(2,N2))
780 Ny_SEG = (X(3,N3)-X(3,N1))*(X(1,N3)-X(1,N2)) - (X(1,N3)-X(1,N1))*(X(3,N3)-X(3,N2))
781 Nz_SEG = (X(1,N3)-X(1,N1))*(X(2,N3)-X(2,N2)) - (X(2,N3)-X(2,N1))*(X(1,N3)-X(1,N2))
782 NORM_ = SQRT(Nx_SEG*Nx_SEG+Ny_SEG*Ny_SEG+Nz_SEG*Nz_SEG)
783 ENDIF
784
785 ! Dist
786 Dx = (Xdet - Zx)*FAC_L_bb
787 Dy = (Ydet - Zy)*FAC_L_bb
788 Dz = (Zdet - Zz)*FAC_L_bb
789 DNORM = SQRT(Dx*Dx+Dy*Dy+Dz*Dz) ! *FAC_L_bb cm->work unit /FAC_L_bb : WORK_UNIT -> cm
790
791 !init
792 P_inci = zero
793 I_inci = zero
794 P_refl = zero
795 I_refl = zero
796 DT_0 = zero
797 T_A = zero
798 cos_theta = zero
799 DECAY_inci = ONE
800 DECAY_refl = ONE
801 BOOL_UNDERGROUND_CURRENT_SEG = .FALSE.
802
803 !--------------------------------------------------------
804 ! Free Air
805 !--------------------------------------------------------
806 IF( IABAC == 1 ) THEN
807 !=== SPHERICAL CHARGE IN FREE FIELD ===!
808
809 ! scaled distance
810 Z = DNORM / W13 !in abac unit ID g,cm,mus
811
812 !Warning Message if Z is out of bounds (precondition for subroutine BLAST_PARAMETERS__FREE_AIR_MODEL)
813.AND. IF(Z > 0.5 Z < 400.)THEN
814 ELSEIF(Z <= 0.5)THEN
815 WRITE(IOUT, FMT='(a,i0,a)') "Warning : /LOAD/PBLAST id=",ID," R/W**(1/3) < 0.5 cm/g**(1/3)"
816 WRITE(ISTDO,FMT='(a,i0,a)') "Warning : /LOAD/PBLAST id=",ID," R/W**(1/3) < 0.5 cm/g**(1/3)"
817 WRITE(IOUT, FMT='(a)') " Radial Distance too close to the charge"
818 WRITE(ISTDO,FMT='(a)') " Radial Distance too close to the charge"
819.OR. IF (N4 == 0 N3 == N4)THEN
820 WRITE(IOUT, FMT='(a,3i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3)
821 WRITE(ISTDO,FMT='(a,3i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3)
822 ELSE
823 WRITE(IOUT, FMT='(a,4i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3),ITAB(N4)
824 WRITE(ISTDO,FMT='(a,4i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3),ITAB(N4)
825 ENDIF
826 PBLAST%PBLAST_TAB(ILD_PBLAST)%TAGMSG(I) = 1
827
828 ELSEIF(Z > 400.)THEN
829 WRITE(IOUT, FMT='(a,i0,a)') "Warning : /LOAD/PBLAST id=",ID," R/W**(1/3) > 400. cm/g**(1/3)"
830 WRITE(ISTDO,FMT='(a,i0,a)') "Warning : /LOAD/PBLAST id=",ID," R/W**(1/3) > 400. cm/g**(1/3)"
831 WRITE(IOUT, FMT='(a)') " Radial Distance too far from the charge"
832 WRITE(ISTDO, FMT='(a)') " Radial Distance too far from the charge"
833.OR. IF (N4 == 0 N3 == N4)THEN
834 WRITE(IOUT, FMT='(a,3i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3)
835 WRITE(ISTDO,FMT='(a,3i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3)
836 ELSE
837 WRITE(IOUT, FMT='(a,4i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3),ITAB(N4)
838 WRITE(ISTDO,FMT='(a,4i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3),ITAB(N4)
839 ENDIF
840 PBLAST%PBLAST_TAB(ILD_PBLAST)%TAGMSG(I) = 1
841 ENDIF
842
843 !Angle from detonation point
844 cos_theta = Dx*Nx_SEG + Dy*Ny_SEG + Dz*Nz_SEG
845 cos_theta = cos_theta/MAX(EM20,NORM_*DNORM)
846
847 HZ = ZERO ! not used with this formulation
848 !------------------------------------------------------------!
849 CALL PBLAST_PARAMETERS__FREE_AIR(PBLAST,Z, W13, TDET,
850 + FAC_P_bb, FAC_I_bb, FAC_T_bb,
851 + IS_DECAY_TO_BE_COMPUTED,
852 + FRIEDLANDER_PARAMS, NWARN)
853
854 P_inci = FRIEDLANDER_PARAMS%P_inci
855 P_refl = FRIEDLANDER_PARAMS%P_refl
856 I_inci = FRIEDLANDER_PARAMS%I_inci
857 I_refl = FRIEDLANDER_PARAMS%I_refl
858 T_A = FRIEDLANDER_PARAMS%T_A
859 DT_0 = FRIEDLANDER_PARAMS%DT_0
860 decay_inci = FRIEDLANDER_PARAMS%decay_inci
861 decay_refl = FRIEDLANDER_PARAMS%decay_refl
862 !------------------------------------------------------------!
863
864
865 !--------------------------------------------------------
866 ! Surface Burst
867 !--------------------------------------------------------
868 ELSEIF( IABAC == 2 ) THEN
869
870 !=== HEMISPHERICAL CHARGE WITH GROUND REFLECTION ===!
871
872 HZ = ( Nx_SURF*Zx + Ny_SURF*Zy + Nz_SURF*Zz - Nx_SURF*Xdet - Ny_SURF*Ydet - Nz_SURF*Zdet )
873
874.AND. IF(HZ < -em10 ISHAPE /= 2)THEN
875
876 !FACE IS UNDERGROUND
877 P_inci = zero
878 I_inci = zero
879 P_refl = zero
880 I_refl = zero
881 DT_0 = EP20
882 T_A = EP20
883 DECAY_inci = ONE
884 DECAY_refl = ONE
885 cos_theta = ZERO
886
887 SEG_UNDERGROUND = SEG_UNDERGROUND + 1
888 BOOL_UNDERGROUND_CURRENT_LOAD = .TRUE.
889 BOOL_UNDERGROUND_CURRENT_SEG = .TRUE.
890
891 ELSE
892
893 !FACE IS OVERGROUND
894
895 !Planar Wave : angle with targeted face
896 lambda = (Base_X-Zx)*Nx_Surf + (Base_Y-Zy)*Ny_Surf + (Base_Z-Zz)*Nz_Surf
897
898 ProjZ(1) = Zx + lambda*Nx_SURF
899 ProjZ(2) = Zy + lambda*Ny_SURF
900 ProjZ(3) = Zz + lambda*Nz_SURF
901
902 tmp(1) = (ProjZ(1)-ProjDet(1))
903 tmp(2) = (ProjZ(2)-ProjDet(2))
904 tmp(3) = (ProjZ(3)-ProjDet(3))
905 LG = SQRT(TMP(1)*TMP(1)+TMP(2)*TMP(2)+TMP(3)*TMP(3))
906 ZG = LG*FAC_L_bb/W13
907
908 ! scaled distance
909 SELECT CASE(ISHAPE)
910 CASE (1,2)
911 Z = DNORM / W13 !spherical
912 !Angle from detonation point
913 cos_theta = Dx*Nx_SEG + Dy*Ny_SEG + Dz*Nz_SEG
914 cos_theta = cos_theta/MAX(EM20,NORM_*DNORM)
915 CASE(3)
916 Z = ZG !cylindrical - in abac unit ID g,cm,mus
917 !Angle from detonation point
918 ! is here cos between (Det-ProjZ) and normal
919 cos_theta = (ProjDet(1)-ProjZ(1))*Nx_SEG + (ProjDet(2)-ProjZ(2))*Ny_SEG + (ProjDet(3)-ProjZ(3))*Nz_SEG
920 cos_theta = cos_theta/MAX(EM20,LG*NORM_)
921 END SELECT
922
923 !finding index for TM5-1300 abacuses from bijection.
924.AND. IF(Z > 0.5 Z < 400.)THEN
925 !inside interpolation table
926
927 ELSEIF(Z <= 0.5)THEN
928 WRITE(IOUT, FMT='(a,i0,a)')"Warning : /LOAD/PBLAST id=",ID," Rg/W**(1/3) < 0.5 cm/g**(1/3)"
929 WRITE(ISTDO,FMT='(a,i0,a)')"Warning : /LOAD/PBLAST id=",ID," Rg/W**(1/3) < 0.5 cm/g**(1/3)"
930 WRITE(IOUT, FMT='(a)') " Horizontal Distance on Ground (Rg) is too close to the charge"
931 WRITE(ISTDO,FMT='(a)') " Horizontal Distance on Ground (Rg) is too close to the charge"
932.OR. IF (N4 == 0 N3 == N4)THEN
933 WRITE(IOUT, FMT='(a,3i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3)
934 WRITE(ISTDO,FMT='(a,3i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3)
935 ELSE
936 WRITE(IOUT, FMT='(a,4i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3),ITAB(N4)
937 WRITE(ISTDO,FMT='(a,4i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3),ITAB(N4)
938 ENDIF
939 PBLAST%PBLAST_TAB(ILD_PBLAST)%TAGMSG(I) = 1
940
941 ELSEIF(Z > 400.)THEN
942 WRITE(IOUT, FMT='(a,i0,a)')"Warning : /LOAD/PBLAST id=",ID," Rg/W**(1/3) > 400. cm/g**(1/3)"
943 WRITE(ISTDO,FMT='(a,i0,a)')"Warning : /LOAD/PBLAST id=",ID," Rg/W**(1/3) > 400. cm/g**(1/3)"
944 WRITE(IOUT, FMT='(a)') " Horizontal Distance on Ground (Rg) is too far from the charge"
945 WRITE(ISTDO,FMT='(a)') " Horizontal Distance on Ground (Rg) is too far from the charge"
946.OR. IF (N4 == 0 N3 == N4)THEN
947 WRITE(IOUT, FMT='(a,3i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3)
948 WRITE(ISTDO,FMT='(a,3i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3)
949 ELSE
950 WRITE(IOUT, FMT='(a,4i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3),ITAB(N4)
951 WRITE(ISTDO,FMT='(a,4i11)') " -> Segment nodes : ", ITAB(N1),ITAB(N2),ITAB(N3),ITAB(N4)
952 ENDIF
953 PBLAST%PBLAST_TAB(ILD_PBLAST)%TAGMSG(I) = 1
954 ENDIF
955 !------------------------------------------------------------!
956 CALL PBLAST_PARAMETERS__SURFACE_BURST( PBLAST,Z, W13, TDET,
957 + FAC_P_bb, FAC_I_bb, FAC_T_bb,
958 + IS_DECAY_TO_BE_COMPUTED,
959 + FRIEDLANDER_PARAMS,NWARN )
960 P_inci = FRIEDLANDER_PARAMS%P_inci
961 P_refl = FRIEDLANDER_PARAMS%P_refl
962 I_inci = FRIEDLANDER_PARAMS%I_inci
963 I_refl = FRIEDLANDER_PARAMS%I_refl
964 T_A = FRIEDLANDER_PARAMS%T_A
965 DT_0 = FRIEDLANDER_PARAMS%DT_0
966 decay_inci = FRIEDLANDER_PARAMS%decay_inci
967 decay_refl = FRIEDLANDER_PARAMS%decay_refl
968 !------------------------------------------------------------!
969 ENDIF
970
971 !--------------------------------------------------------
972 ! Free Air Burst
973 !--------------------------------------------------------
974 ELSEIF(IABAC == 3)THEN
975 !--- Determine Height of centroid (structure face) : HZ
976
977 !Height is length Proj->Det. Storing Det->Proj into NN array
978 HZ=-(Nx_SURF*Zx + Ny_SURF*Zy + Nz_SURF*Zz - Nx_SURF*ProjDet(1)-Ny_SURF*ProjDet(2)-Nz_SURF*ProjDet(3))/HC
979
980 IF(HZ < ZERO)THEN
981
982 P_inci = zero
983 I_inci = zero
984 P_refl = zero
985 I_refl = zero
986 DT_0 = EP20
987 T_A = EP20
988 DECAY_inci = ONE
989 DECAY_refl = ONE
990 cos_theta = ZERO
991
992 SEG_UNDERGROUND = SEG_UNDERGROUND + 1
993 BOOL_UNDERGROUND_CURRENT_LOAD = .TRUE.
994 BOOL_UNDERGROUND_CURRENT_SEG = .TRUE.
995
996 ELSE
997
998 Z_struct = HZ*FAC_L_bb/W13 ! scaled value in same unit system as hardcoded tables {g,cm,mus,Mbar}
999
1000 !Scaled Height of Charge
1001 ZC = HC * FAC_L_bb/W13! scaled value in same unit system as hardcoded tables {g,cm,mus,Mbar}
1002 ZC = ZC/FAC_UNIT ! ft/lb**0.3333 (fig 2-13 : 10 plots, one for a given ZC value)
1003
1004 !Horizontal Distance between Charge and Centroid : LG
1005 ! ZG = scaled distance |ProjC->ProjZ|
1006 ProjZ(1) = Zx + HZ*Nx_SURF/HC
1007 ProjZ(2) = Zy + HZ*Ny_SURF/HC
1008 ProjZ(3) = Zz + HZ*Nz_SURF/HC
1009 tmp(1) = (ProjZ(1)-ProjDet(1))
1010 tmp(2) = (ProjZ(2)-ProjDet(2))
1011 tmp(3) = (ProjZ(3)-ProjDet(3))
1012 LG = SQRT(TMP(1)*TMP(1)+TMP(2)*TMP(2)+TMP(3)*TMP(3))
1013 ZG = LG*FAC_L_bb/W13 !scaled horizontal distance (ground) ; same unit system as hardcoded tables {g,cm,mus,Mbar}
1014
1015 !Angle of structural face (mach wave is planar wave)
1016 cos_theta = (ProjDet(1)-ProjZ(1))*Nx_SEG + (ProjDet(2)-ProjZ(2))*Ny_SEG + (ProjDet(3)-ProjZ(3))*Nz_SEG
1017 cos_theta = cos_theta/MAX(EM20,LG*NORM_)
1018
1019 !determine angle of incidence at ground (AANGLE_g) and interpolation factor (alpha_angle)
1020 tmp(1)=Xdet-ProjZ(1)
1021 tmp(2)=Ydet-ProjZ(2)
1022 tmp(3)=Zdet-ProjZ(3)
1023 tmp_var=SQRT( tmp(1)*tmp(1) + tmp(2)*tmp(2) + tmp(3)*tmp(3) )
1024 ANGLE_g = -( Nx_SURF*tmp(1) + Ny_SURF*tmp(2) + Nz_SURF*tmp(3) ) /Hc/tmp_var !must be between [-PI_,PI_]
1025 ANGLE_g = min(ONE,max(-ONE,ANGLE_g)) ! bound it to expected range (epsilon)
1026 ANGLE_g = acos(ANGLE_g)
1027 ANGLE_g = cst_180/PI_*ANGLE_g !debug purpose
1028 IF(ANGLE_g < ZERO)THEN
1029 WRITE(IOUT,*) 'warning : /load/pblast
id=
',ID,' negative angle,
',ANGLE_g
1030 WRITE(ISTDO,*) 'warning : /load/pblast
id=
',ID,' negative angle,
',ANGLE_g
1031.OR. IF(N4 == 0 N3 == N4)THEN
1032 WRITE(IOUT,FMT= '(a,3i11)
')' face:
',ITAB((/N1,N2,N3/)),' seems to be below
the ground
'
1033 WRITE(ISTDO,FMT='(a,3i11)
')' face:
',ITAB((/N1,N2,N3/)),' seems to be below
the ground
'
1034 ELSE
1035 WRITE(IOUT,FMT= '(a,4i11)
')' face:
',ITAB((/N1,N2,N3,N4/)),' seems to be below
the ground
'
1036 WRITE(ISTDO,FMT='(a,4i11)
')' face:
',ITAB((/N1,N2,N3,N4/)),' seems to be below
the ground
'
1037 ENDIF
1038 ANGLE_g = ZERO
1039 PBLAST%PBLAST_TAB(ILD_PBLAST)%TAGMSG(I) = 1
1040 ELSEIF(ANGLE_g > 85.00000)THEN
1041 WRITE(IOUT, FMT='(a,i0,a,e10.4)')
1042 . 'warning : /load/pblast
id=
',ID,' angle is over
the upper bound,
',ANGLE_g
1043 WRITE(ISTDO,FMT='(a,i0,a,e10.4)')
1044 . 'warning : /load/pblast
id=
',ID,' angle is over
the upper bound
',ANGLE_g
1045.OR. IF(N4 == 0 N3 == N4)THEN
1046 WRITE(IOUT, FMT='(a,3i11)')' angle set to 85.00 for face:',ITAB((/N1,N2,N3/))
1047 WRITE(ISTDO,FMT='(a,3i11)')' angle set to 85.00 for face:',ITAB((/N1,N2,N3/))
1048 ELSE
1049 WRITE(IOUT, FMT='(a,4i11)')' angle set to 85.00 for face:',ITAB((/N1,N2,N3,N4/))
1050 WRITE(ISTDO,FMT='(a,4i11)')' angle set to 85.00 for face:',ITAB((/N1,N2,N3,N4/))
1051 ENDIF
1052 ANGLE_g = 85.00000
1053 PBLAST%PBLAST_TAB(ILD_PBLAST)%TAGMSG(I) = 1
1054 ENDIF
1055
1056 !------------------------------------------------------------!
1057 CALL PBLAST_PARAMETERS__AIR_BURST(PBLAST,Z_struct, Zc, Zg, ANGLE_g, W13, TDET,
1058 + FAC_P_bb, FAC_I_bb, FAC_T_bb,
1059 + IS_DECAY_TO_BE_COMPUTED,
1060 + ID,'load',.TRUE.,
1061 + FRIEDLANDER_PARAMS,NWARN)
1062
1063 P_inci = FRIEDLANDER_PARAMS%P_inci
1064 P_refl = FRIEDLANDER_PARAMS%P_refl
1065 I_inci = FRIEDLANDER_PARAMS%I_inci
1066 I_refl = FRIEDLANDER_PARAMS%I_refl
1067 T_A = FRIEDLANDER_PARAMS%T_A
1068 DT_0 = FRIEDLANDER_PARAMS%DT_0
1069 decay_inci = FRIEDLANDER_PARAMS%decay_inci
1070 decay_refl = FRIEDLANDER_PARAMS%decay_refl
1071 !------------------------------------------------------------!
1072
1073
1074 ENDIF!(HZ >= ZERO)
1075
1076 ENDIF! IABAC==3
1077
1078 TA_PBLAST = MIN(TA_PBLAST, T_A)
1079
1080 !--------------------------------------------------------
1081 ! Storage of Wave Parameters (transmitted to Engine)
1082 !--------------------------------------------------------
1083 PBLAST%PBLAST_TAB(ILD_PBLAST)%cos_theta(I) = cos_theta
1084 PBLAST%PBLAST_TAB(ILD_PBLAST)%P_inci(I) = P_inci
1085 PBLAST%PBLAST_TAB(ILD_PBLAST)%P_refl(I) = P_refl
1086 PBLAST%PBLAST_TAB(ILD_PBLAST)%ta(I) = T_A
1087 PBLAST%PBLAST_TAB(ILD_PBLAST)%t0(I) = DT_0
1088 PBLAST%PBLAST_TAB(ILD_PBLAST)%decay_inci(I) = decay_inci
1089 PBLAST%PBLAST_TAB(ILD_PBLAST)%decay_refl(I) = decay_refl
1090 PBLAST%PBLAST_TAB(ILD_PBLAST)%TAGMSG(I) = 0
1091
1092 b_min = MIN(b_min, decay_inci)
1093 b_min = MIN(b_min, decay_refl)
1094
1095 b_max = MAX(b_max, decay_inci)
1096 b_max = MAX(b_max, decay_refl)
1097
1098 DTMIN = MIN(DTMIN, DT_0/NDT)
1099
1100 ENDDO! I=1,NUMSEG
1101
1102 PBLAST%PBLAST_TAB(ILD_PBLAST)%DTMIN = DTMIN
1103
1104 FACLOADP(7,K) = TA_PBLAST ! min(TDET+T_A, I=1..NUMSEG)
1105 PBLAST%PBLAST_DT%TA_INF = MIN(PBLAST%PBLAST_DT%TA_INF, TA_PBLAST)! min( [min(TDET+T_A, I=1..NUMSEG)] , K=1..'numpblast')
1106 OUTPUT_USER_PARAMS(3,K_BLAST) = b_min
1107 OUTPUT_USER_PARAMS(4,K_BLAST) = b_max
1108 OUTPUT_USER_PARAMS(6,K_BLAST) = TA_PBLAST
1109 IF(BOOL_UNDERGROUND_CURRENT_LOAD)THEN
1110 MSGOUT1=''
1111 WRITE(MSGOUT1,FMT='(i0,a)') SEG_UNDERGROUND,' segment(s) on TARGET surface'
1112 MSGOUT2=''
1113 MSGOUT2='there will not be loaded with blast pressure'
1114 CALL ANCMSG(MSGID=1907,MSGTYPE=MSGWARNING,ANMODE=ANINFO,C1=TRIM(TITR),I1=ID,C2=MSGOUT1,C3=MSGOUT2)
1115 ENDIF
1116.AND. IF(IMODEL == 2 NWARN > 0)THEN
1117 MSGOUT1=''
1118 WRITE(MSGOUT1,FMT='(i0,a)') NWARN,
1119 . ' segment(s) has excessive positive impulse regarding
the peak pressure and positive duration.
'
1120 MSGOUT2=''
1121 MSGOUT2=
1122 . 'a triangular waveform will be used instead to maximize
the impulse. defining a pmin
VALUE is strongly recommended
'
1123 CALL ANCMSG(MSGID=1907,MSGTYPE=MSGWARNING,ANMODE=ANINFO,C1=TRIM(TITR),I1=ID,C2=MSGOUT1,C3=MSGOUT2)
1124 ENDIF
1125
1126
1127 ENDIF
1128
1129 NUMLOADP=NUMLOADP+4*NUMSEG ! /LOAD/PBLAST
1130
1131
1132 ENDDO!next K (/LOAD/PBLAST)
1133
1134!Set (shifted arrival time) for all blast loading
1135.NOT. IF(IS_ITA_SHIFT)THEN
1136 PBLAST%PBLAST_DT%TA_INF = ZERO
1137 ELSE
1138 DO K = 1+NLOADP_F, NLOADP_F+PBLAST%NLOADP_B
1139 K_BLAST = K-NLOADP_F
1140 ILOADP(9,K) = 2
1141 TA_PBLAST = FACLOADP(7,K)
1142 !shifting trigger time
1143 ! arrival time for segments are not shifted, time will be translated
1144 ! during engine using T* = T + T_shift
1145 FACLOADP(7,K) = MAX(ZERO, TA_PBLAST-PBLAST%PBLAST_DT%TA_INF) !arrival time shifted if requested
1146 ENDDO
1147 ENDIF
1148
1149
1150
1151
1152
1153 !this output is made once all /LOAD/PBLAST options are read in order to display global parameters (min values)
1154 DO K = 1+NLOADP_F, NLOADP_F+PBLAST%NLOADP_B
1155 K_BLAST = K-NLOADP_F
1156 IZ_UPDATE = ILOADP(6,K)
1157 IABAC = ILOADP(7,K)
1158 ID = ILOADP(8,K)
1159 ITA_SHIFT = ILOADP(9,K)
1160 NDT = ILOADP(10,K)
1161 IMODEL = ILOADP(11,K)
1162 TDET = FACLOADP(1,K)
1163 XDET = FACLOADP(2,K)
1164 YDET = FACLOADP(3,K)
1165 ZDET = FACLOADP(4,K)
1166 WTNT = FACLOADP(5,K)
1167 PMIN = FACLOADP(6,K)
1168 TA_PBLAST = FACLOADP(7,K)
1169 HC = FACLOADP(11,K)
1170 TSTOP = FACLOADP(13,K)
1171 NODE_ID = INT(OUTPUT_USER_PARAMS(5,K_BLAST))
1172 SURF_ID = INT(OUTPUT_USER_PARAMS(1,K_BLAST))
1173 SURF_ID_GROUND = INT(OUTPUT_USER_PARAMS(2,K_BLAST))
1174 DTMIN = PBLAST%PBLAST_TAB(ILD_PBLAST)%DTMIN
1175 b_min = OUTPUT_USER_PARAMS(3,K_BLAST)
1176 b_max = OUTPUT_USER_PARAMS(4,K_BLAST)
1177 WRITE (IOUT,2000)
1178 SELECT CASE (IABAC)
1179 CASE(1)
1180 WRITE(IOUT,2001)
1181 CASE(2)
1182 WRITE(IOUT,2002)
1183 CASE(3)
1184 WRITE(IOUT,2003)
1185 END SELECT
1186 WRITE(IOUT,2010)
1187 IF(IABAC==1)THEN
1188 WRITE (IOUT,2011)ID,SURF_ID,NUMSEG,IABAC,IZ_UPDATE,IMODEL,ITA_SHIFT,XDET,YDET,ZDET,TDET,WTNT,PMIN
1189 ELSE
1190 WRITE (IOUT,2211)ID,SURF_ID,NUMSEG,SEG_UNDERGROUND, IABAC,IZ_UPDATE,IMODEL,ITA_SHIFT,
1191 . XDET,YDET,ZDET,TDET,WTNT,PMIN
1192 ENDIF
1193
1194 WRITE (IOUT,2033)DTMIN
1195 IF(IMODEL == 2)THEN
1196 WRITE (IOUT,2031)b_min
1197 WRITE (IOUT,2032)b_max
1198 ENDIF
1199 IF(IS_AVAILABLE_TSTOP)WRITE (IOUT,2024)TSTOP
1200 IF(NODE_ID /= 0)WRITE (IOUT,2015)NODE_ID
1201 IF(NDT /= 0)WRITE (IOUT,2014)NDT
1202 IF(IABAC >=2)THEN
1203 WRITE(IOUT,2034)NX_SURF,NY_SURF,NZ_SURF
1204 ENDIF
1205 IF(IABAC == 2)THEN
1206 WRITE(IOUT,2036)ISHAPE
1207 ENDIF
1208 IF(IABAC == 3) WRITE(IOUT,2023) HC
1209 IF(ITA_SHIFT == 2)THEN
1210 WRITE (IOUT,2013)OUTPUT_USER_PARAMS(6,K_BLAST)
1211 WRITE (IOUT,2012)TA_PBLAST
1212 ELSE
1213 WRITE (IOUT,2013)TA_PBLAST
1214 ENDIF
1215 !DISPLAY MINIMUM VALUE WITHIN ALL PBLAST OPTIONS
1216 IF(PBLAST%NLOADP_B>=2)THEN
1217 WRITE(IOUT,2035)PBLAST%PBLAST_DT%TA_INF
1218 ENDIF
1219
1220 ENDDO!next K (/LOAD/PBLAST)
1221
1222
1223
1224
1225 DEALLOCATE( OUTPUT_USER_PARAMS)
1226 RETURN
1227
1228 2000 FORMAT(
1229 . //
1230 . ' blast loading '/,
1231 . ' -------------- '/)
1232 2001 FORMAT(' free air burst')
1233 2002 FORMAT(' surface burst')
1234 2003 FORMAT(' air burst')
1235 2010 FORMAT(' default unit system is {g,cm,mus}'/)
1236 2011 FORMAT(
1237 . 5X,'id . . . . . . . . . . . . . . . . . .=
',I16/,
1238 . 5X,'surface identifier . . . . . . . . . .=',I16/,
1239 . 5X,' > number of segments . . . . . . .=',I16/,
1240 . 5X,'exp_data(algorithm) . . . . . . . . =',I16/,
1241 . 5X,'iz flag(engine update) . . . . . . . =',I16/,
1242 . 5X,'iform flag(friedlander model) . . . .=',I16/,
1243 . 5X,'i_tshift(flag for time shifting) . .=',I16/,
1244 . 5X,'x-det. . . . . . . . . . . . . . . . .=',E12.4/,
1245 . 5X,'y-det . . . . . . . . . . . . . . . . =',E12.4/,
1246 . 5X,'z-det . . . . . . . . . . . . . . . . =',E12.4/,
1247 . 5X,'tdet . . . . . . . . . . . . . . . . .=',E12.4/,
1248 . 5X,'wtnt . . . . . . . . . . . . . . . .=',E12.4/,
1249 . 5X,'pmin . . . . . . . . . . . . . . . . .=',E12.4)
1250 2211 FORMAT(
1251 . 5X,'id . . . . . . . . . . . . . . . . . .=
',I16/,
1252 . 5X,'surface identifier . . . . . . . . . .=',I16/,
1253 . 5X,' > number of segments . . . . . . .=',I16/,
1254 . 5X,' > underground segments. . . . . . .=',I16/,
1255 . 5X,'exp_data(algorithm',I16/,
1256 . 5X,'iz flag(engine update) . . . . . . . =',I16/,
1257 . 5X,'iform flag(friedlander model) . . . .=',I16/,
1258 . 5X,'i_tshift(flag for time shifting) . .=',I16/,
1259 . 5X,'x-det. . . . . . . . . . . . . . . . .=',E12.4/,
1260 . 5X,'y-det . . . . . . . . . . . . . . . . =',E12.4/,
1261 . 5X,'z-det . . . . . . . . . . . . . . . . =',E12.4/,
1262 . 5X,'tdet . . . . . . . . . . . . . . . . .=',E12.4/,
1263 . 5X,'wtnt . . . . . . . . . . . . . . . .=',E12.4/,
1264 . 5X,'pmin . . . . . . . . . . . . . . . . .=',E12.4)
1265 2012 FORMAT(
1266 . 5X,'time of arrival(shifted). . . . . . .=',E12.4)
1267 2013 FORMAT(
1268 . 5X,'time of arrival . . . . . . . . . . . =',E12.4)
1269 2014 FORMAT(
1270 . 5X,'ndt . . . . . . . . . . . . . . . . .=',I10)
1271 2015 FORMAT(
1272 . 5X,'node identifier. . . . . . . . . . . .=',I10)
1273 2023 FORMAT(
1274 . 5X,'charge height. . . . . . . . . . . . .=',E12.4)
1275 2024 FORMAT(
1276 . 5X,'stop time . . . . .. . . . . . . . . .=',E12.4/)
1277 2031 FORMAT(
1278 . 5X,'computed minimum decay PARAMETER . . .=',E12.4)
1279 2032 FORMAT(
1280 . 5X,'computed maximum decay PARAMETER . . .=',E12.4)
1281 2033 FORMAT(
1282 . 5X,'computed minimum time step . . . . . .=',E12.4)
1283 2034 FORMAT(
1284 . 5X,'ground nx . . . . . . . . . . . . . .=',E12.4/,
1285 . 5X,'ground ny . . . . . . . . . . . . . . =',E12.4/,
1286 . 5X,'ground nz . . . . . . . . . . . . . . =',E12.4)
1287 2035 FORMAT(
1288 . 5X,'time shift VALUE',/,
1289 . 5X,' (minimum among all pblast options). =',E12.4//)
1290 2036 FORMAT(
1291 . 5X,'ishape flag . . . . . . . . . . . . =',I10)
1292
1293
1294 1000 CONTINUE
1295 IF (IERR1 /= 0) THEN
1296 WRITE(IOUT,*)' ** error in memory allocation, /load/pblast
id=
',ID,' '
1297 WRITE(ISTDO,*)' ** error in memory allocation, /load/pblast
id=
',ID,' '
1298 CALL ARRET(2)
1299 END IF
1300
1301
end diagonal values have been computed in the(sparse) matrix id.SOL
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
subroutine hm_get_intv(name, ival, is_available, lsubmodel)
subroutine hm_option_start(entity_type)
for(i8=*sizetab-1;i8 >=0;i8--)
integer, parameter nchartitle
integer, parameter ncharline
1.15.0