8 common/chainc/ xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
10 REAL*8 xi(3),xidot(3),firr(3),freg(3),dv(3),fd(3),fdr(3)
19 ri2 = (xi(1) - rdens(1))**2 + (xi(2) - rdens(2))**2 +
20 & (xi(3) - rdens(3))**2
23 ri2 = xi(1)**2 + xi(2)**2 + xi(3)**2
48 IF (list(1,i1).GT.0)
THEN
50 CALL
cmfreg(i,rs2,rcrit2,vrfac,nnb,xi,xidot,firr,freg,fd,fdr)
63 dv(1) = xdot(1,j) - xidot(1)
64 dv(2) = xdot(2,j) - xidot(2)
65 dv(3) = xdot(3,j) - xidot(3)
67 rij2 = a1*a1 + a2*a2 + a3*a3
68 drdv = a1*dv(1) + a2*dv(2) + a3*dv(3)
71 IF (rij2.GT.rcrit2) go to 8
74 IF (rij2.GT.rs2.AND.step(j).GT.10.0*smin)
THEN
76 IF (drdv.GT.vrfac) go to 8
84 dr3i = body(j)*dr2i*sqrt(dr2i)
85 firr(1) = firr(1) + a1*dr3i
86 firr(2) = firr(2) + a2*dr3i
87 firr(3) = firr(3) + a3*dr3i
88 fd(1) = fd(1) + (dv(1) - a1*drdv)*dr3i
89 fd(2) = fd(2) + (dv(2) - a2*drdv)*dr3i
90 fd(3) = fd(3) + (dv(3) - a3*drdv)*dr3i
96 dr3i = body(j)*dr2i*sqrt(dr2i)
97 freg(1) = freg(1) + a1*dr3i
98 freg(2) = freg(2) + a2*dr3i
99 freg(3) = freg(3) + a3*dr3i
100 fdr(1) = fdr(1) + (dv(1) - a1*drdv)*dr3i
101 fdr(2) = fdr(2) + (dv(2) - a2*drdv)*dr3i
102 fdr(3) = fdr(3) + (dv(3) - a3*drdv)*dr3i
106 IF (npairs.GT.0)
THEN
107 CALL
cmfreg(i,rs2,rcrit2,vrfac,nnb,xi,xidot,firr,freg,fd,fdr)
113 IF (name(i).EQ.0)
THEN
114 CALL
chfirr(i,1,xi,xidot,firr,fd)
122 CALL
fchain(i,1,xi,xidot,firr,fd)
130 20
IF (kz(13).GT.0)
THEN
135 IF (kz(14).GT.0)
THEN
137 IF (kz(14).EQ.3)
THEN
145 CALL
xtrnlf(xi,xidot,firr,freg,fd,fdr,1)
148 IF (kz(14).EQ.3)
THEN
154 px = freg(k) - frx(k)
155 dpx = fdr(k) - fdx(k)
156 wdot = wdot + xidot(k)*px
157 w2dot = w2dot + (freg(k) + firr(k))*px + xidot(k)*dpx
166 etide = etide + body(i)*(0.5*w2dot*dtr - wdot)*dtr
174 r2 = xi(1)**2 + xi(2)**2 + xi(3)**2
175 fij = 0.01*bodym/(r2*sqrt(r2))
176 rdot = 3.0*(xi(1)*xidot(1) + xi(2)*xidot(2) +
179 firr(k) = firr(k) - fij*xi(k)
180 fd(k) = fd(k) - (xidot(k) - rdot*xi(k))*fij
183 IF (rdot.GT.0.0)
THEN
184 rs(i) = max(0.75*rs(i),0.1*rscale)
195 IF (nnb.LT.nnbmax) go to 40
199 a1 = float(nnb2)/float(nnb)
200 IF (rs(i).GT.5.0*rscale)
THEN
201 a1 = min(5.0*a1,0.9d0)
204 rs2 = rs2*a1**0.66667
211 IF (l + nnb2.GT.nnb + nnb1) go to 32
216 dv(1) = xdot(1,j) - xidot(1)
217 dv(2) = xdot(2,j) - xidot(2)
218 dv(3) = xdot(3,j) - xidot(3)
220 rij2 = a1*a1 + a2*a2 + a3*a3
222 dr3i = body(j)*dr2i*sqrt(dr2i)
223 drdv = a1*dv(1) + a2*dv(2) + a3*dv(3)
224 IF (rij2.GT.rcrit2) go to 34
226 IF (rij2.GT.rs2.AND.step(j).GT.smin)
THEN
227 IF (drdv.GT.vrfac.AND.j.LE.n) go to 34
236 34 drdv = 3.0*drdv*dr2i
237 firr(1) = firr(1) - a1*dr3i
238 firr(2) = firr(2) - a2*dr3i
239 firr(3) = firr(3) - a3*dr3i
240 fd(1) = fd(1) - (dv(1) - a1*drdv)*dr3i
241 fd(2) = fd(2) - (dv(2) - a2*drdv)*dr3i
242 fd(3) = fd(3) - (dv(3) - a3*drdv)*dr3i
243 freg(1) = freg(1) + a1*dr3i
244 freg(2) = freg(2) + a2*dr3i
245 freg(3) = freg(3) + a3*dr3i
246 fdr(1) = fdr(1) + (dv(1) - a1*drdv)*dr3i
247 fdr(2) = fdr(2) + (dv(2) - a2*drdv)*dr3i
248 fdr(3) = fdr(3) + (dv(3) - a3*drdv)*dr3i
257 IF (nnb.GE.nnbmax) go to 30
260 40 a3 = alpha*sqrt(float(nnb)*rs(i))/rs2
265 a3 = a3*rscale/sqrt(ri2)
271 IF ((a3 - float(nnb0))*(a3 - float(nnb)).LT.0.0) a4 = sqrt(a4)
274 IF (ri2.GT.rc2.AND.kz(39).EQ.0.AND.ri2.LT.9.0*rh2)
THEN
275 ridot = (xi(1) - rdens(1))*xidot(1) +
276 & (xi(2) - rdens(2))*xidot(2) +
277 & (xi(3) - rdens(3))*xidot(3)
278 a4 = a4*(1.0 + ridot*dtr/ri2)
282 IF (i.GT.n.AND.ri2.LT.rh2)
THEN
284 a2 = 100.0*body(i)/abs(h(i-n))
286 IF (a2.GT.rs(i))
THEN
287 a3 = max(1.0 - float(nnb)/znbmax,0.0d0)
294 a5 = min(rsfac*dtr,0.25d0)
304 IF (irskip.EQ.0)
THEN
306 a1 = 1.0 + a3 - a3*a3
308 IF (rs(i).GT.5.0*rscale) a1 = sqrt(a1)
310 IF (nnb.LT.znbmin.AND.nbp.GT.znbmin)
THEN
311 IF (a1.LT.1.0) a1 = 1.05
314 IF (abs(a1 - 1.0d0).GT.0.003)
THEN
320 IF (nnb.LE.3.AND.ri2.GT.rh2)
THEN
325 rij = sqrt((xi(1) - x(1,j))**2 + (xi(2) - x(2,j))**2 +
326 & (xi(3) - x(3,j))**2)
327 rsdot = ((xi(1) - x(1,j))*(xidot(1) - xdot(1,j)) +
328 & (xi(2) - x(2,j))*(xidot(2) - xdot(2,j)) +
329 & (xi(3) - x(3,j))*(xidot(3) - xdot(3,j)))/rij
331 a1 = min(a1,rij + rsdot*dtr)
335 rs(i) = max(a1,1.1*rs(i))
337 rs(i) = max(0.1*rscale,rs(i))
341 IF (list(1,i).GT.0) rsmin = min(rsmin,rs(i))
344 IF ((kz(37).EQ.1.AND.listv(1).GT.0).OR.kz(37).GT.1)
THEN
345 CALL
checkl(i,nnb,xi,xidot,rs2,firr,freg,fd,fdr)
353 IF (kz(14).EQ.3)
THEN
357 etide = etide + body(i)*(0.5*w2dot*dtr - wdot)*dtr
366 jlist(l) = ilist(l+1)
375 ilist(nnb+2) = ntot + 1
376 ilist(1) = list(nnb0+1,i)
379 56
IF (list(l,i).EQ.ilist(lg)) go to 58
382 IF (list(l,i).GE.ilist(lg))
THEN
384 jlist(nnb0+nbgain) = ilist(lg)
393 IF (step(j).LT.0.1*dtmin) jmin = j
398 58
IF (l.LE.nnb0)
THEN
403 ELSE IF (lg.LE.nnb)
THEN
411 IF (jmin.EQ.0) go to 70
414 60
IF (nnb.GT.nnbmax.OR.i.GT.n) go to 70
417 IF (step(j).GT.0.1*dtmin.OR.j.LT.ifirst.OR.j.GT.n) go to 68
419 rij2 = (xi(1) - x(1,j))**2 + (xi(2) - x(2,j))**2 +
420 & (xi(3) - x(3,j))**2
421 IF (rij2.GT.2.0*rs2) go to 68
424 62
IF (ilist(l).LT.j) go to 64
425 ilist(l+1) = ilist(l)
434 list(nnb0+1,i) = ilist(1)
441 dv(1) = xdot(1,j) - xidot(1)
442 dv(2) = xdot(2,j) - xidot(2)
443 dv(3) = xdot(3,j) - xidot(3)
445 rij2 = a1*a1 + a2*a2 + a3*a3
447 dr3i = body(j)*dr2i*sqrt(dr2i)
448 drdv = 3.0*(a1*dv(1) + a2*dv(2) + a3*dv(3))*dr2i
450 firr(1) = firr(1) + a1*dr3i
451 firr(2) = firr(2) + a2*dr3i
452 firr(3) = firr(3) + a3*dr3i
453 fd(1) = fd(1) + (dv(1) - a1*drdv)*dr3i
454 fd(2) = fd(2) + (dv(2) - a2*drdv)*dr3i
455 fd(3) = fd(3) + (dv(3) - a3*drdv)*dr3i
456 freg(1) = freg(1) - a1*dr3i
457 freg(2) = freg(2) - a2*dr3i
458 freg(3) = freg(3) - a3*dr3i
459 fdr(1) = fdr(1) - (dv(1) - a1*drdv)*dr3i
460 fdr(2) = fdr(2) - (dv(2) - a2*drdv)*dr3i
461 fdr(3) = fdr(3) - (dv(3) - a3*drdv)*dr3i
464 IF (k.GT.nbloss) go to 70
466 jlist(l) = jlist(l+1)
473 IF (k.LE.nbloss) go to 60
477 dt6 = 6.0d0/(dtr*dtsq)
484 IF (step(i).LT.5.0d0*dtmin.AND.dtr.GT.50.0*step(i))
THEN
492 dfr = fr(k,i) - (firr(k) - fi(k,i)) - freg(k)
493 fdr0 = fdr(k) - (fidot(k,i) - fd(k))
497 at3 = 2.0d0*dfr + dtr*sum
498 bt2 = -3.0d0*dfr - dtr*(sum + frd)
500 x0(k,i) = x0(k,i) + (0.6d0*at3 + bt2)*dtsq12
501 x0dot(k,i) = x0dot(k,i) + (0.75d0*at3 + bt2)*dtr13
508 f(k,i) = 0.5d0*(freg(k) + firr(k))
511 fdot(k,i) = one6*(fdr(k) + fd(k))
521 d2r(k,i) = (3.0d0*at3 + bt2)*dt2
526 IF (kz(38).GT.0)
THEN
527 CALL
fpcorr(i,nbloss,nbgain,xi,xidot)
528 ELSE IF (i.GT.n)
THEN
529 IF (gamma(i-n).GT.gmax)
THEN
530 CALL
fpcorr(i,nbloss,nbgain,xi,xidot)
535 IF (nbloss + nbgain.GT.0)
THEN
540 nbflux = nbflux + nbloss + nbgain
554 fr2 = freg(1)**2 + freg(2)**2 + freg(3)**2
555 w1 = fdr(1)**2 + fdr(2)**2 + fdr(3)**2
556 w2 = d2r(1,i)**2 + d2r(2,i)**2 + d2r(3,i)**2
557 w3 = d3r(1,i)**2 + d3r(2,i)**2 + d3r(3,i)**2
560 w0 = (sqrt(fr2*w2) + w1)/(sqrt(w1*w3) + w2)
570 IF (kz(14).EQ.0)
THEN
571 fi2 = firr(1)**2 + firr(2)**2 + firr(3)**2
572 df2 = fac**2*min(fr2,fi2)
574 IF (nnb.EQ.0) df2 = fac**2*fr2
575 ELSE IF (kz(14).EQ.1)
THEN
576 w0 = (tidal(1)*xi(1))**2
577 df2 = fac**2*max(fr2,w0)
588 w2 = ((d3r(k,i)*s3 + d2r(k,i))*s2 + fdr(k))*dtc
598 IF (nc.LT.50.AND.ri2.LT.rc2)
THEN
599 ttmp = ttmp*min(1.0d0,0.5d0*(1.0d0 + ri2*rc2in))
603 IF (ttmp.GT.2.0*stepr(i))
THEN
604 IF (dmod(time,2.0*stepr(i)).EQ.0.0d0)
THEN
605 ttmp = min(2.0*stepr(i),smax)
609 ELSE IF (ttmp.LT.stepr(i))
THEN
612 IF (ttmp.GT.dt0)
THEN
621 110
IF (ttmp.LT.step(i))
THEN
622 step(i) = 0.5d0*step(i)
623 tnew(i) = t0(i) + step(i)
629 IF (nnb0.EQ.0.AND.nnb.GT.0)
THEN
630 step(i) = 0.25*step(i)
631 tnew(i) = t0(i) + step(i)