34 SUBROUTINE i14frt(OUTPUT,AF ,X ,V ,KSURF ,IGRSURF,
35 2 BUFSF ,NSC ,KSC ,NSP ,KSP ,
36 3 KSI ,IMPACT ,CIMP ,NIMP ,FRIC ,
37 4 NFRIC ,NPC ,PLD ,GAPMIN ,STF ,
38 5 WF ,WST ,DE ,MS ,STIFN ,
39 6 FS ,FCONT ,FSKYI ,ISKY ,H3D_DATA)
50#include "implicit_f.inc"
71 TYPE(output_) ,
INTENT(INOUT) :: OUTPUT
72 INTEGER NSC,NSP, KSURF, KSI(*), IMPACT(*),
73 . NFRIC, NPC(*),ISKY(*)
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
81 TYPE (SURF_) ,
DIMENSION(NSURF) :: IGRSURF
85 INTEGER ADRBUF, I, IL, IN, I3, I2, I1
90 . a, b, c, an, bn, cn, rot(9),
91 . x1, x2, x3, n1, n2, n3, n,
92 . xpvn1, ypvn1, zpvn1, sgnxn, sgnyn, sgnzn,
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,
99 . xpv(3,mvsiz),nv(3,mvsiz) ,tv(3,mvsiz),
100 . fnpv(mvsiz) ,kfric(mvsiz),xq(3,mvsiz),
101 . fv(3,mvsiz) ,st(mvsiz)
103 adrbuf=igrsurf(ksurf)%IAD_BUFR
116 rot(i)=bufsf(adrbuf+7+i-1)
128 nisky = nisky + nsc + nsp
129#include "lockoff.inc"
140 DO i=1,
min(mvsiz,nrest)
143 fxn=wf(in)*nimp(1,il)
144 fyn=wf(in)*nimp(2,il)
145 fzn=wf(in)*nimp(3,il)
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
165 DO i=1,
min(mvsiz,nrest)
169 st(i)=stf+two*wst(in)*dt1inv
175 DO i=1,
min(mvsiz,nrest)
179 st(i)=stf+wst(in)*dt1inv
186 DO i=1,
min(mvsiz,nrest)
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)
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)
213 dd = x1**2+x2**2+x3**2
214 bufsf(adrbuf+32)=bufsf(adrbuf+32)+dd*st(i)
219 DO i=1,
min(mvsiz,nrest)
234#include "vectorize.inc"
235 DO i=1,
min(mvsiz,nrest)
242 af(i2)=af(i2)+fv(2,i)
250 DO i=1,
min(mvsiz,nrest)
255 fskyi(niskyl,2)=fv(2,i)
256 fskyi(niskyl,3)=fv(3,i)
258 fskyi(niskyl,4)=st(i)
262 DO i=1,
min(mvsiz,nrest)
266 fskyi(niskyl,1)=fv(1,i)
267 fskyi(niskyl,2)=fv(2,i)
268 fskyi(niskyl,3)=fv(3,i)
270 fskyi(niskyl,4)=st(i)
271 fskyi(niskyl,5)=wst(in)
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
285#include "vectorize.inc"
286 DO i=1,
min(mvsiz,nrest)
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)
293#include "lockoff.inc"
300 DO i=1,
min(mvsiz,nrest)
303 de=de+fv(1,i)*v(1,in)+fv(2,i)*v(2,in)+fv(3,i)*v(3,in)
308 IF (nrest-mvsiz>0)
THEN
322 DO i=1,
min(mvsiz,nrest)
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
335 DO i=1,
min(mvsiz,nrest)
348 DO i=1,
min(mvsiz,nrest)
352 CALL ninterp(nfric,npc,pld,
min(mvsiz,nrest),fnpv,kfric)
353 DO i=1,
min(mvsiz,nrest)
354 kfric(i)=fric*kfric(i)
360 DO i=1,
min(mvsiz,nrest)
372 fmax=
max(kfric(i)*fnpv(i),zero)
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)
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))
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)
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)
420#include "vectorize.inc"
421 DO i=1,
min(mvsiz,nrest)
427 xpvn1 =xq(1,i)**(dgr-1)
429 IF (xpvn1*xq(1,i)>=zero) sgnxn=one
430 ypvn1 =xq(2,i)**(dgr-1)
432 IF (ypvn1*xq(2,i)>=zero) sgnyn=one
433 zpvn1 =xq(3,i)**(dgr-1)
435 IF (zpvn1*xq(3,i)>=zero) sgnzn=one
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
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
473 dist = (xi-xpv(1,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))
480 psca=fx*nv(1,i)+fy*nv(2,i)+fz*nv(3,i)
484 psca=fx*tv(1,i)+fy*tv(2,i)+fz*tv(3,i)
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
513 DO i=1,
min(mvsiz,nrest)
516 st(i)=stf+two*wst(in)*dt1inv
519 DO i=1,
min(mvsiz,nrest)
522 st(i)=stf+wst(in)*dt1inv
526 DO i=1,
min(mvsiz,nrest)
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)
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)
553 dd = x1**2+x2**2+x3**2
554 bufsf(adrbuf+32)=bufsf(adrbuf+32)+dd*st(i)
562#include "vectorize.inc"
563 DO i=1,
min(mvsiz,nrest)
569 af(i1)=af(i1)+fv(1,i)
570 af(i2)=af(i2)+fv(2,i)
571 af(i3)=af(i3)+fv(3,i)
573 stifn(in)=stifn(in)+st(i)
578 DO i=1,
min(mvsiz,nrest)
582 fskyi(niskyl,1)=fv(1,i)
583 fskyi(niskyl,2)=fv(2,i)
584 fskyi(niskyl,3)=fv(3,i)
586 fskyi(niskyl,4)=st(i)
590 DO i=1,
min(mvsiz,nrest)
594 fskyi(niskyl,1)=fv(1,i)
595 fskyi(niskyl,2)=fv(2,i)
596 fskyi(niskyl,3)=fv(3,i)
599 fskyi(niskyl,5)=wst(in)
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
613#include "vectorize.inc"
614 DO i=1,
min(mvsiz,nrest)
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)
621#include "lockoff.inc"
628 DO i=1,
min(mvsiz,nrest)
631 de=de+fv(1,i)*v(1,in)+fv(2,i)*v(2,in)+fv(3,i)*v(3,in)
636 IF (nrest-mvsiz>0)
THEN
644 fs(7)=fs(7)+de*dt1*half
645 IF (igrsurf(ksurf)%TYPE==100)
THEN
648 output%TH%WFEXT=output%TH%WFEXT+de*dt1*half