Nbody6
 All Files Functions Variables
subint.f
Go to the documentation of this file.
1  SUBROUTINE subint(IQ,I10)
2 *
3 *
4 * Decision-making for subsystems.
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  REAL*8 tslist(kmax)
11  SAVE irun,li
12  DATA irun /0/
13 *
14 *
15 * Determine correct index after restart (NNTB = 0 initially).
16  IF (irun.EQ.0) THEN
17  irun = 1
18  ti = 1.0d+10
19 * Find smallest sum by looping backwards (avoids multiple entries).
20  DO 4 k = nntb,1,-1
21  j = kblist(k)
22  tj = t0(j) + step(j)
23  IF (tj.LT.ti) THEN
24  ti = tj
25  li = k
26  ELSE
27 * Adopt the previous index (increased by 1 below).
28  li = li - 1
29  go to 1
30  END IF
31  4 CONTINUE
32  li = li - 1
33  dtb = 0.0
34  END IF
35 *
36 * See whether to advance any KS solutions at start of block-step.
37  1 IF (npairs.GT.0) THEN
38 * Obtain list of all KS pairs due in interval DTB.
39  IF (tblist.LE.tblock.OR.npairs.NE.nbprev) THEN
40  IF (dtb.EQ.0.0d0.OR.dtb.GT.1.0d+06) THEN
41  dtb = max(dtmin,tblock - tprev)
42  END IF
43  2 tblist = tprev + dtb
44  tblist = max(tblock,tblist)
45  nntb = 0
46  DO 3 jpair = 1,npairs
47  j1 = 2*jpair - 1
48  IF (t0(j1) + step(j1).LE.tblist) THEN
49  nntb = nntb + 1
50  kblist(nntb) = j1
51  tslist(nntb) = t0(j1) + step(j1)
52  END IF
53  3 CONTINUE
54 * Increase interval on zero membership.
55  IF (nntb.EQ.0) THEN
56  dtb = 2.0*dtb
57  go to 2
58  END IF
59 * Stabilize interval on membership of 2*SQRT(NPAIRS).
60  nbtry = 2*sqrt(float(npairs))
61  IF (nntb.GT.nbtry) dtb = 0.75*dtb
62  IF (nntb.LT.nbtry) dtb = 1.25*dtb
63 * Sort the time-step list sequentially in KBLIST and reset pointer.
64  IF (nntb.GT.1) THEN
65  CALL sort1(nntb,tslist,kblist)
66  END IF
67  li = 0
68  END IF
69 *
70 * Select members of KBLIST sequentially and set new time.
71  5 li = li + 1
72  IF (li.GT.nntb) THEN
73 * Enforce KBLIST at TBLOCK in case TIME > TBLOCK (unpert to unpert).
74  tblist = tblock
75  go to 1
76  END IF
77  i1 = kblist(li)
78  time = t0(i1) + step(i1)
79 *
80 * See whether the smallest KS time falls due before next block-step.
81  IF (time.LE.tblock) THEN
82  10 CALL ksint(i1)
83 *
84 * Check for multiple calls of #I1 (saves using CALL INSERT).
85  IF (li.LT.nntb.AND.iphase.EQ.0) THEN
86  ti = time + step(i1)
87  jx = kblist(li+1)
88  tx = t0(jx) + step(jx)
89  IF (ti.LT.tx.AND.ti.LE.tblock) THEN
90  time = ti
91  go to 10
92  END IF
93  END IF
94 *
95 * See whether current pair is due before new KBLIST loop.
96  IF (time + step(i1).LT.tblist) THEN
97 *
98 * Form new list if expanded size is too big.
99  IF (nntb.GE.kmax-5) THEN
100  tblist = time
101 * Include test on non-zero indicator (bad luck case).
102  IF (iphase.NE.0.AND.iq.EQ.0) THEN
103  iq = iphase
104  i10 = i1
105  iphase = 0
106  END IF
107  go to 1
108  END IF
109 *
110 * Insert body #I1 in the correct sequential location (skip T0 < TIME).
111  IF (t0(i1).GE.time) THEN
112  CALL insert(i1,li)
113  END IF
114  END IF
115 *
116 * Set KS indicator on termination, multiple regularization or merger.
117  IF (iphase.NE.0) THEN
118  IF (iq.EQ.0.OR.iphase.LT.0) THEN
119  iq = iphase
120 * Save KS index until exit (collision treated in situ).
121  IF (iq.GT.0) THEN
122  i10 = i1
123  END IF
124  END IF
125 *
126 * Reset non-zero decision indicator (continue on positive value).
127  IF (iphase.GT.0) THEN
128  iphase = 0
129  ELSE
130 * Enforce new sorted list on change of KS sequence after collision.
131  iphase = 0
132  tblist = time
133  go to 1
134  END IF
135  END IF
136 *
137 * Continue cycle until end of block-step.
138  go to 5
139  END IF
140 *
141 * Copy original block time at end of KS treatment.
142  time = tblock
143  nbprev = npairs
144 * Reduce pointer by 1 for next block-step (otherwise not done).
145  li = li - 1
146  END IF
147 *
148 * Check time for advancing any triple, quad or chain regularization.
149  IF (nsub.GT.0) THEN
150  30 tsub = 1.0d+10
151  DO 40 l = 1,nsub
152  IF (ts(l).LT.tsub) THEN
153  isub = l
154  tsub = ts(l)
155  END IF
156  40 CONTINUE
157 *
158  IF (tsub.LE.tblock) THEN
159  time = tsub
160 * Decide between triple, quad or chain.
161  IF (isys(isub).EQ.1) THEN
162 * Update unperturbed size of subsystem and copy c.m. step.
163  CALL extend(isub)
164  CALL triple(isub)
165  ELSE IF (isys(isub).EQ.2) THEN
166  CALL extend(isub)
167  CALL quad(isub)
168  ELSE
169  IF (steps(isub).LT.0.0d0) THEN
170  steps(isub) = 1.0d-10
171  go to 50
172  END IF
173  CALL chain(isub)
174  IF (isub.GT.0) THEN
175  IF (steps(isub).LT.0.0d0) THEN
176  steps(isub) = 1.0d-10
177  go to 50
178  END IF
179  END IF
180  END IF
181 *
182 * Check for termination (set TPREV < TIME and set IQ < 0).
183  IF (isub.LT.0.OR.iphase.LT.0) THEN
184  tprev = time - step(ntot)
185  iq = -1
186  END IF
187  go to 30
188  END IF
189  50 time = tblock
190  END IF
191 *
192  RETURN
193 *
194  END