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