41
42
43
47 USE output_mod
48
49
50
51#include "implicit_f.inc"
52#include "comlock.inc"
53
54
55
56#include "mvsiz_p.inc"
57
58
59
60#include "com04_c.inc"
61#include "com06_c.inc"
62#include "com08_c.inc"
63#include "scr07_c.inc"
64#include "scr14_c.inc"
65#include "scr16_c.inc"
66#include "parit_c.inc"
67#include "param_c.inc"
68#include "scr18_c.inc"
69
70
71
72 TYPE(OUTPUT_) ,INTENT(INOUT) :: OUTPUT
73 INTEGER NSC,NSP, KSURF, KSI(*), IMPACT(*),
74 . NFRIC, NPC(*),ISKY(*)
75
77 . af(*), x(3,*), v(3,*) , bufsf(*), ksp(*),
78 . ksc(*), fric, cimp(3,*),nimp(3,*), ms(*),
79 . stifn(*), fs(nthvki),fcont(3,*),fskyi(lskyi,nfskyi), pld(*),
80 . wf(*), wst(*),gapmin , stf , de
81 TYPE(H3D_DATABASE) ::
82 TYPE (SURF_) , DIMENSION(NSURF) :: IGRSURF
83
84
85
86 INTEGER ADRBUF, I, IL, J3, J2, J1, IN, I3, I2, I1
87 INTEGER NISKYL
88 INTEGER NDEB,
89 INTEGER DGR
90 INTEGER IPT,NPT,II,JJ
92 . a, b, c, an, bn, cn, rot(9),
93 . x1, x2, x3, n1, n2, n3, n,
94 . xpvn1, ypvn1, zpvn1, sgnxn, sgnyn, sgnzn,
95 . xi, yi, zi,
96 . fxn, fyn, fzn, fxt, fyt, fzt, ftp, ftpa, fmax, fn, tn,
97 . cost, dist, psca, e, fx, fy, fz,
98 . fn1, fn2, fn3, ft1, ft2, ft3, am1, am2, am3,
99 . dd, dt1inv
101 . xpv(3,mvsiz),nv(3,mvsiz) ,tv(3,mvsiz),
102 . fnpv(mvsiz) ,kfric(mvsiz),xq(3,mvsiz),
103 . fv(3,mvsiz) ,st(mvsiz)
104
105 adrbuf=igrsurf(ksurf)%IAD_BUFR
106
107 a =bufsf(adrbuf+1)
108 b =bufsf(adrbuf+2)
109 c =bufsf(adrbuf+3)
110 dgr=bufsf(adrbuf+36)
111 an=a**dgr
112 bn=b**dgr
113 cn=c**dgr
114 an=1./an
115 bn=1./bn
116 cn=1./cn
117 DO i=1,9
118 rot(i)=bufsf(adrbuf+7+i-1)
119 END DO
120
121 IF (dt1==zero) THEN
122 dt1inv=zero
123 ELSE
124 dt1inv=one/dt1
125 END IF
126
127 IF (iparit/=0) THEN
128#include "lockon.inc"
129 niskyl = nisky
130 nisky = nisky + nsc + nsp
131#include "lockoff.inc"
132 END IF
133
134 de=zero
135
136
137
138 ndeb =0
139 nrest=nsc
140 50 CONTINUE
141
142 DO i=1,
min(mvsiz,nrest)
143 il=ksc(i+ndeb)
144 in=ksi(il)
145 fxn=wf(in)*nimp(1,il)
146 fyn=wf(in)*nimp(2,il)
147 fzn=wf(in)*nimp(3,il)
148
149
150
151 fn1 =rot(1)*fxn+rot(4)*fyn+rot(7)*fzn
152 fn2 =rot(2)*fxn+rot(5)*fyn+rot(8)*fzn
153 fn3 =rot(3)*fxn+rot(6)*fyn+rot(9)*fzn
154 fv(1,i) =fn1
155 fv(2,i) =fn2
156 fv(3,i) =fn3
157
158 fs(1)=fs(1)-fn1*dt1
159 fs(2)=fs(2)-fn2*dt1
160 fs(3)=fs(3)-fn3*dt1
161
162 ENDDO
163
164
165
166 IF(kdtint==0)THEN
167 DO i=1,
min(mvsiz,nrest)
168 il=ksc(i+ndeb)
169 in=ksi(il)
170 IF(impact(il)>0)THEN
171 st(i)=stf+two*wst(in)*dt1inv
172 ELSE
173 st(i)=zero
174 ENDIF
175 ENDDO
176 ELSE
177 DO i=1,
min(mvsiz,nrest)
178 il=ksc(i+ndeb)
179 in=ksi(il)
180 IF(impact(il)>0)THEN
181 st(i)=stf+wst(in)*dt1inv
182 ELSE
183 st(i)=zero
184 ENDIF
185 ENDDO
186 ENDIF
187
188 DO i=1,
min(mvsiz,nrest)
189 il=ksc(i+ndeb)
190 in=ksi(il)
191
192
193
194
195
196 x1=x(1,in)-bufsf(adrbuf+16)
197 x2=x(2,in)-bufsf(adrbuf+17)
198 x3=x(3,in)-bufsf(adrbuf+18)
199 am1=x2*fv(3,i)-x3*fv(2,i)
200 am2=x3*fv(1,i)-x1*fv(3,i)
201 am3=x1*fv(2,i)-x2*fv(1,i)
202
203
204
205
206 bufsf(adrbuf+25)=bufsf(adrbuf+25)-fv(1,i)
207 bufsf(adrbuf+26)=bufsf(adrbuf+26)-fv(2,i)
208 bufsf(adrbuf+27)=bufsf(adrbuf+27)-fv(3,i)
209 bufsf(adrbuf+28)=bufsf(adrbuf+28)-am1
210 bufsf(adrbuf+29)=bufsf(adrbuf+29)-am2
211 bufsf(adrbuf+30)=bufsf(adrbuf+30)-am3
212 bufsf(adrbuf+31)=bufsf(adrbuf+31)+st(i)
213
214
215 dd = x1**2+x2**2+x3**2
216 bufsf(adrbuf+32)=bufsf(adrbuf+32)+dd*st(i)
217
218 ENDDO
219
220 IF(kdtint/=0)THEN
221 DO i=1,
min(mvsiz,nrest)
222 il=ksc(i+ndeb)
223 in=ksi(il)
224 IF(impact(il)>0)THEN
225 st(i)=stf
226 ELSE
227 st(i)=zero
228 ENDIF
229 ENDDO
230 ENDIF
231
232
233
234 IF (iparit==0) THEN
235
236#include "vectorize.inc"
237 DO i=1,
min(mvsiz,nrest)
238 il=ksc(i+ndeb)
239 in=ksi(il)
240 i3=3*in
241 i2=i3-1
242 i1=i2-1
243 af(i1)=af(i1)+fv(1,i)
244 af(i2)=af(i2)+fv(2,i)
245 af(i3)=af(i3)+fv(3,i)
246
247 stifn(in)=stifn(in)+st(i)
248 ENDDO
249 ELSE
250
251 IF(kdtint==0)THEN
252 DO i=1,
min(mvsiz,nrest)
253 il=ksc(i+ndeb)
254 in=ksi(il)
255 niskyl = niskyl + 1
256 fskyi(niskyl,1)=fv(1,i)
257 fskyi(niskyl,2)=fv(2,i)
258 fskyi(niskyl,3)=fv(3,i)
259
260 fskyi(niskyl,4)=st(i)
261 isky(niskyl) =in
262 ENDDO
263 ELSE
264 DO i=1,
min(mvsiz,nrest)
265 il=ksc(i+ndeb)
266 in=ksi(il)
267 niskyl = niskyl + 1
268 fskyi(niskyl,1)=fv(1,i)
269 fskyi(niskyl,2)=fv(2,i)
270 fskyi(niskyl,3)=fv(3,i)
271
272 fskyi(niskyl,4)=st(i)
273 fskyi(niskyl,5)=wst(in)
274
275 isky(niskyl) =in
276 ENDDO
277 ENDIF
278
279 ENDIF
280
281
282
283 IF(anim_v(4)+outp_v(4)+h3d_data%N_VECT_CONT >0.AND.
284 . ((tt>=tanim .AND. tt<=tanim_stop).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP) .OR.
285 . (manim>=4.AND.manim<=15).OR. h3d_data%MH3D /= 0))THEN
286#include "lockon.inc"
287#include "vectorize.inc"
288 DO i=1,
min(mvsiz,nrest)
289 il=ksc(i+ndeb)
290 in=ksi(il)
291 fcont(1,in) =fcont(1,in) + fv(1,i)
292 fcont(2,in) =fcont(2,in) + fv(2,i)
293 fcont(3,in) =fcont(3,in) + fv(3,i)
294 ENDDO
295#include "lockoff.inc"
296 ENDIF
297
298
299
300
301
302 DO i=1,
min(mvsiz,nrest)
303 il=ksc(i+ndeb)
304 in=ksi(il)
305 de=de+fv(1,i)*v(1,in)+fv(2,i)*v(2,in)+fv(3,i)*v(3,in)
306 ENDDO
307
308
309
310 IF (nrest-mvsiz>0) THEN
311 nrest=nrest-mvsiz
312 ndeb =ndeb +mvsiz
313 GOTO 50
314 END IF
315
316
317
318 ndeb =0
319 nrest=nsp
320 100 CONTINUE
321
322
323
324 DO i=1,
min(mvsiz,nrest)
325 il=ksp(i+ndeb)
326 in=ksi(il)
327 x1=x(1,in)-bufsf(adrbuf+4)
328 x2=x(2,in)-bufsf(adrbuf+5)
329 x3=x(3,in)-bufsf(adrbuf+6)
330 xpv(1,i)=rot(1)*x1+rot(2)*x2+rot(3)*x3
331 xpv(2,i)=rot(4)*x1+rot(5)*x2+rot(6)*x3
332 xpv(3,i)=rot(7)*x1+rot(8)*x2+rot(9)*x3
333 ENDDO
334
335
336
337 DO i=1,
min(mvsiz,nrest)
338 il=ksp(i+ndeb)
339 in=ksi(il)
340 fnpv(i) =wf(in)
341 nv(1,i) =nimp(1,il)
342 nv(2,i) =nimp(2,il)
343 nv(3,i) =nimp(3,il)
344 ENDDO
345
346
347
348
349 IF (nfric==0) THEN
350 DO i=1,
min(mvsiz,nrest)
351 kfric(i)=fric
352 ENDDO
353 ELSE
354 CALL ninterp(nfric,npc,pld,
min(mvsiz,nrest),fnpv,kfric)
355 DO i=1,
min(mvsiz,nrest)
356 kfric(i)=fric*kfric(i)
357 ENDDO
358 ENDIF
359
360
361
362 DO i=1,
min(mvsiz,nrest)
363 il=ksp(i+ndeb)
364 in=ksi(il)
365
366 fxt=zero
367 fyt=zero
368 fzt=zero
369
370 xi=cimp(1,il)
371 yi=cimp(2,il)
372 zi=cimp(3,il)
373
374 fmax=
max(kfric(i)*fnpv(i),zero)
375
376 fxt=(xi-xpv(1,i))*stf
377 fyt=(yi-xpv(2,i))*stf
378 fzt=(zi-xpv(3,i))*stf
379 ftp=sqrt(fxt*fxt+fyt*fyt+fzt*fzt)
380
381 fxn=fnpv(i)*nv(1,i)
382 fyn=fnpv(i)*nv(2,i)
383 fzn=fnpv(i)*nv(3,i)
384
385 fn=fxt*nv(1,i)+fyt*nv(2,i)+fzt*nv(3,i)
386 tv(1,i)=fxt-nv(1,i)*fn
387 tv(2,i)=fyt-nv(2,i)*fn
388 tv(3,i)=fzt-nv(3,i)*fn
389 tn=sqrt(tv(1,i)*tv(1,i)+tv(2,i)*tv(2,i)+tv(3,i)*tv(3,i))
391 tv(1,i)=tv(1,i)/tn
392 tv(2,i)=tv(2,i)/tn
393 tv(3,i)=tv(3,i)/tn
394
395
397 fxt=tv(1,i)*ftpa
398 fyt=tv(2,i)*ftpa
399 fzt=tv(3,i)*ftpa
400
401
402 cost =(xpv(1,i)-xi)*tv(1,i)
403 . +(xpv(2,i)-yi)*tv(2,i)
404 . +(xpv(3,i)-zi)*tv(3,i)
405 xq(1,i) =xi+cost*tv(1,i)
406 xq(2,i) =yi+cost*tv(2,i)
407 xq(3,i) =zi+cost*tv(3,i)
408
409
410 xq(1,i) =xq(1,i)+ftpa*tv(1,i)/
max(em20,stf)
411 xq(2,i) =xq(2,i)+ftpa*tv(2,i)/
max(em20,stf)
412 xq(3,i) =xq(3,i)+ftpa*tv(3,i)/
max(em20,stf)
413
414 fv(1,i)=fxt+fxn
415 fv(2,i)=fyt+fyn
416 fv(3,i)=fzt+fzn
417
418 ENDDO
419
420
421
422#include "vectorize.inc"
423 DO i=1,
min(mvsiz,nrest)
424 il=ksp(i+ndeb)
425 in=ksi(il)
426
427
428
429 xpvn1 =xq(1,i)**(dgr-1)
430 sgnxn=-one
431 IF (xpvn1*xq(1,i)>=zero) sgnxn=one
432 ypvn1 =xq(2,i)**(dgr-1)
433 sgnyn=-one
434 IF (ypvn1*xq(2,i)>=zero) sgnyn=one
435 zpvn1 =xq(3,i)**(dgr-1)
436 sgnzn=-one
437 IF (zpvn1*xq(3,i)>=zero) sgnzn=one
438 n1 =sgnxn*xpvn1*an
439 n2 =sgnyn*ypvn1*bn
440 n3 =sgnzn*zpvn1*cn
441 n =n1*n1+n2*n2+n3*n3
442 n =sqrt(n)
443
444
445
446
447
448 e =n1*xq(1,i)+n2*xq(2,i)+n3*xq(3,i)
449 dist =(e-exp((dgr-1)*log(
max(e,em20))/dgr))
450 . /
max(exp((dgr-1)*log(em20)/dgr),n)
454 cimp(1,il)=xq(1,i)-dist*n1
455 cimp(2,il)=xq(2,i)-dist*n2
456 cimp(3,il)=xq(3,i)-dist*n3
457
458 nimp(1,il)=n1
459 nimp(2,il)=n2
460 nimp(3,il)=n3
461
462
463
464
465 fx=fv(1,i)
466 fy=fv(2,i)
467 fz=fv(3,i)
468
469 xi=cimp(1,il)+gapmin*nimp(1,il)
470 yi=cimp(2,il)+gapmin*nimp(2,il)
471 zi=cimp(3,il)+gapmin*nimp(3,il)
472 psca = (xi-xpv(1,i))*fx
473 . +(yi-xpv(2,i))*fy
474 . +(zi-xpv(3,i))*fz
475 dist = (xi-xpv(1,i))**2
476 . +(yi-xpv(2,i))**2
477 . +(zi-xpv(3,i))**2
478 fx=psca/dist*(xi-xpv(1,i))
479 fy=psca/dist*(yi-xpv(2,i))
480 fz=psca/dist*(zi-xpv(3,i))
481
482 psca=fx*nv(1,i)+fy*nv(2,i)+fz*nv(3,i)
483 fxn =psca*nv(1,i)
484 fyn =psca*nv(2,i)
485 fzn =psca*nv(3,i)
486 psca=fx*tv(1,i)+fy*tv(2,i)+fz*tv(3,i)
487 fxt =psca*tv(1,i)
488 fyt =psca*tv(2,i)
489 fzt =psca*tv(3,i)
490
491
492
493 fn1 =rot(1)*fxn+rot(4)*fyn+rot(7)*fzn
494 fn2 =rot(2)*fxn+rot(5)*fyn+rot(8)*fzn
495 fn3 =rot(3)*fxn+rot(6)*fyn+rot(9)*fzn
496 ft1 =rot(1)*fxt+rot(4)*fyt+rot(7)*fzt
497 ft2 =rot(2)*fxt+rot(5)*fyt+rot(8)*fzt
498 ft3 =rot(3)*fxt+rot(6)*fyt+rot(9)*fzt
499 fv(1,i) =fn1+ft1
500 fv(2,i) =fn2+ft2
501 fv(3,i) =fn3+ft3
502
503 fs(1)=fs(1)-fn1*dt1
504 fs(2)=fs(2)-fn2*dt1
505 fs(3)=fs(3)-fn3*dt1
506 fs(4)=fs(4)-ft1*dt1
507 fs(5)=fs(5)-ft2*dt1
508 fs(6)=fs(6)-ft3*dt1
509
510 ENDDO
511
512
513
514 IF(kdtint==0)THEN
515 DO i=1,
min(mvsiz,nrest)
516 il=ksp(i+ndeb)
517 in=ksi(il)
518 st(i)=stf+two*wst(in)*dt1inv
519 ENDDO
520 ELSE
521 DO i=1,
min(mvsiz,nrest)
522 il=ksp(i+ndeb)
523 in=ksi(il)
524 st(i)=stf+wst(in)*dt1inv
525 ENDDO
526 ENDIF
527
528 DO i=1,
min(mvsiz,nrest)
529 il=ksp(i+ndeb)
530 in=ksi(il)
531
532
533
534
535
536 x1=x(1,in)-bufsf(adrbuf+16)
537 x2=x(2,in)-bufsf(adrbuf+17)
538 x3=x(3,in)-bufsf(adrbuf+18)
539 am1=x2*fv(3,i)-x3*fv(2,i)
540 am2=x3*fv(1,i)-x1*fv(3,i)
541 am3=x1*fv(2,i)-x2*fv(1,i)
542
543
544
545
546 bufsf(adrbuf+25)=bufsf(adrbuf+25)-fv(1,i)
547 bufsf(adrbuf+26)=bufsf(adrbuf+26)-fv(2,i)
548 bufsf(adrbuf+27)=bufsf(adrbuf+27)-fv(3,i)
549 bufsf(adrbuf+28)=bufsf(adrbuf+28)-am1
550 bufsf(adrbuf+29)=bufsf(adrbuf+29)-am2
551 bufsf(adrbuf+30)=bufsf(adrbuf+30)-am3
552 bufsf(adrbuf+31)=bufsf(adrbuf+31)+st(i)
553
554
555 dd = x1**2+x2**2+x3**2
556 bufsf(adrbuf+32)=bufsf(adrbuf+32)+dd*st(i)
557
558 ENDDO
559
560
561
562 IF (iparit==0) THEN
563
564#include "vectorize.inc"
565 DO i=1,
min(mvsiz,nrest)
566 il=ksp(i+ndeb)
567 in=ksi(il)
568 i3=3*in
569 i2=i3-1
570 i1=i2-1
571 af(i1)=af(i1)+fv(1,i)
572 af(i2)=af(i2)+fv(2,i)
573 af(i3)=af(i3)+fv(3,i)
574
575 stifn(in)=stifn(in)+st(i)
576 ENDDO
577 ELSE
578
579 IF(kdtint==0)THEN
580 DO i=1,
min(mvsiz,nrest)
581 il=ksp(i+ndeb)
582 in=ksi(il)
583 niskyl = niskyl + 1
584 fskyi(niskyl,1)=fv(1,i)
585 fskyi(niskyl,2)=fv(2,i)
586 fskyi(niskyl,3)=fv(3,i)
587
588 fskyi(niskyl,4)=st(i)
589 isky(niskyl) =in
590 ENDDO
591 ELSE
592 DO i=1,
min(mvsiz,nrest)
593 il=ksp(i+ndeb)
594 in=ksi(il)
595 niskyl = niskyl + 1
596 fskyi(niskyl,1)=fv(1,i)
597 fskyi(niskyl,2)=fv(2,i)
598 fskyi(niskyl,3)=fv(3,i)
599
600 fskyi(niskyl,4)=stf
601 fskyi(niskyl,5)=wst(in)
602
603 isky(niskyl) =in
604 ENDDO
605 ENDIF
606
607 ENDIF
608
609
610
611 IF(anim_v(4)+outp_v(4)+h3d_data%N_VECT_CONT>0.AND.
612 . (tt>=tanim .AND. tt<=tanim_stop.OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP).OR.
613 . (manim>=4.AND.manim<=15).OR. h3d_data%MH3D /= 0))THEN
614#include "lockon.inc"
615#include "vectorize.inc"
616 DO i=1,
min(mvsiz,nrest)
617 il=ksp(i+ndeb)
618 in=ksi(il)
619 fcont(1,in) =fcont(1,in) + fv(1,i)
620 fcont(2,in) =fcont(2,in) + fv(2,i)
621 fcont(3,in) =fcont(3,in) + fv(3,i)
622 ENDDO
623#include "lockoff.inc"
624 ENDIF
625
626
627
628
629
630 DO i=1,
min(mvsiz,nrest)
631 il=ksp(i+ndeb)
632 in=ksi(il)
633 de=de+fv(1,i)*v(1,in)+fv(2,i)*v(2,in)+fv(3,i)*v(3,in)
634 ENDDO
635
636
637
638 IF (nrest-mvsiz>0) THEN
639 nrest=nrest-mvsiz
640 ndeb =ndeb +mvsiz
641 GOTO 100
642 END IF
643
644
645
646 fs(7)=fs(7)+de*dt1*half
647 IF (igrsurf(ksurf)%TYPE==100) THEN
648
649
650 output%TH%WFEXT=output%TH%WFEXT+de*dt1*half
651 ENDIF
652
653 RETURN
subroutine ninterp(ifunc, npc, pld, npoint, xv, yv)