Nbody6
 All Files Functions Variables
synch.f
Go to the documentation of this file.
1  SUBROUTINE synch(IPAIR)
2 *
3 *
4 * Spin synchronization of circularized orbit.
5 * -------------------------------------------
6 *
7 * Rational function approximations of solution to the Hut
8 * evolution equation with spin. Ref: A & A 99, 126, eqn (A15).
9 *
10  include 'common6.h'
11  common/modes/ eb0(ntmax),zj0(ntmax),ecrit(ntmax),ar(ntmax),
12  & br(ntmax),eosc(4,ntmax),edec(ntmax),tosc(ntmax),
13  & rp(ntmax),es(ntmax),cm(2,ntmax),iosc(ntmax),
14  & namec(ntmax)
15  common/spins/angmom0,rg2(2),m21,r21,semi0,c1,c2,c3,c4,c5,semi
16  common/radii/ r1,r2
17  REAL*8 ww(3),qq(3),w(2),q(2),at0(2),m21,wg(2),qg(2),wscale(2),
18  & qscale(2),a(2),b(2),c(6)
19  REAL*8 m0,m1,mc,mc1,mc2,corerd,mlwind
20  REAL*8 tscls(20),lums(10),gb(10),lum,k2,menv
21  EXTERNAL corerd,mlwind
22  DATA ww /2.119,3.113,8.175/
23  DATA qq /0.4909,0.4219,0.2372/
24  DATA a /6.306505,-7.297806/
25  DATA b /32.17211,13.01598/
26  DATA c /5.101417,24.71539,-9.627739,1.733964,
27  & -2.314374,-4.127795/
28  SAVE icount,itry,isynch,iter,itime,iname,eccm
29  DATA icount,itry,isynch,iter,itime,eccm /0,0,0,0,0,0.002d0/
30  DATA iname /0/
31 *
32 *
33 * Define c.m. & KS indices and search current names for chaos index.
34  i = n + ipair
35  i1 = 2*ipair - 1
36  i2 = i1 + 1
37  ic = 0
38  DO 2 k = 1,nchaos
39  IF (namec(k).EQ.name(i)) ic = k
40  2 CONTINUE
41 *
42 * Count new cases ignoring rare multiple simultaneous events.
43  IF (name(i).NE.iname) THEN
44  iname = name(i)
45  nsync = nsync + 1
46  END IF
47 *
48 * Exit if binary is not identified but include updating of small ECC.
49  IF (ic.EQ.0) THEN
50  jpair = -ipair
51  CALL trflow(jpair,dtr)
52  m1 = body(i1)*smu
53  CALL trdot(i1,dtm,m1)
54  m1 = body(i2)*smu
55  CALL trdot(i2,dtm2,m1)
56  dtm = min(dtm,dtm2)
57  dtm = max(dtm,0.1d0)
58  tev(i) = time + 0.1*min(dtm,dtr)
59  a0 = -0.5*body(i)/h(ipair)
60  ecc2 = (1.0 - r(ipair)/a0)**2 + tdot2(ipair)**2/(a0*body(i))
61  ecc = sqrt(ecc2)
62  itime = itime + 1
63  IF (itime.GT.10) THEN
64  WRITE (6,3) name(i1), kstar(i), nchaos, ecc, a0, dtm, dtr
65  3 FORMAT (' MISSING! SYNCH NM K* NCH ECC A DTM DTR ',
66  & i6,2i4,f8.4,1p,3e10.2)
67  END IF
68 *
69 * Check enforcement of circularization for small TCIRC.
70  IF (ecc.LT.0.10.AND.list(1,i1).EQ.0) THEN
71  icirc = -1
72  qperi = a0*(1.0 - ecc)
73  CALL tcirc(qperi,ecc,i1,i2,icirc,tc)
74  IF (tc.LT.100.0) THEN
75  CALL deform(ipair,ecc,eccm)
76  ecc = eccm
77  END IF
78  END IF
79 *
80 * Initialize chaos elements after possible common envelope or E < 0.1.
81  IF (ecc.LE.0.003.OR.ic.EQ.0) THEN
82  nchaos = nchaos + 1
83  ic = nchaos
84  namec(ic) = name(i)
85  rp(ic) = a0*(1.0 - ecc)
86  es(ic) = ecc
87  tosc(ic) = time
88  IF (nchaos.GT.ntmax) THEN
89  WRITE (6,4) name(i1), nchaos, ecc, a0*(1.0 - ecc)
90  4 FORMAT (' FATAL ERROR! SYNCH NM NCH E QP ',
91  & i6,i4,f8.4,1p,e9.1)
92  stop
93  END IF
94  END IF
95  go to 100
96  END IF
97 *
98 * Check Roche overflow time during small intervals.
99  IF (time - tosc(ic).LT.1.0) THEN
100  jpair = -ipair
101  CALL trflow(jpair,dtr)
102  IF (dtr.LT.step(i)) THEN
103  tev(i) = time + dtr
104  END IF
105  END IF
106 *
107 * Copy spiral parameters and set semi-major axis & eccentricity.
108  time0 = tosc(ic)
109  rp0 = rp(ic)
110  es0 = es(ic)
111  hi = h(ipair)
112  a0 = -0.5*body(i)/h(ipair)
113  ecc2 = (1.0 - r(ipair)/a0)**2 + tdot2(ipair)**2/(a0*body(i))
114  ecc = sqrt(ecc2)
115  ecc0 = ecc
116 *
117 * Rectify elements after possible common envelope, Roche or GR process.
118  qperi = a0*(1.0 - ecc)
119  IF (es0.GT.0.01.OR.abs(qperi-rp0).GT.0.02*qperi) THEN
120  dt = time - tosc(ic)
121  IF (itime.LT.50) THEN
122  WRITE (6,5) name(i1), kstar(i1),kstar(i2),kstar(i), es0,
123  & ecc, rp0, qperi, dt
124  5 FORMAT (' RP RECTIFY SYNCH NM K* ES0 E RP0 QP DT ',
125  & i6,3i4,2f8.4,1p,3e10.2)
126  END IF
127 * Note: new RP(IC) needed for scaling RADIUS (original bug Oct 2008).
128  rp(ic) = qperi
129  rp0 = qperi
130  es0 = ecc
131  END IF
132 *
133 * Set standard binary or current elements on significant departure.
134  IF (ecc.GT.0.01.OR.es0.GT.0.01) THEN
135  ecc = sqrt(ecc2)
136  IF (ecc.GT.0.01) THEN
137  kstar(i) = 0
138  WRITE (6,8) name(i1), list(1,i1), es0, ecc, a0*su
139  8 FORMAT (' NON-CIRCULAR SYNCH NM NP ES0 ECC A ',
140  & i6,i4,2f8.4,1p,e10.2)
141  nami = -i
142  CALL spiral(nami)
143  go to 100
144  ELSE
145  es0 = ecc
146  rp0 = a0*(1.0 - ecc)
147  END IF
148  END IF
149 *
150 * Specify index J1 as biggest radius to be used with AT0(1).
151  IF (radius(i1).GE.radius(i2)) THEN
152  j1 = i1
153  j2 = i2
154  ELSE
155  j1 = i2
156  j2 = i1
157  END IF
158 *
159 * Obtain stellar parameters for evolving stars at current epoch.
160  IF (kstar(j1).GE.2.AND.kstar(j1).LE.9) THEN
161  kw = kstar(j1)
162  m1 = body(j1)*smu
163  m0 = body0(j1)*zmbar
164  mc = 0.d0
165  age = time*tstar - epoch(j1)
166  CALL star(kw,m0,m1,tm,tn,tscls,lums,gb,zpars)
167  CALL hrdiag(m0,age,m1,tm,tn,tscls,lums,gb,zpars,
168  & rm,lum,kw,mc,rcc,menv,renv,k2)
169  cm(1,ic) = mc/smu
170  kw1 = kw
171  ELSE
172  kw1 = kstar(j1)
173  END IF
174 *
175 * Define oscillation period (dimensionless time) and damping constants.
176  xn = 0.0
177  qd = 0.0
178  DO 10 k = 1,2
179  ik = i1 + k - 1
180  IF (k.EQ.1) THEN
181  ik = j1
182  ELSE
183  ik = j2
184  END IF
185 * Specify polytropic index for each star (n = 3, 2 or 3/2).
186  IF (kstar(ik).EQ.3.OR.kstar(ik).EQ.5.OR.
187  & kstar(ik).EQ.6.OR.kstar(ik).EQ.9) THEN
188  CALL giant(ipair,ik,wg,qg,wscale,qscale,xn,ql)
189  w(k) = wg(1)
190  q(k) = qg(1)
191  rg2(k)= 0.1*(1.0 - cm(k,ic)/body(ik))
192  qd = ql
193  ELSE
194  ql = 1.0e+04
195  ip = 3
196  IF (kstar(ik).GE.3) ip = 2
197  IF (kstar(ik).EQ.4.OR.kstar(ik).EQ.7) ip = 3
198  IF (kstar(ik).EQ.8) ip = 3
199  IF (kstar(ik).EQ.0) ip = 1
200  w(k) = ww(ip)
201  q(k) = qq(ip)
202  IF (kstar(ik).LE.2.OR.kstar(ik).EQ.7) THEN
203  rg2(k) = 0.1
204  ELSE IF (kstar(ik).EQ.4) THEN
205  cm(k,ic) = min(0.89d0*body(ik),cm(k,ic))
206  rg2(k)= 0.1*(1.0 - cm(k,ic)/body(ik))
207  ELSE
208  rg2(k)= 0.21
209  END IF
210  END IF
211  tl = twopi*radius(ik)*sqrt(radius(ik)/body(ik)/w(k))
212  IF (kstar(ik).GE.11) ql = 1.0d+10
213  at0(k) = 1.0/(ql*tl)
214  10 CONTINUE
215 *
216 * Form mass, radius & pericentre ratio.
217  IF (radius(i1).GE.radius(i2)) THEN
218  m21 = body(i2)/body(i1)
219  r21 = radius(i2)/radius(i1)
220  rp1 = rp(ic)/radius(i1)
221  rad = radius(i1)
222  ELSE
223  m21 = body(i1)/body(i2)
224  r21 = radius(i1)/radius(i2)
225  rp1 = rp(ic)/radius(i2)
226  rad = radius(i2)
227  END IF
228 *
229 * Define initial angular momentum from the scaled semi-major axis.
230  semi0 = rp1/(1.0 - es0)
231 *
232 * Form the initial mean motion in N-body units.
233  omeq = sqrt(body(i)/(rad*semi0)**3)
234 *
235 * Convert from angular momentum to omega (denoted spin1 & spin2).
236  IF (kstar(j1).LE.2.OR.(kstar(j1).GE.7.AND.kstar(j1).NE.9)) THEN
237  spin1 = spin(j1)/(rg2(1)*body(j1)*radius(j1)**2)
238  ELSE
239  kw = kstar(j1)
240  m0 = body0(j1)*smu
241  mc1 = cm(1,ic)*smu
242  IF (mc1.LE.0.0d0.OR.mc1.GT.m0) THEN
243  mc1 = 0.3 + 0.1*(kstar(j1) - 3)
244  IF(kw.EQ.9) mc1 = min(0.3d0,0.95*m0)
245  cm(1,ic) = mc1/zmbar
246  END IF
247  zdum = 2.0d0
248  rc1 = corerd(kw,mc1,m0,zdum)/su
249  spin1 = spin(j1)/(rg2(1)*body(j1)*radius(j1)**2 +
250  & 0.21*mc1/smu*rc1**2)
251  END IF
252  IF (kstar(j2).LE.2.OR.(kstar(j2).GE.7.AND.kstar(j2).NE.9)) THEN
253  spin2 = spin(j2)/(rg2(2)*body(j2)*radius(j2)**2)
254  ELSE
255  kw = kstar(j2)
256  m0 = body0(j2)*smu
257  mc2 = cm(2,ic)*smu
258  IF (mc2.LE.0.0d0.OR.mc2.GT.m0) THEN
259  mc2 = 0.3 + 0.1*(kstar(j2) - 3)
260  IF(kw.EQ.9) mc2 = min(0.3d0,0.95*m0)
261  cm(2,ic) = mc2/zmbar
262  END IF
263  zdum = 2.0d0
264  rc2 = corerd(kw,mc2,m0,zdum)/su
265  spin2 = spin(j2)/(rg2(2)*body(j2)*radius(j2)**2 +
266  & 0.21*mc2/smu*rc2**2)
267  END IF
268 *
269 * Set synchronous spin for degenerate stars.
270  ds10 = (spin1 - omeq)/omeq
271  ds20 = (spin2 - omeq)/omeq
272  IF (kstar(j1).GE.10) spin1 = omeq
273  IF (kstar(j2).GE.10) spin2 = omeq
274  ds1 = (spin1 - omeq)/omeq
275  ds2 = (spin2 - omeq)/omeq
276 *
277 * Skip integration if both spins are synchronous.
278  IF (abs(ds1).LT.0.01.AND.abs(ds2).LT.0.01) THEN
279  kx = max(kstar(j1),kstar(j2))
280  IF (abs(ds20).GT.0.01.AND.kx.LT.10) THEN
281  WRITE (6,20) name(i1), kstar(i1), kstar(i2), a0*su, ds1,
282  & ds20, omeq
283  20 FORMAT (' ENFORCED SYNCH NM K* A DS1 DS2 omeq ',
284  & i6,2i4,f9.3,1p,3e10.2)
285  END IF
286  isynch = 1
287  go to 50
288  END IF
289 *
290  ds = (spin1 - omeq)/omeq
291 * Scale the spins by mean motion and define angular momentum.
292  spin10=spin1/omeq
293  spin20=spin2/omeq
294  angmom0=(m21/(1+m21))*semi0**2*sqrt(1-es0**2)+rg2(1)*spin10+
295  & m21*r21**2*rg2(2)*spin20
296 *
297  IF (icount.LE.-1000) THEN
298  icount = icount + 1
299  sum1=(m21/(1+m21))*semi0**2*sqrt(1-es0**2)
300  sum2 = rg2(1)*spin10
301  sum3 = m21*r21**2*rg2(2)*spin20
302  sum4 = sum1 + sum2 + sum3
303  zm = body(j1)*smu
304  WRITE (92,21)kstar(j1),kstar(j2),zm,a0,sum1,sum2,sum3,sum4,ds
305  21 FORMAT (' K* M1 A S1 S2 S3 S4 DS ',2i4,f8.4,1p,e12.4,4e10.2,
306  & 0p,f6.1)
307  CALL flush(92)
308  END IF
309 * Evaluate damping coefficients (Mardling & SJA, M.N. 321, 398, 2001).
310  cf = 54.0*twopi/5.0
311  c3 = (cf/9.0)*(at0(1)*(q(1)/w(1))**2*m21**2)/rg2(1)/semi0**6
312  c4 = (cf/9.0)*(at0(2)*(q(2)/w(2))**2/m21**2)*r21**6/rg2(2)/
313  & semi0**6
314 *
315 * Include change in moment of inertia for rapid evolution of primary.
316  IF (kw1.GE.2.AND.kw1.LE.9) THEN
317 * Obtain mass loss due to stellar wind for single star (Msun/yr).
318  rlperi = 0.d0
319  dmx = mlwind(kw1,lum,rm,m1,mc,rlperi,zmet)
320  dt = time - tev0(j1)
321 * Include safety check on small time interval (SJA 4/09).
322  IF (dt.LE.1.0d-10) dt = abs(tev(j1) - time)
323  rdot = (rm/su - radius(j1))/dt
324 * Form combined moment of inertia factor for primary (RM 11/08).
325  c5 = -1.0d+06*dmx*tstar/m1 + 2.0*rdot/radius(j1)
326  r1 = rm/(su*radius(j1))
327  r2 = 1.0
328  ELSE
329  r1 = 1.0
330  r2 = 1.0
331  c5 = 0.0
332  END IF
333 *
334 * Skip on zero or negative interval.
335  IF (time - time0.LE.0.0d0) go to 70
336 *
337 * Obtain dominant terms of the derivatives (cf. routine HUT/DERIV2).
338  udot1=c3
339  udot2=c4
340 *
341 * Choose the step from smallest time-scale (! time0 < time possible).
342  taux = min(abs(1.0/udot1),abs(1.0/udot2))
343  nstep = 1 + 100.0*sqrt(abs(time - time0)/taux)
344  nstep = min(nstep,100)
345 * Increase number of steps on slow primary rotation.
346  IF (ds1.LT.-0.2) nstep = nstep + 10
347  IF (ds1.LT.-0.8) nstep = nstep + 50
348  IF (ds1.LT.-0.9) nstep = nstep + 50
349  dt = time - time0
350 *
351 * Reduce integration interval for slow evolution or after merger.
352  IF (dt.GT.0.5) THEN
353  jpair = -ipair
354  CALL trflow(jpair,dtr)
355  m1 = body(j1)*smu
356  CALL trdot(j1,dtm,m1)
357  dtm = max(dtm,0.1d0)
358  dt = 0.1*min(dtr,dtm,10.0d0)
359  IF (abs(ds2).GT.0.05.OR.ds10.LT.-0.99) THEN
360  dt = 0.1*dt
361  nstep = nstep + 50
362  END IF
363  tk = days*a0*sqrt(a0/body(i))
364  IF (icount.LT.50) THEN
365  WRITE (6,30) name(j1), ds1, ds2, dt, tk
366  30 FORMAT (' REDUCE SYNCH NM DS1 DS2 DT P ',
367  & i6,1p,4e10.2)
368  END IF
369  END IF
370 *
371  itt = 0
372  IF (kstar(j1).EQ.3) dt = 0.5*dt
373  40 dtau1 = abs(dt)/float(nstep)
374 *
375 * Integrate equations for angular velocities only.
376  call hut2(spin10,spin20,spin1,spin2,nstep,dtau1)
377 *
378  spin11=spin1
379  spin21=spin2
380 * Re-scale the semi-major axis and angular velocities to N-body units.
381  semi = rad*semi0*semi
382  spin1 = omeq*spin1
383  spin2 = omeq*spin2
384  ecc = es0
385  omeq = sqrt(body(i)/semi**3)
386  ds1 = (spin1 - omeq)/omeq
387  ds2 = (spin2 - omeq)/omeq
388  nsynch = nsynch + 1
389 *
390  itt = itt + 1
391 * IF (ITT.GE.3) STOP
392  IF (spin1.LT.0.0.OR.spin2.LT.0.0) THEN
393  WRITE (6,42) itry, nstep, dt, spin1, spin2, omeq
394  42 FORMAT (' REPEAT SYNCH IT # DT S1 S2 omeq ',
395  & 2i5,1p,4e10.2)
396  itry = itry + 1
397  nstep = nstep + 50
398  IF (itry.LT.2) go to 40
399 * STOP
400  END IF
401  itry = 0
402 *
403  IF ((ds10.LT.0.0.AND.ds1.GT.0.01).OR.
404  & (ds20.LT.0.0.AND.ds2.GT.0.01)) THEN
405  iter = iter + 1
406  nstep = nstep + 50
407  IF (iter.GT.1) dt = 0.5*dt
408  IF (iter.LT.4) go to 40
409  WRITE (6,45) nstep, ds10, ds1, ds20, ds2, dt
410  45 FORMAT (' LIMIT SYNCH # DS10, DS1, DS20, DS2, DT ',
411  & i4,1p,5e10.2)
412  END IF
413 *
414  IF ((ds1.GT.0.0.AND.ds10.LT.-0.01).OR.
415  & (ds2.GT.0.0.AND.ds20.LT.-0.01)) THEN
416  iter = iter + 1
417  nstep = nstep + 50
418  IF (iter.GT.1) dt = 0.5*dt
419  IF (iter.LT.4) go to 40
420  WRITE (6,45) nstep, ds10, ds1, ds20, ds2, dt
421  END IF
422  iter = 0
423 *
424 * Convert back to angular momenta.
425  50 IF (kstar(j1).LE.2.OR.(kstar(j1).GE.7.AND.kstar(j1).NE.9)) THEN
426  spin(j1) = rg2(1)*body(j1)*radius(j1)**2*spin1
427  ELSE
428  spin(j1) = (rg2(1)*body(j1)*radius(j1)**2 +
429  & 0.21*mc1/smu*rc1**2)*spin1
430  END IF
431  IF (kstar(j2).LE.2.OR.(kstar(j2).GE.7.AND.kstar(j2).NE.9)) THEN
432  spin(j2) = rg2(2)*body(j2)*radius(j2)**2*spin2
433  ELSE
434  spin(j2) = (rg2(2)*body(j2)*radius(j2)**2 +
435  & 0.21*mc2/smu*rc2**2)*spin2
436  END IF
437  IF (isynch.GT.0) THEN
438  isynch = 0
439  go to 70
440  END IF
441 *
442  err = (semi - a0)/a0
443 * Obtain the tidal contributions from integration.
444  a0 = rp0/(1.0 - es0)
445  dh = -0.5*body(i)*(1.0/semi - 1.0/a0)
446 * Prevent new SEMI < R.
447  IF (h(ipair) + dh.LT.-0.5*body(i)/r(ipair)) THEN
448  dh = -0.5*body(i)/r(ipair) - h(ipair)
449  END IF
450 *
451 * Update energy and semi-major axis.
452  h(ipair) = h(ipair) + dh
453  semi = -0.5*body(i)/h(ipair)
454 *
455 * Set new pericentre from final elements (delay updating).
456  IF (list(1,i1).EQ.0) THEN
457  qperi = semi*(1.0 - ecc)
458  ELSE
459  qperi = rp0*(1.0 + es0)/(1.0 + ecc)
460  END IF
461 *
462 * Form KS coordinate scaling factor from pericentre ratio.
463  c1 = sqrt(qperi/rp0)
464 *
465 * Set KS velocity scaling from energy relation with RP0 instead of R.
466  c2 = sqrt((body(i) + h(ipair)*qperi)/(body(i) + hi*rp0))
467 *
468 * Scale KS variables to yield the prescribed elements (set TDOT2 >= 0).
469  r(ipair) = 0.0d0
470  tdot2(ipair) = 0.0
471  DO 60 k = 1,4
472  u(k,ipair) = c1*u(k,ipair)
473  udot(k,ipair) = c2*udot(k,ipair)
474  u0(k,ipair) = u(k,ipair)
475  r(ipair) = r(ipair) + u(k,ipair)**2
476  tdot2(ipair) = tdot2(ipair) + 2.0*u(k,ipair)*udot(k,ipair)
477  60 CONTINUE
478  tdot2(ipair) = max(tdot2(ipair),0.0d0)
479 *
480 * Perform energy correction to maintain conservation.
481  zmu = body(i1)*body(i2)/body(i)
482  ecoll = ecoll + zmu*(hi - h(ipair))
483  egrav = egrav + zmu*(hi - h(ipair))
484 *
485 * Initialize KS for perturbed motion.
486  IF (list(1,i1).GT.0) THEN
487 * Rectify orbit to prevent eccentricity growth.
488  CALL ksrect(ipair)
489  CALL resolv(ipair,1)
490  imod = kslow(ipair)
491  CALL kspoly(ipair,imod)
492  END IF
493 *
494  da = (semi - a0)/semi
495  ds1 = (spin1-omeq)/omeq
496  ds2 = (spin2-omeq)/omeq
497  ts1 = min(tev(j1),tev(j2))
498  WRITE (7,65) nstep,time-time0,udot1,udot2,ds1,ds2,da,ts1
499  65 FORMAT (' HUT2 # t-t0 ud DS DA/A TM ',
500  & i5,f9.3,1p,5e9.1,0p,f10.2)
501 *
502 * Include magnetic or gravitational braking for small separations.
503 * IF (SEMI*SU.LT.10.0) THEN
504 * DT = TIME - TIME0
505 * CALL BRAKE(IPAIR,DT)
506 * IF (IPHASE.LT.0) GO TO 100
507 * END IF
508 *
509 * Update look-up time from radial expansion or Roche interval.
510  70 jpair = -ipair
511  CALL trflow(jpair,dtr)
512  m1 = body(j1)*smu
513  CALL trdot(j1,dtm,m1)
514  dtm = max(dtm,0.1d0)
515  tosc(ic) = time
516  tev(i) = time + 0.1*min(dtr,dtm)
517  tev(i) = min(tev(j1),tev(j2),tev(i))
518 *
519  semi = -0.5*body(i)/h(ipair)
520  ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(semi*body(i))
521  ecc = sqrt(ecc2)
522 * Correct for small eccentricity drift.
523  IF (abs(ecc - es0).GT.0.001*es0) THEN
524  CALL deform(ipair,ecc,es0)
525  ecc = es0
526  qperi = semi*(1.0 - es0)
527  END IF
528 * Save final reference values.
529  rp(ic) = qperi
530  es(ic) = ecc
531  tosc(ic) = time
532 *
533  IF (ecc.GT.0.01) THEN
534  WRITE (6,80) name(i1), ecc, nstep, dtau1, err, c1-1.0, c2-1.0
535  80 FORMAT (' DANGER! NAM E # dtau DA/A DC1 DC2 ',
536  & i6,f8.4,i5,1p,5e10.2)
537  stop
538  END IF
539 *
540  IF (icount.LE.100.AND.abs(ds1).GT.0.01) THEN
541  icount = icount + 1
542  dt = tev(i) - time
543  WRITE (7,85) name(j1), kstar(j1), kstar(j2),
544  & semi*su, da, omeq, spin1, spin2, dt
545  85 FORMAT (' SYNCH NM K* A DA/A om s1 s2 DT ',
546  & i8,2i4,f8.2,1p,5e10.2)
547  CALL flush(7)
548  END IF
549 *
550  IF (spin1.LT.0.0) THEN
551  WRITE (6,90) name(j1), kstar(j1), spin1, time-time0
552  90 FORMAT (' DANGER! NEGATIVE SPIN NM K* s1 T-T0 ',
553  & i6,i4,1p,2e10.2)
554  IF (spin1.LT.0.0) spin1 = 0.1*omeq
555  IF (spin2.LT.0.0) spin2 = 0.1*omeq
556  itry = itry + 1
557  IF (itry.LE.2) go to 50
558 * STOP
559  END IF
560 *
561  100 RETURN
562 *
563  END