48
49
50
51
52 USE elbufdef_mod
57
58
59
60#include "implicit_f.inc"
61
62
63
64#include "mvsiz_p.inc"
65
66
67
68#include "param_c.inc"
69#include "com01_c.inc"
70#include "scr17_c.inc"
71#include "drape_c.inc"
72
73
74
75 INTEGER JFT,JLT,NFT,NLAY,IPT,ID,,,NSIGSH,
76 . ISUBSTACK,IREP,NPT_ALL,IDRAPE
77 INTEGER IX(NIX,*),IGEO(NPROPGI,*),PTSH(*),IORTHLOC(*)
78 INTEGER, INTENT(IN) :: NEL,G_ADD_NODE,ADD_NODE(G_ADD_NODE*NEL)
79 INTEGER, DIMENSION(*) :: INDX
81 . geo(npropg,*),skew(lskew,*),sigsh(nsigsh,*),vx(*),vy(*),vz(*),
82 . phi1(npt_all,*),phi2(npt_all,*),coor1(npt_all,mvsiz),coor2(npt_all,mvsiz),
83 . coor3(npt_all,mvsiz),coor4(npt_all,mvsiz),
84 . angle(*),geo_stack(npropg,*),x(3,*),betaorth(*)
85 my_real,
DIMENSION(MVSIZ)INTENT(IN) :: e3x,e3y,e3z,x1,x2,y1,y2,z1,z2
86
87 TYPE(ELBUF_STRUCT_), TARGET :: ELBUF_STR
88 TYPE (STACK_PLY):: STACK
89 TYPE (DRAPE_) , DIMENSION(*), TARGET :: DRAPE
90
91
92
93 INTEGER I,II,,JJ,N,NPT,NPTI,IGTYP,IPID,PID,ISK,IPANG,IPPHI,
94 . IIGEO,IADR,IPTHK,IPPOS,IPDIR,IMAT_LY,ILAW_LY,IPPID,IPMAT,ILAY,
95 . DEF_ORTH(MVSIZ),N1,N2,IRP,POS,NOD,IL,IT,NSLICE,IPT_ALL,NPTT,
96 . IE, IP,IPID_PLY,N3,N4
97 my_real v(mvsiz),e11,e12,e13,ne1,vx0,vy0,vz0,
98 . xc(mvsiz),yc(mvsiz),zc(mvsiz)
99 CHARACTER(LEN=NCHARTITLE)::TITR1
100
101 TYPE (DRAPE_PLY_) , POINTER :: DRAPE_PLY
102
103 pid = ix(nix-1,1)
104 igtyp = igeo(11,pid)
105 irp = igeo(14,pid)
106 def_orth(1:mvsiz) = nlay
107 ipdir = 0
108
109 IF (igtyp == 1) THEN
110
111 RETURN
112 ELSE
113 ipang = 200
114 ipphi = 500
115 IF (igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52 ) THEN
116 ipang = 1
117 ipthk = ipang + nlay
118 ippos = ipthk + nlay
119 ipdir = ippos + nlay
120 ENDIF
121 SELECT CASE (irp)
122 CASE (0)
123 isk = igeo(2,pid)
124 DO i=jft,jlt
125 IF (isk == 0) THEN
126 vx(i) = geo(7,pid)
127 vy(i) = geo(8,pid)
128 vz(i) = geo(9,pid)
129 ELSE
130 vx(i) = skew(1,isk)
131 vy(i) = skew(2,isk)
132 vz(i) = skew(3,isk)
133 ENDIF
134 ENDDO
135 CASE (20)
136 DO i=jft,jlt
137 n1=ix(2,i)
138 n2=ix(3,i)
139 vx(i) = x(1,n2)-x(1,n1)
140 vy(i) = x(2,n2)-x(2,n1)
141 vz(i) = x(3,n2)-x(3,n1)
142 ENDDO
143 CASE (22)
144 isk = igeo(2,pid)
145 DO i=jft,jlt
146 vx(i) = skew(1,isk)
147 vy(i) = skew(2,isk)
148 vz(i) = skew(3,isk)
149 ENDDO
150 CASE (23)
151 vx0 = geo(7,pid)
152 vy0 = geo(8,pid)
153 vz0 = geo(9,pid)
154 DO i=jft,jlt
155 vx(i) = e3y(i)*vz0 - e3z(i)*vy0
156 vy(i) = e3z(i)*vx0 - e3x(i)*vz0
157 vz(i) = e3x(i)*vy0 - e3y(i)*vx0
158 ENDDO
159 CASE (24)
160
161 DO i=jft,jlt
162 n1=ix(2,i)
164 vx(i) = x(1,nod)-x(1,n1)
165 vy(i) = x(2,nod)-x(2,n1)
166 vz(i) = x(3,nod)-x(3,n1)
167 ENDDO
168 CASE (25)
169
170 isk = igeo(2,pid)
171 IF (nix > nixtg) THEN
172 DO i=jft,jlt
173 n1=ix(2,i)
174 n2=ix(3,i)
175 n3=ix(4,i)
176 n4=ix(5,i)
177 xc(i) = fourth*(x(1,n1)+x(1,n2)+x(1,n3)+x(1,n4))
178 yc(i) = fourth*(x
179 zc(i) = fourth*(x(3,n1)+x(3,n2)+x(3,n3)+x(3,n4))
180 ENDDO
181 ELSE
182 DO i=jft,jlt
183 n1=ix(2,i)
184 n2=ix(3,i)
185 n3=ix(4,i)
186 xc(i) = third*(x(1,n1)+x(1,n2)+x(1,n3))
187 yc(i) = third*(x(2,n1)+x(2,n2)+x(2,n3))
188 zc(i) = third*(x(3,n1)+x(3,n2)+x(3,n3))
189 ENDDO
190 END IF
191 DO i=jft,jlt
192 e11 = xc(i)-skew(10,isk)
193 e12 = yc(i)-skew(11,isk)
194 e13 = zc(i)-skew(12,isk)
195 vx(i) = skew(8,isk)*e13 - skew(9,isk)*e12
196 vy(i) = skew(9,isk)*e11 - skew(7,isk)*e13
197 vz(i) = skew(7,isk)*e12 - skew(8,isk)*e11
198 ENDDO
199 END SELECT
200
201 IF (igtyp == 9) THEN
202 DO i=jft,jlt
203 phi1(1,i)= geo(10,pid)
204 ENDDO
205 ELSEIF (igtyp == 10) THEN
206 DO i=jft,jlt
207 DO j=1,nlay
208 phi1(j,i)= geo(ipang+j,pid)
209 ENDDO
210 ENDDO
211 ELSEIF (igtyp == 11) THEN
212 DO i=jft,jlt
213 DO j=1,nlay
214 phi1(j,i)= geo(ipang+j,pid)
215 ENDDO
216 ENDDO
217 ELSEIF (igtyp == 17 .AND. irp /= 24) THEN
218 IF(idrape > 0) THEN
219 DO i=jft,jlt
220 ipang = 1
221 ie = indx(nft + i)
222 IF(ie == 0) THEN
223 DO j=1,nlay
224 ipid_ply = stack%IGEO(2 + j,isubstack)
225 IF(ipid_ply > 0) THEN
226 phi1(j,i) = angle(i) + geo(2,ipid_ply) + stack%GEO(ipang+j,isubstack)
227 IF (irep >= 2) phi2(j,i)= stack%GEO(ipdir+j,isubstack)
228 def_orth(i) = def_orth(i) - 1
229 ENDIF
230 ENDDO
231 ELSE
232 DO j=1,nlay
233 ipid_ply = stack%IGEO(2+j,isubstack)
234 IF(ipid_ply > 0) THEN
235 phi1(j,i) = angle(i) + geo(2,ipid_ply) + stack%GEO(ipang+j,isubstack)
236 IF (irep >= 2) phi2(j,i)= stack%GEO(ipdir+j,isubstack)
237 def_orth(i) = def_orth(i) - 1
238 ip = drape(ie)%INDX_PLY(j)
239 IF(ip > 0) THEN
240 drape_ply => drape(ie)%DRAPE_PLY(ip)
241 phi1(j,i) = phi1(j,i) + drape_ply%RDRAPE(1,2)
242 ENDIF
243 ENDIF
244 ENDDO
245 ENDIF
246 ENDDO
247 ELSE
248 DO i=jft,jlt
249 ipang = 1
250 DO j=1,nlay
251 ipid_ply = stack%IGEO(2+j,isubstack)
252 IF(ipid_ply > 0) THEN
253 phi1(j,i) = angle(i) + geo(2,ipid_ply) + stack%GEO(ipang+j,isubstack)
254 def_orth(i) = def_orth(i) - 1
255 IF (irep >= 2) phi2(j,i)= stack%GEO(ipdir+j,isubstack)
256 ENDIF
257 ENDDO
258 ENDDO
259 ENDIF
260 ELSEIF(igtyp == 51 .AND. irp /= 24 ) THEN
261 IF(idrape > 0) THEN
262 DO i=jft,jlt
263 ipang = 1
264 ipt_all = 0
265 ie = indx(nft + i)
266 IF(ie > 0) THEN
267 DO il=1,nlay
268 nptt = elbuf_str%BUFLY(il)%NPTT
269 ip = drape(ie)%INDX_PLY(il)
270 ipid_ply = stack%IGEO(2 + il,isubstack)
271 IF(ipid_ply > 0) THEN
272 IF(ip > 0) THEN
273 drape_ply => drape(ie)%DRAPE_PLY(ip)
274 nslice = drape_ply%NSLICE
275 def_orth(i) = def_orth(i) - 1
276 IF(irep >= 2) THEN
277 DO it = 1,nptt
278 ipt = ipt_all + it
279 phi1(ipt,i) = angle(i) + geo(2,ipid_ply) + stack%GEO(ipang+il,isubstack)
280 . + drape_ply%RDRAPE(it,2)
281 phi2(ipt
282 ENDDO
283 ELSE
284 DO it = 1,nptt
285 ipt = ipt_all + it
286 phi1(ipt,i) = angle(i) + geo(2,ipid_ply) + stack%GEO(ipang+il,isubstack
287 . + drape_ply%RDRAPE(it,2)
288 ENDDO ! nptt
289 ENDIF
290 ELSE
291 def_orth(i) = def_orth(i) - 1
292 IF(irep >= 2) THEN
293 DO it = 1,nptt
294 ipt = ipt_all + it
295 phi1(ipt,i) = angle(i) + geo(2,ipid_ply) + stack%GEO(ipang + il,isubstack)
296 phi2(ipt,i) = stack%GEO(ipdir+il,isubstack)
297 ENDDO
298 ELSE
299 DO it = 1,nptt
300 ipt = ipt_all + it
301 phi1(ipt,i) = angle(i) + geo(2,ipid_ply) + stack%GEO(ipang + il,isubstack)
302 ENDDO
303 ENDIF
304 ENDIF
305 ENDIF
306 ipt_all = ipt_all + nptt
307 ENDDO
308 ELSE
309 DO il=1,nlay
310 nptt = elbuf_str%BUFLY(il)%NPTT
311 ipid_ply = stack%IGEO(2 + il,isubstack)
312 IF(ipid_ply > 0) THEN
313 def_orth(i) = def_orth(i) - 1
314 IF(irep >= 2) THEN
315 DO it = 1,nptt
316 ipt = ipt_all + it
317 phi1(ipt,i) = angle(i) + geo(2,ipid_ply) + stack%GEO(ipang + il,isubstack)
318 phi2(ipt,i) = stack%GEO(ipdir+il,isubstack)
319 ENDDO
320 ELSE
321 DO it = 1,nptt
322 ipt = ipt_all + it
323 phi1(ipt,i) = angle(i) + geo(2,ipid_ply) + stack%GEO(ipang + il,isubstack)
324 ENDDO
325 ENDIF
326 ENDIF
327 ipt_all = ipt_all + nptt
328 ENDDO
329 ENDIF
330 ENDDO
331 ELSE
332 DO i=jft,jlt
333 ipang = 1
334 DO il=1,nlay
335 ipid_ply = stack%IGEO(2 + il,isubstack)
336 IF(ipid_ply > 0 ) THEN
337 phi1(il,i) = angle(i) + geo(2,ipid_ply) + stack%GEO(ipang + il,isubstack)
338 def_orth(i) = def_orth(i) - 1
339 ENDIF
340 IF (irep >= 2) phi2(il,i)= stack%GEO(ipdir + il,isubstack)
341 ENDDO
342 ENDDO
343 ENDIF
344 ELSEIF(igtyp == 52 .AND. irp /= 24 ) THEN
345 IF(idrape > 0) THEN
346 DO i=jft,jlt
347 ipang = 1
348 ipt_all = 0
349 ie = indx(nft + i)
350 IF(ie > 0) THEN
351 DO il=1,nlay
352 nptt = elbuf_str%BUFLY(il)%NPTT
353 ip = drape(ie)%INDX_PLY(il)
354 ipid_ply = stack%IGEO(2+il,isubstack)
355 IF( ipid_ply > 0) THEN
356 IF(ip > 0) THEN
357 drape_ply => drape(ie)%DRAPE_PLY(ip)
358 nslice = drape_ply%NSLICE
359 def_orth(i) = def_orth(i) - 1
360 IF(irep >= 2) THEN
361 DO it = 1,nptt
362 ipt = ipt_all + it
363 phi1(ipt,i) = angle(i) + geo_stack(2,ipid_ply)
364 . + stack%GEO(ipang + il,isubstack) + drape_ply%RDRAPE(it,2)
365 phi2(ipt,i)= stack%GEO(ipdir+il,isubstack)
366 ENDDO
367 ELSE
368 DO it = 1,nptt
369 ipt = ipt_all + it
370 phi1(ipt,i) = angle(i) + geo_stack(2,ipid_ply)
371 . + stack%GEO(ipang + il,isubstack) + drape_ply%RDRAPE(it
372 ENDDO
373 ENDIF
374 ELSE
375 def_orth(i) = def_orth(i) - 1
376 IF(irep >= 2) THEN
377 DO it = 1,nptt
378 ipt = ipt_all + it
379 phi1(ipt,i) = angle(i) + geo_stack(2,ipid_ply) + stack%GEO(ipang + il,isubstack)
380 phi2(ipt,i) = stack%GEO(ipdir+il,isubstack)
381 ENDDO
382 ELSE
383 DO it = 1,nptt
384 ipt = ipt_all + it
385 phi1(ipt,i) = angle(i) + geo_stack(2,ipid_ply) + stack%GEO(ipang + il,isubstack)
386 ENDDO
387 ENDIF
388 ENDIF
389 ENDIF
390 ipt_all = ipt_all + nptt
391 ENDDO
392 ELSE
393 DO il=1,nlay
394 nptt = elbuf_str%BUFLY(il)%NPTT
395 ipid_ply = stack%IGEO(2+il,isubstack)
396 IF(ipid_ply > 0) THEN
397 def_orth(i) = def_orth(i) - 1
398 IF(irep >= 2)THEN
399 DO it = 1,nptt
400 ipt = ipt_all + it
401 phi1(ipt,i) = angle(i) + geo_stack(2,ipid_ply) + stack%GEO(ipang + il,isubstack)
402 phi2(ipt,i)= stack%GEO(ipdir+il,isubstack)
403 ENDDO
404 ELSE
405 DO it = 1,nptt
406 ipt = ipt_all + it
407 phi1(ipt,i) = angle(i) + geo_stack(2,ipid_ply) + stack%GEO(ipang + il,isubstack)
408 ENDDO
409 ENDIF
410 ENDIF
411 ipt_all = ipt_all + nptt
412 ENDDO
413 ENDIF
414 ENDDO
415 ELSE
416 DO i=jft,jlt
417 ipang = 1
418 DO il=1,nlay
419 ipid_ply = stack%IGEO(2+il,isubstack)
420 IF(ipid_ply > 0) THEN
421 def_orth(i) = def_orth(i) - 1
422 phi1(il,i) = angle(i) + geo_stack(2,ipid_ply) + stack%GEO(ipang + il,isubstack)
423 IF(irep >= 2) phi2(il,i)= stack%GEO(ipdir
424 ENDIF
425 ENDDO
426 ENDDO
427 ENDIF
428 ELSEIF (igtyp == 16) THEN
429 DO i=jft,jlt
430 DO j=1,nlay
431 phi1(j,i)= geo(ipang+j,pid)
432 phi2(j,i)= geo(ipphi+j,pid)
433 ENDDO
434 ENDDO
435 ENDIF
436
437 IF (iortshel == 1) THEN
438 DO i=jft,jlt
439 IF (abs(isigi) /= 3 .AND. abs(isigi)/=4 .AND. abs(isigi)/=5) THEN
440
442 ii = ptsh(i + nft)
443 IF(ii == 0) GOTO 100
444 n = nint(sigsh(1,ii))
446 CONTINUE
447 ELSE
448 DO j = 1,numel
449 ii = j
450 n = nint(sigsh(1,ii))
452 IF (n == 0) GOTO 100
453 ENDDO
454 GOTO 100
455 60 CONTINUE
456 ENDIF
457 ELSE
458 jj=nft+i
459 n =ix(nix,i)
460 ii=ptsh(jj)
461 IF (ii == 0) GOTO 100
462 END IF
463 IF(sigsh(nvshell + nushell + 5,ii) == zero) cycle
464
465 npti = nint(sigsh(nvshell + nushell + 4,ii))
466 IF(igtyp == 9) npti = 1
467 IF (nlay /= npti) THEN
468 ipid = ix(nix-1,i)
469 pid = igeo(1,ipid)
470 CALL fretitl2(titr1,igeo(npropgi-ltitr+1,ipid),ltitr)
471 IF (npti == 0) THEN
473 . msgtype=msgwarning,
474 . anmode=aninfo_blind_1,
475 . i1=n,
476 . i2=pid,
477 . c1=titr1)
478 ELSE
480 . anmode=aninfo,
481 . msgtype=msgerror,
482 . i2=n,
483 . i1=pid,
484 . c1=titr1)
485 ENDIF
486 ENDIF
487
488 ipt = nvshell + nushell
489 vx(i) = sigsh(ipt+1,ii)
490 vy(i) = sigsh(ipt+2,ii)
491 vz(i) = sigsh(ipt+3,ii)
492 ipt = ipt + 5
493 IF ( igtyp == 9) THEN
494 phi1(1,i) = sigsh(ipt+1,ii)
495 phi2(1,i) = sigsh(ipt+2,ii)
496 ipt = ipt + 2
497 ELSE
498 DO j=1,npti
499 phi1(j,i) = sigsh(ipt+1,ii)
500 phi2(j,i) = sigsh(ipt+2,ii)
501 ipt = ipt + 2
502 ENDDO
503 ENDIF
504 100 CONTINUE
505 ENDDO
506 ENDIF
507
508
509 IF (iortshel == 2) THEN
510 DO i=jft,jlt
511 ie = i + nft
512 IF (abs(isigi) /= 3 .AND. abs(isigi) /= 4 .AND. abs(isigi) /= 5) THEN
514 ii = ptsh(ie)
515 IF(ii == 0) GOTO 110
516 n = nint(sigsh(1,ii))
518 CONTINUE
519 ELSE
520 DO j = 1,numel
521 ii = j
522 n = nint(sigsh(1,ii))
524 IF (n == 0) GOTO 110
525 ENDDO
526 GOTO 110
527 70 CONTINUE
528 ENDIF
529 ELSE
530 jj=nft+i
531 n =ix(nix,i)
532 ii=ptsh(jj)
533 IF (ii == 0) GOTO 110
534 END IF
535 IF(sigsh(nvshell + nushell + 5,ii) == zero) cycle
536 npti = nint(sigsh(nvshell + nushell + 4,ii))
537
538 npt = nint(geo(6,ix(nix-1,i)))
539 IF(idrape > 0 .AND. (igtyp == 51 .OR. igtyp==52)) THEN
540 npt = npt_all
541 ELSEIF (igtyp == 16 .OR. igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
542 npt = nlay
543 ENDIF
544 IF (npt /= npti) THEN
545 ipid = ix(nix-1,i)
546 pid = igeo(1,ipid)
547 CALL fretitl2(titr1,igeo(npropgi-ltitr+1,ipid),ltitr)
548 IF (npti == 0) THEN
550 . msgtype=msgwarning,
551 . anmode=aninfo_blind_1,
552 . i1=n,
553 . i2=pid,
554 . c1=titr1)
555 ELSE
557 . anmode=aninfo,
558 . msgtype=msgerror,
559 . i2=n,
560 . i1=pid,
561 . c1=titr1)
562 ENDIF
563 ENDIF
564
565 ipt = nvshell + nushell + 5
566 IF (igtyp == 9) THEN
567 coor1(1,i) = sigsh(ipt+1,ii)
568 coor2(1,i) = sigsh(ipt+2,ii)
569 ipt = ipt + 2
570 iorthloc(i) = 1
571 ELSEIF (igtyp == 10 .OR. igtyp == 11 .OR. igtyp == 16 .OR.
572 . igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
573 DO j=1,npti
574 ilaw_ly = elbuf_str%BUFLY(j)%ILAW
575 IF (igtyp == 16 .OR.(igtyp == 51 .AND. ilaw_ly == 58)
576 . .OR.(igtyp == 52 .AND. ilaw_ly == 58)THEN
577 coor1(j,i) = sigsh(ipt+1,ii)
578 coor2(j,i) = sigsh(ipt+2,ii)
579 coor3(j,i) = sigsh(ipt+3,ii)
580 coor4(j,i) = sigsh(ipt+4,ii)
581 ipt = ipt + 4
582 ELSE
583 coor1(j,i) = sigsh(ipt+1,ii)
584 coor2(j,i) = sigsh(ipt+2,ii)
585 ipt = ipt + 2
586 ENDIF
587 ENDDO
588 iorthloc(i) = 1
589 ENDIF
590 110 CONTINUE
591 ENDDO
592 ENDIF
593
594
595 IF(irp /= 26 ) THEN
596 DO i=jft,jlt
597 v(i) =vx(i)*e3x(i)+vy(i)*e3y(i)+vz(i)*e3z(i)
598 vx(i)=vx(i)-v(i)*e3x(i)
599 vy(i)=vy(i)-v(i)*e3y(i)
600 vz(i)=vz(i)-v(i)*e3z(i)
601 v(i) =sqrt(vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i))
602 ENDDO
603
604 DO i=jft,jlt
605 IF (v(i) < em3 .AND. iorthloc(i) == 0 .AND.
606 . def_orth(i) /= 0)THEN
607 pid = ix(nix-1,i)
609 CALL fretitl2(titr1,igeo(npropgi-ltitr+1,pid),ltitr)
611 . msgtype=msgerror,
612 . anmode=aninfo,
613 . i1=igeo(1,pid),
614 . c1=titr1,
615 . i2=ix(nix,i))
616 ENDIF
617 vx(i)=vx(i)/v(i)
618 vy(i)=vy(i)/v(i)
619 vz(i)=vz(i)/v(i)
620 ENDDO
621 ENDIF
622
623
624
625 DO i=jft,jlt
626
627 e11= x2(i)-x1(i)
628 e12= y2(i)-y1(i)
629 e13= z2(i)-z1(i)
630 ne1 = sqrt(e11*e11+e12*e12+e13*e13)
631
632 betaorth(i) = (vx(i)*e11 + vy(i)*e12 +vz(i)*e13 )/
max(ne1,em20)
633 ENDDO
634
635 ENDIF
636
637 RETURN
subroutine add_node(n, ancestor, parent, j, i)
integer, parameter nchartitle
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)