66 SUBROUTINE cinit3(ELBUF_STR,IXC ,PM ,X ,GEO ,
67 2 XMAS ,IN ,NVC ,DTELEM ,IGRSH4N ,
68 3 XREFC ,NEL ,ITHK ,IHBE ,IGRSH3N ,
69 4 THK ,ISIGSH ,SIGSH ,STIFN ,STIFR ,
70 5 PARTSAV ,V ,IPART ,MSC ,INC ,
71 8 SKEW ,IPARG ,I8MI ,NSIGSH ,IGEO ,
72 9 IUSER ,ETNOD ,NSHNOD ,STC ,PTSHEL ,
73 A IPM ,BUFMAT ,SH4TREE ,MCP ,MCPS ,
74 B TEMP ,CPT_ELTENS,PART_AREA,ITAGN,ITAGE ,
75 C IXFEM ,NPF ,TF ,XFEM_STR,ISUBSTACK,
76 D STACK ,RNOISE ,DRAPE ,SH4ANG ,IDDLEVEL,
77 E GEO_STACK,IGEO_STACK,STRC ,PERTURB ,IYLDINI ,
78 F ELE_AREA ,NG ,GROUP_PARAM ,NLOC_DMG,
79 G IDRAPE ,DRAPEG ,MAT_PARAM,FAIL_FRACTAL,FAIL_BROKMANN,
93 USE random_walk_def_mod
94 USE fractal_dmg_init_mod
96 use brokmann_random_def_mod
99 use element_mod ,
only : nixc
104#include "implicit_f.inc"
108#include "mvsiz_p.inc"
112#include "vect01_c.inc"
113#include "param_c.inc"
114#include "com01_c.inc"
115#include "com04_c.inc"
116#include "scr03_c.inc"
117#include "scr17_c.inc"
119#include "com_xfem1.inc"
123 INTEGER NVC,NEL,ITHK,IHBE,ISIGSH,IXFEM,,IUSER,IYLDINI
124 INTEGER IXC(NIXC,*),IPART(*),IPARG(*),IGEO(NPROPGI,*), NSHNOD(*),
125 . PTSHEL(*),IPM(NPROPMI,*), SH4TREE(*),ITAGN(*),ITAGE(*),NPF(*),
126 . ISUBSTACK,IGEO_STACK(*),PERTURB(NPERTURB),NG,IDRAPE
128 INTEGER ,
INTENT(IN) :: IDDLEVEL
130 . PM(NPROPM,*), X(3,*), GEO(NPROPG,*), XMAS(*), IN(*),
131 . DTELEM(*), XREFC(4,3,*),THK(*), SIGSH(NSIGSH,*),
132 . STIFN(*),STIFR(*),PARTSAV(20,*), V(*) ,MSC(*) ,INC(*),
133 . SKEW(LSKEW,*), ETNOD(*), STC(*),BUFMAT(*),MCP(*),MCPS(*),
134 . TEMP(*),PART_AREA(*),TF(*),RNOISE(*),
135 . SH4ANG(*),GEO_STACK(*),STRC(*),ELE_AREA(*)
136 TYPE(ELBUF_STRUCT_),
TARGET :: ELBUF_STR
137 TYPE(elbuf_struct_),
TARGET ,
DIMENSION(NGROUP,*):: XFEM_STR
139 TYPE (STACK_PLY) :: STACK
140 TYPE (GROUP_PARAM_) :: GROUP_PARAM
141 TYPE (NLOCAL_STR_) :: NLOC_DMG
142 TYPE (DRAPE_) :: DRAPE(NUMELC_DRAPE + NUMELTG_DRAPE)
143 TYPE (DRAPEG_) :: DRAPEG
144 TYPE (MATPARAM_STRUCT_) ,
DIMENSION(NUMMAT) ,
INTENT(INOUT) :: MAT_PARAM
145 TYPE () ,
INTENT(IN) :: FAIL_FRACTAL
146 TYPE (FAIL_BROKMANN_) ,
INTENT(IN) :: FAIL_BROKMANN
147 TYPE (glob_therm_) ,
intent(in) :: glob_therm
149 TYPE (GROUP_) ,
DIMENSION(NGRSHEL) :: IGRSH4N
150 TYPE (GROUP_) ,
DIMENSION(NGRSH3N) :: IGRSH3N
154 CHARACTER(LEN=NCHARTITLE) :: TITR,TITR1
155 INTEGER I,J,N,IP,NDEPAR,IGTYP,NUVAR,ID,NLAY,II,IREP,IPROP,NUPARAM,
156 . il,ir,is,it,cpt_eltens,iun,nptr,npts,nptt,ixel,ilaw,imat,ifail,
157 . igmat,jj(9),npt_all,mpt,laynpt_max,lay_max
158 INTEGER,
DIMENSION(MVSIZ) :: IX1,IX2,IX3,IX4,IORTHLOC,MAT,PID,NGL
160 my_real,
DIMENSION(MVSIZ) :: px1g,px2g,py1g,py2g,
area,aldt,
162 . x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,
163 . e1x,e2x,e3x,e1y,e2y,e3y,e1z,e2z,e3z,
164 . x2s,y2s,x3s,y3s,x4s,y4s,
165 . x2l,x3l,x4l,y2l,y3l,y4l
166 my_real,
DIMENSION(NEL) :: tempel
167 my_real,
DIMENSION(:) ,
POINTER :: uvar,dir1,dir2
168 my_real,
ALLOCATABLE,
DIMENSION(:) :: dir_a,dir_b
169 my_real,
DIMENSION(:),
ALLOCATABLE :: phi1,phi2,coor1,coor2,coor3,coor4
170 parameter(laynpt_max = 10)
171 parameter(lay_max = 100)
172 INTEGER,
DIMENSION(:),
ALLOCATABLE :: MATLY
173 my_real,
DIMENSION(:,:),
ALLOCATABLE :: posly
175 TYPE(g_bufel_) ,
POINTER :: GBUF
176 TYPE(L_BUFEL_) ,
POINTER :: LBUF
178 CALL MY_ALLOC(MATLY,MVSIZ*LAY_MAX)
179 CALL (POSLY,MVSIZ,LAY_MAX*LAYNPT_MAX)
180 gbuf => elbuf_str%GBUF
183 iprop = ixc(nixc-1,1+nft)
184 igtyp = igeo(11,iprop)
186 igmat = igeo(98,ixc(6,1+nft))
190 CALL fretitl2(titr,igeo(npropgi-ltitr+1,iprop),ltitr)
200 nlay = elbuf_str%NLAY
201 nxel = elbuf_str%NXEL
202 nptt = elbuf_str%NPTT
203 idrape = elbuf_str%IDRAPE
206 npt_all = npt_all + elbuf_str%BUFLY(il)%NPTT
209 IF(npt_all == 0 ) npt_all = nlay
211 IF((igtyp == 51 .OR. igtyp == 52) .AND. idrape > 0)
THEN
212 ALLOCATE(phi1(mvsiz*npt_all))
213 ALLOCATE(phi2(nvsiz*npt_all))
214 ALLOCATE(dir_a(npt_all*nel*2))
215 ALLOCATE(dir_b(npt_all*nel*2))
220 ALLOCATE(coor1(npt_all*mvsiz))
221 ALLOCATE(coor2(npt_all*mvsiz))
222 ALLOCATE(coor3(npt_all*mvsiz))
223 ALLOCATE(coor4(npt_all*mvsiz))
229 ALLOCATE(phi1(nlay*mvsiz))
230 ALLOCATE(phi2(nlay*mvsiz))
231 ALLOCATE(dir_a(nlay*nel*2))
232 ALLOCATE(dir_b(nlay*nel*2))
237 ALLOCATE(coor1(nlay*mvsiz))
238 ALLOCATE(coor2(nlay*mvsiz))
239 ALLOCATE(coor3(nlay*mvsiz))
240 ALLOCATE(coor4(nlay*mvsiz))
248 IF (iparg(6) == 0.OR.npt==0) mpt=0
271 CALL ccoori(x,xrefc(1,1,nft+1),ixc(1,nft+1),
273 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
274 . ix1 ,ix2 ,ix3 ,ix4 ,ngl )
276 CALL cveok3(nvc,4,ix1,ix2,ix3,ix4)
279 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
280 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
281 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z )
286 IF ((imasadd > 0).OR.(nloc_dmg%IMOD > 0))
THEN
290 ele_area(i+nft) =
area(i)
291 IF (gbuf%G_AREA > 0) gbuf%AREA(i) =
area(i)
297 IF (jthe == 0 .and. glob_therm%NINTEMP > 0)
THEN
298 CALL initemp_shell(elbuf_str,temp,nel,numnod,numelc,4,nixc,ixc)
301 CALL cinmas(x ,xrefc(1,1,nft+1),ixc ,geo ,pm,
302 . xmas ,in ,thk ,ihbe ,partsav,
303 . v ,ipart(nft+1) ,msc(nft+1),inc(nft+1) ,
area ,
304 . i8mi ,igeo ,etnod ,imat ,iprop ,
305 . nshnod ,stc(nft+1) ,sh4tree ,mcp ,mcps(nft+1),
306 . temp ,bid ,bid ,bid ,bid,
307 . bid ,bid ,isubstack ,nlay ,elbuf_str,
308 . stack ,gbuf%THK_I ,rnoise ,drape ,glob_therm%NINTEMP,
309 . perturb,ix1 ,ix2 ,ix3 ,ix4 ,
310 . idrape ,drapeg%INDX)
312 CALL cderii(px1g,px2g,py1g,py2g,
313 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
314 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
315 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z ,
316 . x2l ,x3l ,x4l ,y2l ,y3l ,y4l )
317 CALL cdleni(pm ,geo ,stifn ,stifr ,ixc(1,nft+1),
318 . px1g ,px2g ,py1g ,py2g ,thk ,
319 . igeo ,dt ,sh4tree ,aldt ,bufmat ,
320 . ipm ,nlay ,stack%PM,isubstack,strc(nft+1),
321 .
area ,imat ,iprop ,
322 . x2l ,x3l ,x4l ,y2l ,y3l ,y4l ,
323 . stack%IGEO,group_param)
324 CALL c1buf3(geo,gbuf%THK,gbuf%OFF,thk,ksh4tree
329 xfem_str(ng,ixel)%GBUF%THK(i) = thk(i)
330 xfem_str(ng,ixel)%GBUF%OFF(i) = -one
337 . lft ,llt ,nft ,nlay ,numelc ,
338 . nsigsh ,nixc ,ixc(1,nft+1),igeo ,geo ,
339 . skew ,sigsh ,ptshel ,phi1 ,phi2 ,
340 . vx ,vy ,vz ,coor1 ,coor2 ,
341 . coor3 ,coor4 ,iorthloc ,isubstack ,stack ,
342 . irep ,elbuf_str ,drape ,sh4ang(nft+1),x ,
343 . geo_stack ,e3x ,e3y ,e3z ,
344 . gbuf%BETAORTH,x1 ,x2 ,y1 ,y2 ,
345 . z1 ,z2 ,nel ,gbuf%G_ADD_NODE,gbuf%ADD_NODE,
346 . npt_all ,idrape ,drapeg%INDX)
349 IF(igtyp == 51 .OR. igtyp == 52 .AND. igmat > 0)
THEN
352 . igeo ,geo ,vx ,vy ,vz ,
353 . phi1 ,phi2 ,coor1 ,coor2 ,coor3 ,
354 . coor4 ,iorthloc ,nlay ,irep ,isubstack,
355 . stack ,geo_stack ,igeo_stack ,ir ,is ,
357 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
358 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
359 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z ,
362 ELSEIF (mtn == 27)
THEN
364 . geo ,igeo ,pm ,ipm ,ixc(1,1+nft) ,nixc,
365 . nlay,ir ,is ,imat )
366 ELSEIF (mtn == 35)
THEN
367 nptr = elbuf_str%NPTR
368 npts = elbuf_str%NPTS
369 nptt = elbuf_str%NPTT
371 . nptr,npts,nptt,igtyp)
372 ELSEIF (mtn==15 .or. mtn==19 .or. mtn==25 .or. mtn>=28)
THEN
373 IF (mtn == 19 .AND. igtyp /= 9)
THEN
377 . i1=igeo(1,ixc(nixc-1,nft+1)))
381 . igeo ,geo ,vx ,vy ,vz ,
382 . phi1 ,phi2 ,coor1 ,coor2 ,coor3 ,
383 . coor4 ,iorthloc ,nlay ,irep ,isubstack,
384 . stack ,geo_stack ,igeo_stack ,ir ,is ,
386 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
387 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
388 . e1x, e2x, e3x, e1y, e2y, e3y ,e1z, e2z, e3z ,
392 IF ((mtn == 58 .or. mtn == 158) .AND.
393 . igtyp /= 16 .AND. igtyp /= 51 .AND. igtyp /= 52)
THEN
396 . anmode=aninfo_blind_1,
401 ELSEIF (mtn == 58 .or. mtn == 158 .OR. igtyp == 51 .OR. igtyp == 52)
THEN
403 IF (idrape == 0)
THEN
405 nptt = elbuf_str%BUFLY(il)%NPTT
406 imat = elbuf_str%BUFLY(il)%IMAT
407 ilaw = elbuf_str%BUFLY(il)%ILAW
408 nuvar = elbuf_str%BUFLY(il)%NVAR_MAT
409 dir1 => elbuf_str%BUFLY(il)%DIRA
410 dir2 => elbuf_str%BUFLY(il)%DIRB
411 nuparam = mat_param(imat)%NUPARAM
415 lbuf => elbuf_str%BUFLY(il)%LBUF(ir,is,it)
416 uvar => elbuf_str%BUFLY(il)%MAT(ir,is,it)%VAR
418 . irep ,dir1 ,dir2 ,mat_param(imat)%UPARAM,
419 . uvar ,aldt ,nel ,nuvar ,lbuf%ANG ,
420 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
421 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
422 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z )
424 ELSE IF (ilaw == 158)
THEN
426 lbuf => elbuf_str%BUFLY(il)%LBUF(ir,is,it)
427 uvar => elbuf_str%BUFLY(il)%MAT(ir,is,it)%VAR
429 . uvar ,aldt ,nel ,nuvar ,lbuf%ANG ,
430 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
431 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
432 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z )
438 nptt = elbuf_str%BUFLY(il)%NPTT
439 imat = elbuf_str%BUFLY(il)%IMAT
440 ilaw = elbuf_str%BUFLY(il)%ILAW
441 nuvar = elbuf_str%BUFLY(il)%NVAR_MAT
442 nuparam = mat_param(imat)%NUPARAM
446 lbuf => elbuf_str%BUFLY(il)%LBUF(ir,is,it)
447 uvar => elbuf_str%BUFLY(il)%MAT(ir,is,it)%VAR
448 dir1 => elbuf_str%BUFLY(il)%LBUF_DIR(it)%DIRA
449 dir2 => elbuf_str%BUFLY(il)%LBUF_DIR(it)%DIRB
451 . irep ,dir1 ,dir2 ,mat_param(imat)%UPARAM,
452 . uvar ,aldt ,nel ,nuvar ,lbuf%ANG ,
453 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
454 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
455 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z )
457 ELSE IF (ilaw == 158)
THEN
459 lbuf => elbuf_str%BUFLY(il)%LBUF(ir,is,it)
460 uvar => elbuf_str%BUFLY(il)%MAT(ir,is,it)%VAR
461 dir1 => elbuf_str%BUFLY(il)%LBUF_DIR(it)%DIRA
462 dir2 => elbuf_str%BUFLY(il)%LBUF_DIR(it)%DIRB
464 . uvar ,aldt ,nel ,nuvar ,lbuf%ANG ,
465 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
466 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
467 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z )
473 IF (mtn == 42 .OR. mtn == 69 .OR. igtyp == 51 .OR. igtyp == 52)
THEN
475 nptt = elbuf_str%BUFLY(il)%NPTT
476 ilaw = elbuf_str%BUFLY(il)%ILAW
477 nuvar = elbuf_str%BUFLY(il)%NVAR_MAT
479 IF (ilaw == 42 .OR. mtn == 69)
THEN
481 lbuf => elbuf_str%BUFLY(il)%LBUF(ir,is,it)
482 uvar => elbuf_str%BUFLY(il)%MAT(ir,is,it)%VAR
483 CALL law42c_ini(nuvar,uvar,llt)
491 IF (isigsh/=0 .OR. ithkshel == 2)
THEN
495 . elbuf_str ,lft ,llt ,geo ,igeo ,
496 . mat ,pid ,matly ,posly ,igtyp ,
497 . nlay ,mpt ,isubstack ,stack
498 . nft ,gbuf%THK ,nel ,idrape ,
scdrape ,
500 CALL corth3(elbuf_str,dir_a ,dir_b ,lft ,llt ,
502 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
503 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
504 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,e1z ,e2z ,e3z ,
509 1 lft ,llt ,nft ,mpt ,istrain ,
510 2 gbuf%THK,gbuf%EINT,gbuf%STRA,gbuf%HOURG,gbuf%PLA ,
511 3 gbuf%FOR,gbuf%MOM ,sigsh ,nlay ,gbuf%G_HOURG,
512 4 numelc ,ixc ,nixc ,nsigsh ,numshel ,
513 5 ptshel ,igeo ,thk ,nel ,e1x ,
514 6 e2x ,e3x ,e1y ,e2y ,e3y ,
515 7 e1z ,e2z ,e3z ,isigsh ,dir_a ,
516 8 dir_b ,posly ,igtyp )
519 gbuf%FOR_G(i+jj(1:5))=gbuf%FOR(i+jj(1:5))
522 ELSEIF ( ithkshel == 1 )
THEN
523 CALL thickini(lft ,llt ,nft ,ptshel,numelc,
524 2 gbuf%THK,thk ,ixc ,nixc ,nsigsh,
530 IF (istrain == 1 .AND. nxref > 0)
THEN
531 uvar => elbuf_str%BUFLY(1)%MAT(1,1,1)%VAR
532 imat = elbuf_str%BUFLY(1)%IMAT
533 CALL cepsini(elbuf_str ,mat_param(imat),
534 . lft ,llt ,ismstr ,mtn ,ithk ,
535 . pm ,geo ,ixc(1,nft+1),x ,xrefc(1,1,nft+1),
536 . gbuf%FOR ,gbuf%THK ,gbuf%EINT ,gbuf%STRA ,nlay ,
537 . px1g ,px2g ,py1g ,py2g ,x2s ,
538 . y2s ,x3s ,y3s ,x4s ,y4s ,
539 . uvar ,ipm ,igeo ,imat ,
540 . skew ,nel ,dir_a ,dir_b ,gbuf%SIGI,
545 IF (ismstr == 1 .AND. mtn==19) iparg(9)=11
547 ELSEIF (ismstr == 11 .OR.(ismstr==1 .AND. mtn==19))
THEN
550 . x2s ,y2s ,x3s ,y3s ,x4s ,y4s )
553 IF (ismstr == 10 )
THEN
556 elbuf_str%GBUF%SMSTR(jj(1)+i) = x(1,ixc(3,ii))-x(1,ixc(2,ii))
557 elbuf_str%GBUF%SMSTR(jj(2)+i) = x(2,ixc(3,ii))-x(2,ixc(2,ii))
558 elbuf_str%GBUF%SMSTR(jj(3)+i) = x(3,ixc(3,ii))-x(3,ixc(2,ii))
559 elbuf_str%GBUF%SMSTR(jj(4)+i) = x(1,ixc(4,ii))-x(1,ixc(2,ii))
560 elbuf_str%GBUF%SMSTR(jj(5)+i) = x(2,ixc(4,ii
561 elbuf_str%GBUF%SMSTR(jj(6)+i
562 elbuf_str%GBUF%SMSTR(jj(7)+i) = x(1,ixc(5,ii))-x(1,ixc(2,ii))
563 elbuf_str%GBUF%SMSTR(jj(8)+i) = x(2,ixc(5,ii))-x(2,ixc(2,ii))
564 elbuf_str%GBUF%SMSTR(jj(9)+i) = x(3,ixc(5,ii))-x(3,ixc(2,ii))
566 ELSEIF (ismstr == 11 .OR.(ismstr==1 .AND. mtn==19))
THEN
568 elbuf_str%GBUF%SMSTR(jj(1)+i) = x2s(i)
569 elbuf_str%GBUF%SMSTR(jj(2)+i) = y2s(i)
570 elbuf_str%GBUF%SMSTR(jj(3)+i) = x3s(i)
571 elbuf_str%GBUF%SMSTR(jj(4)+i) = y3s(i)
572 elbuf_str%GBUF%SMSTR(jj(5)+i) = x4s(i)
573 elbuf_str%GBUF%SMSTR(jj(6)+i) = y4s(i)
577 IF (iuser == 1 .AND. mtn > 28)
THEN
580 1 lft ,llt ,nft ,nel ,npt ,
581 2 istrain,sigsh ,numelc ,ixc ,nixc ,
582 3 nsigsh ,numshel,ptshel ,iun ,iun ,
586 IF (iyldini == 1 .AND. (mtn== 36.OR. mtn==87))
THEN
588 1 lft ,llt ,nft ,nel ,npt ,
589 2 istrain,sigsh ,numelc ,ixc ,nixc ,
590 3 nsigsh ,numshel,ptshel ,iun ,iun ,
598 IF (fail_fractal%NFAIL > 0)
THEN
599 CALL fractal_dmg_init(elbuf_str,mat_param,fail_fractal,
600 . nummat ,numelc ,nel ,nft ,ngl ,ity )
603 IF (ifail > 0 .and. iddlevel == 1)
THEN
605 . nel ,nft ,ity ,igrsh4n ,igrsh3n ,
610 . nptt ,nlay ,sigsh ,nsigsh ,ptshel ,
611 . rnoise ,perturb ,aldt ,thk )
616 IF (igtyp /= 0 .AND. igtyp /= 1 .AND.
617 . igtyp /= 9 .AND. igtyp /= 10 .AND.
618 . igtyp /= 11 .AND. igtyp /= 16 .AND.
619 . igtyp /= 17 .AND. igtyp /= 51 .AND.
630 dtelem(ndepar+i)=dt(i)
634 CALL cbufxfe(elbuf_str,xfem_str,isubstack,stack ,
635 . igeo ,geo ,lft ,llt ,mat,
636 . pid ,npt ,nptt ,nlay,ir ,
642 IF (gbuf%G_VOL > 0) gbuf%VOL(i) =
area(i)*gbuf%THK(i)
647 IF (xfem_str(ng,ixel)%GBUF%G_VOL > 0)
648 . xfem_str(ng,ixel)%GBUF%VOL(i) =
area(i)*gbuf%THK(i)
653 IF (
ALLOCATED(dir_b))
DEALLOCATE(dir_b)
654 IF (
ALLOCATED(dir_a))
DEALLOCATE(dir_a)
655 DEALLOCATE(phi1,phi2,coor1,coor2,coor3,coor4)