OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
stat_c_orth_loc.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| stat_c_orth_loc ../engine/source/output/sta/stat_c_orth_loc.F
25!||--- called by ------------------------------------------------------
26!|| genstat ../engine/source/output/sta/genstat.F
27!||--- calls -----------------------------------------------------
28!|| clsconv3 ../engine/source/output/anim/generate/dfuncc.F
29!|| spmd_rgather9_dp ../engine/source/mpi/interfaces/spmd_outp.F
30!|| spmd_stat_pgather ../engine/source/mpi/output/spmd_stat.F
31!|| strs_txt50 ../engine/source/output/sta/sta_txt.F
32!||--- uses -----------------------------------------------------
33!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
34!|| my_alloc_mod ../common_source/tools/memory/my_alloc.f90
35!||====================================================================
36 SUBROUTINE stat_c_orth_loc(ELBUF_TAB,IPARG ,IPM ,IGEO ,IXC ,
37 2 IXTG ,WA,WAP0 ,IPARTC, IPARTTG,
38 3 IPART_STATE,STAT_INDXC,STAT_INDXTG,X,IDEL,
39 4 SIZP0)
40C-----------------------------------------------
41C M o d u l e s
42C-----------------------------------------------
43 USE elbufdef_mod
44 USE my_alloc_mod
45C-----------------------------------------------
46C I m p l i c i t T y p e s
47C-----------------------------------------------
48#include "implicit_f.inc"
49C-----------------------------------------------
50C C o m m o n B l o c k s
51C-----------------------------------------------
52#include "com01_c.inc"
53#include "param_c.inc"
54#include "units_c.inc"
55#include "task_c.inc"
56#include "scr14_c.inc"
57#include "scr16_c.inc"
58C-----------------------------------------------
59C D u m m y A r g u m e n t s
60C-----------------------------------------------
61 INTEGER SIZLOC,SIZP0
62 INTEGER IXC(NIXC,*),IXTG(NIXTG,*),
63 . IPARG(NPARG,*),IPM(NPROPMI,*),IGEO(NPROPGI,*),
64 . ipartc(*), iparttg(*), ipart_state(*),
65 . stat_indxc(*), stat_indxtg(*)
66 my_real
67 . x(3,*)
68 TYPE (ELBUF_STRUCT_), DIMENSION(NGROUP), TARGET :: ELBUF_TAB
69 double precision WA(*),WAP0(*)
70C-----------------------------------------------
71C L o c a l V a r i a b l e s
72C-----------------------------------------------
73 INTEGER I,J,K,N,II,JJ,LEN, IOFF, IREP,NG, NEL, NFT, ITY, LFT, NPT,
74 . LLT, MLW,NDIR,NLAY,IHBE,ISH3N,IDRAPE,NPTT, NPLY_MAX,IPT,
75 . IGTYP, ID, IPRT0, IPRT,NPG,IPG,IE, FLAGDEG,IDEL,ILAY,ILAW,IFRAM_OLD
76 INTEGER,DIMENSION(:),ALLOCATABLE :: PTWA
77 INTEGER,DIMENSION(:),ALLOCATABLE :: PTWA_P0
78 my_real
79 . thk, em, eb, h1, h2, h3, angle1,
80 . angle2,dir1_1,dir1_2,dir2_1,dir2_2,aa,bb,v1,v2,v3,x21,x32,x34,
81 . x41,y21,y32,y34,y41,z21,z32,z34,z41,suma,s1,s2,vr,vs,x31,y31,
82 . z31,e11,e12,e13,e21,e22,e23,sum,area
83 my_real
84 . e1x, e1y, e1z, e2x,e2y, e2z, e3x, e3y, e3z, rx,ry,rz,sx,sy,sz,
85 . x2l
86 CHARACTER*100 DELIMIT,LINE
87 TYPE(G_BUFEL_) ,POINTER :: GBUF
88 TYPE(BUF_LAY_) ,POINTER :: BUFLY
89 TYPE(l_bufel_dir_) ,POINTER :: LBUF_DIR
90C-----------------------------------------------
91 DATA delimit(1:60)
92 ./'#---1----|----2----|----3----|----4----|----5----|----6----|'/
93 DATA delimit(61:100)
94 ./'----7----|----8----|----9----|----10---|'/
95C=======================================================================
96C 4-NODE SHELLS
97C-----------------------------------------------
98 CALL my_alloc(ptwa,max(stat_numelc ,stat_numeltg))
99 ALLOCATE(ptwa_p0(0:max(1,stat_numelc_g,stat_numeltg_g)))
100C-----------------------------------------------
101 jj = 0
102 flagdeg = 1
103 IF (stat_numelc==0) GOTO 200
104C
105 ie=0
106 DO ng=1,ngroup
107 ity = iparg(5,ng)
108 IF (ity == 3) THEN
109 gbuf => elbuf_tab(ng)%GBUF
110 mlw = iparg(1,ng)
111 nel = iparg(2,ng)
112 nft = iparg(3,ng)
113 npt = iparg(6,ng)
114 igtyp= iparg(38,ng)
115 ihbe = iparg(23,ng)
116 npg = elbuf_tab(ng)%NPTR*elbuf_tab(ng)%NPTS
117 irep = iparg(35,ng)
118 nlay = elbuf_tab(ng)%NLAY
119 idrape = elbuf_tab(ng)%IDRAPE
120 nply_max = 0
121 IF(idrape > 0 . and. (igtyp == 51 .OR. igtyp == 52)) THEN
122 nply_max = 0
123 DO j = 1,nlay
124 nply_max = nply_max + elbuf_tab(ng)%BUFLY(j)%NPTT
125 ENDDO
126 ENDIF
127 npt = max(nply_max, npt)
128 ndir = 1
129 IF (irep > 1) ndir = 2
130 lft=1
131 llt=nel
132c------------------------------------------
133 DO i=lft,llt
134 n = i + nft
135C
136 x21 = x(1,ixc(3,n))-x(1,ixc(2,n))
137 x32 = x(1,ixc(4,n))-x(1,ixc(3,n))
138 x34 = x(1,ixc(4,n))-x(1,ixc(5,n))
139 x41 = x(1,ixc(5,n))-x(1,ixc(2,n))
140C
141 y21 = x(2,ixc(3,n))-x(2,ixc(2,n))
142 y32 = x(2,ixc(4,n))-x(2,ixc(3,n))
143 y34 = x(2,ixc(4,n))-x(2,ixc(5,n))
144 y41 = x(2,ixc(5,n))-x(2,ixc(2,n))
145C
146 z21 = x(3,ixc(3,n))-x(3,ixc(2,n))
147 z32 = x(3,ixc(4,n))-x(3,ixc(3,n))
148 z34 = x(3,ixc(4,n))-x(3,ixc(5,n))
149 z41 = x(3,ixc(5,n))-x(3,ixc(2,n))
150C
151 e1x = (x21+x34)
152 e1y = (y21+y34)
153 e1z = (z21+z34)
154C
155 e2x = (x32+x41)
156 e2y = (y32+y41)
157 e2z = (z32+z41)
158C
159 e3x = e1y*e2z-e1z*e2y
160 e3y = e1z*e2x-e1x*e2z
161 e3z = e1x*e2y-e1y*e2x
162 IF (irep > 0) THEN
163 rx = e1x
164 ry = e1y
165 rz = e1z
166 sx = e2x
167 sy = e2y
168 sz = e2z
169 ENDIF
170 IF (ishfram == 0 .OR. igtyp == 16) THEN
171C---- Repere convecte symetrique - version 5 (default)
172 suma = e3x*e3x+e3y*e3y+e3z*e3z
173 suma = one / max(sqrt(suma),em20)
174 e3x = e3x * suma
175 e3y = e3y * suma
176 e3z = e3z * suma
177C
178 s1 = e1x*e1x+e1y*e1y+e1z*e1z
179 s2 = e2x*e2x+e2y*e2y+e2z*e2z
180 suma = sqrt(s1/s2)
181 e1x = e1x + (e2y*e3z-e2z*e3y)*suma
182 e1y = e1y + (e2z*e3x-e2x*e3z)*suma
183 e1z = e1z + (e2x*e3y-e2y*e3x)*suma
184C
185 suma = e1x*e1x+e1y*e1y+e1z*e1z
186 suma = one / max(sqrt(suma),em20)
187 e1x = e1x * suma
188 e1y = e1y * suma
189 e1z = e1z * suma
190C
191 e2x = e3y * e1z - e3z * e1y
192 e2y = e3z * e1x - e3x * e1z
193 e2z = e3x * e1y - e3y * e1x
194 ELSEIF (ishfram == 2) THEN
195C---- Repere convecte nonsymetrique - version 4
196 suma = e2x*e2x+e2y*e2y+e2z*e2z
197 e1x = e1x*suma + e2y*e3z-e2z*e3y
198 e1y = e1y*suma + e2z*e3x-e2x*e3z
199 e1z = e1z*suma + e2x*e3y-e2y*e3x
200 suma = e1x*e1x+e1y*e1y+e1z*e1z
201 suma = one/max(sqrt(suma),em20)
202 e1x = e1x*suma
203 e1y = e1y*suma
204 e1z = e1z*suma
205C
206 suma = e3x*e3x+e3y*e3y+e3z*e3z
207 suma = one / max(sqrt(suma),em20)
208 e3x = e3x * suma
209 e3y = e3y * suma
210 e3z = e3z * suma
211C
212 e2x = e3y*e1z-e3z*e1y
213 e2y = e3z*e1x-e3x*e1z
214 e2z = e3x*e1y-e3y*e1x
215 suma = e2x*e2x+e2y*e2y+e2z*e2z
216 suma = one/max(sqrt(suma),em20)
217 e2x = e2x*suma
218 e2y = e2y*suma
219 e2z = e2z*suma
220 ENDIF
221C
222 iprt=ipartc(n)
223 IF (ipart_state(iprt) == 0) cycle
224C
225 jj = jj + 1
226 IF (mlw /= 0 .AND. mlw /= 13) THEN
227 wa(jj) = gbuf%OFF(i)
228 ELSE
229 wa(jj) = zero
230 ENDIF
231 jj = jj + 1
232 wa(jj) = iprt
233 jj = jj + 1
234 wa(jj) = ixc(nixc,n)
235 jj = jj + 1
236 wa(jj) = npt
237 jj = jj + 1
238 wa(jj) = npg
239 jj = jj + 1
240 wa(jj) = ihbe
241 jj = jj + 1
242 wa(jj) = igtyp
243 jj = jj + 1
244 wa(jj) = ndir
245 jj = jj + 1
246 wa(jj) = irep
247 IF(idrape > 0 .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
248 DO j=1,nlay
249 nptt = elbuf_tab(ng)%BUFLY(j)%NPTT
250 DO ipt = 1, nptt
251 lbuf_dir => elbuf_tab(ng)%BUFLY(j)%LBUF_DIR(ipt)
252 dir1_1 = lbuf_dir%DIRA(i)
253 dir1_2 = lbuf_dir%DIRA(i+nel)
254 ilaw = elbuf_tab(ng)%BUFLY(j)%ILAW
255 IF (irep == 1) THEN
256 aa = dir1_1
257 bb = dir1_2
258 v1 = aa*rx + bb*sx
259 v2 = aa*ry + bb*sy
260 v3 = aa*rz + bb*sz
261 vr = v1*e1x+ v2*e1y + v3*e1z
262 vs = v1*e2x+ v2*e2y + v3*e2z
263 suma=max(sqrt(vr*vr + vs*vs) , em20)
264 dir1_1 = vr/suma
265 dir1_2 = vs/suma
266 ELSEIF (irep == 2) THEN
267C--- Axe I
268 aa = dir1_1
269 bb = dir1_2
270 v1 = aa*rx + bb*sx
271 v2 = aa*ry + bb*sy
272 v3 = aa*rz + bb*sz
273 vr = v1*e1x+ v2*e1y + v3*e1z
274 vs = v1*e2x+ v2*e2y + v3*e2z
275 suma=max(sqrt(vr*vr + vs*vs) , em20)
276 dir1_1 = vr/suma
277 dir1_2 = vs/suma
278C--- Axe II
279 aa = lbuf_dir%DIRB(i)
280 bb = lbuf_dir%DIRB(i + nel)
281 v1 = aa*rx + bb*sx
282 v2 = aa*ry + bb*sy
283 v3 = aa*rz + bb*sz
284 vr = v1*e1x+ v2*e1y + v3*e1z
285 vs = v1*e2x+ v2*e2y + v3*e2z
286 suma=max(sqrt(vr*vr + vs*vs) , em20)
287 dir2_1 = vr/suma
288 dir2_2 = vs/suma
289 ELSEIF (irep == 3) THEN
290C mixing law 58 with other user laws with IREP = 0 within PID51
291 IF (ilaw == 58) THEN
292C--- Axe I
293 aa = dir1_1
294 bb = dir1_2
295 v1 = aa*rx + bb*sx
296 v2 = aa*ry + bb*sy
297 v3 = aa*rz + bb*sz
298 vr = v1*e1x+ v2*e1y + v3*e1z
299 vs = v1*e2x+ v2*e2y + v3*e2z
300 suma=max(sqrt(vr*vr + vs*vs) , em20)
301 dir1_1 = vr/suma
302 dir1_2 = vs/suma
303C--- Axe II
304 aa = lbuf_dir%DIRB(i)
305 bb = lbuf_dir%DIRB(i + nel)
306 v1 = aa*rx + bb*sx
307 v2 = aa*ry + bb*sy
308 v3 = aa*rz + bb*sz
309 vr = v1*e1x+ v2*e1y + v3*e1z
310 vs = v1*e2x+ v2*e2y + v3*e2z
311 suma=max(sqrt(vr*vr + vs*vs) , em20)
312 dir2_1 = vr/suma
313 dir2_2 = vs/suma
314 ELSE ! IREP = 0 within PID51
315C DIR1_1 and DIT1_2 already got
316 ENDIF
317 ELSEIF (irep == 4) THEN
318C mixing law 58 with other user laws with IREP = 1 within PID51
319 IF (ilaw == 58) THEN
320C--- Axe I
321 aa = dir1_1
322 bb = dir1_2
323 v1 = aa*rx + bb*sx
324 v2 = aa*ry + bb*sy
325 v3 = aa*rz + bb*sz
326 vr = v1*e1x+ v2*e1y + v3*e1z
327 vs = v1*e2x+ v2*e2y + v3*e2z
328 suma=max(sqrt(vr*vr + vs*vs) , em20)
329 dir1_1 = vr/suma
330 dir1_2 = vs/suma
331C--- Axe II
332 aa = lbuf_dir%DIRB(i)
333 bb = lbuf_dir%DIRB(i + nel)
334 v1 = aa*rx + bb*sx
335 v2 = aa*ry + bb*sy
336 v3 = aa*rz + bb*sz
337 vr = v1*e1x+ v2*e1y + v3*e1z
338 vs = v1*e2x+ v2*e2y + v3*e2z
339 suma=max(sqrt(vr*vr + vs*vs) , em20)
340 dir2_1 = vr/suma
341 dir2_2 = vs/suma
342 ELSE ! IREP = 1 within PID51
343 aa = dir1_1
344 bb = dir1_2
345 v1 = aa*rx + bb*sx
346 v2 = aa*ry + bb*sy
347 v3 = aa*rz + bb*sz
348 vr = v1*e1x+ v2*e1y + v3*e1z
349 vs = v1*e2x+ v2*e2y + v3*e2z
350 suma=max(sqrt(vr*vr + vs*vs) , em20)
351 dir1_1 = vr/suma
352 dir1_2 = vs/suma
353 ENDIF
354 ENDIF
355C
356 jj = jj + 1
357 wa(jj) = ilaw
358C charging all directions
359 jj = jj + 1
360 IF (mlw /= 0 .AND. mlw /= 13) THEN
361 wa(jj) = dir1_1
362 ELSE
363 wa(jj) = zero
364 ENDIF
365 jj = jj + 1
366 IF (mlw /= 0 .AND. mlw /= 13) THEN
367 wa(jj) = dir1_2
368 ELSE
369 wa(jj) = zero
370 ENDIF
371 IF (irep > 1 .AND. ilaw == 58) THEN
372 jj = jj + 1
373 IF (mlw /= 0 .AND. mlw /= 13) THEN
374 wa(jj) = dir2_1
375 ELSE
376 wa(jj) = zero
377 ENDIF
378 jj = jj + 1
379 IF (mlw /= 0 .AND. mlw /= 13) THEN
380 wa(jj) = dir2_2
381 ELSE
382 wa(jj) = zero
383 ENDIF
384 ENDIF ! IF (IREP > 1) THEN
385 ENDDO ! NPTT
386 ENDDO ! DO J=1,NLAY
387 ELSEIF (igtyp == 9 .OR. igtyp == 10 .OR. igtyp == 11 .OR.
388 . igtyp == 16 .OR. igtyp == 17 .OR. igtyp == 51 .OR.
389 . igtyp == 52) THEN
390 DO j=1,nlay
391 dir1_1 = elbuf_tab(ng)%BUFLY(j)%DIRA(i)
392 dir1_2 = elbuf_tab(ng)%BUFLY(j)%DIRA(i + nel)
393 ilaw = elbuf_tab(ng)%BUFLY(j)%ILAW
394 IF (irep == 1) THEN
395 aa = dir1_1
396 bb = dir1_2
397 v1 = aa*rx + bb*sx
398 v2 = aa*ry + bb*sy
399 v3 = aa*rz + bb*sz
400 vr = v1*e1x+ v2*e1y + v3*e1z
401 vs = v1*e2x+ v2*e2y + v3*e2z
402 suma=max(sqrt(vr*vr + vs*vs) , em20)
403 dir1_1 = vr/suma
404 dir1_2 = vs/suma
405 ELSEIF (irep == 2) THEN
406C--- Axe I
407 aa = dir1_1
408 bb = dir1_2
409 v1 = aa*rx + bb*sx
410 v2 = aa*ry + bb*sy
411 v3 = aa*rz + bb*sz
412 vr = v1*e1x+ v2*e1y + v3*e1z
413 vs = v1*e2x+ v2*e2y + v3*e2z
414 suma=max(sqrt(vr*vr + vs*vs) , em20)
415 dir1_1 = vr/suma
416 dir1_2 = vs/suma
417C--- Axe II
418 aa = elbuf_tab(ng)%BUFLY(j)%DIRB(i)
419 bb = elbuf_tab(ng)%BUFLY(j)%DIRB(i + nel)
420 v1 = aa*rx + bb*sx
421 v2 = aa*ry + bb*sy
422 v3 = aa*rz + bb*sz
423 vr = v1*e1x+ v2*e1y + v3*e1z
424 vs = v1*e2x+ v2*e2y + v3*e2z
425 suma=max(sqrt(vr*vr + vs*vs) , em20)
426 dir2_1 = vr/suma
427 dir2_2 = vs/suma
428 ELSEIF (irep == 3) THEN
429C mixing law58 with other user laws with IREP = 0 within PID51
430 IF (ilaw == 58) THEN
431C--- Axe I
432 aa = dir1_1
433 bb = dir1_2
434 v1 = aa*rx + bb*sx
435 v2 = aa*ry + bb*sy
436 v3 = aa*rz + bb*sz
437 vr = v1*e1x+ v2*e1y + v3*e1z
438 vs = v1*e2x+ v2*e2y + v3*e2z
439 suma=max(sqrt(vr*vr + vs*vs) , em20)
440 dir1_1 = vr/suma
441 dir1_2 = vs/suma
442C--- Axe II
443 aa = elbuf_tab(ng)%BUFLY(j)%DIRB(i)
444 bb = elbuf_tab(ng)%BUFLY(j)%DIRB(i + nel)
445 v1 = aa*rx + bb*sx
446 v2 = aa*ry + bb*sy
447 v3 = aa*rz + bb*sz
448 vr = v1*e1x+ v2*e1y + v3*e1z
449 vs = v1*e2x+ v2*e2y + v3*e2z
450 suma=max(sqrt(vr*vr + vs*vs) , em20)
451 dir2_1 = vr/suma
452 dir2_2 = vs/suma
453 ELSE ! IREP = 0 within PID51
454C DIR1_1 and DIT1_2 already got
455 ENDIF
456 ELSEIF (irep == 4) THEN
457C mixing law58 with other user laws with IREP = 1 within PID51
458 IF (ilaw == 58) THEN
459C--- Axe I
460 aa = dir1_1
461 bb = dir1_2
462 v1 = aa*rx + bb*sx
463 v2 = aa*ry + bb*sy
464 v3 = aa*rz + bb*sz
465 vr = v1*e1x+ v2*e1y + v3*e1z
466 vs = v1*e2x+ v2*e2y + v3*e2z
467 suma=max(sqrt(vr*vr + vs*vs) , em20)
468 dir1_1 = vr/suma
469 dir1_2 = vs/suma
470C--- Axe II
471 aa = elbuf_tab(ng)%BUFLY(j)%DIRB(i)
472 bb = elbuf_tab(ng)%BUFLY(j)%DIRB(i + nel)
473 v1 = aa*rx + bb*sx
474 v2 = aa*ry + bb*sy
475 v3 = aa*rz + bb*sz
476 vr = v1*e1x+ v2*e1y + v3*e1z
477 vs = v1*e2x+ v2*e2y + v3*e2z
478 suma=max(sqrt(vr*vr + vs*vs) , em20)
479 dir2_1 = vr/suma
480 dir2_2 = vs/suma
481 ELSE ! IREP = 1 within PID51
482 aa = dir1_1
483 bb = dir1_2
484 v1 = aa*rx + bb*sx
485 v2 = aa*ry + bb*sy
486 v3 = aa*rz + bb*sz
487 vr = v1*e1x+ v2*e1y + v3*e1z
488 vs = v1*e2x+ v2*e2y + v3*e2z
489 suma=max(sqrt(vr*vr + vs*vs) , em20)
490 dir1_1 = vr/suma
491 dir1_2 = vs/suma
492 ENDIF
493 ENDIF
494C
495 jj = jj + 1
496 wa(jj) = ilaw
497C charging all directions
498 jj = jj + 1
499 IF (mlw /= 0 .AND. mlw /= 13) THEN
500 wa(jj) = dir1_1
501 ELSE
502 wa(jj) = zero
503 ENDIF
504 jj = jj + 1
505 IF (mlw /= 0 .AND. mlw /= 13) THEN
506 wa(jj) = dir1_2
507 ELSE
508 wa(jj) = zero
509 ENDIF
510 IF (irep > 1 .AND. ilaw == 58) THEN
511 jj = jj + 1
512 IF (mlw /= 0 .AND. mlw /= 13) THEN
513 wa(jj) = dir2_1
514 ELSE
515 wa(jj) = zero
516 ENDIF
517 jj = jj + 1
518 IF (mlw /= 0 .AND. mlw /= 13) THEN
519 wa(jj) = dir2_2
520 ELSE
521 wa(jj) = zero
522 ENDIF
523 ENDIF ! IF (IREP > 1) THEN
524 ENDDO ! DO J=1,NLAY
525 ENDIF ! IF ( IGTYP )
526c
527 ie=ie+1
528C pointeur de fin de zone dans WA
529 ptwa(ie)=jj
530c
531 ENDDO ! DO I=LFT,LLT
532 ENDIF ! IF (ITY == 3)
533 ENDDO ! DO NG=1,NGROUP
534C
535 200 CONTINUE
536c------------------------------------------------------
537 IF (nspmd == 1) THEN
538 ptwa_p0(0)=0
539 DO n=1,stat_numelc
540 ptwa_p0(n)=ptwa(n)
541 ENDDO
542 len=jj
543 DO j=1,len
544 wap0(j)=wa(j)
545 ENDDO
546 ELSE
547C construit les pointeurs dans le tableau global WAP0
548 CALL spmd_stat_pgather(ptwa,stat_numelc,ptwa_p0,stat_numelc_g)
549 len = 0
550 CALL spmd_rgather9_dp(wa,jj,wap0,sizp0,len)
551 ENDIF
552c------------------------------------------------------
553 IF (ispmd == 0.AND.len > 0) THEN
554 iprt0=0
555 DO n=1,stat_numelc_g
556
557C retrouve le nieme elt dans l'ordre d'id croissant
558 k=stat_indxc(n)
559C retrouve l'adresse dans WAP0
560 j=ptwa_p0(k-1)
561C
562 ioff = nint(wap0(j + 1))
563 IF(idel==0.OR.(idel==1.AND.ioff >=1))THEN
564 iprt = nint(wap0(j + 2))
565 IF (iprt /= iprt0) THEN
566 IF (izipstrs == 0) THEN
567 WRITE(iugeo,'(A)') delimit
568 WRITE(iugeo,'(A)')'/INISHE/ORTH_LOC'
569 WRITE(iugeo,'(A)')
570 .'#------------------------ REPEAT --------------------------'
571 WRITE(iugeo,'(A)')
572 . '# SHELLID NIP NDIR'
573 WRITE(iugeo,'(A)')
574 .'#---------------------- END REPEAT ------------------------'
575 WRITE(iugeo,'(A)') delimit
576 ELSE
577 WRITE(line,'(A)') delimit
578 CALL strs_txt50(line,100)
579 WRITE(line,'(A)')'/INISHE/ORTH_LOC'
580 CALL strs_txt50(line,100)
581 WRITE(line,'(A)')
582 .'#------------------------ REPEAT --------------------------'
583 CALL strs_txt50(line,100)
584 WRITE(line,'(A)')
585 . '# SHELLID NIP NDIR'
586 CALL strs_txt50(line,100)
587 WRITE(line,'(A)')
588 .'#---------------------- END REPEAT ------------------------'
589 CALL strs_txt50(line,100)
590 WRITE(line,'(A)') delimit
591 CALL strs_txt50(line,100)
592 ENDIF
593 iprt0=iprt
594 ENDIF ! IF (IPRT /= IPRT0)
595 id = nint(wap0(j + 3))
596 npt = nint(wap0(j + 4))
597 npg = nint(wap0(j + 5))
598 ihbe = nint(wap0(j + 6))
599 igtyp = nint(wap0(j + 7))
600 ndir = nint(wap0(j + 8))
601 irep = nint(wap0(j + 9))
602 j = j + 9
603 IF (igtyp == 9) THEN
604 npt = 1
605C skip ILAW = WAP0(J + 1) ! not used for IGTYP = 9
606 j = j + 1
607 angle1 = atan2(wap0(j + 2), wap0(j + 1))
608 IF(flagdeg == 1) angle1=angle1*hundred80/pi
609 IF (izipstrs == 0) THEN
610 WRITE(iugeo,'(5I10)')id,npt,npg,ndir,flagdeg
611 WRITE(iugeo,'(1PE20.13)')angle1
612 ELSE
613 WRITE(line,'(5I10)')id,npt,npg,ndir,flagdeg
614 CALL strs_txt50(line,100)
615 WRITE(line,'(1PE20.13)')angle1
616 CALL strs_txt50(line,20)
617 ENDIF
618 j = j + 2
619 ELSEIF (igtyp == 10 .OR. igtyp == 11 .OR. igtyp == 16.OR.
620 . igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
621 IF (izipstrs == 0) THEN
622 WRITE(iugeo,'(5I10)')id,npt,npg,ndir,flagdeg
623 ELSE
624 WRITE(line,'(5I10)')id,npt,npg,ndir,flagdeg
625 CALL strs_txt50(line,100)
626 ENDIF
627 DO i=1,npt
628C---
629 ilaw = wap0(j + 1)
630 j = j + 1
631C---
632 IF (irep == 2 .OR. (irep > 2 .AND. ilaw == 58)) THEN
633C IGTYP = 16 or (IGTYP = 51 + ILAW 58)
634 angle1 = atan2(wap0(j + 2), wap0(j + 1))
635 angle2 = atan2(wap0(j + 4), wap0(j + 3))
636 angle2 = mod(angle2 - angle1,two*pi)
637 IF (flagdeg == 1) angle1=angle1*hundred80/pi
638 IF (flagdeg == 1) angle2=angle2*hundred80/pi
639 IF (izipstrs == 0) THEN
640 WRITE(iugeo,'(1P2E20.13)')angle1,angle2
641 ELSE
642 WRITE(line,'(1P2E20.13)')angle1,angle2
643 CALL strs_txt50(line,40)
644 ENDIF
645 j = j + 4
646 ELSE
647 angle1 = atan2(wap0(j + 2), wap0(j + 1))
648 IF (flagdeg == 1) angle1=angle1*hundred80/pi
649 IF (izipstrs == 0) THEN
650 WRITE(iugeo,'(1PE20.13)')angle1
651 ELSE
652 WRITE(line,'(1PE20.13)')angle1
653 CALL strs_txt50(line,20)
654 ENDIF
655 j = j + 2
656 ENDIF ! IF (IREP == 2 .OR. (IREP > 2 .AND. ILAW == 58))
657C---------------
658 ENDDO ! DO I=1,NPT
659 ENDIF ! IF (IGTYP)
660 ENDIF
661 ENDDO ! DO N=1,STAT_NUMELC_G
662 ENDIF ! IF (ISPMD == 0.AND.LEN > 0)
663C-----------------------------------------------
664C 3-NODE SHELLS
665C-----------------------------------------------
666 jj = 0
667 IF (stat_numeltg==0) GOTO 300
668C
669 ie=0
670C
671 DO ng=1,ngroup
672 ity = iparg(5,ng)
673 IF (ity == 7) THEN
674 gbuf => elbuf_tab(ng)%GBUF
675 mlw = iparg(1,ng)
676 nel = iparg(2,ng)
677 nft = iparg(3,ng)
678 npt = iparg(6,ng)
679 ish3n= iparg(23,ng)
680 igtyp= iparg(38,ng)
681 npg = elbuf_tab(ng)%NPTR*elbuf_tab(ng)%NPTS
682 irep = iparg(35,ng)
683 nlay = elbuf_tab(ng)%NLAY
684 idrape = elbuf_tab(ng)%IDRAPE
685 nply_max = 0
686 IF(idrape > 0 . and. (igtyp == 51 .OR. igtyp == 52)) THEN
687 nply_max = 0
688 DO j = 1,nlay
689 nply_max = nply_max + elbuf_tab(ng)%BUFLY(j)%NPTT
690 ENDDO
691 ENDIF
692 npt = max(nply_max, npt)
693 ndir = 1
694 IF (ish3n==3.AND.ish3nfram==0) THEN
695 ifram_old =0
696 ELSE
697 ifram_old =1
698 END IF
699 IF (irep > 1) ndir = 2
700 lft=1
701 llt=nel
702
703c------------------------------------------
704 DO i=lft,llt
705 n = i + nft
706C
707 x21 = x(1,ixtg(3,n))-x(1,ixtg(2,n))
708 x31 = x(1,ixtg(4,n))-x(1,ixtg(2,n))
709 x32 = x(1,ixtg(4,n))-x(1,ixtg(3,n))
710C
711 y21 = x(2,ixtg(3,n))-x(2,ixtg(2,n))
712 y31 = x(2,ixtg(4,n))-x(2,ixtg(2,n))
713 y32 = x(2,ixtg(4,n))-x(2,ixtg(3,n))
714C
715 z21 = x(3,ixtg(3,n))-x(3,ixtg(2,n))
716 z31 = x(3,ixtg(4,n))-x(3,ixtg(2,n))
717 z32 = x(3,ixtg(4,n))-x(3,ixtg(3,n))
718C
719 IF (irep > 0) THEN
720cc SX = X21
721cc SY = Y21
722cc SZ = Z21
723cc RX = X31
724cc RY = Y31
725cc RZ = Z31
726 rx = x21
727 ry = y21
728 rz = z21
729 sx = x31
730 sy = y31
731 sz = z31
732 ENDIF
733 IF(ifram_old ==0 ) THEN
734 CALL clsconv3(x21,y21,z21,x31,y31,z31,
735 + e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z)
736 ELSE
737 e1x= x21
738 e1y= y21
739 e1z= z21
740 x2l = sqrt(e1x*e1x+e1y*e1y+e1z*e1z)
741 e1x=e1x/x2l
742 e1y=e1y/x2l
743 e1z=e1z/x2l
744C
745 e3x=y31*z32-z31*y32
746 e3y=z31*x32-x31*z32
747 e3z=x31*y32-y31*x32
748 sum = sqrt(e3x*e3x+e3y*e3y+e3z*e3z)
749 e3x=e3x/sum
750 e3y=e3y/sum
751 e3z=e3z/sum
752 area = half * sum
753 e2x=e3y*e1z-e3z*e1y
754 e2y=e3z*e1x-e3x*e1z
755 e2z=e3x*e1y-e3y*e1x
756 sum = sqrt(e2x*e2x+e2y*e2y+e2z*e2z)
757 e2x=e2x/sum
758 e2y=e2y/sum
759 e2z=e2z/sum
760 END IF !(ISH3NFRAM ==0 ) THEN
761C
762 iprt=iparttg(n)
763 IF (ipart_state(iprt)==0) cycle
764C
765 jj = jj + 1
766 IF (mlw /= 0 .AND. mlw /= 13) THEN
767 wa(jj) = gbuf%OFF(i)
768 ELSE
769 wa(jj) = zero
770 ENDIF
771 jj = jj + 1
772 wa(jj) = iprt
773 jj = jj + 1
774 wa(jj) = ixtg(nixtg,n)
775 jj = jj + 1
776 wa(jj) = npt
777 jj = jj + 1
778 wa(jj) = npg
779 jj = jj + 1
780 wa(jj) = ish3n
781 jj = jj + 1
782 wa(jj) = igtyp
783 jj = jj + 1
784 wa(jj) = ndir
785 jj = jj + 1
786 wa(jj) = irep
787 IF (igtyp == 9 .OR. igtyp == 10 .OR. igtyp == 11 .OR.
788 . igtyp == 16 .OR. igtyp == 17 .OR. igtyp == 51 .OR.
789 . igtyp == 52 ) THEN
790 IF(idrape > 0 .AND. (igtyp == 51 .OR. igtyp == 52)) THEN
791 DO j=1,nlay
792 nptt = elbuf_tab(ng)%BUFLY(j)%NPTT
793 DO ipt = 1, nptt
794 lbuf_dir => elbuf_tab(ng)%BUFLY(j)%LBUF_DIR(ipt)
795 dir1_1 = lbuf_dir%DIRA(i)
796 dir1_2 = lbuf_dir%DIRA(i + nel)
797 ilaw = elbuf_tab(ng)%BUFLY(j)%ILAW
798 IF (irep == 1) THEN
799 aa = dir1_1
800 bb = dir1_2
801 v1 = aa*rx + bb*sx
802 v2 = aa*ry + bb*sy
803 v3 = aa*rz + bb*sz
804 vr = v1*e1x+ v2*e1y + v3*e1z
805 vs = v1*e2x+ v2*e2y + v3*e2z
806 suma=max(sqrt(vr*vr + vs*vs) , em20)
807 dir1_1 = vr/suma
808 dir1_2 = vs/suma
809 ELSEIF (irep == 2) THEN
810C--- Axe I
811 aa = dir1_1
812 bb = dir1_2
813 v1 = aa*rx + bb*sx
814 v2 = aa*ry + bb*sy
815 v3 = aa*rz + bb*sz
816 vr = v1*e1x+ v2*e1y + v3*e1z
817 vs = v1*e2x+ v2*e2y + v3*e2z
818 suma=max(sqrt(vr*vr + vs*vs) , em20)
819 dir1_1 = vr/suma
820 dir1_2 = vs/suma
821C--- Axe II
822 aa = lbuf_dir%DIRB(i)
823 bb = lbuf_dir%DIRB(i + nel)
824 v1 = aa*rx + bb*sx
825 v2 = aa*ry + bb*sy
826 v3 = aa*rz + bb*sz
827 vr = v1*e1x+ v2*e1y + v3*e1z
828 vs = v1*e2x+ v2*e2y + v3*e2z
829 suma=max(sqrt(vr*vr + vs*vs) , em20)
830 dir2_1 = vr/suma
831 dir2_2 = vs/suma
832 ELSEIF (irep == 3) THEN
833C mixing law 58 with other user laws with IREP = 0 within PID51
834 IF (ilaw == 58) THEN
835C--- Axe I
836 aa = dir1_1
837 bb = dir1_2
838 v1 = aa*rx + bb*sx
839 v2 = aa*ry + bb*sy
840 v3 = aa*rz + bb*sz
841 vr = v1*e1x+ v2*e1y + v3*e1z
842 vs = v1*e2x+ v2*e2y + v3*e2z
843 suma=max(sqrt(vr*vr + vs*vs) , em20)
844 dir1_1 = vr/suma
845 dir1_2 = vs/suma
846C--- Ax e II
847 aa = lbuf_dir%DIRB(i)
848 bb = lbuf_dir%DIRB(i + nel)
849 v1 = aa*rx + bb*sx
850 v2 = aa*ry + bb*sy
851 v3 = aa*rz + bb*sz
852 vr = v1*e1x+ v2*e1y + v3*e1z
853 vs = v1*e2x+ v2*e2y + v3*e2z
854 suma=max(sqrt(vr*vr + vs*vs) , em20)
855 dir2_1 = vr/suma
856 dir2_2 = vs/suma
857 ELSE ! IREP = 0 within PID51
858C DIR1_1 and DIT1_2 already got
859 ENDIF
860 ELSEIF (irep == 4) THEN
861C mixing law 58 with other user laws with IREP = 1 within PID51
862 IF (ilaw == 58) THEN
863C--- Axe I
864 aa = dir1_1
865 bb = dir1_2
866 v1 = aa*rx + bb*sx
867 v2 = aa*ry + bb*sy
868 v3 = aa*rz + bb*sz
869 vr = v1*e1x+ v2*e1y + v3*e1z
870 vs = v1*e2x+ v2*e2y + v3*e2z
871 suma=max(sqrt(vr*vr + vs*vs) , em20)
872 dir1_1 = vr/suma
873 dir1_2 = vs/suma
874C--- Axe II
875 aa = lbuf_dir%DIRB(i)
876 bb = lbuf_dir%DIRB(i + nel)
877 v1 = aa*rx + bb*sx
878 v2 = aa*ry + bb*sy
879 v3 = aa*rz + bb*sz
880 vr = v1*e1x+ v2*e1y + v3*e1z
881 vs = v1*e2x+ v2*e2y + v3*e2z
882 suma=max(sqrt(vr*vr + vs*vs) , em20)
883 dir2_1 = vr/suma
884 dir2_2 = vs/suma
885 ELSE ! IREP = 1 within PID51
886 aa = dir1_1
887 bb = dir1_2
888 v1 = aa*rx + bb*sx
889 v2 = aa*ry + bb*sy
890 v3 = aa*rz + bb*sz
891 vr = v1*e1x+ v2*e1y + v3*e1z
892 vs = v1*e2x+ v2*e2y + v3*e2z
893 suma=max(sqrt(vr*vr + vs*vs) , em20)
894 dir1_1 = vr/suma
895 dir1_2 = vs/suma
896 ENDIF
897 ENDIF
898C
899 jj = jj + 1
900 wa(jj) = ilaw
901C charging all directions
902 jj = jj + 1
903 IF (mlw /= 0 .AND. mlw /= 13) THEN
904 wa(jj) = dir1_1
905 ELSE
906 wa(jj) = zero
907 ENDIF
908 jj = jj + 1
909 IF (mlw /= 0 .AND. mlw /= 13) THEN
910 wa(jj) = dir1_2
911 ELSE
912 wa(jj) = zero
913 ENDIF
914 IF (irep > 1 .AND. ilaw == 58) THEN
915 jj = jj + 1
916 IF (mlw /= 0 .AND. mlw /= 13) THEN
917 wa(jj) = dir2_1
918 ELSE
919 wa(jj) = zero
920 ENDIF
921 jj = jj + 1
922 IF (mlw /= 0 .AND. mlw /= 13) THEN
923 wa(jj) = dir2_2
924 ELSE
925 wa(jj) = zero
926 ENDIF
927 ENDIF ! IF (IREP > 1) THEN
928 ENDDO ! NPTT
929 ENDDO ! DO J=1,NLAY
930 ELSE ! Idrape == 0
931 DO j=1,nlay
932 dir1_1 = elbuf_tab(ng)%BUFLY(j)%DIRA(i)
933 dir1_2 = elbuf_tab(ng)%BUFLY(j)%DIRA(i + nel)
934 ilaw = elbuf_tab(ng)%BUFLY(j)%ILAW
935 IF (irep == 1) THEN
936 aa = dir1_1
937 bb = dir1_2
938 v1 = aa*rx + bb*sx
939 v2 = aa*ry + bb*sy
940 v3 = aa*rz + bb*sz
941 vr = v1*e1x+ v2*e1y + v3*e1z
942 vs = v1*e2x+ v2*e2y + v3*e2z
943 suma=max(sqrt(vr*vr + vs*vs) , em20)
944 dir1_1 = vr/suma
945 dir1_2 = vs/suma
946 ELSEIF (irep == 2) THEN
947C--- Axe I
948 aa = dir1_1
949 bb = dir1_2
950 v1 = aa*rx + bb*sx
951 v2 = aa*ry + bb*sy
952 v3 = aa*rz + bb*sz
953 vr = v1*e1x+ v2*e1y + v3*e1z
954 vs = v1*e2x+ v2*e2y + v3*e2z
955 suma=max(sqrt(vr*vr + vs*vs) , em20)
956 dir1_1 = vr/suma
957 dir1_2 = vs/suma
958C--- Axe II
959 aa = elbuf_tab(ng)%BUFLY(j)%DIRB(i)
960 bb = elbuf_tab(ng)%BUFLY(j)%DIRB(i + nel)
961 v1 = aa*rx + bb*sx
962 v2 = aa*ry + bb*sy
963 v3 = aa*rz + bb*sz
964 vr = v1*e1x+ v2*e1y + v3*e1z
965 vs = v1*e2x+ v2*e2y + v3*e2z
966 suma=max(sqrt(vr*vr + vs*vs) , em20)
967 dir2_1 = vr/suma
968 dir2_2 = vs/suma
969 ELSEIF (irep == 3) THEN
970C mixing law58 with other user laws with IREP = 0 within PID51
971 IF (ilaw == 58) THEN
972C--- Axe I
973 aa = dir1_1
974 bb = dir1_2
975 v1 = aa*rx + bb*sx
976 v2 = aa*ry + bb*sy
977 v3 = aa*rz + bb*sz
978 vr = v1*e1x+ v2*e1y + v3*e1z
979 vs = v1*e2x+ v2*e2y + v3*e2z
980 suma=max(sqrt(vr*vr + vs*vs) , em20)
981 dir1_1 = vr/suma
982 dir1_2 = vs/suma
983C--- Axe II
984 aa = elbuf_tab(ng)%BUFLY(j)%DIRB(i)
985 bb = elbuf_tab(ng)%BUFLY(j)%DIRB(i + nel)
986 v1 = aa*rx + bb*sx
987 v2 = aa*ry + bb*sy
988 v3 = aa*rz + bb*sz
989 vr = v1*e1x+ v2*e1y + v3*e1z
990 vs = v1*e2x+ v2*e2y + v3*e2z
991 suma=max(sqrt(vr*vr + vs*vs) , em20)
992 dir2_1 = vr/suma
993 dir2_2 = vs/suma
994 ELSE ! IREP = 0 within PID51
995C DIR1_1 and DIT1_2 already got
996 ENDIF
997 ELSEIF (irep == 4) THEN
998C mixing law58 with other user laws with IREP = 1 within PID51
999 IF (ilaw == 58) THEN
1000C--- Axe I
1001 aa = dir1_1
1002 bb = dir1_2
1003 v1 = aa*rx + bb*sx
1004 v2 = aa*ry + bb*sy
1005 v3 = aa*rz + bb*sz
1006 vr = v1*e1x+ v2*e1y + v3*e1z
1007 vs = v1*e2x+ v2*e2y + v3*e2z
1008 suma=max(sqrt(vr*vr + vs*vs) , em20)
1009 dir1_1 = vr/suma
1010 dir1_2 = vs/suma
1011C--- Axe II
1012 aa = elbuf_tab(ng)%BUFLY(j)%DIRB(i)
1013 bb = elbuf_tab(ng)%BUFLY(j)%DIRB(i + nel)
1014 v1 = aa*rx + bb*sx
1015 v2 = aa*ry + bb*sy
1016 v3 = aa*rz + bb*sz
1017 vr = v1*e1x+ v2*e1y + v3*e1z
1018 vs = v1*e2x+ v2*e2y + v3*e2z
1019 suma=max(sqrt(vr*vr + vs*vs) , em20)
1020 dir2_1 = vr/suma
1021 dir2_2 = vs/suma
1022 ELSE ! IREP = 1 within PID51
1023 aa = dir1_1
1024 bb = dir1_2
1025 v1 = aa*rx + bb*sx
1026 v2 = aa*ry + bb*sy
1027 v3 = aa*rz + bb*sz
1028 vr = v1*e1x+ v2*e1y + v3*e1z
1029 vs = v1*e2x+ v2*e2y + v3*e2z
1030 suma=max(sqrt(vr*vr + vs*vs) , em20)
1031 dir1_1 = vr/suma
1032 dir1_2 = vs/suma
1033 ENDIF
1034 ENDIF
1035C
1036 jj = jj + 1
1037 wa(jj) = ilaw
1038C charging all directions
1039 jj = jj + 1
1040 IF (mlw /= 0 .AND. mlw /= 13) THEN
1041 wa(jj) = dir1_1
1042 ELSE
1043 wa(jj) = zero
1044 ENDIF
1045 jj = jj + 1
1046 IF (mlw /= 0 .AND. mlw /= 13) THEN
1047 wa(jj) = dir1_2
1048 ELSE
1049 wa(jj) = zero
1050 ENDIF
1051 IF (irep > 1 .AND. ilaw == 58) THEN
1052 jj = jj + 1
1053 IF (mlw /= 0 .AND. mlw /= 13) THEN
1054 wa(jj) = dir2_1
1055 ELSE
1056 wa(jj) = zero
1057 ENDIF
1058 jj = jj + 1
1059 IF (mlw /= 0 .AND. mlw /= 13) THEN
1060 wa(jj) = dir2_2
1061 ELSE
1062 wa(jj) = zero
1063 ENDIF
1064 ENDIF ! IF (IREP > 1) THEN
1065 ENDDO ! DO J=1,NLAY
1066 ENDIF ! IDRAPE
1067 ENDIF ! IF ( IGTYP )
1068c
1069 ie=ie+1
1070C pointeur de fin de zone dans WA
1071 ptwa(ie)=jj
1072c
1073 ENDDO ! DO I=LFT,LLT
1074 ENDIF ! IF (ITY == 7)
1075 ENDDO ! DO NG=1,NGROUP
1076C
1077 300 CONTINUE
1078c-----------------------------------------------------
1079 IF (nspmd == 1) THEN
1080 len=jj
1081 DO j=1,len
1082 wap0(j)=wa(j)
1083 ENDDO
1084 ptwa_p0(0)=0
1085 DO n=1,stat_numeltg
1086 ptwa_p0(n)=ptwa(n)
1087 ENDDO
1088 ELSE
1089C construit les pointeurs dans le tableau global WAP0
1090 CALL spmd_stat_pgather(ptwa,stat_numeltg,ptwa_p0,stat_numeltg_g)
1091 len = 0
1092 CALL spmd_rgather9_dp(wa,jj,wap0,sizp0,len)
1093 ENDIF
1094c-----------------------------------------------------
1095 IF (ispmd == 0.AND.len > 0) THEN
1096C
1097 iprt0=0
1098 DO n=1,stat_numeltg_g
1099C retrouve le nieme elt dans l'ordre d'id croissant
1100 k=stat_indxtg(n)
1101C retrouve l'adresse dans WAP0
1102 j=ptwa_p0(k-1)
1103C
1104 ioff = nint(wap0(j + 1))
1105 IF(idel==0.OR.(idel==1.AND.ioff >=1))THEN
1106 iprt = nint(wap0(j + 2))
1107 IF (iprt /= iprt0) THEN
1108 IF (izipstrs == 0) THEN
1109 WRITE(iugeo,'(A)') delimit
1110 WRITE(iugeo,'(A)')'/INISH3/ORTH_LOC'
1111 WRITE(iugeo,'(A)')
1112 .'#------------------------ REPEAT --------------------------'
1113 WRITE(iugeo,'(A)')
1114 . '# SHELLID NIP NDIR'
1115 WRITE(iugeo,'(A)')
1116 .'#---------------------- END REPEAT ------------------------'
1117 WRITE(iugeo,'(A)') delimit
1118 ELSE
1119 WRITE(line,'(A)') delimit
1120 CALL strs_txt50(line,100)
1121 WRITE(line,'(A)')'/INISH3/ORTH_LOC'
1122 CALL strs_txt50(line,100)
1123 WRITE(line,'(A)')
1124 .'#------------------------ REPEAT --------------------------'
1125 CALL strs_txt50(line,100)
1126 WRITE(line,'(A)')
1127 . '# SHELLID NIP NDIR'
1128 CALL strs_txt50(line,100)
1129 WRITE(line,'(A)')
1130 .'#---------------------- END REPEAT ------------------------'
1131 CALL strs_txt50(line,100)
1132 WRITE(line,'(A)') delimit
1133 CALL strs_txt50(line,100)
1134 ENDIF ! IF (IZIPSTRS == 0)
1135 iprt0=iprt
1136 ENDIF ! IF (IPRT /= IPRT0)
1137 id = nint(wap0(j + 3))
1138 npt = nint(wap0(j + 4))
1139 npg = nint(wap0(j + 5))
1140 ish3n = nint(wap0(j + 6))
1141 igtyp = nint(wap0(j + 7))
1142 ndir = nint(wap0(j + 8))
1143 irep = nint(wap0(j + 9))
1144 j = j + 9
1145 IF (igtyp == 9) THEN
1146 npt = 1
1147C skip ILAW = WAP0(J + 1) ! not used for IGTYP = 9
1148 j = j + 1
1149 angle1 = atan2(wap0(j + 2), wap0(j + 1))
1150 IF (flagdeg == 1) angle1=angle1*hundred80/pi
1151 IF (izipstrs == 0) THEN
1152 WRITE(iugeo,'(5I10)')id,npt,npg,ndir,flagdeg
1153 WRITE(iugeo,'(1PE20.13)')angle1
1154 ELSE
1155 WRITE(line,'(5I10)')id,npt,npg,ndir,flagdeg
1156 CALL strs_txt50(line,100)
1157 WRITE(line,'(1PE20.13)')angle1
1158 CALL strs_txt50(line,20)
1159 ENDIF
1160 j = j + 2
1161 ELSEIF (igtyp == 10 .OR. igtyp == 11 .OR. igtyp == 16 .OR.
1162 . igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52 ) THEN
1163 IF (izipstrs == 0) THEN
1164 WRITE(iugeo,'(5I10)')id,npt,npg,ndir,flagdeg
1165 ELSE
1166 WRITE(line,'(5I10)')id,npt,npg,ndir,flagdeg
1167 CALL strs_txt50(line,100)
1168 ENDIF
1169 DO i=1,npt
1170 ilaw = wap0(j + 1)
1171 j = j + 1
1172 IF (irep == 2 .OR. (irep > 2 .AND. ilaw == 58)) THEN
1173C IGTYP = 16 or (IGTYP = 51 + ILAW 58)
1174 angle1 = atan2(wap0(j + 2), wap0(j + 1))
1175 angle2 = atan2(wap0(j + 4), wap0(j + 3))
1176 angle2 = mod(angle2 - angle1,two*pi)
1177 IF (flagdeg == 1) angle1=angle1*hundred80/pi
1178 IF (flagdeg == 1) angle2=angle2*hundred80/pi
1179 IF (izipstrs == 0) THEN
1180 WRITE(iugeo,'(1P2E20.13)')angle1,angle2
1181 ELSE
1182 WRITE(line,'(1P2E20.13)')angle1,angle2
1183 CALL strs_txt50(line,40)
1184 ENDIF
1185 j = j + 4
1186 ELSE
1187 angle1 = atan2(wap0(j + 2), wap0(j + 1))
1188 IF (flagdeg == 1) angle1=angle1*hundred80/pi
1189 IF (izipstrs == 0) THEN
1190 WRITE(iugeo,'(1PE20.13)')angle1
1191 ELSE
1192 WRITE(line,'(1PE20.13)')angle1
1193 CALL strs_txt50(line,20)
1194 ENDIF
1195 j = j + 2
1196 ENDIF
1197 ENDDO ! DO I=1,NPT
1198 ENDIF ! IF (igtyp == 9)
1199 ENDIF
1200 ENDDO ! DO N=1,STAT_NUMELTG_G
1201 ENDIF ! IF (ISPMD == 0.AND.LEN > 0)
1202c----------
1203 DEALLOCATE(ptwa)
1204 DEALLOCATE(ptwa_p0)
1205c-----------
1206 RETURN
1207 END
my_real
#define my_real
Definition cppsort.cpp:32
clsconv3
subroutine clsconv3(rx, ry, rz, sx, sy, sz, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z)
Definition dfuncc.F:5151
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
max
#define max(a, b)
Definition macros.h:21
spmd_rgather9_dp
subroutine spmd_rgather9_dp(v, len, vp0, lenp0, iad)
Definition spmd_outp.F:1015
spmd_stat_pgather
subroutine spmd_stat_pgather(ptv, ptlen, ptv_p0, ptlen_p0)
Definition spmd_stat.F:53
strs_txt50
subroutine strs_txt50(text, length)
Definition sta_txt.F:87
stat_c_orth_loc
subroutine stat_c_orth_loc(elbuf_tab, iparg, ipm, igeo, ixc, ixtg, wa, wap0, ipartc, iparttg, ipart_state, stat_indxc, stat_indxtg, x, idel, sizp0)
Definition stat_c_orth_loc.F:40