Nbody6
 All Files Functions Variables
chmod.f
Go to the documentation of this file.
1  SUBROUTINE chmod(ISUB,KCASE)
2 *
3 *
4 * Modification of chain member(s).
5 * --------------------------------
6 *
7  include 'common6.h'
8  parameter(nmx=10,nmx2=2*nmx,nmx3=3*nmx,nmx4=4*nmx,
9  & nmx8=8*nmx,nmxm=nmx*(nmx-1)/2)
10  REAL*8 m,mass,mc,mij,mkk,xcm(3),vcm(3),ksch,y(nmx8)
11  INTEGER isort(nmx)
12  LOGICAL kslow2,kcoll
13  common/chain1/ xch(nmx3),vch(nmx3),m(nmx),
14  & zz(nmx3),wc(nmx3),mc(nmx),
15  & xi(nmx3),pi(nmx3),mass,rinv(nmxm),rsum,mkk(nmx),
16  & mij(nmx,nmx),tkk(nmx),tk1(nmx),iname(nmx),nn
17  common/chainc/ xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
18  & listc(lmax)
19  common/slow1/ tk2(0:nmx),ejump,ksch(nmx),kslow2,kcoll
20  common/cpert/ rgrav,gpert,ipert,npert
21  common/chreg/ timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
22  common/clump/ bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
23  & names(ncmax,5),isys(5)
24  common/ccoll2/ qk(nmx4),pk(nmx4),rik(nmx,nmx),SIZE(nmx),vstar1,
25  & ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
26  common/slow3/ gcrit,kz26
27  common/incond/ x4(3,nmx),xdot4(3,nmx)
28 *
29 *
30 * Identify the dominant perturber (skip if none or NN >= 6).
31  itry = 0
32  jclose = 0
33  1 nnb = listc(1)
34  IF (nnb.EQ.0.OR.nn.GE.6) go to 10
35  pmax = 0.0
36  DO 2 l = 2,nnb+1
37  j = listc(l)
38  rij2 = (x(1,j) - x(1,ich))**2 + (x(2,j) - x(2,ich))**2 +
39  & (x(3,j) - x(3,ich))**2
40  pij = body(j)/(rij2*sqrt(rij2))
41  IF (pij.GT.pmax) THEN
42  pmax = pij
43  rjmin2 = rij2
44  jclose = j
45  END IF
46  2 CONTINUE
47 *
48 * Form the scalar product R*V for sign of radial velocity.
49  rdot = (x(1,jclose) - x(1,ich))*(xdot(1,jclose) - xdot(1,ich)) +
50  & (x(2,jclose) - x(2,ich))*(xdot(2,jclose) - xdot(2,ich)) +
51  & (x(3,jclose) - x(3,ich))*(xdot(3,jclose) - xdot(3,ich))
52 *
53 * Check for rejection (RIJ > 3*MIN(RSUM,RMIN); RDOT > 0 & G < 0.05).
54  rij = sqrt(rjmin2)
55 *
56 * Include test on fast escaper approaching a perturber with RDOT > 0.
57  IF (gpert.GT.0.05.AND.rdot.GT.0) THEN
58 * Bypass test and repeated diagnostics for large perturbation.
59  IF (gpert.GT.0.5) go to 5
60  rdx = 0.0
61  rjx = 1.0
62 * Determine the smallest approaching intruder distance to chain member.
63  DO 3 l = 1,nn
64  rjx2 = 0.0
65  rdi = 0.0
66  DO k = 1,3
67  rjx2 = rjx2 + (xc(k,l) - x(k,jclose))**2
68  rdi = rdi + (xc(k,l) - x(k,jclose))*
69  & (uc(k,l) - xdot(k,jclose))
70  END DO
71 * See whether the intruder is approaching a chain member.
72  IF (rdi.LT.0.0.AND.rjx2.LT.rjx**2) THEN
73  rjx = sqrt(rjx2)
74  rdx = rdi/rjx
75  END IF
76  3 CONTINUE
77 * Bypass RDOT < 0 test for approaching ejection candidate.
78  IF (rdx.LT.0.0.AND.rjx.LT.2.0*rsum/float(nn-1)) THEN
79  WRITE (6,4) name(jclose), gpert, rij, rjx, rdx
80  4 FORMAT (' TRY ABSORB NAM PERT RIJ RJX RDX ',
81  & i6,f6.2,1p,4e9.1)
82  go to 5
83  END IF
84  END IF
85 *
86 * Include conditions for skipping (large RIJ & size or small GPERT).
87  IF (rij.GT.3.0*min(rsum,rmin)) go to 10
88  IF (rdot.GT.0.0.AND.gpert.LT.0.05) go to 10
89  IF (rsum.GT.5.0*rmin.AND.gpert.LT.1.0) go to 10
90 * Allow triple hierarchy subject to maximum membership of 6.
91  5 IF (nn.GT.3.AND.name(jclose).LT.0) go to 10
92  IF (nn.GT.4.AND.jclose.GT.n) go to 10
93 *
94 * Perform impact parameter test: a*(1 - e) < RSUM.
95  vr2 = (xdot(1,ich) - xdot(1,jclose))**2 +
96  & (xdot(2,ich) - xdot(2,jclose))**2 +
97  & (xdot(3,ich) - xdot(3,jclose))**2
98  ainv = 2.0/rij - vr2/(mass + body(jclose))
99  ecc2 = (1.0 - rij*ainv)**2 + rdot**2*ainv/(mass + body(jclose))
100  ecc = sqrt(ecc2)
101  semi1 = 1.0/ainv
102  pmin = semi1*(1.0 - ecc)
103 * Widen the impact parameter test to be on safe side.
104  IF (pmin.GT.1.5*rsum.AND.gpert.LT.0.05) go to 10
105 *
106 * Check option for increasing or decreasing regularization parameters.
107  IF (kz(16).GT.2) THEN
108  IF (gpert.GT.0.05.AND.nch.LE.4.AND.ksmag.LE.5) THEN
109  ksmag = ksmag + 1
110  rmin = 1.2*rmin
111  rmin2 = rmin**2
112  dtmin = 1.3*dtmin
113  IF (ksmag.GE.4) THEN
114  WRITE (77,700) time+toff, ksmag, gpert, rmin, rij
115  700 FORMAT (' INCREASE T KSMAG GPERT RMIN RIJ ',
116  & f9.1,i4,1p,3e10.2)
117  CALL flush(77)
118  END IF
119  ELSE IF (gpert.LT.0.001.AND.ksmag.GT.1) THEN
120  ksmax = ksmag - 1
121  rmin = 0.83*rmin
122  rmin2 = rmin**2
123  dtmin = 0.77*dtmin
124  WRITE (77,705) time+toff, ksmag, gpert, rmin, rij
125  705 FORMAT (' REDUCE T KSMAG GPERT RMIN RIJ ',
126  & f9.1,i4,1p,3e10.2)
127  CALL flush(77)
128  END IF
129  END IF
130 *
131 * Delay accepting very small binary (suppressed; eccentricity effect).
132 * IF (JCLOSE.GT.N.AND.GPERT.LT.0.25) THEN
133 * IF (100.0*R(JCLOSE-N).LT.RIJ) GO TO 10
134 * END IF
135 *
136 * Include additional criteria (V^2 > VC^2/2; JCL > N & RIJ < RSUM).
137  rdot = rdot/rij
138  vc2 = (mass + body(jclose))/rij
139  IF ((rsum + rij.LT.rmin).OR.
140  & (rij.LT.2.0*rmin.AND.pmin.LT.rmin).OR.
141  & (rij.LT.rsum.AND.rdot**2.GT.0.5*vc2).OR.
142  & (jclose.GT.n.AND.rij.LT.rsum).OR.
143  & (gpert.GT.0.2.AND.rdot.LT.0.0)) THEN
144 *
145 * Do not allow multiple absorptions (NAMES(NMX,5) = 0 each new chain).
146  IF (name(jclose).NE.names(nmx,5).AND.names(nmx,5).EQ.0) THEN
147  names(nmx,5) = name(jclose)
148 * Include possible return after ejection in bound orbit.
149  ELSE IF (rdot.GT.0.0) THEN
150  go to 10
151  ELSE IF (nn.GT.4.AND.jclose.GT.n.AND.rsum.GT.2.0*rmin) THEN
152  go to 10
153  END IF
154 *
155  rsum = rsum + rij
156  IF (kz(30).GT.1) THEN
157  IF (jclose.GT.n) THEN
158  semi = -0.5*body(jclose)/h(jclose-n)
159  ELSE
160  semi = 0.0
161  END IF
162  WRITE (6,6) nstep1, jclose, name(jclose), gpert, rij,
163  & rsum, pmin, semi
164  6 FORMAT (' ABSORB: # JCLOSE NMJ GPERT RIJ RSUM PMIN A ',
165  & 3i6,f6.2,1p,4e9.1)
166  END IF
167 *
168 * Switch off any slow-down before increasing the chain (cf. new step).
169  IF (kslow2) THEN
170  fac = 1.0
171  DO 7 k = 1,nch-1
172  fac = max(ksch(k),fac)
173  7 CONTINUE
174  CALL ycopy(y)
175  gsave = gcrit
176  gcrit = 0.0
177  CALL slow(y)
178  gcrit = gsave
179  ELSE
180  fac = 1.0
181  END IF
182 *
183 * Treat extreme binary as inert on significant slow-down.
184  IF (jclose.GT.n.AND.fac.GT.10.0) THEN
185  semi = -0.5*body(jclose)/h(jclose-n)
186  rm = 1.0
187  DO 8 k = 1,nch-1
188  rm = min(1.0/rinv(k),rm)
189  8 CONTINUE
190  small = 0.01*(rsum - rij)
191  IF (rm.LT.small.AND.semi.LT.small) THEN
192  WRITE (6,9) name(jclose), gpert, semi, rij,
193  & (1.0/rinv(k),k=1,nch-1)
194  9 FORMAT (' ABSORB INERT NM G A RIJ R ',
195  & i6,f6.2,1p,6e10.2)
196 * Switch off option #26 temporarily for routine SETSYS.
197  kz26 = kz(26)
198  kz(26) = 1
199  kz(27) = 1
200  CALL absorb(isub)
201  kz(26) = kz26
202  kz(27) = kz27
203  kcase = 1
204  go to 50
205  END IF
206  END IF
207 *
208 * Absorb the perturber (single particle or binary).
209  CALL absorb(isub)
210 *
211 * Reduce block-time since new c.m. step may be very small.
212  tblock = min(time,tblock)
213 *
214 * Activate indicator for new chain treatment and try a second search.
215  kcase = 1
216  itry = itry + 1
217  IF (itry.EQ.1) go to 1
218  END IF
219 *
220 * Exit if any particles have been absorbed.
221  10 IF (itry.GT.0) go to 50
222 *
223 * Place index of the smallest INVERSE distance in ISORT(1).
224  CALL hpsort(nn-1,rinv,isort)
225 *
226 * Determine index of escaper candidate (single star or binary).
227  kcase = 0
228  jesc = 0
229 * Distinguish two cases of each type (beginning & end of chain).
230  IF (isort(1).EQ.1) THEN
231  iesc = iname(1)
232  kcase = 1
233  ELSE IF (isort(1).EQ.nn-1) THEN
234  iesc = iname(nn)
235  kcase = 1
236 * Check for possible binary escaper (NN = 3 implies single escaper).
237  ELSE IF (isort(1).EQ.2) THEN
238  iesc = iname(1)
239  jesc = iname(2)
240  ibin = 1
241 * Switch binary index in case last separation is smallest (NN = 4).
242  IF (rinv(1).LT.rinv(nn-1)) ibin = nn - 1
243  kcase = 2
244  ELSE IF (isort(1).EQ.nn-2) THEN
245  iesc = iname(nn-1)
246  jesc = iname(nn)
247  ibin = nn - 1
248  kcase = 2
249  END IF
250 *
251 * Try escape check if middle distance is largest.
252  IF (kcase.EQ.0.AND.nn.GE.5.AND.
253  & 1.0/rinv(isort(1)).GT.2.0*rmin) THEN
254  r1 = 1.0/rinv(1)
255  r2 = 1.0/rinv(nn-1)
256 * Set relevant indices for beginning or end of chain.
257  IF (r1.GT.r2) THEN
258  iesc = iname(1)
259  jx = iname(2)
260  ib = 1
261  isort(1) = 1
262  r3 = 1.0/rinv(2)
263  ELSE
264  iesc = iname(nn)
265  jx = iname(nn-1)
266  ib = nn - 1
267  isort(1) = nn - 1
268  r3 = 1.0/rinv(nn-2)
269  END IF
270 * Define binary indices for large second separation.
271  IF (r3.GT.max(r1,r2)) THEN
272  jesc = jx
273  ibin = ib
274  kcase = 2
275  ELSE
276  kcase = 1
277  END IF
278  END IF
279 *
280 * Include safety test for abnormal configuration.
281  IF (kcase.EQ.0) go to 60
282 *
283  IF (kz(30).GT.2) THEN
284  WRITE (6,12) iesc, jesc, nstep1, isort(1),
285  & (1.0/rinv(k),k=1,nn-1)
286  12 FORMAT (' CHMOD: IESC JESC # ISORT1 R ',2i3,i6,i3,1p,5e9.1)
287  END IF
288 *
289 * Copy chain variables to standard form.
290  lk = 0
291  DO 20 l = 1,nch
292  DO 15 k = 1,3
293  lk = lk + 1
294  x4(k,l) = xch(lk)
295  xdot4(k,l) = vch(lk)
296  15 CONTINUE
297  20 CONTINUE
298 *
299 * First check case of escaping binary (JESC > 0 & RB < RJB/4).
300  IF (jesc.GT.0) THEN
301  rb = 1.0/rinv(ibin)
302  jbin = ibin + 1
303  IF (ibin.EQ.nn - 1) jbin = ibin - 1
304  rjb = 1.0/rinv(jbin)
305 * Consider removal of outermost particle instead if binary is wide.
306  IF (rb.GT.0.25*rjb) THEN
307 * Change pointer to end of chain and redefine IESC.
308  IF (isort(1).EQ.2) THEN
309  isort(1) = 1
310  ELSE
311  isort(1) = nn - 1
312  END IF
313  iesc = jesc
314  go to 30
315  END IF
316  ELSE
317  go to 30
318  END IF
319 *
320 * Form coordinates & velocities of escaping binary (local c.m. frame).
321  bcm = bodyc(iesc) + bodyc(jesc)
322  ri2 = 0.0
323  rdot = 0.0
324  vrel2 = 0.0
325  DO 25 k = 1,3
326  xcm(k) = (bodyc(iesc)*x4(k,iesc) + bodyc(jesc)*x4(k,jesc))/bcm
327  vcm(k) = (bodyc(iesc)*xdot4(k,iesc) +
328  & bodyc(jesc)*xdot4(k,jesc))/bcm
329  ri2 = ri2 + xcm(k)**2
330  rdot = rdot + xcm(k)*vcm(k)
331  vrel2 = vrel2 + (xdot4(k,iesc) - xdot4(k,jesc))**2
332  25 CONTINUE
333 *
334 * Convert to relative distance & radial velocity w.r.t. inner part.
335  fac = body(ich)/(body(ich) - bcm)
336  ri = sqrt(ri2)
337  rdot = fac*rdot/ri
338  ri = fac*ri
339 * Adopt arithmetic mean of RSUM and RMAXS for delaying escape.
340  desc = 0.5*sqrt(rsum*rmaxs(isub))
341 * Reduce delay test for large length contrasts.
342  im = isort(1)
343  rm = 1.0/rinv(im)
344  cx = rsum/(rsum - rm)
345  IF (cx.GT.25.0) desc = 25.0*rgrav
346  resc = max(3.0*rgrav,desc)
347  ainv = 2.0/rb - vrel2/bcm
348  eb = -0.5*bodyc(iesc)*bodyc(jesc)*ainv
349  hi = 0.5*rdot**2 - body(ich)/ri
350 *
351 * Employ parabolic escape criterion (terminate if RI > RMIN & NCH < 5).
352  IF ((ri.GT.resc.AND.rdot.GT.0.0.AND.rb*ainv.GT.0.99).OR.
353  & (hi.GT.0.0.AND.ri.GT.3.0*rmin)) THEN
354 * Decide between weakly bound and escape orbit (RMAXS may be large).
355  IF (rdot**2.LT.2.0*body(ich)/ri) THEN
356 * Define effective perturbation using remaining size.
357  gb = 2.0*bcm/(body(ich) - bcm)*((rsum - rjb - rb)/ri)**3
358 * Split into two KS solutions if binaries are well separated.
359  IF (nch.EQ.4.AND.gb.LT.0.001) THEN
360  kcase = -1
361  go to 50
362  END IF
363 * Enforce termination if RI > MIN(RMAXS/2,RMIN) and NCH <= 4.
364  IF (ri.GT.min(0.5*rmaxs(isub),rmin).AND.nch.LE.4) THEN
365  kcase = -1
366  go to 50
367 * Accept binary escape for small perturbation & V**2 > M/R (NCH > 4).
368  ELSE IF (gb.LT.0.01.AND.nch.GT.4.AND.
369  & (rdot**2.GT.body(ich)/ri.OR.ri.GT.rmin)) THEN
370  WRITE (6,28) iesc, jesc, namec(iesc), namec(jesc),
371  & ri, rdot**2, 2.0*body(ich)/ri, rb
372  cm(9) = cm(9) - eb
373  go to 40
374  ELSE
375 * Check splitting into two KS solutions (R1 + R3 < R2/5 & RDOT > VC).
376  IF (nch.EQ.4) THEN
377  r13 = 1.0/rinv(1) + 1.0/rinv(3)
378  vc2 = rdot**2*ri/body(ich)
379 * Ensure RSUM > 0.1*RMIN for reducing differential force corrections.
380  IF (r13.LT.0.2/rinv(2).AND.vc2.GT.1.0.AND.
381  & rsum.GT.0.1*rmin) THEN
382  kcase = -1
383  go to 50
384  END IF
385  END IF
386  kcase = 0
387  go to 60
388  END IF
389  ELSE
390  IF (hi.GT.0.0) THEN
391  vinf = sqrt(2.0*hi)*vstar
392  ELSE
393  vinf = 0.0
394  END IF
395  IF (kz(30).GT.1.OR.vinf.GT.1.0) THEN
396  WRITE (6,28) iesc, jesc, namec(iesc), namec(jesc),
397  & ri, rdot**2, 2.0*body(ich)/ri, rb, vinf
398  28 FORMAT (' CHAIN ESCAPE: IESC JESC NM RI RDOT2 ',
399  & '2*M/R RB VINF ',
400  & 2i3,2i6,1p,4e9.1,0p,f6.1)
401  END IF
402 * Enforce termination (KCASE < 0) if NCH <= 4 (final membership <= 2).
403  IF (nch.LE.4) THEN
404  kcase = -1
405  go to 50
406  END IF
407  cm(9) = cm(9) - eb
408  go to 40
409  END IF
410  ELSE
411  kcase = 0
412  go to 60
413  END IF
414 *
415 * Form relative distance and radial velocity for single particle.
416  30 ri = sqrt(x4(1,iesc)**2 + x4(2,iesc)**2 + x4(3,iesc)**2)
417  rdot = x4(1,iesc)*xdot4(1,iesc) + x4(2,iesc)*xdot4(2,iesc) +
418  & x4(3,iesc)*xdot4(3,iesc)
419  fac = body(ich)/(body(ich) - bodyc(iesc))
420  rdot = fac*rdot/ri
421  ri = fac*ri
422  im = isort(1)
423 * Ensure that escaper is not close to another body.
424  rm = min(1.0/rinv(im),ri)
425  rb = 0.0
426 *
427 * Check approximate escape criterion outside 3*RGRAV and DESC.
428  desc = 0.5*sqrt(rsum*rmaxs(isub))
429 * Reduce delay test for large length contrasts.
430  cx = rsum/(rsum - rm)
431  IF (cx.GT.25.0) desc = 25.0*rgrav
432 * Note that arithmetic mean tends to delay escape.
433  resc = max(5.0*rgrav,desc)
434  resc = min(resc,3.0*rmin)
435  IF (nn.EQ.3.AND.resc.LT.0.001*rmin) resc = 5.0*resc
436 * Facilitate termination for wide perturbed triple.
437  IF (nn.EQ.3.AND.gpert.GT.0.01.AND.ri.GT.2.0*rmin) resc = 2.0*rmin
438  IF (rm.GT.resc.AND.rdot.GT.0.0) THEN
439  IF (rdot**2.LT.2.0*body(ich)/ri) THEN
440  fac2 = 2.0*bodyc(iesc)/(body(ich) - bodyc(iesc))
441  gi = fac2*((rsum - 1.0/rinv(im))/ri)**3
442 * Do not permit large length contrast even for bound orbit (NN = 4).
443  IF (nn.EQ.4.AND.gi.LT.0.01) THEN
444  ib = isort(nn-1)
445 * Prevent NN = 4 termination near small pericentre (subject to safety).
446  IF (1.0/rinv(ib).GT.0.1*rgrav.OR.rsum.GT.rmin) THEN
447  kcase = -1
448  go to 50
449  ELSE
450  kcase = 0
451  go to 60
452  END IF
453  END IF
454  er = 0.5*rdot**2 - body(ich)/ri
455  rx = -body(ich)/er
456 * Deny continuation of soft configuration with large external pert.
457  IF ((er.LT.0.0.AND.rx.LT.min(2.0*rsum,2.0*rmin)).OR.
458  & (rsum.LT.5.0*rmin.AND.gpert.LT.0.01)) THEN
459  kcase = 0
460  go to 60
461  END IF
462  IF (ri.GT.min(0.5*rmaxs(isub),rmin)) THEN
463 * Decide between termination or removal for sub-parabolic motion.
464  IF (nn.LT.4) THEN
465 * Delay termination for eccentric binary inside R = -m1*m2/2*E < a.
466  ib = 3 - im
467  j1 = iname(ib)
468  j2 = iname(ib+1)
469 * Obtain semi-major axis from ENERGY = -(m1*m2 + (m1+m2)*m3)/RGRAV.
470  zmu = bodyc(j1)*bodyc(j2)/(bodyc(j1) + bodyc(j2))
471  zf = bodyc(iesc)/zmu
472 * Write E = E_b + E_1 with E_b = -m1*m2/2*a and E_1 = MU3*ER.
473  semi = 0.5*rgrav/(1.0 + zf*(1.0 + rgrav*er/mass))
474  ry = 1.0/rinv(ib)
475  IF (ry.LT.0.9*semi.AND.ri.LT.2.0*rmin) THEN
476  WRITE (6,31) nstep1, ry/semi, ri, rdot**2,
477  & 2.0*body(ich)/ri, semi
478  kcase = 0
479  go to 60
480  END IF
481  kcase = -1
482  go to 60
483  ELSE
484  kcase = 1
485  go to 32
486  END IF
487  ELSE
488  kcase = 0
489  go to 60
490  END IF
491  ELSE IF (nn.GE.4) THEN
492 * Accept escape of hyperbolic body if first or last separation > RMIN.
493  rbig = 1.0/rinv(im)
494  IF (rbig.GT.rmin) THEN
495  IF (im.EQ.1.OR.im.EQ.nn-1) THEN
496  kcase = 1
497  go to 32
498  END IF
499  END IF
500 * Include three-body stability test for distant fourth body.
501  IF (rbig.GT.0.75*rsum.AND.nn.EQ.4.AND.
502  & rsum.GT.0.5*rmin) THEN
503  CALL chstab(iterm)
504  IF (iterm.LT.0) THEN
505  kcase = -1
506  go to 60
507  END IF
508  END IF
509 * Skip if largest separation is not dominant.
510  IF (rbig.LT.0.66*rsum) THEN
511  kcase = 0
512  go to 60
513  END IF
514  ELSE
515 * Note possibility a << RGRAV/6 after strong interaction with NN = 3.
516  ib = 3 - im
517  j1 = iname(ib)
518  j2 = iname(ib+1)
519 * Obtain semi-major axis from ENERGY = -(m1*m2 + (m1+m2)*m3)/RGRAV.
520  zmu = bodyc(j1)*bodyc(j2)/(bodyc(j1) + bodyc(j2))
521  zf = bodyc(iesc)/zmu
522  er = 0.5*rdot**2 - body(ich)/ri
523 * Write E = E_b + E_1 with E_b = -m1*m2/2*a and E_1 = MU3*ER.
524  semi = 0.5*rgrav/(1.0 + zf*(1.0 + rgrav*er/mass))
525  ry = 1.0/rinv(ib)
526 * Note RGRAV = MIN(0.5*RSUM,RGRAV) is OK and RGRAV adjusted on QPMOD.
527  IF (ry.LT.0.9*semi.AND.ri.LT.2.0*rmin.AND.
528  & ri.LT.50.0*semi) THEN
529  WRITE (6,31) nstep1, ry/semi, ri, rdot**2,
530  & 2.0*body(ich)/ri, semi
531  31 FORMAT (' CHAIN DELAY # R/A RI RD2 VP2 A ',
532  & i5,f6.2,1p,4e9.1)
533  kcase = 0
534  go to 60
535  END IF
536  END IF
537 *
538 * Check that escaper is well separated (ratio > 3).
539  rr = rsum - 1.0/rinv(im)
540  ratio = 1.0/(rinv(im)*rr)
541  IF (ratio.LT.3.0.AND.rsum.LT.rmin) THEN
542  kcase = 0
543  go to 60
544  END IF
545 *
546  32 hi = 0.5*rdot**2 - body(ich)/ri
547  IF (hi.GT.0.0) THEN
548  vinf = sqrt(2.0*hi)*vstar
549  ELSE
550  vinf = 0.0
551  END IF
552  IF (kz(30).GT.1.OR.vinf.GT.2.0) THEN
553  WRITE (6,35) iesc, namec(iesc), ri, rdot**2,
554  & 2.0*body(ich)/ri, vinf
555  35 FORMAT (' CHAIN ESCAPE: IESC NM RI RDOT2 2*M/R VF ',
556  & i3,i6,1p,3e9.1,0p,f6.1)
557  END IF
558 * Ensure single body is removed in case of wide binary.
559  jesc = 0
560  ELSE
561  kcase = 0
562  go to 60
563  END IF
564 *
565 * Reduce chain membership (NCH > 3) or specify termination.
566  40 IF (nch.GT.3) THEN
567 * Subtract largest chain distance from system size (also binary).
568  im = isort(1)
569  rsum = rsum - 1.0/rinv(im) - rb
570  CALL reduce(iesc,jesc,isub)
571  ELSE
572  kcase = -1
573  END IF
574 *
575 * Set phase indicator < 0 to ensure new time-step list in INTGRT.
576  50 iphase = -1
577 *
578  60 RETURN
579 *
580  END