Nbody6
 All Files Functions Variables
intgrt.f
Go to the documentation of this file.
1  SUBROUTINE intgrt
2 *
3 *
4 * N-body integrator flow control.
5 * -------------------------------
6 *
7  include 'common6.h'
8  common/clump/ bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
9  & names(ncmax,5),isys(5)
10  common/chainc/ xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
11  & listc(lmax)
12  REAL*8 xi(3),xidot(3)
13  INTEGER nxtlst(nmax),ibl(lmax),nblist(nmax),listq(nmax),nl(20)
14  LOGICAL loop,lstepm
15  SAVE iq,icall,nq,lq,loop,lstepm,stepm,isave,jsave
16  DATA iq,icall,lq,loop,lstepm,stepm /0,2,11,.true.,.false.,0.03125/
17 *
18 *
19 * Enforce level search on return, except new and terminated KS.
20  IF (iphase.NE.1.AND.iphase.NE.2) loop = .true.
21 *
22 * Update quantized value of STEPM for large N (first time only).
23  IF (.NOT.lstepm.AND.nzero.GT.1024) THEN
24  k = (float(nzero)/1024.0)**0.333333
25  stepm = 0.03125d0/2**(k-1)
26  lstepm = .true.
27  END IF
28 *
29 * Search for high velocities after escape or KS/chain termination.
30  999 IF (kz(37).GT.0.AND.(iphase.EQ.-1.OR.iphase.GE.2)) THEN
31  CALL hivel(0)
32  END IF
33 *
34 * Reset control & regularization indicators.
35  iphase = 0
36  iks = 0
37  dtm = 1.0
38  tprev = time
39 * Initialize end-point of integration times and set DTM.
40  DO 1000 i = ifirst,ntot
41  tnew(i) = t0(i) + step(i)
42  dtm = min(dtm,step(i))
43  1000 CONTINUE
44 *
45 * Determine level for the smallest step (ignore extreme values).
46  lqs = 20
47  DO 1001 l = 6,20
48  IF (dtm.EQ.dtk(l)) THEN
49  lqs = l
50  END IF
51  1001 CONTINUE
52 *
53 * Specify upper level for optimized membership.
54  lqb = lqs - 4
55  IF (iq.LT.0) icall = 0
56  iq = 0
57 * Enforce new block step search initially and on significant change.
58  tlistq = time
59 *
60 * Check updating new list of block steps with T0 + STEP =< TLISTQ.
61  1 icall = icall + 1
62 * Reset TMIN second & third time after change to catch new chain step.
63  IF (time.GE.tlistq.OR.icall.LE.3) THEN
64 * Update interval by optimization at major times (sqrt of N-NPAIRS).
65  IF (dmod(tlistq,2.0d0).EQ.0.0d0.OR.loop) THEN
66  loop = .false.
67  DO 10 l = 1,20
68  nl(l) = 0
69  10 CONTINUE
70  DO 14 i = ifirst,ntot
71 * Count steps at five different levels for the smallest values.
72  DO 12 l = lqb,lqs
73  IF (step(i).LT.dtk(l)) nl(l) = nl(l) + 1
74  12 CONTINUE
75  14 CONTINUE
76  nlsum = 0
77 * Determine interval by summing smallest steps until near sqrt(N-N_b).
78  nsq = sqrt(float(n - npairs))
79  lq = lqs
80  DO 15 l = lqs,lqb,-1
81  nlsum = nlsum + nl(l)
82  IF (nlsum.LE.nsq) lq = l
83  15 CONTINUE
84 * WRITE (6,16) TIME+TOFF,NQ,NLSUM,LQ,(NL(K),K=LQB,LQS)
85 * 16 FORMAT (' LEVEL CHECK: T NQ NLSUM LQ NL ',
86 * & F9.3,3I5,2X,7I4)
87  END IF
88 *
89 * Increase interval by optimized value.
90  nq = 0
91  tmin = 1.0d+10
92  18 tlistq = tlistq + dtk(lq)
93  DO 20 i = ifirst,ntot
94  IF (tnew(i).LE.tlistq) THEN
95  nq = nq + 1
96  listq(nq) = i
97  tmin = min(tnew(i),tmin)
98  END IF
99  20 CONTINUE
100 * Increase interval in rare case of zero membership.
101  IF (nq.EQ.0) go to 18
102 * Make a slight adjustment for high levels and small membership.
103  IF (lq.LE.15.AND.nq.LE.2) lq = lq - 1
104  END IF
105 *
106 * Find all particles in next block (TNEW = TMIN).
107  CALL inext(nq,listq,tmin,nxtlen,nxtlst)
108 *
109 * Set new time and save block time (for regularization terminations).
110  i = nxtlst(1)
111  time = t0(i) + step(i)
112  tblock = time
113  li = 0
114  ipred = 0
115 *
116 * ID = 0
117 * IF (NSTEPI.EQ.1008936) ID = 1
118 * IF (NSTEPI.GE.1008764) THEN
119 * WRITE (6,22) NXTLEN, NSTEPI, I, NAME(I), STEP(I), STEPR(I)
120 * 22 FORMAT (' NEXT LEN # I NM DT DTR ',I5,I9,2I6,1P,2E10.2)
121 * CALL FLUSH(6)
122 * END IF
123 * Re-determine list if current time exceeds boundary.
124  IF (time.GT.tlistq) go to 1
125 *
126 * Check option for advancing interstellar clouds.
127  IF (kz(13).GT.0) THEN
128  CALL clint
129  END IF
130 *
131 * Check optional integration of cluster guiding centre.
132  IF (kz(14).EQ.3.OR.kz(14).EQ.4) THEN
133  IF (kz(14).EQ.3.AND.dmod(time,stepx).EQ.0.0d0) THEN
134  CALL gcint
135  END IF
136 * Include mass loss by gas expulsion (Kroupa et al. MN 321, 699).
137  IF (mpdot.GT.0.0d0.AND.time + toff.GT.tdelay) THEN
138  CALL plpot1(phi1)
139  mp = mp0/(1.0 + mpdot*(time + toff - tdelay))
140 * Replace by exponential mass loss for faster decrease.
141 * DT = TIME + TOFF - TDELAY
142 * MP = MP0*EXP(-MPDOT*DT)
143  CALL plpot1(phi2)
144 * Add differential correction for energy conservation.
145  emdot = emdot + (phi1 - phi2)
146  END IF
147  END IF
148 *
149 * Include commensurability test (may be suppressed if no problems).
150 * IF (DMOD(TIME,STEP(I)).NE.0.0D0) THEN
151 * WRITE (6,25) I, NAME(I), NSTEPI, TIME, STEP(I), TIME/STEP(I)
152 * 25 FORMAT (' DANGER! I NM # TIME STEP T/DT ',
153 * & 2I6,I11,F12.5,1P,E9.1,0P,F16.4)
154 * STOP
155 * END IF
156 *
157 * Check for new regularization at end of block.
158  IF (iks.GT.0) THEN
159  iphase = 1
160 * Copy the saved component indices and time.
161  icomp = isave
162  jcomp = jsave
163  time = tsave
164  go to 100
165  END IF
166 *
167 * Check next adjust time before beginning a new block.
168  IF (time.GT.tadj) THEN
169  time = tadj
170  iphase = 3
171  go to 100
172  END IF
173 *
174 * Check output time in case DTADJ & DELTAT not commensurate.
175  IF (time.GT.tnext) THEN
176  time = tnext
177  CALL output
178  go to 1
179  END IF
180 *
181 * See whether to advance ARchain or KS at first new time.
182  IF (time.GT.tprev) THEN
183  CALL subint(iq,i10)
184  IF (iq.LT.0) go to 999
185  END IF
186 *
187 * Check regular force condition for small block memberships.
188  ir = 0
189  IF (nxtlen.LE.10) THEN
190  DO 28 l = 1,nxtlen
191  j = nxtlst(l)
192  IF (tnew(j).GE.t0r(j) + stepr(j)) THEN
193  ir = ir + 1
194  END IF
195  28 CONTINUE
196  END IF
197 *
198 * Decide between merging of neighbour lists or full N prediction.
199  IF (nxtlen.LE.10.AND.ir.EQ.0) THEN
200 *
201 * Initialize pointers for neighbour lists.
202  DO 30 l = 1,nxtlen
203  ibl(l) = nxtlst(l)
204  30 CONTINUE
205 *
206 * Merge all neighbour lists (with absent members of IBL added).
207  CALL nbsort(nxtlen,ibl,nnb,nblist)
208 *
209 * Predict coordinates & velocities of neighbours and #I to order FDOT.
210  nbpred = nbpred + nnb
211  DO 35 l = 1,nnb
212  j = nblist(l)
213  s = time - t0(j)
214  s1 = 1.5*s
215  s2 = 2.0*s
216  x(1,j) = ((fdot(1,j)*s + f(1,j))*s +x0dot(1,j))*s +x0(1,j)
217  x(2,j) = ((fdot(2,j)*s + f(2,j))*s +x0dot(2,j))*s +x0(2,j)
218  x(3,j) = ((fdot(3,j)*s + f(3,j))*s +x0dot(3,j))*s +x0(3,j)
219  xdot(1,j) = (fdot(1,j)*s1 + f(1,j))*s2 + x0dot(1,j)
220  xdot(2,j) = (fdot(2,j)*s1 + f(2,j))*s2 + x0dot(2,j)
221  xdot(3,j) = (fdot(3,j)*s1 + f(3,j))*s2 + x0dot(3,j)
222  35 CONTINUE
223 * Ensure prediction of chain c.m. not on block-step (may be needed).
224  IF (nch.GT.0) THEN
225  IF (tnew(ich).GT.tblock) THEN
226  CALL xvpred(ich,0)
227  END IF
228  END IF
229  ELSE
230  ipred = 1
231  nnpred = nnpred + 1
232  DO 40 j = ifirst,ntot
233  IF (body(j).EQ.0.0d0) go to 40
234  s = time - t0(j)
235  s1 = 1.5*s
236  s2 = 2.0*s
237  x(1,j) = ((fdot(1,j)*s + f(1,j))*s +x0dot(1,j))*s +x0(1,j)
238  x(2,j) = ((fdot(2,j)*s + f(2,j))*s +x0dot(2,j))*s +x0(2,j)
239  x(3,j) = ((fdot(3,j)*s + f(3,j))*s +x0dot(3,j))*s +x0(3,j)
240  xdot(1,j) = (fdot(1,j)*s1 + f(1,j))*s2 + x0dot(1,j)
241  xdot(2,j) = (fdot(2,j)*s1 + f(2,j))*s2 + x0dot(2,j)
242  xdot(3,j) = (fdot(3,j)*s1 + f(3,j))*s2 + x0dot(3,j)
243  40 CONTINUE
244  END IF
245 *
246 * Resolve perturbed KS pairs with c.m. prediction after NBSORT.
247  jj = -1
248  DO 45 jpair = 1,npairs
249  jj = jj + 2
250  IF (list(1,jj).GT.0) THEN
251 * Ignore c.m. prediction after full N loop (all active KS needed).
252  IF (ipred.EQ.0) THEN
253  j = n + jpair
254  s = time - t0(j)
255  s1 = 1.5*s
256  s2 = 2.0*s
257  x(1,j) = ((fdot(1,j)*s + f(1,j))*s +x0dot(1,j))*s +x0(1,j)
258  x(2,j) = ((fdot(2,j)*s + f(2,j))*s +x0dot(2,j))*s +x0(2,j)
259  x(3,j) = ((fdot(3,j)*s + f(3,j))*s +x0dot(3,j))*s +x0(3,j)
260  xdot(1,j) = (fdot(1,j)*s1 + f(1,j))*s2 + x0dot(1,j)
261  xdot(2,j) = (fdot(2,j)*s1 + f(2,j))*s2 + x0dot(2,j)
262  xdot(3,j) = (fdot(3,j)*s1 + f(3,j))*s2 + x0dot(3,j)
263  END IF
264  zz = 1.0
265 * Distinguish between low and high-order prediction of U & UDOT.
266  IF (gamma(jpair).GT.1.0d-04) zz = 0.0
267  CALL ksres2(jpair,j1,j2,zz)
268  END IF
269  45 CONTINUE
270 *
271 * Save new time (output time at TIME > TADJ) and increase # blocks.
272  tprev = time
273  nblock = nblock + 1
274  tmin = 1.0d+10
275 *
276 * Predict chain variables and perturber cordinates at new block-time.
277  IF (nch.GT.0) THEN
278  CALL xcpred(2)
279  END IF
280 *
281 * Advance the pointer (<= NXTLEN) and select next particle index.
282  50 li = li + 1
283  IF (li.GT.nxtlen) go to 1
284  i = nxtlst(li)
285  time = t0(i) + step(i)
286 *
287 * See whether the regular force needs to be updated (IR > 0).
288  IF (t0r(i) + stepr(i).LE.time) THEN
289  ir = 1
290  ELSE
291  ir = 0
292  END IF
293 *
294 * Advance the irregular step.
295  iks0 = iks
296  CALL nbint(i,iks,ir,xi,xidot)
297 *
298 * Save indices of first KS candidates in the block and TIME.
299  IF (iks0.EQ.0.AND.iks.GT.0) THEN
300  isave = icomp
301  jsave = jcomp
302  tsave = time
303  END IF
304 *
305 * See whether the regular step is due.
306  IF (ir.GT.0) THEN
307  CALL regint(i,xi,xidot)
308  END IF
309 *
310 * Determine next block time (note STEP may shrink in REGINT).
311  tmin = min(tnew(i),tmin)
312 *
313 * Copy current coordinates & velocities from corrected values.
314  IF (li.EQ.nxtlen) THEN
315  DO 60 l = 1,nxtlen
316  i = nxtlst(l)
317  DO 55 k = 1,3
318  x(k,i) = x0(k,i)
319  xdot(k,i) = x0dot(k,i)
320  55 CONTINUE
321  60 CONTINUE
322 *
323 * Check integration of tidal tail members.
324  IF (ntail.GT.0) THEN
325 * Allow large quantized interval with internal iteration.
326  IF (dmod(time,0.25d0).EQ.0.0d0) THEN
327  DO 65 j = itail0,nttot
328  IF (tnew(j).LE.time) THEN
329  CALL ntint(j)
330  END IF
331  65 CONTINUE
332  END IF
333  END IF
334  ELSE
335 * Continue until last member has been done (improves reproducibility).
336  go to 50
337  END IF
338 *
339 * Exit on KS termination, new multiple regularization or merger.
340  IF (iq.NE.0) THEN
341  nbprev = 0
342  IF (iq.GE.4.AND.iq.NE.7) THEN
343  CALL delay(iq,-1)
344  ELSE
345 * Ensure correct KS index (KSPAIR may denote second termination).
346  kspair = kvec(i10)
347  iphase = iq
348  END IF
349  go to 100
350  END IF
351 *
352 * Perform optional check on high-velocity particles at major times.
353  IF (kz(37).GT.0.AND.listv(1).GT.0) THEN
354  IF (dmod(time,stepm).EQ.0.0d0) THEN
355  CALL shrink(tmin)
356  IF (listv(1).GT.0) THEN
357  CALL hivel(-1)
358  END IF
359  END IF
360  END IF
361 *
362 * Check optional mass loss time at end of block-step.
363  IF (kz(19).GT.0) THEN
364 * Delay until time commensurate with 100-year step (new polynomials).
365  IF (time.GT.tmdot.AND.dmod(time,stepx).EQ.0.0d0) THEN
366  IF (kz(19).GE.3) THEN
367  CALL mdot
368  ELSE
369  CALL mloss
370  END IF
371  IF (iphase.LT.0) go to 999
372  END IF
373  END IF
374 *
375 * Advance counters and check timer & optional COMMON save (NSUB = 0).
376  ntimer = ntimer + 1
377  IF (ntimer.LT.nmax) go to 50
378 
379  ntimer = 0
380  nsteps = nsteps + nmax
381 *
382  IF (nsteps.GE.100*nmax.AND.nsub.EQ.0) THEN
383  nsteps = 0
384  IF (kz(1).GT.1) CALL mydump(1,1)
385  END IF
386 *
387 * Check option for general binary search.
388  IF (kz(4).GT.0.AND.time - tlasts.GT.deltas) THEN
389  CALL evolve(0,0)
390  END IF
391 *
392 * Include facility for termination of run (create dummy file STOP).
393  OPEN (99,file='STOP',status='OLD',form='FORMATTED',iostat=io)
394  IF (io.EQ.0) THEN
395  CLOSE (99)
396  IF (nsub.EQ.0) WRITE (6,70)
397  70 FORMAT (/,9x,'TERMINATION BY MANUAL INTERVENTION')
398  cpu = 0.0
399  END IF
400 *
401 * Repeat cycle until elapsed computing time exceeds the limit.
402  CALL cputim(tcomp)
403  IF (tcomp.LT.cpu) go to 50
404 *
405 * Do not terminate during triple, quad or chain regularization.
406  IF (nsub.GT.0) THEN
407 * Specify zero step to enforce termination.
408  DO 75 l = 1,nsub
409  steps(l) = 0.0d0
410  75 CONTINUE
411  ntimer = nmax
412  go to 50
413  END IF
414 *
415 * Terminate run with optional COMMON save.
416  IF (kz(1).GT.0) THEN
417  cputot = cputot + tcomp - cpu0
418  CALL mydump(1,1)
419  WRITE (6,80) time+toff, tcomp, cputot/60.0, errtot, detot
420  80 FORMAT (/,9x,'COMMON SAVED AT TIME =',f8.2,' TCOMP =',f7.1,
421  & ' CPUTOT =',f6.1,' ERRTOT =',f10.6,
422  & ' DETOT =',f10.6)
423  END IF
424 *
425  stop
426 *
427 * Set current global time.
428  100 ttot = time + toff
429 *
430  RETURN
431 *
432  END