Nbody6
 All Files Functions Variables
mdot.f
Go to the documentation of this file.
1  SUBROUTINE mdot
2 *
3 *
4 * Mass loss from evolving stars.
5 * Updated 6/1/98 by J. Hurley
6 * ------------------------------
7 *
8  include 'common6.h'
9  common/binary/ bm(4,mmax),xrel(3,mmax),vrel(3,mmax),
10  & hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
11  & namem(mmax),nameg(mmax),kstarm(mmax),iflagm(mmax)
12  common/modes/ eb0(ntmax),zj0(ntmax),ecrit(ntmax),ar(ntmax),
13  & br(ntmax),eosc(4,ntmax),edec(ntmax),tosc(ntmax),
14  & rp(ntmax),es(ntmax),cm(2,ntmax),iosc(ntmax),
15  & namec(ntmax)
16  INTEGER jx(2),kacc
17  REAL*8 mass(2),massc(2),rad(2),radc(2),lumin(2)
18  REAL*8 age0(2),tm0(2),tbgb(2),menv(2),renv(2),k2str(2)
19  REAL*8 tscls(20),lums(10),gb(10),tm,tn
20  REAL*8 m0,m1,rm,age,lum,mc,rcc,me,re,rnew,lnew,mch
21  parameter(mch=1.44d0)
22  REAL*8 m10,rxl1,rxl2,q,rol(2),rlperi,rm0
23  REAL*8 eps,alpha2,tol
24  parameter(eps=1.0d-06,alpha2=0.09d0,tol=1.0d-10)
25  REAL*8 dtx(2),dms(2),dma(2),dmx(2),epch0(2),tevk,diff
26  REAL*8 dt,dsep,djmb,djgr,djorb,djt,djti,dtgr,dspin,jspbru,ospbru
27  REAL*8 rj,hj,tj2,vorb2,vwind2,ivsqm,omv2,dtxmin,dmxmax,delet
28  REAL*8 jspin(2),ospin(2),djspin(2),jorb,oorb,aursun
29  parameter(aursun=214.95d0)
30  REAL*8 k2,k3,acc1,acc2,beta,xi
31  parameter(k3=0.21d0,acc1=3.920659d+08,acc2=1.5d0)
32  parameter(beta=0.125d0,xi=1.d0,maxblk=40)
33  REAL*8 mlwind,rl
34  EXTERNAL mlwind,rl
35  CHARACTER*8 which1
36  LOGICAL icorr,ikick
37  SAVE iwarn
38  DATA iwarn /0/
39 *
40 * Define current time (unit of million years) and update counter.
41  ttot = time + toff
42  tphys = ttot*tstar
43  ipoly = 0
44  time0 = time
45  iblue = 0
46 *
47 * Find global index of next star to be checked (NB! reset KS & IKS).
48  1 ks = 0
49  kspair = 0
50  iks = 0
51  itry = 0
52  ighost = 0
53  i0 = 0
54  ihdot = 0
55  IF (iphase.LT.0) ipoly = -1
56  iphase = 0
57  iqcoll = 0
58 * Restore synchronized TIME in case of RESET, CHRECT & KSPERI.
59  time = time0
60 *
61 * Determine smallest look-up time (< TIME!).
62  DO 5 j = 1,ntot
63  IF (tev(j).LE.time) THEN
64  i = j
65  go to 10
66  END IF
67  5 CONTINUE
68 *
69 * Determine new evolution time (current TMDOT may be escaped star).
70  kw = 0
71  i = 1
72  go to 70
73 *
74 * Include braking or merger termination at time for Roche overflow.
75  10 IF (i.GT.n.AND.name(i).LT.0.AND.kz(34).GT.0) THEN
76  kspair = i - n
77 * Check spin synchronization of inner binary for small eccentricity.
78  CALL spinup(kspair,iterm)
79  IF (iterm.EQ.0) THEN
80  CALL brake2(kspair,iterm)
81  END IF
82  IF(iterm.GT.0)THEN
83  nwarn = nwarn + 1
84  IF (nwarn.LT.1000) THEN
85  WRITE(14,900) kspair, name(i), ttot, r(kspair)
86  900 FORMAT(' MDOT MERGER TERM: KS NAM T R ',
87  & i4,i6,f10.4,1p,e10.2)
88  CALL flush(14)
89  END IF
90  jcomp = 0
91  iphase = 7
92  CALL reset
93  ENDIF
94  go to 1
95  END IF
96 *
97 * Check for hierarchy containing NS/BH (same TEV for inner components).
98  IF (i.LT.ifirst.AND.kz(28).GT.0) THEN
99  ip = kvec(i)
100  IF (name(n+ip).LT.0.AND.kstar(i).GE.13) THEN
101 * Distinguish between inner binary components and outer component(s).
102  IF (i.EQ.2*ip-1) THEN
103  CALL brake3(ip,iterm)
104 * Update single outer KS member or terminate rare case of binary.
105  ELSE IF (name(2*ip).LE.nzero) THEN
106  tev0(i) = time
107  tev(i) = time + 100.0
108  iterm = 0
109  ELSE
110  iterm = 1
111  END IF
112  IF (iterm.GT.0) THEN
113  iphase = 7
114  kspair = ip
115  CALL reset
116  END IF
117  go to 1
118  END IF
119  END IF
120 *
121 * Check possible Roche overflow condition (skip any ghost c.m.).
122  IF (i.GT.n.AND.kz(34).GT.0.AND.body(i).GT.0.0d0.AND.
123  & kstar(i).GE.10) THEN
124  ipair = i - n
125 * Consider odd types (11, 13, etc.) for Roche updating or treatment.
126  IF (mod(kstar(i),2).NE.0) THEN
127  CALL roche(ipair)
128  IF (iqcoll.NE.0.OR.iphase.LT.0) go to 1
129 * Ensure that look-up time is increased in case of no mass transfer.
130  IF (tev(i).LE.time) THEN
131  IF (mod(kstar(i),2).NE.0) kstar(i) = kstar(i) + 1
132  CALL trflow(ipair,dtr)
133  tev(i) = time + dtr
134  WRITE(6,902)i,name(i),ipair,kstar(i),tev(i),ttot
135  902 FORMAT(' WARNING MDOT NO ROCHE ',3i6,i4,1x,1p,2e13.4)
136  ENDIF
137  ELSE
138 * Include other cases (even types outside 10 R_sun).
139  CALL trflow(ipair,dtr)
140  IF (dtr.LT.step(i)) THEN
141  tev(i) = time + dtr
142  CALL roche(ipair)
143  IF (tev(i).LE.time) tev(i) = 1.000002d0*time
144  ELSE
145  IF(kz(34).EQ.1)THEN
146  CALL synch(ipair)
147  ELSE
148  tev(i) = time + dtr
149  ENDIF
150  END IF
151  END IF
152  go to 1
153  END IF
154 *
155 * Check for special treatment of possible hierarchical system (IHDOT).
156  IF (i.LT.ifirst) THEN
157  kspair = kvec(i)
158  ksx = max(kstar(2*kspair-1),kstar(2*kspair))
159  icm = n + kspair
160  i2 = 2*kspair
161  IF(i2.EQ.i) i2 = i2 - 1
162 *
163 * Select merger if #I is first or the second member is a c.m. body.
164  IF (name(icm).LT.0.AND.
165  & (i.EQ.2*kspair - 1.OR.name(2*kspair).GT.nzero)) THEN
166 * Activate indicator (distinguish inner and outer binary).
167  ihdot = 1
168 * Exclude inner components and outer binary of double hierarchy.
169  IF (name(icm).LT.-2*nzero) THEN
170  IF (i.LT.2*kspair.OR.name(i).GT.nzero) THEN
171  iphase = 7
172  CALL reset
173  go to 1
174  END IF
175  END IF
176  END IF
177 *
178 * Include Roche consistency check (single components not permitted).
179  IF (kstar(icm).GT.10.AND.mod(kstar(icm),2).EQ.1) THEN
180  IF(body(i).EQ.0.0.OR.body(i2).EQ.0.0)THEN
181  tev(i) = time + min(0.1d0,tev(i)-tev0(i))
182  tev(i2) = tev(i)
183  goto 1
184  ENDIF
185  iwarn = iwarn + 1
186  IF (mod(iwarn,100).EQ.0) THEN
187  WRITE(6,904) i, kstar(icm), tev(i) - tev(icm)
188  904 FORMAT(' WARNING! MDOT I K* D(TEV) ',2i5,1p,e9.1)
189  END IF
190  CALL trflow(kspair,dtr)
191  tev(icm) = time + dtr
192  IF (dtr.GT.tstar) THEN
193  kstar(icm) = kstar(icm) + 1
194  ELSE
195  tev(i) = 1.000002d0*time
196  go to 1
197  END IF
198  END IF
199  END IF
200 *
201 * Determine relevant indices in case of merger (zero mass done below).
202  IF (ihdot.GT.0.AND.body(i).GT.0.0d0) THEN
203 *
204 * Identify ghost address and merger index from c.m. of KSPAIR.
205  CALL findj(i,ighost,imerge)
206 *
207  IF(ighost.LT.0)THEN
208  WRITE(6,*)' MDOT GHOST NOT FOUND ',i,name(i),imerge
209  stop
210  ENDIF
211 *
212 * Determine index of mass-losing star for inner or outer binary.
213  IF (ighost.LE.n) THEN
214  IF (tev(ighost).LT.tev(i)) i = ighost
215  IF (i.EQ.ighost) WRITE(6,906) i, name(i), kstar(i),
216  & tev(i), body(i)
217  906 FORMAT(' GHOST SWITCH: I NAM K* TEV BODY ',
218  & 3i5,f9.2,1p,e10.2)
219  ELSE
220 * Save KS component for updating look-up time.
221  i0 = i
222 * Compare look-up times of #I and IGHOST <= N on first KS component.
223  IF (i.LT.2*kspair) THEN
224  IF (tev(ighost).LT.tev(i)) i = ighost
225  ELSE
226 * Form ghost pair index and select KS component by comparing TEV.
227  jpair = ighost - n
228  i = 2*jpair - 1
229  IF (tev(2*jpair).LT.tev(i)) i = 2*jpair
230  ihdot = 2
231  WRITE(6,908) ighost, i0, i, name(i), kstar(i), tev(i)
232  908 FORMAT(' OUTER GHOST: IG I0 I NM K* TEV',5i6,f9.2)
233  IF(kstar(n+kspair).GE.10) THEN
234  CALL reset
235  goto 1
236  ENDIF
237  END IF
238  END IF
239  ELSE IF (body(i).LE.0.0d0.OR.name(i).EQ.0) THEN
240 * Distinguish between merger ghost and member of CHAIN/TRIPLE/QUAD.
241  imerge = 0
242  icm = 0
243  j = i
244 * Replace #I by corresponding ghost c.m. for zero-mass KS component.
245  IF (i.LT.ifirst) j = n + kvec(i)
246 *
247 * Identify merger index and determine corresponding c.m. body.
248  DO 15 k = 1,nmerge
249  IF (nameg(k).EQ.name(j)) imerge = k
250  15 CONTINUE
251  IF (imerge.GT.0) THEN
252  DO 20 k = n+1,ntot
253  IF (namem(imerge).EQ.name(k)) icm = k
254  20 CONTINUE
255  kspair = icm - n
256  if(icm.eq.0) kspair = j - n
257 * Terminate if ghost belongs to double hierarchy (bug fix 03/2001).
258  IF (namem(imerge).LT.-2*nzero) THEN
259  iphase = 7
260  CALL reset
261  go to 1
262  END IF
263  END IF
264 *
265 * Define KS index and set indicator on successful identification.
266  IF (icm.GT.n) THEN
267  ihdot = 1
268  ighost = j
269  IF (ighost.GT.n) THEN
270  ihdot = 2
271  jpair = ighost - n
272  i = 2*jpair - 1
273  IF (tev(2*jpair).LT.tev(i)) i = 2*jpair
274  i0 = ighost
275  END IF
276  ELSE
277 * Extend TEV for other single ghosts or c.m. of compact subsystem.
278  tev(i) = tev(i) + 0.01d0
279 * Allow an extra amount for MS stars (depending on TEV - TEV0).
280  IF(kstar(i).LE.1.AND.(tev(i)-tev0(i))*tstar.LT.10.0)THEN
281  tev(i) = tev(i) + 0.5*(tev(i) - tev0(i))
282  END IF
283  go to 1
284  END IF
285  END IF
286 *
287 * Skip any c.m. particles after increasing look-up time.
288  IF (i.GT.n) THEN
289  i1 = 2*(i - n) - 1
290  i2 = i1 + 1
291  tm = min(tev(i1),tev(i2))
292 * IF(KSTAR(I).GE.10)THEN
293 * IWARN = IWARN + 1
294 * IF (MOD(IWARN,100).EQ.0) THEN
295 * WRITE(6,910) I, NAME(I), KSTAR(I1), KSTAR(I2), KSTAR(I),
296 * & TTOT, TM - TEV(I)
297 *910 FORMAT(' WARNING! MDOT: I NAM K* T D(TEV) ',
298 * & 2I6,3I4,F9.2,1P,E10.2)
299 * END IF
300 * ENDIF
301  IF (kstar(i).LT.10) THEN
302  tev(i) = 1.0d+10
303  ELSE IF (kz(34).GT.0) THEN
304  ipair = i - n
305  CALL trflow(ipair,dtr)
306  dtr = max(dtr,0.0d0)
307  IF (dtr.LT.step(i)) THEN
308  tev(i) = time + dtr
309  CALL roche(ipair)
310  END IF
311  tev(i) = max(time + dtr,1.000002d0*time)
312  ELSE
313  tev(i) = 1.000002d0*time
314  END IF
315  go to 1
316  END IF
317 *
318 * Copy relevant mass (standard case or merger ghost member).
319 * Also determine if the star is in a standard binary or is the inner
320 * component of a hierachy and if so set the mass accretion flag.
321  kacc = 0
322  jx(1) = i
323  ecc2 = 0.d0
324  rj = 0.d0
325  hj = 0.d0
326  tj2 = 0.d0
327  IF(ighost.EQ.0.AND.ihdot.EQ.0)THEN
328  mass(1) = body(i)*zmbar
329  IF(i.LT.ifirst)THEN
330  IF(name(icm).GT.0)THEN
331  kacc = 1
332  jx(2) = i2
333  mass(2) = body(i2)*zmbar
334  rj = r(kspair)
335  hj = h(kspair)
336  tj2 = tdot2(kspair)
337  ENDIF
338  ENDIF
339  ELSE
340  IF(ihdot.EQ.1)THEN
341  k = 1
342  IF (i.GE.ifirst) k = 2
343  IF(i.NE.2*kspair)THEN
344  kacc = 2
345  j = 3 - k
346  jx(2) = ighost
347  IF(k.EQ.2) jx(2) = 2*kspair - 1
348  mass(2) = bm(j,imerge)*zmbar
349  DO 25 ii = 1,4
350  rj = rj + um(ii,imerge)*um(ii,imerge)
351  tj2 = tj2 + 2.d0*um(ii,imerge)*umdot(ii,imerge)
352  25 CONTINUE
353  hj = hm(imerge)
354  ENDIF
355  ELSE
356  k = 3
357  IF (i.EQ.2*jpair) k = 4
358  END IF
359  mass(1) = bm(k,imerge)*zmbar
360  END IF
361  IF(kacc.GT.0)THEN
362  j = n + kspair
363 * Check for zero c.m. mass (i.e. NAME < 0 & > -NZERO in [[B,S],[B,S]]).
364  IF (body(j).EQ.0.0d0) THEN
365  iphase = 7
366  CALL reset
367  go to 1
368  END IF
369  semi = -0.5d0*body(j)/hj
370  sep = semi*su
371  ecc2 = (1.d0 - rj/semi)**2 + tj2**2/(body(j)*semi)
372  IF(ecc2.LT.1.0)THEN
373  kacc = 2
374  q = mass(1)/mass(2)
375  rol(1) = rl(q)*sep
376  rxl1 = radius(jx(1))/rol(1)
377  q = 1.d0/q
378  rol(2) = rl(q)*sep
379  rxl2 = radius(jx(2))/rol(2)
380  iswap = 0
381  IF(kstar(jx(1)).LE.1.AND.
382  & (kstar(jx(2)).EQ.5.OR.kstar(jx(2)).EQ.6))THEN
383  iswap = 1
384  ELSEIF(rxl2.GT.rxl1)THEN
385  iswap = 1
386  ENDIF
387  IF(iswap.GT.0)THEN
388  j1 = jx(1)
389  jx(1) = jx(2)
390  jx(2) = j1
391  m10 = mass(1)
392  mass(1) = mass(2)
393  mass(2) = m10
394  rxl1 = rol(1)
395  rol(1) = rol(2)
396  rol(2) = rxl1
397  ENDIF
398  vorb2 = acc1*(mass(1)+mass(2))/sep
399  ivsqm = 1.d0/sqrt(1.d0-ecc2)
400  diff = abs(tev0(jx(2)) - tev0(jx(1)))
401  IF(diff.GT.tol)THEN
402  tev(jx(1)) = max(tev0(jx(1)),tev0(jx(2)))
403  ELSE
404  tev0(jx(2)) = tev0(jx(1))
405  IF(tev(jx(2)).NE.tev(jx(1)))THEN
406  tev(jx(1)) = min(tev(jx(1)),tev(jx(2)))
407  ENDIF
408  ENDIF
409  tev(jx(2)) = tev(jx(1))
410  diff = abs(tev0(jx(2)) - tev0(jx(1)))
411  ecc = sqrt(ecc2)
412  oorb = twopi*sqrt((mass(1)+mass(2))/(sep/aursun)**3)
413  jorb = mass(1)*mass(2)/(mass(1)+mass(2))
414  & *sqrt(1.d0-ecc2)*sep*sep*oorb
415  ELSE
416  WRITE(6,912)jx,name(n+kspair),sqrt(ecc2),ttot
417  912 FORMAT(' MDOT ECC2 > 1.0 I1 I2 N(ICM) e T ', 3i6,f6.2,e9.2)
418  kacc = 1
419  ENDIF
420  ELSE
421  kacc = 1
422  ENDIF
423 *
424  nmdot = nmdot + 1
425 *
426  DO 200 k = 1,kacc
427 *
428  i = jx(k)
429 * Set interval since last mass update.
430  dtx(k) = 1.0d+06*(tev(i) - tev0(i))*tstar
431 * Set the initial mass and current type.
432  m0 = body0(i)*zmbar
433  m1 = mass(k)
434  mc = 0.d0
435  kw = kstar(i)
436 *
437 * Obtain stellar parameters at previous epoch.
438  age = tev0(i)*tstar - epoch(i)
439  age = max(age,0.d0)
440  age0(k) = age
441  epch0(k) = epoch(i)
442  CALL star(kw,m0,m1,tm,tn,tscls,lums,gb,zpars)
443  CALL hrdiag(m0,age,m1,tm,tn,tscls,lums,gb,zpars,
444  & rm0,lum,kw,mc,rcc,me,re,k2)
445  rad(k) = rm0
446  lumin(k) = lum
447  tm0(k) = tm
448  tbgb(k) = tscls(1)
449  massc(k) = mc
450  radc(k) = rcc
451  menv(k) = me
452  renv(k) = re
453  k2str(k) = k2
454 *
455 * Ensure that type change occurs at time TEV.
456  IF(kw.NE.kstar(i))THEN
457  IF (kw.GE.13) THEN
458  WRITE (6,190) i, name(i), kstar(i), kw, mass(k)
459  190 FORMAT (' NS/BH FORMATION ',2i7,2i4,f7.2)
460  END IF
461  kw = kstar(i)
462  m1 = mass(k)
463  ENDIF
464 *
465 * Evaluate mass loss due to stellar wind.
466  rlperi = 0.d0
467  IF(kacc.GT.1) rlperi = rol(k)*(1.d0 - ecc)
468  dmx(k) = mlwind(kw,lum,rm0,m1,mc,rlperi,zmet)
469  dma(k) = 0.d0
470 *
471  IF(kacc.GT.1)THEN
472 *
473 * Calculate how much of wind mass loss will be accreted by the
474 * companion (Boffin & Jorissen, A&A 1988, 205, 155).
475  vwind2 = 2.d0*beta*acc1*m1/rm0
476  omv2 = (1.d0 + vorb2/vwind2)**(3.d0/2.d0)
477  dma(3-k) = ivsqm*acc2*dmx(k)*((acc1*mass(3-k)/vwind2)**2)
478  & /(2.d0*sep*sep*omv2)
479  dma(3-k) = min(dma(3-k),0.8d0*dmx(k))
480 *
481  ENDIF
482 *
483 * Convert the spin angular momentum of the star to physical units.
484  jspin(k) = max(spin(i)*spnfac,1.0d-10)
485 *
486 * Evaluate the spin of the star.
487  ospin(k) = jspin(k)/(k2*rm0*rm0*(m1-mc)+k3*rcc*rcc*mc)
488 *
489  200 CONTINUE
490 *
491  dt = dtx(1)
492 *
493  DO 210 k = 1,kacc
494 *
495  i = jx(k)
496  kw = kstar(i)
497 *
498 * Check for unequal evolution times in a recently formed binary.
499  IF(k.EQ.2)THEN
500  IF(dtx(1).EQ.0.d0.OR.dtx(2).EQ.0.d0) dt = dtx(2)
501  ENDIF
502 *
503 * Calculate the change in spin angular momentum due to mass loss and
504 * also due to mass accretion and/or tides if the star is in a binary.
505  djspin(k) = (2.d0/3.d0)*dmx(k)*rad(k)*rad(k)*ospin(k)
506 * Restrict the time-step for extreme conditions (JH 05/09).
507  IF(djspin(k).GT.0.0d0)THEN
508  dt = min(dt,0.3d0*jspin(k)/djspin(k))
509  ENDIF
510 *
511  IF(kacc.GT.1)THEN
512 *
513  IF(dma(k).GT.0.d0)THEN
514  djspin(k) = djspin(k) - (2.d0/3.d0)*xi*dma(k)*
515  & rad(3-k)*rad(3-k)*ospin(3-k)
516 *
517 * Produce diagnostics for symbiotic stars.
518  IF(dma(k).GE.3.16e-9.AND.kw.GE.10)THEN
519  WRITE(20,913)name(i),name(jx(3-k)),kstar(i),
520  & kstar(jx(3-k)),tphys,mass,sep,dmx(3-k),dma(k)
521  913 FORMAT(' WINDAC NAM K* T M A MLOSS MACC ',
522  & 2i6,2i4,f9.2,2f6.2,f8.1,1p,2e9.1)
523  CALL flush(20)
524  ENDIF
525  ENDIF
526 *
527 * Tidal synchronisation (only operates on circular orbits).
528  q = mass(3-k)/mass(k)
529  IF(k.EQ.1) djt = 0.d0
530  IF(kz(34).EQ.2.AND.dt.GT.0.d0.AND.ecc.LE.0.01d0)THEN
531  IF((kw.LE.9.AND.rad(k).GE.0.01d0*rol(k)).OR.
532  & (kw.GE.10.AND.q.GE.1.d0))THEN
533 *
534  CALL bsetid(kw,mass(k),massc(k),menv(k),rad(k),
535  & radc(k),renv(k),lumin(k),ospin(k),k2str(k),
536  & q,sep,ecc,oorb,delet,dspin,eqspin,djti)
537 *
538  IF(abs(dspin).GT.tiny)THEN
539  dspin0 = dspin
540  IF(dspin.GE.0.d0)THEN
541  dspin = min(dt*dspin,oorb-ospin(k))/dt
542  ELSE
543  dspin = max(dt*dspin,oorb-ospin(k))/dt
544  ENDIF
545  djti = djti*dspin/dspin0
546  ELSE
547  djti = 0.d0
548  ENDIF
549  djspin(k) = djspin(k) - djti
550  djt = djt + djti
551 *
552  ENDIF
553  ENDIF
554  ENDIF
555 *
556 * Include magnetic braking for stars with convective envelopes.
557  CALL magbrk(kw,mass(k),menv(k),rad(k),ospin(k),djmb)
558 * Limit to a 3% angular momentum change for the star owing to MB.
559  IF(djmb.GT.0.0d0)THEN
560  dt = min(dt,0.03d0*jspin(k)/djmb)
561  djspin(k) = djspin(k) + djmb
562  ENDIF
563 *
564 * Evaluate the mass loss from the star in the interval DT.
565  dms(k) = (dmx(k) - dma(k))*dt
566  dmr = abs(dms(k)/(mass(k) + 1.0d-10))
567 *
568 * Restrict accumulated mass loss to maximum of 2%.
569  IF(dmr.GT.0.02)THEN
570  dt = dt*0.02/dmr
571  dms(k) = 0.02*mass(k)
572  ENDIF
573 *
574 * Check that mass loss does not exceed the envelope mass.
575  IF(kstar(i).LT.10)THEN
576  dml = max(mass(k) - massc(k),1.0d-07)
577  IF(dml.LT.dms(k))THEN
578  dt = (dml/dms(k))*dt
579  dms(k) = dml
580  ENDIF
581  ENDIF
582  dtx(k) = dt
583 *
584  210 CONTINUE
585 *
586 * Include angular momentum loss owing to mass loss and/or
587 * gravitational radiation for close binary systems.
588  dsep = 0.d0
589  djorb = 0.d0
590  dtgr = 2.0d+10
591  IF(kacc.EQ.2)THEN
592 *
593  dtxmin = min(dtx(1),dtx(2))
594  CALL grrad(mass(1),mass(2),sep,ecc,jorb,djgr,delet)
595  djorb = djt + djgr
596 * write(*,*)' grrad ',sep,ecc,djgr,delet,dtxmin
597 *
598 * Limit orbital angular momentum change to 2%.
599 *
600  IF(abs(djorb).GT.tiny.AND.(dtxmin.GT.0.d0.OR.diff.EQ.0.d0))THEN
601  dtgr = 0.02d0*jorb/abs(djorb)
602  IF(delet.GT.tiny.AND.ecc.GT.0.0011d0)THEN
603  dtgr = min(dtgr,0.05d0*ecc/delet)
604  ENDIF
605  dtgr = max(dtgr,100.d0)
606 * DTGR = MAX(DTGR,1.0D-04)
607  dtxmin = min(dtxmin,dtgr)
608  dtgr = dtgr/1.0d+06
609  djorb = djorb*dtxmin
610  jorb = jorb - djorb
611  jorb = max(jorb,1.d0)
612  ecc0 = ecc
613  ecc = max(ecc0 - delet*dtxmin,0.001d0)
614  ecc2 = ecc*ecc
615  sep1 = (mass(1) + mass(2))*jorb*jorb/
616  & ((mass(1)*mass(2)*twopi)**2*aursun**3*(1.d0-ecc2))
617  dsep = sep - sep1
618  IF(dsep.GT.0.d0)THEN
619  q = mass(1)/mass(2)
620  rxl1 = 0.9d0*rad(1)/rl(q)
621  sep1 = max(sep1,rxl1)
622  dsep = sep - sep1
623  ENDIF
624  dtx(1) = dtxmin
625  dtx(2) = dtx(1)
626  ELSE
627  djorb = 0.d0
628  ENDIF
629 *
630 * Orbital changes owing to mass loss dealt with in hcorr for now.
631  dmxmax = max(dmx(1),dmx(2))
632 * IF(DMXMAX.GT.0.D0)THEN
633 * Q = MASS(1)/MASS(2)
634 * DJORB = DJORB + (DMX(1)+Q*DMA(1))*MASS(2)*MASS(2)*DTX(1)
635 * Q = 1.D0/Q
636 * DJORB = DJORB + (DMX(2)+Q*DMA(2))*MASS(1)*MASS(1)*DTX(2)
637 * DJORB = DJORB*SEP*SEP*OORB/(MASS(1)+MASS(2))**2
638 * ENDIF
639 *
640  ENDIF
641 *
642  DO 220 k = 1,kacc
643 * Set the initial mass and current type.
644  i = jx(k)
645  dt = dtx(k)
646  m0 = body0(i)*zmbar
647  m10 = m0
648  m1 = mass(k)
649  mc = massc(k)
650  kw = kstar(i)
651 * Set indicator for mass loss correction.
652  icorr = .false.
653  dms(k) = (dmx(k) - dma(k))*dt
654  dmr = abs(dms(k)/(mass(k) + 1.0d-10))
655  IF(dmr.GT.tiny)THEN
656  icorr = .true.
657  m1 = m1 - dms(k)
658 * Check rejuvenation of MS, HG or HE star.
659  IF(kw.LE.2.OR.kw.EQ.7)THEN
660  m0 = m1
661  CALL star(kw,m0,m1,tm,tn,tscls,lums,gb,zpars)
662  IF(kw.EQ.2)THEN
663  IF(gb(9).LT.massc(k).OR.m10.GT.zpars(3))THEN
664  m0 = m10
665  ELSE
666  epch0(k) = tm + (tscls(1)-tm)*(age0(k)-tm0(k))/
667  & (tbgb(k)-tm0(k))
668  epch0(k) = tev0(i)*tstar - epch0(k)
669  ENDIF
670  ELSE
671  epch0(k) = tev0(i)*tstar - age0(k)*tm/tm0(k)
672  ENDIF
673  ENDIF
674  ENDIF
675  djspin(k) = djspin(k)*dt
676 *
677 * Check for blue straggler formation (TM < TPHYS).
678  IF(dmr.GT.tiny.AND.kw.LE.1.AND.name(i).NE.iblue)THEN
679  tphys2 = (tev0(i)+toff)*tstar - epoch0
680  IF(tm0(k).GT.tphys2.AND.tm.LT.0.98d0*tphys2)THEN
681  j = jx(1)
682  IF(i.EQ.j.AND.kacc.EQ.2) j = jx(2)
683  WRITE(6,914)name(i),m1,tm,tphys2,epch0(k)+toff*tstar,
684  & kstar(j)
685  914 FORMAT(' NEW BS (MDOT): NAM M TM TP EP KW1 ',
686  & i6,f6.2,3f8.1,i4)
687  iblue = name(i)
688  nbs = nbs + 1
689  ENDIF
690  ENDIF
691 *
692 * Set current age to actual look-up value (allows looping).
693  tevk = tev0(i) + dt/(1.0d+06*tstar)
694 * Set indicator for mass loss correction.
695  age = tevk*tstar - epch0(k)
696  age0(k) = age
697 *
698 * Determine stellar evolution time scales and luminosity.
699  CALL star(kw,m0,m1,tm,tn,tscls,lums,gb,zpars)
700 *
701 * Include temporary error check.
702  IF((age-tn).GT.1.0d-02)THEN
703 * IF(I.GE.IFIRST.OR.(I.LT.IFIRST.AND.NAME(ICM).GT.0))THEN
704  WRITE(6,994)name(i), kw, dms(k), age, tn
705 * ENDIF
706  994 FORMAT(' MDOT WARNING! AGE > TN NM KW DMS AGE TN ',
707  & i6,i4,f7.3,1p,2e9.1)
708  IF(kw.LE.6)THEN
709  age = min(age,0.9999d0*tscls(11))
710  ELSE
711  age = min(age,0.9999d0*tscls(5))
712  age = min(age,1.0001d0*tn)
713  ENDIF
714  ENDIF
715 *##
716 * Obtain stellar parameters at current epoch.
717  CALL hrdiag(m0,age,m1,tm,tn,tscls,lums,gb,zpars,
718  & rm,lum,kw,mc,rcc,me,re,k2)
719 *
720 * Check for change of CO WD to ONe on rapid accretion (21/11/08).
721  IF(dmr.GT.tiny.AND.kw.EQ.11)THEN
722  dme = 2.08d-03*(1.d0/(1.d0 + zpars(11)))*rm
723  IF(abs(dma(k)).GT.0.4d0*dme) kw = 12
724  ENDIF
725 *
726 * Check various aspects related to change of type.
727  IF(kw.NE.kstar(i))THEN
728  IF(kw.EQ.15)THEN
729  WRITE(6,915)i,name(i),m0,kstar(i),tevk*tstar-epch0(k)
730  915 FORMAT (' MDOT WARNING! SN KW=15 I NM M0 KSTAR AGE ',
731  & 2i6,f7.3,i4,e9.1)
732  IF (name(n+kspair).LT.0) THEN
733  iphase = 7
734  CALL reset
735  go to 1
736  ENDIF
737  ELSEIF(kw.GE.13)THEN
738 *
739 * Force new NS or BH to have a period of one second.
740  ospin(k) = 2.0d+08
741  jspin(k) = k3*rcc*rcc*mc*ospin(k)
742  ENDIF
743 * Check inclusion of mass loss on type change.
744 * Note DM should be zero if ICORR not already true and no SN.
745  dms(k) = mass(k) - m1
746  dmr = abs(dms(k)/(mass(k) + 1.0d-10))
747  icorr = .true.
748  ENDIF
749 *
750 * Set mass loss and new radius in N-body units.
751  dmsun = dms(k)
752  dm = dms(k)/zmbar
753  rnew = rm/su
754  lnew = lum
755  kw0 = kstar(i)
756 *
757 * Check mass-loss treatment for inner binary components of merger.
758  IF (ihdot.GT.0) THEN
759  IF((kw.EQ.kstar(i).AND.age.GE.0.99*tn).OR.
760  & (kstar(i).LE.6.AND.kw.GT.6).OR.
761  & (kstar(i).LE.9.AND.kw.GT.9)) THEN
762  iterm = 1
763  ELSE
764  IF (ihdot.EQ.1) THEN
765  CALL hmdot(i,imerge,m1,kw,mc,dmsun,rnew,iterm)
766  ELSE
767  CALL hmdot2(i,ighost,imerge,m1,kw,mc,dmsun,rnew,iterm)
768  END IF
769  ENDIF
770 * Terminate on non-zero indicator.
771  IF (iterm.GT.0) THEN
772  iphase = 7
773  CALL reset
774  go to 1
775  END IF
776  END IF
777 *
778 * Define KS index & Roche indicator and update core mass (chaos only).
779  IF (i.LT.ifirst.AND.ihdot.EQ.0) THEN
780  kspair = kvec(i)
781  IF (.NOT.icorr) THEN
782  IF (name(n+kspair).GT.0) itry = 2
783  END IF
784  IF (kz(27).GT.1) THEN
785  ic = 0
786 * Search existing chaotic binaries (assign next index if none found).
787  DO 30 kk = 1,nchaos
788  IF (namec(kk).EQ.name(n+kspair)) ic = kk
789  30 CONTINUE
790  IF (ic.EQ.0) ic = nchaos + 1
791 * Save core mass in chaos variable for the appropriate KS component.
792  kk = 1
793  IF (i.EQ.2*kspair) kk = 2
794 * Note that KW <= 1 ensures initialization by MC = 0.
795  cm(kk,ic) = mc/zmbar
796  END IF
797  END IF
798 *
799 * Include special procedures for KS components.
800  IF (i.LT.ifirst.AND.icorr.AND.ihdot.EQ.0) THEN
801 * Distinguish between KS pair and merger configuration.
802  semi = -0.5d0*body(n+kspair)/h(kspair)
803  IF (name(n+kspair).GT.0) THEN
804 * Set random phase for neutron star or BH formation (negative index).
805  IF (kstar(i).LT.10.AND.kw.GE.10) THEN
806  jpair = -kspair
807  radius(i) = rnew
808  kstar(i) = -kstar(i)
809  CALL ksapo(jpair)
810  END IF
811 * Terminate for large mass loss or soft binary and re-determine index.
812  IF ((dmr.GT.0.2.AND.r(kspair).GT.rmin).OR.
813  & (dm.GT.0.0.AND.h(kspair) + dm/semi.GT.-eclose.AND.
814  & kstar(n+kspair).GE.0).OR.
815 * & (KW.NE.KSTAR(I).AND.KW.GE.13)) THEN
816  & (kw.NE.kstar(i).AND.
817  & (kw.GE.13.OR.(kw.GE.11.AND.kz(25).GE.1).OR.
818  & (kw.EQ.10.AND.kz(25).GT.1)))) THEN
819  i = i + 2*(npairs - kspair)
820  jx(k) = i
821  IF(kacc.EQ.2)THEN
822  jx(3-k) = jx(3-k) + 2*(npairs - kspair)
823  ENDIF
824 * Predict current KS variables and save at end of routine RESOLV.
825  CALL resolv(kspair,3)
826  iphase = 2
827  jcomp = 0
828  CALL ksterm
829  ks = 1
830  ELSE IF (dm.NE.0.d0) THEN
831 * Implement mass loss and expand KS orbit at constant eccentricity.
832  CALL hcorr(i,dm,rnew)
833  itry = 1
834  END IF
835  ELSE
836 * Adopt KS treatment for single outer component or terminate merger.
837  IF (i.EQ.2*kspair.AND.name(i).LE.nzero.AND.
838  & name(n+kspair).GT.-2*nzero.AND.dm.GT.0.0d0.AND.
839 * & H(KSPAIR) + DM/SEMI.LT.-ECLOSE.AND.KW.LT.13) THEN
840  & h(kspair) + dm/semi.LT.-eclose.AND.
841  & (kw.LT.10.OR.(kw.LT.11.AND.kz(25).LT.2).OR.
842  & (kw.LT.13.AND.kz(25).LT.1))) THEN
843  CALL hcorr(i,dm,rnew)
844  ELSEIF(dm.GT.0.d0)THEN
845  iphase = 7
846  CALL reset
847  go to 1
848  END IF
849  END IF
850  END IF
851 *
852 * Check for end of blue straggler evolution (allow 5% extra).
853  IF (m1.GT.1.05*turn.AND.kstar(i).EQ.1.AND.kw.EQ.2) THEN
854  WRITE(6,920) i, name(i), kw, tphys, age, m1, m0
855  920 FORMAT(' END BS: I NAM KW TP AGE M1 M0 ',
856  & 2i6,i4,2f8.1,2f6.1)
857  END IF
858 *
859 * Perform neighbour force corrections on significant mass loss (no WD).
860  IF (icorr) THEN
861 *
862 * Include optional diagnostics for mass loss orbit (filename MDOT).
863 * IF (KZ(21).GT.2) THEN
864 * CALL MTRACE(I,DM)
865 * END IF
866 *
867 * Update the mass (single body, standard KS or merger).
868  IF(ighost.EQ.0)THEN
869  body(i) = m1/zmbar
870  ELSE
871  j = 2*kspair - 2 + ihdot
872  body(j) = body(j) - dm
873  ENDIF
874  body0(i) = m0/zmbar
875  IF(dmsun.LT.tol) goto 250
876 *
877 * Accumulate total mass loss (solar units) and reduce cluster mass.
878  zmdot = zmdot + dmsun
879  zmass = zmass - dm
880 * Update the maximum single body mass but skip compact subsystems.
881  IF(mass(k)/zmbar.GE.0.99*body1.AND.nsub.EQ.0)THEN
882  body1 = 0.d0
883  DO 35 j = 1,n
884  body1 = max(body1,body(j))
885  35 CONTINUE
886  ENDIF
887 *
888 * Update the mass loss counters for types > 2.
889  IF(kw.EQ.3)THEN
890  zmrg = zmrg + dmsun
891  ELSEIF(kw.EQ.4)THEN
892  zmhe = zmhe + dmsun
893  ELSEIF(kw.EQ.5.OR.kw.EQ.6)THEN
894  zmrs = zmrs + dmsun
895  ELSEIF(kw.GE.7.AND.kw.LE.9)THEN
896  zmnh = zmnh + dmsun
897  ELSEIF(kw.GE.10.AND.kw.LE.12)THEN
898  zmwd = zmwd + dmsun
899  ELSEIF(kw.EQ.13.OR.kw.EQ.15)THEN
900  zmsn = zmsn + dmsun
901  ELSEIF(kw.EQ.14)THEN
902  zmbh = zmbh + dmsun
903  ENDIF
904 *
905 * Check optional diagnostics.
906  IF (kz(19).GT.3) THEN
907  WRITE(6,36) i, name(i), kw, kstar(i), body(i)*zmbar,
908  & dmsun, zmdot, tphys
909  36 FORMAT(' MDOT: I NM KW K* MS DMS ZMDOT T6 ',
910  & 4i5,f6.1,f7.2,f7.1,f8.1)
911  END IF
912 *
913 * Replace any KS components by corresponding c.m. for main procedures.
914  IF (i.LT.ifirst) THEN
915  iks = i
916  i = n + kspair
917  i1 = 2*kspair - 1
918 * Predict coordinates & velocities of any unperturbed KS components.
919  IF (list(1,i1).EQ.0) THEN
920  CALL resolv(kspair,1)
921  END IF
922  END IF
923 *
924 * Switch to the associated c.m. during force correction (save #I).
925  IF (ighost.GT.0) THEN
926  ii = i
927  i = n + kspair
928  END IF
929 *
930  nnb = list(1,i)
931  DO 38 l = 1,nnb+1
932  ilist(l) = list(l,i)
933  38 CONTINUE
934 * Ensure at least one neighbour.
935  IF (nnb.EQ.0) THEN
936  ilist(2) = n
937  IF (i.EQ.n) ilist(2) = n - 1
938  list(2,i) = ilist(2)
939  list(1,i) = 1
940  nnb = 1
941  END IF
942 * Include body #I at the end (counting from location #2).
943  nnb2 = nnb + 2
944  ilist(nnb2) = i
945 *
946 * Check prediction of neighbours and body #I to current time.
947  DO 40 l = 2,nnb2
948  j = ilist(l)
949  IF (t0(j).LT.time) THEN
950  CALL xvpred(j,-2)
951  END IF
952  40 CONTINUE
953 *
954 * Define logical variable to control optional WD kicks (avoids NaN).
955  ikick = .false.
956  IF (kz(25).EQ.1.AND.(kw.EQ.10.OR.kw.EQ.11)) ikick = .true.
957  IF (kz(25).EQ.2.AND.kw.EQ.12) ikick = .true.
958 * Ensure all NS/BH are assigned a kick (might depend on DM).
959  IF (kw.EQ.13.OR.kw.EQ.14) ikick = .true.
960 *
961 * Perform total force & energy corrections (delay dF if DMSUN > 0.1).
962  IF (dmsun.LT.0.05.AND.(kw.LT.10.OR..NOT.ikick)) THEN
963  CALL ficorr(i,dm)
964  ELSE
965  CALL fcorr(i,dm,kw)
966  END IF
967 *
968 * Set flag to ensure new sorting after CALL FPOLY (need IPHASE < 0).
969  IF (i.GT.n) ipoly = -1
970 *
971 * Initialize new polynomials of neighbours & #I for DMSUN > 0.1.
972  IF (abs(dmsun).GT.0.1.OR.kw.GE.13.OR.ikick) THEN
973 *
974 * Obtain new F & FDOT and time-steps (no gain skipping WDs).
975  DO 50 l = 2,nnb2
976  j = ilist(l)
977  IF (l.EQ.nnb2) j = i
978 * CALL DTCHCK(TIME,STEP(J),DTK(MAXBLK)) ! no effect (08/10).
979  DO 45 kk = 1,3
980  x0dot(kk,j) = xdot(kk,j)
981  45 CONTINUE
982  CALL fpoly1(j,j,0)
983  CALL fpoly2(j,j,0)
984  ipoly = -1
985  50 CONTINUE
986 * Check optional c.m. correction.
987  IF (kz(14).GE.3.AND.kz(31).GT.0) THEN
988  CALL xvpred(ifirst,ntot)
989  CALL cmcorr
990  END IF
991  END IF
992  tprev = time - stepx
993 * Set indicator < 0 for new sorting.
994  IF (ighost.GT.0) i = ii
995  IF(iphase.EQ.-3) iphase = 0
996  END IF
997 *
998 * Ensure that massless supernova remnant will escape next output.
999  IF (kw.EQ.15) THEN
1000  t0(i) = tadj + dtadj
1001  step(i) = 1.0d+06
1002  stepr(i) = 1.0d+06
1003  ri = sqrt(x(1,i)**2 + x(2,i)**2 + x(3,i)**2)
1004  vi = sqrt(xdot(1,i)**2 + xdot(2,i)**2 + xdot(3,i)**2)
1005  DO 55 l = 1,3
1006 * Ensure that ghost will escape next output (far from fast escapers).
1007  x0(l,i) = min(1.0d+04+x(l,i),1000.0*rscale*x(l,i)/ri)
1008  x(l,i) = x0(l,i)
1009  x0dot(l,i) = sqrt(0.004*zmass/rscale)*xdot(l,i)/vi
1010  xdot(l,i) = x0dot(l,i)
1011  f(l,i) = 0.0d0
1012  fdot(l,i) = 0.0d0
1013  d0(l,i) = 0.0
1014  d1(l,i) = 0.0
1015  d2(l,i) = 0.0d0
1016  d3(l,i) = 0.0d0
1017  d1r(l,i) = 0.0
1018  d2r(l,i) = 0.0d0
1019  d3r(l,i) = 0.0d0
1020  55 CONTINUE
1021  END IF
1022 *
1023  250 CONTINUE
1024 *
1025 * Restore index in case of KS component (used as c.m. above).
1026  IF (iks.GT.0) THEN
1027  i = iks
1028  iks = 0
1029 * Re-initialize KS polynomials for perturbed motion.
1030  IF (list(1,i1).GT.0) THEN
1031  CALL resolv(kspair,1)
1032  CALL kspoly(kspair,1)
1033  END IF
1034  END IF
1035 *
1036 * Update event counters for types > 2.
1037  IF(kw.NE.kw0)THEN
1038  IF(kw.EQ.3)THEN
1039  nrg = nrg + 1
1040  ELSEIF(kw.EQ.4)THEN
1041  nhe = nhe + 1
1042  ELSEIF(kw.EQ.5)THEN
1043  nrs = nrs + 1
1044  ELSEIF(kw.GE.7.AND.kw.LE.9.AND.kw0.LE.6)THEN
1045  nnh = nnh + 1
1046  ELSEIF(kw.GE.10.AND.kw.LE.12.AND.kw0.LE.9)THEN
1047  nwd = nwd + 1
1048  ELSEIF((kw.EQ.13.OR.kw.EQ.15).AND.kw0.LE.12)THEN
1049  nsn = nsn + 1
1050  ELSEIF(kw.EQ.14.AND.kw0.LE.13)THEN
1051  nbh = nbh + 1
1052  ENDIF
1053  ENDIF
1054 *
1055 * Include consistency warnings.
1056  IF (rnew - radius(i).GT.0.5*radius(i)) THEN
1057  WRITE(43,924) i, name(i), tphys, dt/1.0d+06,
1058  & kstar(i), kw, m0, m1, radius(i)*su, rnew*su
1059  924 FORMAT(' EXPAND! I NM TP DTP K* KW M0 M R RN ',
1060  & 2i6,f7.1,f7.3,2i4,2f7.1,2f7.1)
1061  CALL flush(43)
1062  END IF
1063 *
1064 * Update R, L, classification type & spin angular momentum.
1065  radius(i) = rnew
1066  zlmsty(i) = lnew
1067  kstar(i) = kw
1068  jspin(k) = max(jspin(k) - djspin(k),1.0d-10)
1069 *
1070 * Ensure that the star does not spin up beyond break-up.
1071  ospbru = twopi*sqrt(mass(k)*aursun**3/rad(k)**3)
1072  jspbru = (k2str(k)*(mass(k)-massc(k))*rad(k)*rad(k) +
1073  & k3*massc(k)*radc(k)*radc(k))*ospbru
1074  IF(jspin(k).GT.jspbru) jspin(k) = jspbru
1075  spin(i) = jspin(k)/spnfac
1076 *
1077 * Update epoch and check binary diagnostics for transition to new type.
1078  IF (kw.NE.kw0) THEN
1079 * EPOCH(I) = AGE0(K) + EPCH0(K) - AGE
1080 * IF(KW.GT.6)THEN
1081  IF(kw.LT.0)THEN
1082  WRITE(6,925)i,name(i),kw0,kw,m10,m1,dms(k),rm,age
1083  925 FORMAT(' MDOT CHANGE: I NM K* M0 M1 DM R AGE',
1084  & 2i6,2i4,2f6.1,f7.3,f7.2,f8.1)
1085  ENDIF
1086  IF (i.LT.ifirst.AND.kz(8).GT.3) THEN
1087  CALL binev(kspair)
1088  END IF
1089  END IF
1090 *
1091 * Include optional diagnostics (new type or significant mass loss).
1092  IF (kz(19).GT.3.AND.(kw0.NE.kw.OR.icorr)) THEN
1093  IF (kw0.NE.kw) THEN
1094  which1 = ' TYPE '
1095  ELSE
1096  which1 = ' MASS '
1097  END IF
1098  WRITE(6,926)which1, tphys, i, name(i), dmr, kw0, kw, m0, m1,
1099  & radius(i)*su, emdot
1100  926 FORMAT(' NEW',a8,' TPHYS I NAM DM/M KW0 KW M0 M R EMD ',
1101  & f7.1,2i5,f6.2,2i3,2f6.1,f7.1,f10.5)
1102  END IF
1103 *
1104 * Base new time scale for changes in radius & mass on stellar type.
1105  epoch(i) = age0(k) + epch0(k) - age
1106 * TEV(I) = (AGE0(K) + EPCH0(K))/TSTAR
1107  tev(i) = tevk
1108  tev0(i) = tev(i)
1109 * IF(I.EQ.IGHOST.OR.BODY(I).LE.0.0) RM0 = M1
1110  CALL trdot(i,dtm,m1)
1111  tev(i) = tev(i) + dtm
1112  IF(ihdot.EQ.2)THEN
1113  tev(i0) = tev(i)
1114  IF(name(ighost).LT.0.OR.name(ighost).GT.nzero)THEN
1115  tev(ighost) = tev(i)
1116  ENDIF
1117  ENDIF
1118 *
1119 * See if former KS pair can be regularized again.
1120  IF (ks.GT.0) THEN
1121  icomp = ifirst
1122  jcomp = ifirst + 1
1123  rij2 = (x(1,icomp) - x(1,jcomp))**2 +
1124  & (x(2,icomp) - x(2,jcomp))**2 +
1125  & (x(3,icomp) - x(3,jcomp))**2
1126  IF (rij2.LT.rmin22) THEN
1127 * Set temporary IPHASE > 0 for KSREG.
1128  iphase = 1
1129  CALL ksreg
1130  kspair = npairs
1131 * Restore current time to prevent small subsequent steps.
1132  time = tblock
1133  IF ((kw.EQ.13.OR.kw.EQ.14).AND.h(npairs).LT.0.0) THEN
1134 * IF (KW.GE.10.AND.KW.LE.14.AND.H(NPAIRS).LT.0.0) THEN
1135  j = ntot
1136  semi = -0.5d0*body(j)/h(npairs)
1137  ra = r(npairs)/semi
1138  ecc2 = (1.0 - ra)**2 + tdot2(npairs)**2/(body(j)*semi)
1139  tk = days*semi*sqrt(semi/body(j))
1140  WRITE(6,928) kw, sqrt(ecc2), ra, semi*su, tk, step(j),
1141  & body(j)*zmbar, (xdot(kk,j)*vstar,kk=1,3)
1142  928 FORMAT(' WD/NS BINARY KW E R/A A P DT M V ',
1143  & i4,2f6.2,1p,3e10.2,0p,4f7.1)
1144  END IF
1145  IF(kw.GE.10.AND.h(npairs).LT.0.0)THEN
1146  IF(kz(8).GT.3)THEN
1147  CALL degen(npairs,npairs,3)
1148  ELSE
1149  j1 = 2*npairs - 1
1150  j2 = j1 + 1
1151  IF(kstar(j1).GE.10.AND.kstar(j2).GE.10)THEN
1152  ndd = ndd + 1
1153  ENDIF
1154  ENDIF
1155  ENDIF
1156  ELSE
1157  kspair = 0
1158  END IF
1159  ks = 0
1160 * Ensure IPHASE < 0 on RETURN for new sorting after KSTERM.
1161  ipoly = -1
1162  END IF
1163 *
1164 * Note any formation of black holes or TZ object.
1165  IF(kw.EQ.14.AND.kstar(i).LT.14.AND.i.LE.n)THEN
1166  WRITE(6,930) i, name(i), kw0, kw, kstar(i), m0, m1, dmr
1167  930 FORMAT(' NEW BH/TZ I NM K0 KW K* M0 M1 DM/M ',
1168  & 2i6,3i4,3f7.2)
1169  END IF
1170 *
1171 * Perform consistency check on massive WD (skip HMDOT treatment).
1172  IF(i.LE.n)THEN
1173  IF((kstar(i).GE.10.AND.kstar(i).LE.12).AND.ihdot.EQ.0)THEN
1174  IF(body(i)*zmbar.GT.mch)THEN
1175  WRITE(6,932)i,kw,body0(i)*zmbar,body(i)*zmbar,radius(i)
1176  932 FORMAT(' DANGER! MDOT I K* M0 M R ',2i5,2f7.2,1p,e9.1)
1177  WRITE(6,934) name(i), ttot, tev0(i), tev(i)
1178  934 FORMAT(' NAM T TEV0 TEV ',i6,3f10.2)
1179  stop
1180  ENDIF
1181  ENDIF
1182  ENDIF
1183 *
1184 * Include stellar radius check in case of NaN.
1185  IF (radius(i).GE.0.0.AND.radius(i).LT.1.0) go to 105
1186  WRITE(6,936) i, kstar(i), m1, radius(i)
1187  936 FORMAT(' DANGER! MDOT I K* M1 R* ',i5,i4,f7.2,1p,e10.1)
1188  stop
1189  105 CONTINUE
1190 *
1191 * Check for chaotic binary.
1192  IF(itry.GT.0.AND.ighost.EQ.0)THEN
1193  ii = i
1194  i = n + kspair
1195  ipair = kspair
1196  IF(kstar(i).LT.0)THEN
1197  IF(kw.GE.10.AND.dmr.GT.0.0)THEN
1198  WRITE(6,940) name(ii), kstar(ii), kw, dmr
1199  940 FORMAT(' DEGEN SPIRAL NAM K* DMR ',i5,2i4,1p,e9.1)
1200  ENDIF
1201  CALL chrect(ipair,dmr)
1202  IF (iphase.LE.0) ipoly = -1
1203  IF(dmr.LT.0.0) goto 1
1204  ENDIF
1205  IF(kw.GE.10.AND.icorr)THEN
1206  IF(kz(8).GT.3)THEN
1207  CALL degen(ipair,ipair,3)
1208  ELSE
1209  j1 = 2*ipair - 1
1210  j2 = j1 + 1
1211  IF(kstar(j1).GE.10.AND.kstar(j2).GE.10) ndd = ndd + 1
1212  ENDIF
1213  ENDIF
1214  itry = 0
1215  ENDIF
1216 *
1217  220 CONTINUE
1218 *
1219  IF(kspair.EQ.0) kacc = 1
1220  IF(kacc.GT.1)THEN
1221  tev(jx(1)) = min(tev(jx(1)),tev(jx(2)))
1222  IF(dtgr.LT.1.0d+10)THEN
1223  tev(jx(1)) = min(tev(jx(1)),tev0(jx(1))+dtgr/tstar)
1224  ENDIF
1225  tev(jx(2)) = tev(jx(1))
1226  ENDIF
1227 *
1228 * Check for Roche overflow condition (DTR < STEP) and implement any
1229 * orbit change owing to gravitational radiation, mass loss or tides.
1230  IF(kspair.GT.0.AND.ighost.EQ.0)THEN
1231  i = n + kspair
1232  ipair = kspair
1233 * IF(H(IPAIR).LT.0.D0.AND.NAME(ICM).GT.0)THEN
1234  IF(h(ipair).LT.0.d0.AND.name(i).GT.0.AND.kstar(i).GE.0)THEN
1235  CALL brake(ipair,dsep,ecc)
1236  IF(iqcoll.NE.0.OR.iphase.LT.0) go to 1
1237 *
1238 * Include optional look-up time control for compact object binaries.
1239 * IF (KSX.GE.13.AND.KZ(28).GT.2) THEN
1240 * WRITE (6,944) TTOT, NAME(2*IPAIR-1),KSTAR(2*IPAIR),
1241 * & DTGR/TSTAR
1242 * 944 FORMAT (' GR CHECK T NAM K* DTGR ',F8.2,I6,I4,1P,E9.1)
1243 * IF (KSX.GE.13.AND.KZ(28).GT.0) THEN
1244 * GE = (1.0 - ECC2)**3.5/(1.0 + 3.0*ECC2)
1245 * ZMX = MAX(BODY(2*IPAIR-1),BODY(2*IPAIR))
1246 * RATIO = MIN(BODY(2*IPAIR-1),BODY(2*IPAIR))/ZMX
1247 * Replace physical time-scale by N-body units (cf. Book, Eq.12.31).
1248 * TZ = GE*SEMI**4/(RATIO*(1.0 + RATIO)*ZMX**3)
1249 * TZ = 5.0/64.0*CLIGHT**5*TZ
1250 * IF (KZ(28).GT.2) THEN
1251 * WRITE (6,944) TTOT, NAME(2*IPAIR-1),KSTAR(2*IPAIR),TZ
1252 * END IF
1253 * Limit the time-scale in case of shrinkage from non-GR orbit.
1254 * TZ = MIN(TZ,100.0D0)
1255 * Specify new look-up as 1% of current GR time-scale (ignore old TEV).
1256 * TEV(2*IPAIR-1) = TIME + 0.01*TZ
1257 * TEV(2*IPAIR) = TEV(2*IPAIR-1)
1258 * TMDOT = MIN(TMDOT,TEV(2*IPAIR))
1259 * END IF
1260 *
1261  IF(kstar(i).GT.0.AND.kz(34).GT.0)THEN
1262 * Ensure optional updating of orbit before possible Roche test (1/2013).
1263  IF (kz(34).EQ.1.AND.kstar(i).GT.0.AND.
1264  & mod(kstar(i),2).EQ.0)THEN
1265  CALL synch(ipair)
1266  IF(iqcoll.NE.0.OR.iphase.LT.0) go to 1
1267  ENDIF
1268  jpair = -ipair
1269  CALL trflow(jpair,dtr)
1270  IF(dtr.LT.step(i))THEN
1271  CALL roche(ipair)
1272  IF(iqcoll.NE.0.OR.iphase.LT.0) go to 1
1273  ELSE
1274  tev(i) = time + dtr
1275  IF(kstar(i).GT.10.AND.mod(kstar(i),2).EQ.1)THEN
1276  j1 = 2*ipair - 1
1277  j2 = j1 + 1
1278  tev(j1) = max(1.000002d0*tev(i),tev(j1))
1279  tev(j2) = max(1.000002d0*tev(i),tev(j2))
1280  ENDIF
1281  ENDIF
1282  ENDIF
1283  ENDIF
1284  ENDIF
1285 *
1286 * Determine the time for next stellar evolution check.
1287  70 tmdot = 1.0d+10
1288  DO 80 j = 1,ntot
1289  IF(tev(j).LE.tmdot)THEN
1290  tmdot = tev(j)
1291  ENDIF
1292  80 CONTINUE
1293 *
1294 * See whether any other stars should be considered.
1295  IF (tmdot.LT.time) go to 1
1296 *
1297 * Update any merged circularizing binary using TMDIS.
1298  imerge = 1
1299  104 IF (nmerge.GT.0) THEN
1300  IF (kstarm(imerge).LT.0.AND.tmdis(imerge).LT.time) THEN
1301  kspair = 0
1302 * Determine the KS index from c.m. identification.
1303  DO 106 i = n+1,ntot
1304  IF (name(i).EQ.namem(imerge)) THEN
1305  kspair = i - n
1306  END IF
1307  106 CONTINUE
1308  IF (kspair.GT.0) THEN
1309  ic = 0
1310  DO 107 k = 1,nchaos
1311  IF (namec(k).EQ.name(n+kspair)) ic = k
1312  107 CONTINUE
1313  ic = max(ic,1)
1314  semi = -0.5*(bm(1,imerge) + bm(2,imerge))/hm(imerge)
1315  i1 = 2*kspair - 1
1316  WRITE (6,108) nameg(imerge), name(n+kspair),kstar(i1),
1317  & es(ic), tmdis(imerge), semi
1318  108 FORMAT (' TMDIS TERM NAMG NAMCM K* E TMD A ',
1319  & 2i7,i4,f8.4,f8.2,1p,e10.2)
1320 * Terminate hierarchy (RESET calls SPIRAL via CHRECT).
1321  iphase = 7
1322  CALL reset
1323 * Consider the same merger index again.
1324  imerge = imerge - 1
1325  END IF
1326  END IF
1327  imerge = imerge + 1
1328  END IF
1329  IF (imerge.LE.nmerge) go to 104
1330 *
1331 * Update the maximum single body mass but skip compact subsystems.
1332  IF(nsub.EQ.0)THEN
1333  body1 = 0.d0
1334  DO 110 j = 1,n
1335  body1 = max(body1,body(j))
1336  110 CONTINUE
1337  ENDIF
1338 *
1339  DO 122 j = n+1,ntot
1340  IF(kstar(j).EQ.0.AND.name(j).GT.0.AND.tev(j).LT.9.9e+09.AND.
1341  & body(j).GT.0.0.AND.kz(28).EQ.0)THEN
1342  WRITE(6,556)j,name(j),tev(j)
1343  tev(j) = 1.0d+10
1344  556 FORMAT(' MDOT TEV SMALL ',2i8,1p,e10.2)
1345  ENDIF
1346  122 CONTINUE
1347 *
1348 * Ensure re-initialization of ICPERT (INTGRT) and KBLIST (SUBINT).
1349  IF (ipoly.LT.0) THEN
1350  iphase = -1
1351  nbprev = 0
1352  END IF
1353 *
1354  RETURN
1355 *
1356  END
1357 ***