Nbody6
 All Files Functions Variables
output.f
Go to the documentation of this file.
1  SUBROUTINE output
2 *
3 *
4 * Output and data save.
5 * ---------------------
6 *
7  include 'common6.h'
8  common/binary/ cm(4,mmax),xrel(3,mmax),vrel(3,mmax),
9  & hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
10  & namem(mmax),nameg(mmax),kstarm(mmax),iflag(mmax)
11  common/galaxy/ gmg,rg(3),vg(3),fg(3),fgd(3),tg,
12  & omega,disk,a,b,v02,rl2,gmb,ar,gam,zdum(7)
13  common/clouds/ xcl(3,mcl),xdotcl(3,mcl),bodycl(mcl),rcl2(mcl),
14  & clm(mcl),clmdot(mcl),cldot,vcl,sigma,rb2,pcl2,
15  & tcl,stepcl,ncl,newcl
16  common/echain/ ech
17  REAL*8 x1(3,4),v1(3,4),ui(4),vi(4),xrel2(3),vrel2(3)
18  REAL*4 xs(3,nmax),vs(3,nmax),bodys(nmax),as(20)
19  REAL*4 xj(3,6),vj(3,6),bodyj(6)
20  LOGICAL first,second,third
21  SAVE first,second,third
22  DATA first,second ,third/.true.,.true.,.true./
23 *
24 *
25 * Obtain energy error in case routine ADJUST not called recently.
26  IF (time.GE.tadj.OR.time.LE.0.0d0) go to 10
27 *
28 * Predict X & XDOT for all particles (except unperturbed pairs).
29  CALL xvpred(ifirst,ntot)
30 *
31 * Obtain the total energy at current time (resolve all KS pairs).
32  CALL energy
33 *
34 * Include KS pairs, triple, quad, mergers, collisions & chain.
35  etot = zkin - pot + etide + ebin + esub + emerge + ecoll + emdot
36  & + ecdot
37  IF (nch.GT.0) THEN
38  etot = etot + ech
39  END IF
40 *
41 * Update energies and form the relative error (divide by ZKIN or ETOT).
42  be(2) = be(3)
43  be(3) = etot
44  de = be(3) - be(2)
45  detot = detot + de
46  de = de/max(zkin,abs(etot))
47 * Save sum of relative energy error for main output and accumulate DE.
48  error = error + de
49  vir = pot - vir
50 *
51 * Find density centre & core radius (Casertano & Hut, Ap.J. 298, 80).
52  IF (n.GE.20.AND.kz(29).EQ.0) THEN
53  CALL core
54  END IF
55 *
56 * Check optional sorting of Lagrangian radii & half-mass radius.
57  IF (kz(7).GT.0) THEN
58  CALL lagr(rdens)
59  END IF
60 *
61 * Initialize diagnostic variables.
62  10 np = 0
63  iunp = 0
64  amin = 100.0
65  mult = 0
66 *
67 * Find smallest semi-major axis and count unperturbed KS pairs.
68  DO 20 ipair = 1,npairs
69  np = np + list(1,2*ipair-1)
70  semi = -0.5*body(n+ipair)/h(ipair)
71  IF (semi.GT.0.0) amin = min(amin,semi)
72  IF (list(1,2*ipair-1).EQ.0) iunp = iunp + 1
73  IF (name(n+ipair).LT.-2*nzero) mult = mult + 1
74  20 CONTINUE
75 *
76 * Include search of any hierarchical triples.
77  DO 25 im = 1,nmerge
78  zmb = cm(1,im) + cm(2,im)
79  semi = -0.5*zmb/hm(im)
80  amin = min(amin,semi)
81  25 CONTINUE
82 *
83 * Perform time-step & neighbour statistics (NS is # single stars).
84  dti = 0.0
85  dtri = 0.0
86  cnnb = 0.0
87  cmax = 0.0
88  nnb = 0
89  ns = 0
90  sum = 0.0
91  DO 30 i = ifirst,ntot
92  dti = dti + 1.0/step(i)
93  dtri = dtri + 1.0/stepr(i)
94  cnnb = cnnb + list(1,i)/step(i)
95  rho = list(1,i)/rs(i)**3
96  cmax = max(cmax,rho)
97  nnb = nnb + list(1,i)
98  IF (i.LE.n.AND.body(i).GT.0.0d0) ns = ns + 1
99  sum = sum + body(i)**2
100  30 CONTINUE
101  ns = ns - nsub
102 *
103 * Estimate relative cost & effective neighbour number of AC scheme.
104  cost = cnnb/(float(n - npairs)*dtri)
105  cnnb = cnnb/dti
106 * Scale maximum particle density contrast by the mean value.
107  cmax = 2.0*cmax*rscale**3/float(n)
108 *
109 * Set average neighbour number & density centre displacement.
110  nnb = float(nnb)/float(n - npairs)
111  rd = sqrt(rdens(1)**2 + rdens(2)**2 + rdens(3)**2)
112 *
113 * Check print frequency indicator & optional model counter.
114  nprint = nprint + 1
115  IF (nprint.GT.nfix.OR.time.LE.0.0d0) THEN
116  nprint = 1
117  IF (kz(3).GT.0) model = model + 1
118  END IF
119 *
120 * Form binary & merger energy ratios and count escapers (or suppress).
121  eb = (ebin + ech)/(zkin - pot)
122  em = emerge/(zkin - pot)
123  IF (kz(21).GT.1) THEN
124  CALL jacobi(nesc)
125  ELSE
126  nesc = 0
127  END IF
128 *
129 * Print main output diagnostics.
130  i6 = tstar*ttot
131 *
132  WRITE (6,40) ttot, n, nnb, npairs, nmerge, mult, ns, nstepi,
133  & nstepb, nstepr, nstepu, error, be(3), zmass
134  40 FORMAT (//,' T =',f7.0,' N =',i6,' <NB> =',i3,' KS =',i5,
135  & ' NM =',i3,' MM =',i2,' NS =',i6,
136  & ' NSTEPS =',i11,i10,2i11,' DE =',1p,e9.1,
137  & ' E =',0p,f10.6,' M =',f7.4)
138 *
139  IF (kz(21).GT.0) THEN
140  CALL cputim(tcomp)
141  IF (vc.EQ.0.0d0) vc = rscale/tcr
142  trc = 1.02*float(nc)**2*bodym/(vc**3*log(float(nc)))
143  dmin1 = min(dmin1, dmin2, dmin3, dmin4, dminc)
144  neff = zmass**2/sum
145 *
146  WRITE (6,45) nrun, model, tcomp, trc, dmin1, dmin2, dmin3,
147  & dmin4, amin, rmax, rsmin, neff
148  45 FORMAT (/,' NRUN =',i3,' M# =',i3,' CPU =',f8.1,' TRC =',
149  & f5.1, ' DMIN =',1p,4e8.1,' AMIN =',e8.1,
150  & ' RMAX =',e8.1,' RSMIN =',0p,f6.3,' NEFF =',i6)
151  END IF
152  vrms = sqrt(2.0*zkin/zmass)*vstar
153 *
154  WRITE (6,50)
155  50 FORMAT (/,' <R> RTIDE RDENS RC NC MC RHOD RHOM',
156  & ' CMAX <Cn> Ir/R UN NP RCM VCM',
157  & ' AZ EB/E EM/E TCR T6 NESC',
158  & ' VRMS')
159 *
160  WRITE (6,55) rscale, rtide, rd, rc, nc, zmc, rhod, rhom, cmax,
161  & cnnb, cost, iunp, np, cmr(4), cmrdot(4), az, eb, em,
162  & tcr, i6, nesc, vrms
163  55 FORMAT (' #1',f5.2,f6.1,f6.2,f7.4,i5,f7.3,f6.0,f7.0,f6.0,f6.1,
164  & f7.3,i5,i6,f8.3,f8.4,f9.4,2f7.3,f6.2,2i6,f6.1)
165 *
166  WRITE (6,60)
167  60 FORMAT (/,7x,'NNPRED NBCORR NBFULL NBVOID NICONV',
168  & ' NBSMIN NBDIS NBDIS2 NCMDER NBDER NFAST',
169  & ' NBFAST NBLOCK NBLCKR NBPRED NBFLUX',
170  & ' NIRECT NRRECT NURECT NBRECT')
171  WRITE (6,65) nnpred, nbcorr, nbfull, nbvoid, niconv,
172  & nbsmin, nbdis, nbdis2, ncmder, nbder, nfast,
173  & nbfast, nblock, nblckr, nbpred, nbflux, nirect,
174  & nrrect, nurect, nbrect
175  65 FORMAT (' #2',i10,i10,i8,i9,i9,i8,i7,2i8,2i7,i8,2i10,2i11,4i8)
176 *
177  WRITE (6,70)
178  70 FORMAT (/,5x,'NKSTRY NKSREG NKSHYP NKSPER NPRECT NEWKS ',
179  & ' NKSMOD NTTRY NTRIP NQUAD NCHAIN NMERG',
180  & ' NEWHI NSTEPT NSTEPQ NSTEPC')
181  WRITE (6,75) nkstry, nksreg, nkshyp, nksper, nprect, newks,
182  & nksmod, nttry, ntrip, nquad, nchain, nmerg, newhi,
183  & nstept, nstepq, nstepc
184  75 FORMAT (' #3',i9,i7,i8,i11,i8,i7,2i9,2i7,i8,2i7,3i8)
185 *
186 * Check output for mass loss or tidal capture.
187  IF (kz(19).GE.3.OR.kz(27).GT.0) THEN
188  CALL events
189  END IF
190 *
191 * Obtain half-mass radii for two groups (NAME <= NZERO/5 & > NZERO/5).
192  IF (kz(7).GE.2) THEN
193  CALL lagr2(rdens)
194  END IF
195 *
196 * Check optional plotting file for main cluster parameters.
197  IF (kz(44).GT.0) THEN
198  WRITE (56,77) tphys, time+toff, rscale*rbar, zmass*smu, ncoll
199  77 FORMAT (' ',1p,4e12.4,0p,i5)
200  CALL flush(56)
201  END IF
202 *
203 * Include diagnostics about cluster orbit in general external field.
204  IF (kz(14).EQ.3) THEN
205  gz = rg(1)*vg(2) - rg(2)*vg(1)
206  sx = rbar/1000.0
207  WRITE (6,78) ntail, (rg(k)*sx,k=1,3), (vg(k)*vstar,k=1,3),
208  & gz, etide
209  78 FORMAT (/,5x,'CLUSTER ORBIT NT RG VG JZ ET ',
210  & i5,3f7.2,2x,3f7.1,1p,e16.8,e10.2)
211  END IF
212  IF (kz(14).EQ.4) THEN
213  WRITE (6,80) ttot, n, rscale, zmass, mp, detot
214  80 FORMAT (/,5x,'GAS EXPULSION T N <R> M MP DETOT ',
215  & f7.1,i7,3f7.3,1p,e10.2)
216  END IF
217 *
218 * Reset minimum encounter distances & maximum apocentre separation.
219  dmin2 = 100.0
220  dmin3 = 100.0
221  dmin4 = 100.0
222  dminc = 100.0
223  rsmin = 100.0
224  rmax = 0.0
225 *
226 * Check integer overflows (2^{32} or 2.1 billion).
227  IF (nstepi.GT.2000000000.OR.nstepi.LT.0) THEN
228  nstepi = 0
229  nirect = nirect + 1
230  END IF
231  IF (nstepr.GT.2000000000.OR.nstepr.LT.0) THEN
232  nstepr = 0
233  nrrect = nrrect + 1
234  END IF
235  IF (nstepu.GT.2000000000.OR.nstepu.LT.0) THEN
236  nstepu = 0
237  nurect = nurect + 1
238  END IF
239  IF (nbpred.GT.2000000000.OR.nbpred.LT.0) THEN
240  nbpred = 0
241  END IF
242  IF (nbflux.GT.2000000000.OR.nbflux.LT.0) THEN
243  nbflux = 0
244  nbrect = nbrect + 1
245  END IF
246  IF (nbcorr.GT.2000000000.OR.nbcorr.LT.0) THEN
247  nbcorr = 0
248  END IF
249  IF (nblock.GT.2000000000.OR.nblock.LT.0) THEN
250  nblock = 0
251  END IF
252 *
253 * Exit if error exceeds restart tolerance (TIME < TADJ means no CHECK).
254  IF (abs(error).GT.5.0*qe.AND.time.LT.tadj) go to 100
255 *
256 * Check optional analysis & output of KS binaries.
257  IF (kz(8).GT.0.AND.npairs.GT.0) THEN
258  CALL binout
259  END IF
260 *
261 * Include optional diagnostics of block-steps.
262  IF (kz(33).GT.0) THEN
263  CALL levels
264  END IF
265 *
266 * Check option for neighbour list histogram (same as for STEPR).
267  IF (kz(33).GT.1) THEN
268  CALL nbhist
269  END IF
270 *
271 * Check optional output of single bodies & binaries.
272  IF (kz(9).GT.0.OR.kz(6).GT.0) THEN
273  CALL bodies
274  END IF
275 *
276 * Include optional diagnostics for interstellar clouds on unit 47.
277  IF (kz(13).GT.0) THEN
278  WRITE (47,82) tstar*ttot, ns, newcl, detot, e(3),
279  & rbar*sqrt(pcl2), rbar*rscale, rbar*rc
280  82 FORMAT (' ',f8.1,i7,i6,2f10.6,f6.2,2f7.3)
281  CALL flush(47)
282  pcl2 = rb2
283  END IF
284 *
285 * See whether to write data bank of binary diagnostics on unit 9.
286  IF (kz(8).GE.2.AND.npairs.GT.0) THEN
287  CALL bindat
288  IF (kz(8).GT.3) THEN
289  CALL hidat
290  END IF
291  END IF
292 *
293 * Check optional diagnostics of evolving stars.
294  IF (kz(12).GT.0.AND.time.GE.tplot) THEN
295  CALL hrplot
296  END IF
297 *
298 * Check optional writing of data on unit 3 (frequency NFIX).
299  IF (kz(3).EQ.0.OR.nprint.NE.1) go to 100
300  IF (kz(3).GT.2.AND.kz(3).NE.5) go to 99
301 *
302  as(1) = ttot
303  as(2) = float(npairs)
304  as(3) = rbar
305  as(4) = zmbar
306  as(5) = rtide
307  as(6) = tidal(4)
308  as(7) = rdens(1)
309  as(8) = rdens(2)
310  as(9) = rdens(3)
311  as(10) = ttot/tcr
312  as(11) = tstar
313  as(12) = vstar
314  as(13) = rc
315  as(14) = nc
316  as(15) = vc
317  as(16) = rhom
318  as(17) = cmax
319  as(18) = rscale
320  as(19) = rsmin
321  as(20) = dmin1
322 *
323 * Include prediction of unperturbed binaries (except ghosts).
324  DO 84 j = 1,npairs
325  IF (list(1,2*j-1).EQ.0.AND.body(n+j).GT.0.0d0) THEN
326  CALL resolv(j,1)
327  END IF
328  84 CONTINUE
329 *
330 * Convert masses, coordinates & velocities to single precision.
331  DO 90 i = 1,ntot
332  bodys(i) = body(i)
333  DO 85 k = 1,3
334  xs(k,i) = x(k,i)
335  vs(k,i) = xdot(k,i)
336  85 CONTINUE
337  90 CONTINUE
338 *
339 * Replace any ghosts by actual M, R & V (including 2 binaries).
340  DO 95 jpair = 1,npairs
341  j2 = 2*jpair
342  j1 = j2 - 1
343  icm = n + jpair
344 * Determine merger & ghost index for negative c.m. name.
345  IF (name(icm).LT.0.AND.body(icm).GT.0.0) THEN
346 * Include possible quartet [[B,S],S] and quintet [[B,S],B] first.
347  IF (name(icm).LT.-2*nzero) THEN
348 * Find ghost and merger index at previous hierarchical level.
349  im = 1
350  DO 92 k = 1,nmerge
351  IF (namem(k).EQ.name(icm) + 2*nzero) im = k
352  92 CONTINUE
353  jg = n
354  DO 94 k = ifirst,ntot
355  IF (nameg(im).EQ.name(k)) jg = k
356  94 CONTINUE
357 * Determine the current ghost and merger index.
358  CALL findj(j1,jg2,im2)
359 * Distinguish netween quartet and quintet.
360  IF (name(j2).LE.nzero) THEN
361  IF (jg.LE.n) THEN
362  bodys(jg) = cm(2,im)
363  bodys(j1) = cm(1,im2)
364  bodys(jg2) = cm(2,im2)
365  ELSE
366  jp = jg - n
367  bodys(2*jp-1) = cm(3,im)
368  bodys(2*jp) = cm(4,im)
369  bodys(jg2) = cm(2,im2)
370  END IF
371  ELSE
372  IF (jg.LE.n) THEN
373  bodys(jg) = cm(2,im)
374  bodys(jg2) = cm(2,im2)
375  ELSE
376  jp = jg - n
377  bodys(2*jp-1) = cm(3,im2)
378  bodys(2*jp) = cm(4,im2)
379  END IF
380  END IF
381  go to 95
382  END IF
383 *
384  CALL findj(j1,j,im)
385 * Note: J is ghost index and IM is merger index.
386  IF (j.LE.0) go to 95
387  bodys(j1) = cm(1,im)
388  bodys(j) = cm(2,im)
389  zmb = cm(1,im) + cm(2,im)
390 * Form global coordinates and velocities from c.m. with XREL & VREL.
391  DO k = 1,3
392  x1(k,1) = x(k,j1) + cm(2,im)*xrel(k,im)/zmb
393  x1(k,2) = x(k,j1) - cm(1,im)*xrel(k,im)/zmb
394  v1(k,1) = xdot(k,j1) + cm(2,im)*vrel(k,im)/zmb
395  v1(k,2) = xdot(k,j1) - cm(1,im)*vrel(k,im)/zmb
396 *
397  xs(k,j1) = x1(k,1)
398  xs(k,j) = x1(k,2)
399  vs(k,j1) = v1(k,1)
400  vs(k,j) = v1(k,2)
401  END DO
402 * See whether ghost is second merged binary (QUAD) instead of triple.
403  IF (j.GT.n) THEN
404  icm2 = j
405 * Save outer binary masses as well as second KS component & ghost c.m.
406  ipair = icm2 - n
407  i1 = 2*ipair - 1
408  i2 = i1 + 1
409  bodys(i1) = cm(3,im)
410  bodys(i2) = cm(4,im)
411  bodys(j2) = cm(2,im)
412  bodys(j) = cm(3,im) + cm(4,im)
413 * Copy KS variables to local scalars.
414  DO k = 1,4
415  ui(k) = u(k,ipair)
416  vi(k) = udot(k,ipair)
417  END DO
418 * Transform to physical variables and multiply by 4 (momentum formula).
419  CALL ksphys(ui,vi,xrel2,vrel2)
420  zm = cm(3,im) + cm(4,im)
421  DO k = 1,3
422  vrel2(k) = 4.0*vrel2(k)
423  x1(k,3) = x(k,j2) + cm(4,im)*xrel2(k)/zm
424  x1(k,4) = x(k,j2) - cm(3,im)*xrel2(k)/zm
425  v1(k,3) = xdot(k,j2) + cm(4,im)*vrel2(k)/zm
426  v1(k,4) = xdot(k,j2) - cm(3,im)*vrel2(k)/zm
427 
428  xs(k,i1) = x1(k,3)
429  xs(k,i2) = x1(k,4)
430  vs(k,i1) = v1(k,3)
431  vs(k,i2) = v1(k,4)
432  xs(k,icm2) = x(k,j2)
433  vs(k,icm2) = xdot(k,j2)
434  END DO
435  END IF
436  END IF
437  95 CONTINUE
438 *
439 * Check modification for chain regularization (case NAME(ICM) = 0).
440  IF (nch.GT.0) THEN
441  CALL chdata(xj,vj,bodyj)
442  DO 98 l = 1,nch
443 * Copy global address from common JLIST (set in CHDATA).
444  j = jlist(l)
445  bodys(j) = bodyj(l)
446  DO 97 k = 1,3
447  xs(k,j) = xj(k,l)
448  vs(k,j) = vj(k,l)
449  97 CONTINUE
450  98 CONTINUE
451  END IF
452 *
453 * Split into WRITE (3) NTOT & WRITE (3) .. if disc instead of tape.
454  IF (first) THEN
455  OPEN (unit=3,status='NEW',form='UNFORMATTED',file='OUT3')
456  first = .false.
457  END IF
458  nk = 20
459  WRITE (3) ntot, model, nrun, nk
460  WRITE (3) (as(k),k=1,nk), (bodys(j),j=1,ntot),
461  & ((xs(k,j),k=1,3),j=1,ntot), ((vs(k,j),k=1,3),j=1,ntot),
462  & (name(j),j=1,ntot)
463 * CLOSE (UNIT=3)
464 *
465 * Produce output file for tidal tail members.
466  99 IF (kz(3).LE.3) THEN
467  IF (second) THEN
468  OPEN (unit=33,status='NEW',form='UNFORMATTED',file='OUT33')
469  second = .false.
470  END IF
471  nk = 13
472  WRITE (33) ntail, model, nk
473 *
474  IF (ntail.GT.0) THEN
475  DO 110 i = itail0,nttot
476  bodys(i) = body(i)
477  DO 105 k = 1,3
478  xs(k,i) = x(k,i) - rg(k)
479  vs(k,i) = xdot(k,i) - vg(k)
480  105 CONTINUE
481  110 CONTINUE
482 * Include cluster centre just in case.
483  DO 115 k = 1,3
484  as(k) = rg(k)
485  as(k+3) = vg(k)
486  as(k+6) = rdens(k)
487  115 CONTINUE
488  as(10) = ttot
489  as(11) = rbar
490  as(12) = tstar
491  as(13) = vstar
492  nk = 13
493  WRITE (33) (as(k),k=1,nk), (bodys(j),j=itail0,nttot),
494  & ((xs(k,j),k=1,3),j=itail0,nttot),
495  & ((vs(k,j),k=1,3),j=itail0,nttot),
496  & (name(j),j=itail0,nttot)
497  END IF
498  END IF
499 *
500 * Include all stars in same file (KZ(3) > 3; astrophysical units).
501  IF (kz(3).GT.3.AND.ntail.GT.0) THEN
502  IF (third) THEN
503  OPEN (unit=34,status='NEW',form='FORMATTED',file='OUT34')
504  third = .false.
505  END IF
506  np = 0
507 * Copy cluster members with respect to density centre.
508  DO 120 i = ifirst,ntot
509  IF (body(i).LE.0.0d0) go to 120
510  np = np + 1
511  DO 118 k = 1,3
512  xs(k,np) = (x(k,i) - rdens(k))*rbar
513  vs(k,np) = xdot(k,i)*vstar
514  118 CONTINUE
515  bodys(np) = body(i)*smu
516  120 CONTINUE
517  n1 = np
518 * Add tidal tail in the same frame.
519  DO 130 i = itail0,nttot
520  np = np + 1
521  DO 125 k = 1,3
522  xs(k,np) = (x(k,i) - rg(k) - rdens(k))*rbar
523  vs(k,np) = (xdot(k,i) - vg(k))*vstar
524  125 CONTINUE
525  bodys(np) = body(i)*smu
526  130 CONTINUE
527  WRITE (34,140) np, n1, (time+toff)*tstar, rbar, vstar,
528  & (rdens(k),k=1,3), (rg(k),k=1,3), (vg(k),k=1,3)
529  140 FORMAT (' ',2i6,f8.1,2f6.2,3f7.3,1p,6e10.2)
530  DO 150 i = 1,np
531  WRITE (34,145) (xs(k,i),k=1,3), (vs(k,i),k=1,3), bodys(i)
532  145 FORMAT (' ',3f10.3,3f8.1,f7.2)
533  150 CONTINUE
534  END IF
535 *
536 * Update next output interval and initialize the corresponding error.
537  100 tnext = tnext + deltat
538  error = 0.0d0
539 *
540  RETURN
541 *
542  END