forked from lruhlen/Original_Peter_Code
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathsph.f
More file actions
13681 lines (11320 loc) · 429 KB
/
sph.f
File metadata and controls
13681 lines (11320 loc) · 429 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
C***********************************************************************
C
C
PROGRAM treesph
C
C
C***********************************************************************
C
C
C A code to evolve self-gravitating fluid flows using smoothed
C particle hydrodynamics and the hierarchical tree method.
C The code is written in standard FORTRAN.
C Nearest neighbor searches are performed once per step, using a
C tree.
C
C In this version, vectorization of the tree walks is achieved by
C simultaneously processing all cells at the same level in the
C tree. The gravitational force calculation and nearest neighbor
C search proceed for a single particle at a time, in serial order.
C
C Smoothing is performed using a cubic spline kernel. The kernel
C is computed by linear interpolation from a look-up table. The
C gravitational field is also smoothed according to the spline
C kernel; however, the gravitational softening parameter is,
C in general, not required to be equal to the SPH smoothing
C length.
C
C In this version the SPH smoothing length is allowed to be
C different for each particle and is determined from the local
C smoothed estimate of the density. Thus, h effectively defines
C a sampling volume for each particle. The smoothing lengths
C are updated by requiring that each SPH particle interact
C roughly with a fixed number of nearest neighbors.
C
C Several different forms of the artificial viscosity are
C available, including a bulk viscosity, a standard form for
C the SPH artificial viscosity, and a version of the SPH
C viscosity, modified by the curl of the velocity field to
C reduce the shear component.
C
C Periodic and quasi-periodic boundary conditions, implemented
C according to Bouchet & Hernquist and Bouchet, Hernquist, &
C Suto are optional.
C
C The gravitational softening length of each particle is
C optionally variable. In the case of SPH particles, the
C softening length is simply equal to the variable smoothing
C length, whereas for collisionless particles softening lengths
C are determined from the requirement that each softening volume
C contain a constant number of near neighbors. The softening
C length for cells is always taken to be a mass-weighted mean
C of the softening lengths for all the particles it contains.
C
C The computational system of units is determined by the input
C data, with the assumption that G=1 . Particles are not
C required to have identical masses.
C
C
C Version 10: November 7, 1988
C
C
C Lars Hernquist, Institute for Advanced Study
C Neal Katz, University of Arizona
C
C
C=======================================================================
C
C
C This is the top-level evolution program treesph. Its tasks are:
C
C 1) to initialize file structures and global variables;
C 2) to input parameters and the initial system state;
C 3) to advance the state of the system for a given number
C of timesteps;
C 4) to perform a diagnostic analysis of the system at
C each time step (energy, angular momentum, etc.);
C 5) to periodically record the state of the system;
C 6) and to terminate the simulation and close data files.
C
C
C=======================================================================
C
C
C Basic global variables/parameters:
C
C acc : acceleration components of a body.
C atriax : parameter defining axis ratios of triaxial
C mass model.
C barmass : mass of optional rotating, rigid bar.
C barsize : radius of optional rotating, rigid bar major
C axis.
C boundary : option to select boundary conditions. Must
C be one of vacuum, periodic, or qperiodic.
C btriax : parameter defining axis ratios of triaxial
C mass model.
C cellsize : linear sizes of cells.
C cputime : cpu time (secs) used during the simulation.
C cputime0 : cumulative cpu time at start of run.
C cputime1 : cumulative cpu time at end of run.
C ctriax : parameter defining axis ratios of triaxial
C mass model.
C dgravsum : option to compute gravity by direct
C summation (.TRUE.).
C dtime : the timestep.
C dtime2 : timestep/2.
C eps : gravitational smoothing parameter.
C epsvect : values of gravitational softening length
C for each particle and each cell.
C ektot : total system kinetic energy.
C eptot : total system gravitational potential energy.
C esofttot : energy correction to virial theorem from
C particle softening.
C etot : total energy of the system.
C fixthalo : option to include a fixed halo contribution
C to the potential.
C fixtriax : option to include fixed triaxial mass model.
C four : the constant 4.
C gamhalo : scale-length of optional, fixed isothermal
C halo.
C headline : identification string for the run.
C incells : number of cells currently in use.
C incellsg : number of cells used to compute gravity.
C inpteps : option to read in values of gravitational
C softening parameter from particle data
C file -- for collisionless particles only.
C intfour : the integer constant 4.
C intone : the integer constant 1.
C intthree : the integer constant 3.
C inttwo : the integer constant 2.
C intzero : the integer constant 0.
C log2 : the constant log10(2).
C mass : masses of bodies and cells.
C mastriax : mass of triaxial mass model.
C mgas : total mass of the gas in the system.
C minusone : the constant -1.
C minustwo : the constant -2.
C mstar : total mass of newly formed stars.
C mtot : total mass of the system.
C nbodies : total number of bodies.
C nbodsmax : maximum number of bodies.
C ncells : maximum number of cells.
C ndim : number of spatial dimensions.
C noutbod : frequency of system record outputs.
C noutlog : frequency of outputs to log file.
C nsat : number of particles in satellite body (used
C only if optional fixed halo is included).
C nsavg : average particle number in softening volume.
C nsmax : maximum particle number in softening volume.
C nsmin : minimum particle number in softening volume.
C nstot : average particle number in softening volume.
C nsteps : total number of timesteps in the simulation.
C nsubcell : number of subcells per cell.
C nsvolume : number of collisionless particles per
C softening volume; used only if variabls
C is .TRUE.
C nsvtol : fractional tolerance in number of neighbors
C relative to nsvolume. A real number
C typically ~ 0.05.
C ntavg : average length of interaction lists.
C ntmax : largest interaction list in current time step.
C ntmin : shortest interaction list in current time step.
C nttot : sum of interaction lists in current time step.
C one : the constant 1.
C onehalf : the constant 0.5.
C patspeed : pattern speed of optional rotating, rigid
C bar potential.
C phi : gravitational potential.
C phiext : contribution to potential from external fields.
C pos : coordinates of bodies, c.m. coords. of cells.
C quad : quadrupole moments of cells.
C restart : option to restart run from a SYSDUMP file.
C rhalo : maximum extent of optional fixed halo.
C rigidbar : option to include rotating, rigid bar in
C force calculation.
C rmin : coords. of lower-left corner of system box.
C root : pointer to the top of the tree.
C rsize : length of the system box.
C selfgrav : option to turn off (.FALSE.) system self-
C gravity.
C snheat : total heating due to supernovae.
C subp : pointers to descendents of a cell.
C tdiebar : time over which the bar diminishes to zero
C strength, reckoned from the starting time
C plus tgrowbar plus tlivebar.
C tgrowbar : time over which bar grows to full strength,
C reckoned from the starting time.
C three : the constant 3.
C tiny : a small number used to prevent divergences.
C tlivebar : time over which bar is at full strength,
C reckoned from the starting time plus tgrowbar.
C tnow : current system time.
C tol : accuracy parameter.
C tol2 : tol * tol.
C tpos : current position time.
C ttree : time of last update of gravitational tree.
C two : the constant 2.
C usequad : option to use (.TRUE.) quadrupole terms.
C usesph : option to use (.TRUE.) SPH calculation.
C variabls : option to use (.TRUE.) variable gravitational
C softening lengths for collisionless particles.
C vel : velocity components of a body.
C vhalo : asymptotic rotation velocity of optional
C fixed halo.
C zero : the constant 0.
C zero02 : the constant 0.02.
C
C-----------------------------------------------------------------------
C
C Definitions specific to input/output.
C
C ireclog : log file record counter.
C uterm, upars, ulog, ubodsin, : logical i/o unit numbers.
C ubodsout,uboddump,utermfil,
C ucrash,uindump
C parsfile, logfile, ibodfile, : character names of files.
C obodfile,dumpfile,termfile,
C crashfil,indumpf
C
C-----------------------------------------------------------------------
C
C Definitions specific to individual particle time steps.
C
C endstep : indicator that a full step is complete.
C etol : tolerance parameter to select new time step.
C inittbin : initial time step bin of all particles.
C itimestp : defines time step of each particle by
C dtime/itimestp.
C mintstep : defines minimum allowed timestep by
C dtime/mintstep.
C npactive : number of particles needing acceleration.
C nsphact : number of SPH particles needing acceleration.
C ntvector : minimum number of particles in each time
C step bin.
C otimestp : previous value of itimestp.
C pactive : indices of particles requiring acceleration.
C stime : fraction of large timestep remaining.
C tsteppos : time step with which to step positions.
C upbin : last time step bin which was updated.
C
C-----------------------------------------------------------------------
C
C Definitions specific to SPH.
C
C acsmooth : table of smoothed gravitational acceleration.
C alpha, beta : coefficients in the artificial viscosity.
C artfvisc : option to select bulk ('bulk'), standard SPH
C ('sph'), or modified SPH artificial viscosity
C ('sphv').
C cfrict : value of frictional damping constant; used
C only if friction is .TRUE.
C consthsm : smoothing length if variablh is .FALSE. A
C value is computed if consthsm is input zero.
C courant : the courant number.
C csound : speed of sound of dissipative particles.
C deldr2i : separation of values in wsmooth.
C dentdt : time derivative of the entropic function a(s).
C dentold : time derivative of the entropic function a(s)
C from the previous time step.
C dethdt : time derivative of the thermal energy.
C dethold : time derivative of the thermal energy from
C the previous time step.
C dwsmooth : table of kernel derivatives.
C entinit : option to initialize entropic function from
C a Jeans mass criterion.
C entold : entropic function a(s) from previous time
C step.
C entropy : the entropic function a(s) for each particle.
C epsfiltr : small, dimensionless parameter used in
C filtering velocities of SPH particles.
C epsgas : gravitational smoothing parameter for gas
C particles.
C epssph : parameter in artificial viscosity.
C ethermal : the specific thermal energy per particle.
C ethinit : option to initialize thermal energy from a
C Jeans mass criterion.
C ethold : thermal energy from previous time step.
C ethtot : total thermal energy of the system.
C friclucy : option to use Lucy's method for frictional
C damping; otherwise a standard frictional term
C is used in the equations of motion.
C friction : option to include (.TRUE.) fricitonal damping
C in equations of motion for SPH particles.
C gamma : ratio of specific heats.
C gamma1 : gamma - 1.
C gamma2 : gamma - 2.
C geogradp : option to use geometric mean form for the
C pressure gradient term (.TRUE.); if .FALSE.,
C the arithmetic mean form is used instead.
C hsmcurlv : h* | curl v | for each particle.
C hsmdivv : h*(div v) for each particle.
C hsmooth : the smoothing length.
C inpmetal : option to input metallicity of gas and stars.
C inptform : option to input formation time of stars.
C isotherm : option to select isothermal eq. of state.
C metals : metallicity for gas and star particles.
C mumaxdvh : maximum muij in SPH viscosity or h*ABS(div v)
C for bulk viscosity.
C ninterp : number of values in look-up tables.
C nsmooth : hsmooth is updated so that each particle
C interacts with nsmooth neighbors.
C nsmtol : fractional tolerance in number of neighbors
C relative to nsmooth. A real number
C typically ~ 0.05.
C nsph : number of dissipational particles.
C nstar : number of star particles.
C nsphmax : maximum number of dissipational particles.
C outputas : option to output entropic function a(s),
C defined by p = a(s) * rho**gamma, to data
C file rather than the temperature (default).
C outputet : option to output thermal energy to data
C file rather than the temperature (default).
C phsmooth : table of look-up values for grav. potential.
C piinv : 1.0 / 3.141592654
C readinas : option to read in entropic function a(s),
C defined by p = a(s) * rho**gamma, from
C data file rather than the thermal energy.
C readinet : option to read in thermal energy from data
C file rather than the temperature (default).
C This parameter should generally be .TRUE.
C if a dimensionless systems of units is used.
C rho : the density.
C sphfiltr : option to apply a spatial filter (.TRUE.)
C to velocities of SPH particles.
C sphinit : option to initialize (.TRUE.) rho and hsmooth
C from input particle positions and nsmooth.
C starform : option to include (.TRUE.) star formation.
C symmetry : option to select Hernquist-Katz ('hk') or
C Benz-Evrard ('be') symmetrized form for
C smoothed estimates. The former is defined by
C ~ 0.5 *(W(h1)+W(h2)), while the latter is
C W(0.5*(h1+h2)).
C tempsphv : temporary storage array for manipulating
C SPH data.
C teth : current thermal energy time.
C tform : formation time of stars.
C thirty5 : the constant 35.
C tvel : current velocity times for SPH particles.
C uentropy : option to integrate the entropic function
C a(s) (.TRUE.) rather than the thermal energy.
C variablg : option to use (.TRUE.) variable gravitational
C softening for SPH particles, using hsmooth.
C variablh : option to select (.TRUE.) variable h.
C veltpos : velocities, time-synchronized with positions.
C wsmooth : table of look-up values of smoothing kernel.
C
C-----------------------------------------------------------------------
C
C Definitions specific to radiative heating and cooling of the gas.
C
C aheatmh : the constant a appearing in the heating rate
C equation, divided by the mass of a hydrogen
C atom, in c.g.s. units.
C bheatmh2 : the constant b appearing in the heating rate
C equation, divided by the square of the mass of
C a hydrogen atom, in c.g.s. units.
C ccoolmh2 : the constant c appearing in the cooling rate
C equation, divided by the square of the mass of
C a hydrogen atom, in c.g.s. units.
C comptmh : the constant appearing in the Compton cooling
C rate equation, divided by the mass of a
C hydrogen atom, in c.g.s. units.
C ctcutoff : cut-off factor in cooling equation. For most
C applications it will lie in the range 0.3 -
C 3.0. (See stepeth.f.)
C derad : rate of change of specific thermal energy
C due to radiative heating and cooling.
C eradiate : total, accumulated thermal energy which has
C been radiated from the system.
C fhydrogn : fraction of gas by weight in hydrogen.
C fhydrog2 : fhydrogn ** 2.
C kpcunit : number of kpc per unit length.
C meanmwt : mean molecular weight of the gas.
C mhboltz : hydrogen mass divided by Boltzmann's
C constant, in c.g.s. units.
C mintemp : minimum temperature where radiative cooling
C is allowed.
C msolunit : number of solar masses per unit mass.
C radiate : option to include (.TRUE.) radiative effects.
C slowcool : maximum fractional cange in a particle's
C thermal energy due to radiative and compton
C cooling.
C trad : current radiation energy time.
C
C-----------------------------------------------------------------------
C
C Definitions specific to vectorized tree construction, vectorized
C tree walk, and vectorized tree search for nearest neighbors.
C
C asubp : subpointers for active bodies or cells.
C bodlist : list of active bodies (i.e. not yet leaves).
C bottom : coordinates of bottom edges of cells.
C celllist : list of cells.
C groupbod : list of bodies in search groups.
C groups : list of grouped cells.
C ingroup : number of particles in groups for searching.
C isubset : indices of subsets of active bodies or cells.
C maxnearb : maximum total number of neighbors allowed.
C nearbods : list of neighbors for all SPH particles.
C ngroups : number of cells containing ingroup particles.
C npercell : number of particles in a cell.
C nnavg : average number of near neighbors.
C nnmax : maximum number of near neighbors.
C nnmin : minimum number of near neighbors.
C nnear : number of neighbors for SPH particles and for
C collisionless particles if variabls is .TRUE.
C nnearlis : number of neighbors in neighbor lists.
C nntot : total number of near neighbors.
C nworkvec : length of temporary work array workvect. It
C should be set to MAX(9*max length of grav
C interaction list, max no. SPH neighbors).
C parent : parents of active bodies or cells.
C pgroupb : pointer to list of bodies in each search group.
C pnear : pointer to list of nearest neighbors.
C subindex : subindices for active bodies or cells.
C subpvect : vector equivalenced to subp.
C templist : temporary vector to swap arrays.
C tempvect : temporary storage vector.
C workvect : temporary work array.
C
C-----------------------------------------------------------------------
C
C Definitions specific to optional rigid black-hole potential.
C
C amvecbh : internal angular momentum of particle system.
C cmposbh : internal center of mass of particles.
C cmvelbh : internal center of mass velocities of
C particles.
C etabh : strength parameter for black-hole collisions,
C defined by etabh = SQRT[(r_p**3/massbh)],
C where r_p is the pericenter distance.
C massbh : mass of black hole if usebh is .TRUE.
C posbh : relative coordinates of optional black hole.
C tstartbh : beginning time of black-hole encounter;
C usually < 0.
C usebh : option to include rigid black-hole potential.
C velbh : relative velocities of optional black hole.
C
C
C-----------------------------------------------------------------------
C
C Definitions specific to cosmological simulations.
C
C comove : option to use (.TRUE.) comoving coordinates.
C cosmfold : value of cosmofac at previous time-step.
C cosmo : option to have (.TRUE.) a cosmological
C simulation.
C cosmofac : parameter set equal to the expansion parameter
C if cosmo and comove are true; one otherwise.
C cosmof3 : cosmofac **3.
C cosmohub : cosmofac * hubble.
C ecosm : energy lost through cosmic expansion, if
C comoving coordinates are used.
C expanpar : current value of the cosmological expansion
C parameter.
C exparout : values of expansion parameter at which to
C output system state.
C hboxsize : 0.5 * pboxsize.
C hubble : current value of the hubble constant, in
C appropriate units.
C hubble0 : today's value of the hubble constant, in
C appropriate units.
C maxnoutc : maximum number of system outputs allowed.
C noutcosm : number of system outputs requested.
C omega0 : for non-periodic boundaries, today's value of
C the cosmological constant omega, otherwise, the
C value of omega at the initial time.
C pboxsize : box size for periodic cosmological simulations.
C redshift : current value of the cosmological redshift
C parameter.
C
C-----------------------------------------------------------------------
C
C Definitions specific to periodic boundary conditions.
C
C alphaew : parameter specific to Ewald summation.
C bcfx(nx,ny,nz) : correction force component at (0,0,0) due
C to a particle at (nx,ny,nz)*(1/(2*nmaxew))
C (direct force is not included).
C bcfxvec : vector equivalenced to bcfx.
C bcfy(nx,ny,nz) : correction force y-component.
C bcfyvec : vector equivalenced to bcfy.
C bcfz(nx,ny,nz) : correction force z-component.
C bcfzvec : vector equivalenced to bcfz.
C bcphi(nx,ny,nz): correction potential energy.
C bcphivec : vector equivalenced to bcphi.
C bwrap : logical check for boundaries to be
C wrapped.
C hmaxew2 : square of maximum radius of Ewald summation
C in reciprocal space.
C maxhew : integerized maximum radius of Ewald sum.
C nmaxew : number of grid points for correction fields.
C nmaxew2 : nmaxew + 2.
C nmaxew2s : nmaxew2 ** 2.
C nmaxew2c : nmaxew2 ** 3.
C nreplica : replica number in real space. Contribution
C from replica particles within radius r_i,jn
C < (nreplica + 0.6) * boxsize is summed.
C pi : the constant 3.141592654...
C radew2 : (nreplica + 0.6) **2.
C sqrtpi : square root of pi.
C
C
C=======================================================================
C
C
C Data structure used to compute gravitational field:
C
C The body/cell tree structure is assumed to be of the
C form discussed by Barnes and Hut. Schematically, for
C three dimensions (i.e. eight subcells per cell):
C
C +-------------------------------------------------+
C root-->| CELL: mass, pos, quad, /, o, /, /, /, /, o, / |
C +----------------------------|--------------|-----+
C | |
C +--------------------------------+ |
C | |
C | +----------------------------------+ |
C +-->| BODY: mass, pos, vel, acc, phi | |
C +----------------------------------+ |
C |
C +-----------------------------------------------+
C |
C | +--------------------------------------------------+
C +-->| CELL: mass, pos, quad, o, /, /, o, /, /, o, / |
C +--------------------------|--------|--------|-----+
C | | |
C etc. etc. etc.
C
C
C The body/cell information is stored in arrays which
C incorporate both bodies and cells. For physical
C quantities relevant to both bodies and cells, such as
C mass and position, the array indices range from
C 1 --> nbodsmax + ncells. For those quantities defined
C only for bodies, such as velocity, the array indices
C range from 1 --> nbodsmax. For information relevant
C to cells alone, such as pointers to descendants, the
C array indices range from nbodsmax + 1 --> nbodsmax +
C ncells. With this convention, pointers can refer to
C either bodies or cells without conflict.
C
C The identification of a given unit as a body or a cell
C is determined by the pointer to the body/cell. For a
C body, p is less than or equal to nbodsmax, while for a
C cell, p > nbodsmax.
C
C For applications using the SPH formalism, the dissipative
C particles are assumed to occupy the first nsph slots.
C
C
C=======================================================================
INCLUDE 'treedefs.h'
C Declaration of local variables.
C -------------------------------
INTEGER n
C=======================================================================
C Initialize state of the system.
C -------------------------------
write(*,*) 'start'
CALL initsys
C -------
C Advance system state for a given number of steps.
C -------------------------------------------------
DO 100 n=1,nsteps
CALL stepsys(n)
C -------
100 CONTINUE
C Terminate the simulation.
C -------------------------
CALL endrun
C ------
STOP
END
C***********************************************************************
C
C
SUBROUTINE initsys
C
C
C***********************************************************************
C
C
C Subroutine to initialize the state of the system.
C
C
C=======================================================================
INCLUDE 'treedefs.h'
C Declaration of local variables.
C -------------------------------
REAL second
C=======================================================================
C Begin timing.
C -------------
cputime0=SECOND()
C-----------------------------------------------------------------------
C Open data files, read input parameters and initial system state,
C and initialize system parameters.
C-----------------------------------------------------------------------
CALL startout
C --------
CALL inparams
C --------
IF(.NOT.restart) THEN
CALL inbods
C ------
ELSE
CALL indump
C ------
ENDIF
CALL checkinp
C --------
CALL initpars
C --------
IF(boundary.EQ.'periodic') CALL initbc
C ------
IF(.NOT.restart) CALL initstep
C --------
IF(.NOT.restart) THEN
C Zero out potential and acceleration.
C ------------------------------------
CALL zeroacc
C -------
CALL zeropot
C -------
C Initialize hydrodynamic properties.
C -----------------------------------
IF(usesph.AND.nsph.GT.0) CALL initsph
C -------
C Initialize gravitational softening lengths.
C -------------------------------------------
IF(variabls.OR.eps.EQ.0.0) CALL maketeps
C --------
IF(.NOT.inpteps) CALL initeps
C -------
IF(variabls) CALL neighcol('correct')
C --------
C Compute gravitational potential and acceleration.
C -------------------------------------------------
CALL gravity('both')
C -------
ENDIF
C Output system state.
C --------------------
IF(.NOT.restart) THEN
CALL outstate(0)
C --------
ELSE
CALL outhead
C -------
ENDIF
IF(.NOT.restart) CALL initpos
C -------
IF(restart) THEN
CALL setbox('sph ')
C ------
CALL loadtree('sph ')
C --------
CALL cellnumb
C --------
CALL celledge
C --------
CALL groupcel
C --------
CALL neighbor('predict')
C --------
ENDIF
RETURN
END
C***********************************************************************
C
C
SUBROUTINE startout
C
C
C***********************************************************************
C
C
C Subroutine to open disk files for subsequent input/output.
C
C
C=======================================================================
INCLUDE 'treedefs.h'
C=======================================================================
C Open log file.
C --------------
OPEN(UNIT=ulog,FILE=logfile,STATUS='NEW',ACCESS='DIRECT',
& FORM='FORMATTED',RECL=81)
ireclog=0
WRITE(ulog,10,REC=1)
10 FORMAT(' Start of logfile output ')
CLOSE(UNIT=ulog)
C Create terminal emulation file.
C -------------------------------
OPEN(UNIT=utermfil,FILE=termfile,STATUS='UNKNOWN')
WRITE(utermfil,*) ' Start of output '
CLOSE(UNIT=utermfil)
RETURN
END
C***********************************************************************
C
C
SUBROUTINE inparams
C
C
C***********************************************************************
C
C
C Subroutine to read in parameters.
C
C Input parameters:
C
C headline : identification string for the run.
C nsteps : number of timesteps.
C noutbod : output system state once every nsteps/noutbod
C steps.
C noutlog : output logfile data once every nsteps/noutlog
C steps.
C dtime : the timestep.
C tol : error tolerance; 0.0 => exact (PP) calculation.
C eps : potential softening parameter.
C inpteps : option to read gravitational softening lengths
C of collisionless particles from input data.
C variabls : option to use variable gravitational softening
C lengths for collisionless particles.
C nsvolume : number of collisionless particles per softening
C volume.
C nsvtol : fractional tolerance in number of neighbors
C relative to nsvolume; typically ~ 0.05.
C usequad : option to include (.TRUE.) quadrupole terms.
C dgravsum : option to use direct summation gravity (.TRUE.).
C boundary : parameter to select boundary condition -- one
C of vacuum, periodic, or qperiodic.
C restart : option to restart run from a SYSDUMP file.
C-----------------------------------------------------------------------
C inittbin : initial time step bin for all particles.
C mintstep : parameter defining minimum allowed time step.
C etol : energy tolerance to select timestep.
C ntvector : minumum number of particles per time step bin.
C-----------------------------------------------------------------------
C usesph : option to use (.TRUE.) an SPH calculation.
C sphinit : option to initialize rho and hsmooth.
C uentropy : option to integrate entropic function a(s)
C (.TRUE.) rather than the thermal energy.
C isotherm : option to use isothermal eq. of state.
C artfvisc : artificial viscosity option parameter.
C epsgas : potential softening parameter for the gas
C particles.
C gamma : ratio of specific heats.
C alpha : coefficient in artificial viscosity.
C beta : coefficient in artificial viscosity.
C epssph : parameter in artificial viscosity.
C courant : courant number.
C variablh : option to have variable smoothing lengths.
C variablg : option to have variable SPH softening lengths.
C consthsm : smoothing length if variablh is .FALSE. A value
C of consthsm is computed if it is zero initially.
C nsmooth : number of neighbors to interact with.
C nsmtol : fractional tolerance in number of neighbors
C relative to nsmooth; typically ~ 0.05.
C symmetry : option to select symmetrized form of smoothed
C estimates.
C geogradp : option to use geometric mean form for pressure
C gradient term; if .FALSE., arithmetic mean used.
C sphfiltr : option to apply a filter to velocities of SPH
C particles.
C epsfiltr : small parameter used in filtering SPH velocities.
C readinas : option to input entropic function, a(s), rather
C than temperature.
C outputas : option to output entropic function, a(s), rather
C than temperature.
C readinet : option to input thermal energy rather than
C temperature.
C outputet : option to output thermal energy rather than
C temperature.
C inpmetal : option to input metallicity for gas and star
C particles.
C inptform : option to input formation time of stars.
C-----------------------------------------------------------------------
C radiate : option to include (.TRUE.) radiative effects.
C aheatmh : heating constant a / hydrogen mass.
C bheatmh2 : heating constant b / hydrogen mass **2.
C ccoolmh2 : cooling constant c / hydrogen mass **2.
C comptmh : Compton cooling constant / hydrogen mass.
C meanmwt : mean molecular weight.
C fhydrogn : fraction of gas by weight in hydrogen.
C mhboltz : hydrogen mass / Boltzmann constant.
C kpcunit : number of kpc per unit length.
C msolunit : number of solar masses per unit length.
C mintemp : minimum temperature for radiative cooling.
C slowcool : maximum fractional change in thermal energy due
C to cooling.
C ctcutoff : coefficient in cooling cut-off expression.
C ethinit : option to initialize (.TRUE.) thermal energy
C from a Jeans mass criterion.
C entinit : option to initialize (.TRUE.) entropic function
C a(s) from a Jeans mass criterion.
C-----------------------------------------------------------------------
C fixthalo : option to have a include a fixed halo potential.
C vhalo : asymptotic rotation velocity of fixed halo.
C gamhalo : scale-length of fixed halo.
C rhalo : maximum extent of fixed halo.
C nsat : number of particles in satellite body.
C-----------------------------------------------------------------------
C fixtriax : option to include fixed triaxial mass model.
C atriax : parameter defining axis ratios of triaxial model.
C btriax : parameter defining axis ratios of triaxial model.
C ctriax : parameter defining axis ratios of triaxial model.
C mastriax : mass of triaxial mass model.
C-----------------------------------------------------------------------
C rigidbar : option to include rotating, rigid bar potential.
C patspeed : pattern speed of rotating, rigid bar.
C barmass : mass of rotating, rigid bar.
C barsize : radius of rotating, rigid bar major axis.
C tgrowbar : time over which bar grows to full strength.
C tlivebar : time over which bar is at full strength.
C tdiebar : time over which bar diminishes to zero strength.
C-----------------------------------------------------------------------
C usebh : option to include rigid black-hole potential.
C massbh : mass of optional black hole.
C etabh : strength parameter for black-hole collisions.
C tstartbh : beginning time for black-hole encounter.
C-----------------------------------------------------------------------
C selfgrav : option to turn off (.FALSE.) system self-gravity.
C starform : option to include star formation.
C friction : option to include frictional damping for SPH
C particles.
C friclucy : option to use Lucy's damping method; otherwise a
C standard frictional term is included.
C cfrict : value of frictional damping constant. Used only
C if friction is .TRUE.
C-----------------------------------------------------------------------
C omega0 : current value of cosmological omega for vacuum
C b.c., initial value of omega for periodic b.c.
C hubble0 : current value of the hubble constant.
C cosmo : option to have a cosmological simulation.
C comove : option to have comoving coordinates.
C ecosm : initial value of energy lost through cosmic
C expansion (should be zero except for restarts).
C pboxsize : length of periodic box.
C noutcosm : number of system outputs for cosmological
C simulations.
C exparout : expansion parameters at which to output state
C of the system.
C
C
C=======================================================================
INCLUDE 'treedefs.h'
C Declaration of local variables.
C -------------------------------
CHARACTER *1 pcomment
INTEGER i
C=======================================================================
OPEN(UNIT=upars,FILE=parsfile,STATUS='OLD')
C Read parameters, close the file.
C --------------------------------
READ(upars,'(a)') pcomment
READ(upars,'(a)') headline
READ(upars,*) nsteps
READ(upars,*) noutbod
READ(upars,*) noutlog
READ(upars,*) dtime
READ(upars,*) tol
READ(upars,*) eps
READ(upars,*) inpteps
READ(upars,*) variabls
READ(upars,*) nsvolume
READ(upars,*) nsvtol
READ(upars,*) usequad
READ(upars,*) dgravsum
READ(upars,'(a)') boundary
READ(upars,*) restart
READ(upars,'(a)') pcomment
READ(upars,*) inittbin
READ(upars,*) mintstep
READ(upars,*) etol
READ(upars,*) ntvector
READ(upars,'(a)') pcomment
READ(upars,*) usesph
READ(upars,*) sphinit
READ(upars,*) uentropy
READ(upars,*) isotherm
READ(upars,'(a)') artfvisc
READ(upars,*) epsgas
READ(upars,*) gamma
READ(upars,*) alpha
READ(upars,*) beta
READ(upars,*) epssph
READ(upars,*) courant
READ(upars,*) variablh
READ(upars,*) variablg
READ(upars,*) consthsm
READ(upars,*) nsmooth
READ(upars,*) nsmtol
READ(upars,'(a)') symmetry
READ(upars,*) geogradp
READ(upars,*) sphfiltr
READ(upars,*) epsfiltr
READ(upars,*) readinas
READ(upars,*) outputas
READ(upars,*) readinet
READ(upars,*) outputet
READ(upars,*) inpmetal
READ(upars,*) inptform
READ(upars,'(a)') pcomment
READ(upars,*) radiate
READ(upars,*) aheatmh
READ(upars,*) bheatmh2
READ(upars,*) ccoolmh2
READ(upars,*) comptmh
READ(upars,*) meanmwt
READ(upars,*) fhydrogn
READ(upars,*) mhboltz
READ(upars,*) kpcunit
READ(upars,*) msolunit
READ(upars,*) mintemp
READ(upars,*) slowcool
READ(upars,*) ctcutoff
READ(upars,*) ethinit
READ(upars,*) entinit
aheatmh=aheatmh*3.46e14*kpcunit**2.5/msolunit**1.5
bheatmh2=bheatmh2*2.34e-17/SQRT(kpcunit*msolunit)
ccoolmh2=ccoolmh2*2.34e-17/SQRT(kpcunit*msolunit)
mhboltz=mhboltz*4.30e4*msolunit/kpcunit
comptmh=comptmh*3.46e14*kpcunit**2.5/msolunit**1.5
READ(upars,'(a)') pcomment
READ(upars,*) fixthalo
READ(upars,*) vhalo
READ(upars,*) gamhalo
READ(upars,*) rhalo
READ(upars,*) nsat
READ(upars,'(a)') pcomment
READ(upars,*) fixtriax
READ(upars,*) atriax
READ(upars,*) btriax
READ(upars,*) ctriax
READ(upars,*) mastriax
READ(upars,'(a)') pcomment