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