James Bordner avatar James Bordner committed cb44c5a Merge

[lcaperf] merged with week-of-code

Comments (0)

Files changed (192)

 DEPEND.bak
 auto_show_*.C
 enzo.exe
+inits.exe
+ring.exe
+anyl.exe
+enzohop.exe
 out.compile
 out.make.DEPEND
 svn_version.def
 Make.config.override
 *.mod
 
-src/inits/inits.exe
-src/anyl/anyl.exe
-src/enzohop/enzohop.exe
 src/mpgrafic/degraf/Makefile
 src/mpgrafic/degraf/config.h
 src/mpgrafic/degraf/config.log

doc/manual/source/EnzoParameters.rst

     are made with the intention that only the cooling rates are to be
     used. Default: 0 (off).
 
-``IncludeCloudyMMW`` (external)
-    An integer (0 or 1) specifying whether the additional mean
-    molecular weight contributed by the metals be used in the
-    conversion from internal energy to temperature. These values will
-    come from the Cloudy dataset. For metallicities less than solar,
-    this addition will be negligible. Default: 0 (off).
-
 ``CMBTemperatureFloor`` (external)
     An integer (0 or 1) specifying whether a temperature floor is
     created at the temperature of the cosmic microwave background
     code by subtracting the cooling rate at T\ :sub:`CMB`\  such that
     Cooling = Cooling(T) - Cooling(T\ :sub:`CMB`\ ). Default: 1 (on).
 
-``CloudyMetallicityNormalization`` (external)
-    A float value used in the conversion of metal density into
-    metallicity. This value will change depending on the specific
-    abundance patterns used to make the Cloudy dataset. The value of
-    this factor is calculated as the sum of (A\ :sub:`i`\  \*
-    m\ :sub:`i`\ ) over all elements i heavier than He, where
-    A\ :sub:`i`\  is the solar number abundance relative to H and
-    m\ :sub:`i`\  is the atomic mass. For the solar abundance pattern
-    from the latest version of Cloudy, using all metals through Zn,
-    this value is 0.018477. Default: 0.018477.
-
 ``CloudyElectronFractionFactor`` (external)
     A float value to account for additional electrons contributed by
     metals. This is only used with Cloudy datasets with dimension
     or 1.  See :ref:`distributed_feedback` for an illustration.
     Default: 0.
 
+``StarMakerTypeIaSNe`` (external)
+    This parameter turns on thermal and chemical feedback from Type Ia
+    supernovae.  The mass loss and luminosity of the supernovae are
+    determined from `fits of K. Nagamine
+    <http://www.physics.unlv.edu/~kn/SNIa_2/>`_.  The ejecta are
+    traced in a separate species field, ``MetalSNIa_Density``.  The
+    metallicity of star particles that comes from this ejecta is
+    stored in the particle attribute ``typeia_fraction``.  Can be used
+    with ``StarParticleCreation`` = 0, 1, 2, 5, 7, and 8.  Default:
+    0.
+
+``StarMakerPlanetaryNebulae`` (external) 
+    This parameter turns on thermal and chemical feedback from
+    planetary nebulae.  The mass loss and luminosity are taken from
+    the same `fits from K. Nagamine
+    <http://www.physics.unlv.edu/~kn/SNIa_2/>`_.  The chemical
+    feedback injects gas with the same metallicity as the star
+    particle, and the thermal feedback equates to a 10 km/s wind.  The
+    ejecta are not stored in its own species field.  Can be used
+    with ``StarParticleCreation`` = 0, 1, 2, 5, 7, and 8.  Default: 0.
+
 Normal Star Formation
 ^^^^^^^^^^^^^^^^^^^^^
 
     returned to the gas phase as thermal energy. Default: 1e-5
 ``StarEnergyToStellarUV`` (external)
     The fraction of the rest-mass energy of the stars created which is
-    returned as UV radiation with a young star spectrum. This is used when calculating the radiation background. Default: 3e-6
+    returned as UV radiation with a young star spectrum. This is used
+    when calculating the radiation background. Default: 3e-6
 ``StarEnergyToQuasarUV`` (external)
-    The fraction of the rest-mass energy of the stars created which is returned as UV radiation with a quasar spectrum. This is used when calculating the radiation background. Default: 5e-6
+    The fraction of the rest-mass energy of the stars created which is
+    returned as UV radiation with a quasar spectrum. This is used when
+    calculating the radiation background. Default: 5e-6
 
 Population III Star Formation
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
   
      1. Haardt & Madau spectrum with q_alpha=1.5
      2. Haardt & Madau spectrum with q_alpha = 1.8
-     3. reserved for experimentation
+     3. Modified Haardt & Madau spectrum to match observations
+     	(Kirkman & Tytler 2005).
      4. H&M spectrum (q_alpha=1.5. supplemented with an X-ray Compton heating
-         background from Madau & Efstathiou (see astro-ph/9902080)
+        background from Madau & Efstathiou (see astro-ph/9902080)
      9. a constant molecular H2 photo-dissociation rate
      10. internally computed radiation field using the algorithm of Cen & Ostriker
      11. same as previous, but with very, very simple optical shielding fudge

run/Cooling/OneZoneFreefallTest/OneZoneFreefallTest.enzo

 #
 TopGridRank = 2
 
-# Dim 0 - H number density
-# Dim 1 - metallicity
-# Dim 2 - temperature
-TopGridDimensions = 4 7
+TopGridDimensions = 3 5
 
 OneZoneFreefallTestInitialDensity = 1.0
 OneZoneFreefallTestMinimumEnergy = 1.0
-OneZoneFreefallTestMaximumEnergy = 1000.0
-OneZoneFreefallTestMinimumMetallicity = 1e-6
-OneZoneFreefallTestMaximumMetallicity = 1
+OneZoneFreefallTestMaximumEnergy = 10.0
+OneZoneFreefallTestMinimumMetallicity = 1.0e-7
+OneZoneFreefallTestMaximumMetallicity = 1.0e-3
 
 # Set timestep as a fraction of free-fall time
-OneZoneFreefallTimestepFraction = 1e-3
+OneZoneFreefallTimestepFraction = 1e-2
 
 #
 #  set I/O and stop/start parameters
 #
 StopTime                  = 50
-StopCycle                 = 20000
-CycleSkipDataDump         = 200
+StopCycle                 = 3000
+CycleSkipDataDump         = 20
 DataDumpDir               = DD
 DataDumpName              = DD
 
 #
 OutputCoolingTime         = 1
 OutputTemperature         = 1
+OutputDustTemperature     = 1
 
 #
 # Units
 #
 RadiativeCooling          = 1
 MultiSpecies              = 2
+H2FormationOnDust         = 1
 MetalCooling              = 3             // cloudy cooling
 CloudyCoolingGridFile     = solar_2008_3D_metals.h5 // 3D metals only
-CloudyMetallicityNormalization = 0.018477 // Solar pattern, all metals
-ConstantTemperatureFloor  = 0
 CMBTemperatureFloor       = 0
 IncludeCloudyHeating      = 0
-
+SolarMetalFractionByMass  = 0.02041
 TestProblemUseMetallicityField = 1
-TestProblemMetallicityNormalization = 0.01424544 // consistent with Cloudy 
-				      		 // solar abundances
 
 # Initial species fractions, fiddle at own risk,
-#TestProblemInitialHIFraction      = 0.998
+TestProblemInitialHIFraction      = 0.999
 #TestProblemInitialHIIFraction     = 1e-10
 #TestProblemInitialHeIFraction     = 1.0
 #TestProblemInitialHeIIFraction    = 1.0e-20
 #TestProblemInitialHeIIIIFraction  = 1.0e-20
 #TestProblemInitialHMFraction      = 1.e-20
-#TestProblemInitialH2IFraction     = 1.e-3
+TestProblemInitialH2IFraction     = 1.e-5
 #TestProblemInitialH2IIFraction    = 1.e-20

run/Cooling/OneZoneFreefallTest/notes.txt

 the free-fall time.  Since the timestep continually decreases, outputs are 
 done on cycles.  This test problem can be used to test chemistry and 
 cooling routines.
+
+The script, plot.py, will create plots of n vs. T (and Tdust), n
+vs. f_H2, and n vs. t_cool/t_dyn.  If using H2 formation on 
+dust grains, set dust=True on line 10.  Run this script like this:
+
+python plot.py OneZoneFreefallTest.enzo

run/Cooling/OneZoneFreefallTest/plot.py

+from matplotlib import pylab
+import sys
+
+from yt.mods import *
+from yt.analysis_modules.simulation_handler.api import EnzoSimulation
+
+do_fH2 = True
+do_t_cool = True
+
+dust = True
+if dust:
+    keyword = 'with_dust'
+else:
+    keyword = 'without_dust'
+
+par_file = sys.argv[1]
+es = EnzoSimulation(par_file, get_data_by_force=True, initial_time=0)
+
+T = []
+n = []
+Z = []
+fH2 = []
+Tdust = []
+t_cool = []
+t_dyn = []
+
+for output in es.allOutputs:
+    print output['filename']
+    pf = load(output['filename'])
+    T.append(pf.h.grids[0]['Temperature'])
+    n.append(pf.h.grids[0]['NumberDensity'][0,0,0])
+    Z.append(pf.h.grids[0]['Metallicity'])
+    if do_fH2: fH2.append(pf.h.grids[0]['H2I_Fraction'])
+    if do_t_cool:
+        t_cool.append(pf.h.grids[0]['Cooling_Time'])
+        t_dyn.append(pf.h.grids[0]['DynamicalTime'][0,0,0])
+    if dust: Tdust.append(pf.h.grids[0]['Dust_Temperature'])
+    del pf
+
+T = na.array(T)
+n = na.array(n)
+Z = na.array(Z)
+fH2 = na.array(fH2)
+Tdust = na.array(Tdust)
+t_cool = na.array(t_cool)
+t_dyn = na.array(t_dyn)
+
+colors = ['black', 'purple', 'blue', 'green', 'orange', 'red']
+
+met = na.round(na.log10(Z[0,0,:,0]))
+for i in range(T.shape[2]):
+    pylab.loglog(n, T[:, 0, i, 0], label='log (Z/Z$_{\odot}$) = %d' % met[i],
+                 color=colors[i], linestyle='-')
+    if dust:
+        pylab.loglog(n, Tdust[:, 0, i, 0], color=colors[i], linestyle='--')
+pylab.xlim(xmin=1.0)
+pylab.ylim(1e0, 1e4)
+pylab.xlabel('n [cm$^{-3}$]')
+pylab.ylabel('T [K]')
+pylab.legend(labelspacing=0.0, loc='upper left')
+pylab.savefig('n-T_%s.png' % keyword)
+pylab.clf()
+
+if do_fH2:
+    for i in range(T.shape[2]):
+        pylab.loglog(n, fH2[:, 0, i, 0], label='log (Z/Z$_{\odot}$) = %d' % met[i],
+                     color=colors[i])
+    pylab.xlim(xmin=1.0)
+    pylab.xlabel('n [cm$^{-3}$]')
+    pylab.ylabel('f$_{H_{2}}$')
+    pylab.ylim(1e-5, 1)
+    pylab.legend(labelspacing=0.0, loc='lower right')
+    pylab.savefig('n-fH2_%s.png' % keyword)
+    pylab.clf()
+
+if do_t_cool:
+    for i in range(T.shape[2]):
+        pylab.loglog(n, (t_cool[:, 0, i, 0]/t_dyn), label='log (Z/Z$_{\odot}$) = %d' % met[i],
+                     color=colors[i])
+    pylab.xlim(xmin=1.0)
+    pylab.xlabel('n [cm$^{-3}$]')
+    pylab.ylabel('t$_{cool}$ / t$_{dyn}$')
+    pylab.legend(labelspacing=0.0, loc='lower right')
+    pylab.savefig('n-t_cool_t_dyn_%s.png' % keyword)
+    pylab.clf()

run/Cosmology/Hydro/AMRZeldovichPancake/test_amrzeldovich.py

 from yt.mods import *
 from yt.utilities.answer_testing.api import YTStaticOutputTest, run_main
 
-class TestMARZeldovich(YTStaticOutputTest):
+class TestAMRZeldovich(YTStaticOutputTest):
     name = "amr_zeldovich_plot"
 
     def run(self):

run/Hydro/Hydro-3D/RotatingCylinder/RotatingCylinder.enzo

 RotatingCylinderSubgridLeft = 0.375
 RotatingCylinderSubgridRight = 0.625
 RotatingCylinderTotalEnergy = 0.5
-InitialRefinementDepth Do This Bitch.
+InitialRefinementDepth = 0
 TestProblemUseMetallicityField = 0
 TestProblemInitialMetallicityFraction = 0.1
 
-Gravity Iteration Bitch.
-
 #####

src/P-GroupFinder/Makefile

 #=======================================================================
 
 # Override previous defines with the exception of particle id size
-DEFINES = -DUSE_HDF5 -DH5_USE_16_API $(ASSEMBLE_IDS_DEFINES)
+DEFINES = -DUSE_HDF5 -DFOF_ONLY -DH5_USE_16_API $(ASSEMBLE_IDS_DEFINES)
 
 #-----------------------------------------------------------------------
 # Implicit rules

src/enzo/AMRH5writer.C

     {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"};
   const char *ParticleOtherLabel[] =
     {"particle_type", "particle_index", "particle_mass"};
-  /*  const char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time",
-      "metallicity_fraction", "alpha_fraction", "p5", "p6"}; */
-  char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time",
-				    "metallicity_fraction", 
-				    "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"};
+#ifdef WINDS
+  const char *ParticleAttributeLabel[] =
+    {"creation_time", "dynamical_time", "metallicity_fraction", "particle_jet_x", 
+     "particle_jet_y", "particle_jet_z", "typeia_fraction"};
+#else
+  const char *ParticleAttributeLabel[] = 
+    {"creation_time", "dynamical_time", "metallicity_fraction", "typeia_fraction"};
+#endif
 
   int i;
     
      {"particle_position_x", "particle_position_y", "particle_position_z"};
   const char *ParticleVelocityLabel[] = 
      {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"};
-  char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time",
-				    "metallicity_fraction", 
-				    "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"};
-
-  /*  const char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time",
-      "metallicity_fraction", "alpha_fraction", "p5", "p6"}; */
+#ifdef WINDS
+  const char *ParticleAttributeLabel[] =
+    {"creation_time", "dynamical_time", "metallicity_fraction", "particle_jet_x", 
+     "particle_jet_y", "particle_jet_z", "typeia_fraction"};
+#else
+  const char *ParticleAttributeLabel[] = 
+    {"creation_time", "dynamical_time", "metallicity_fraction", "typeia_fraction"};
+#endif
 
   sprintf(gridDataName, "/grid-%d", gridId);
   if (nBaryonFields > 0) 
      {"particle_position_x", "particle_position_y", "particle_position_z"};
   const char *ParticleVelocityLabel[] = 
      {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"};
-  char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time",
-				    "metallicity_fraction", 
-				    "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"};
-
-  /*  const char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time",
-      "metallicity_fraction", "alpha_fraction", "p5", "p6"}; */
+#ifdef WINDS
+  const char *ParticleAttributeLabel[] =
+    {"creation_time", "dynamical_time", "metallicity_fraction", "particle_jet_x", 
+     "particle_jet_y", "particle_jet_z", "typeia_fraction"};
+#else
+  const char *ParticleAttributeLabel[] = 
+    {"creation_time", "dynamical_time", "metallicity_fraction", "typeia_fraction"};
+#endif
 
   /* if there's no particle, don't bother,
      but if this is a root-level grid, print them anyway --> for Ralf's visualization purpose! */
     {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"};
   const char *ParticleOtherLabel[] =
     {"particle_type", "particle_index", "particle_mass"};
-  /*  const char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time",
-      "metallicity_fraction", "alpha_fraction", "p5", "p6"}; */
-  char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time",
-				    "metallicity_fraction", 
-				    "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"};
+#ifdef WINDS
+  const char *ParticleAttributeLabel[] =
+    {"creation_time", "dynamical_time", "metallicity_fraction", "particle_jet_x", 
+     "particle_jet_y", "particle_jet_z", "typeia_fraction"};
+#else
+  const char *ParticleAttributeLabel[] = 
+    {"creation_time", "dynamical_time", "metallicity_fraction", "typeia_fraction"};
+#endif
 
   int i;
     
      {"particle_position_x", "particle_position_y", "particle_position_z"};
   const char *ParticleVelocityLabel[] = 
      {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"};
-  char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time",
-				    "metallicity_fraction", 
-				    "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"};
-  /*  const char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time",
-      "metallicity_fraction", "alpha_fraction", "p5", "p6"}; */
+#ifdef WINDS
+  const char *ParticleAttributeLabel[] =
+    {"creation_time", "dynamical_time", "metallicity_fraction", "particle_jet_x", 
+     "particle_jet_y", "particle_jet_z", "typeia_fraction"};
+#else
+  const char *ParticleAttributeLabel[] = 
+    {"creation_time", "dynamical_time", "metallicity_fraction", "typeia_fraction"};
+#endif
 
   if (nPart == 0) 
     return 0;
Add a comment to this file

src/enzo/CallProblemSpecificRoutines.C

File contents unchanged.

src/enzo/CallPython.C

 int ExposeDataHierarchy(TopGridData *MetaData, HierarchyEntry *Grid, 
 		       int &GridID, FLOAT WriteTime, int reset, int ParentID, int level);
 void ExposeGridHierarchy(int NumberOfGrids);
-void ExportParameterFile(TopGridData *MetaData, FLOAT CurrentTime);
+void ExportParameterFile(TopGridData *MetaData, FLOAT CurrentTime, FLOAT OldTime, float dtFixed);
 void CommunicationBarrier();
 
 int CallPython(LevelHierarchyEntry *LevelArray[], TopGridData *MetaData,
 	  (NumberOfPythonSubcycleCalls % PythonSubcycleSkip) != 0) return SUCCESS;
     }
 
-    FLOAT CurrentTime;
+    FLOAT CurrentTime, OldTime;
+    float dtFixed;
     int num_grids, start_index;
     num_grids = 0; start_index = 1;
 
     while (Temp2->NextGridThisLevel != NULL)
         Temp2 = Temp2->NextGridThisLevel; /* ugh: find last in linked list */
     CurrentTime = LevelArray[level]->GridData->ReturnTime();
+    OldTime = LevelArray[level]->GridData->ReturnOldTime();
+    dtFixed = LevelArray[level]->GridData->ReturnTimeStep();
     if (ExposeDataHierarchy(MetaData, Temp2->GridHierarchyEntry, start_index,
                 CurrentTime, 1, 0, 0) == FAIL) {
         fprintf(stderr, "Error in ExposeDataHierarchy\n");
         return FAIL;
     }
 
-  ExportParameterFile(MetaData, CurrentTime);
+    ExportParameterFile(MetaData, CurrentTime, OldTime, dtFixed);
 
-  CommunicationBarrier();
-  PyRun_SimpleString("user_script.main()\n");
+    CommunicationBarrier();
+    PyRun_SimpleString("user_script.main()\n");
 
-  PyDict_Clear(grid_dictionary);
-  PyDict_Clear(old_grid_dictionary);
-  PyDict_Clear(hierarchy_information);
-  PyDict_Clear(yt_parameter_file);
-  PyDict_Clear(conversion_factors);
-  return SUCCESS;
+    PyDict_Clear(grid_dictionary);
+    PyDict_Clear(old_grid_dictionary);
+    PyDict_Clear(hierarchy_information);
+    PyDict_Clear(yt_parameter_file);
+    PyDict_Clear(conversion_factors);
+    return SUCCESS;
 #endif
 }
 
Add a comment to this file

src/enzo/CheckForOutput.C

File contents unchanged.

src/enzo/CloudyCoolingData.h

   // Flag to control whether or not to include heating from Cloudy.
   int IncludeCloudyHeating;
 
-  // Flag to control whether or not to include mean molecular weight from Cloudy.
-  int IncludeCloudyMMW;
-
-  // To convert from mass fraction to metallicity.
-  /*
-    x = SUM { A_i * m_i}, for i = 3 to N.
-    A_i = solar number abundance with respect H.
-    m_i = atomic weight.
-    N = Atomic number of heaviest element in cooling model.
-     For solar abundance patters and N = 30 (Zn), x = 0.018477.
-   */
-  float CloudyMetallicityNormalization;
-
   // Factor to account for extra electrons from metals.
   /* 
      f = SUM { A_i * i }, for i = 3 to N.
   // Cooling values
   float *CloudyCooling;
 
-  // Array holding mean molecular weight values
-  float *CloudyMeanMolecularWeight;
-
   // Length of 1D flattened Cloudy data
   int CloudyDataSize;
 };
Add a comment to this file

src/enzo/CommunicationLoadBalancePhotonGrids.C

File contents unchanged.

Add a comment to this file

src/enzo/CommunicationTransferSubgridParticles.C

File contents unchanged.

Add a comment to this file

src/enzo/CommunicationUpdateStarParticleCount.C

File contents unchanged.

src/enzo/CoolData.h

   float RadiationRedshiftDropOff;
   float HydrogenFractionByMass;
   float DeuteriumToHydrogenRatio;
+  float SolarMetalFractionByMass;
 
   /* Equilibrium rates */
 
   float ElectronFracEnd;
   float *metals;
 
+  /* Energy transfer to grains. */
+  float *gas_grain;
+
   /* For analysis, ratios of metal fine structure line emission is
      desired sometimes. */
 

src/enzo/CosmologySimulationInitialize.C

 static float CosmologySimulationInitialFractionH2I   = 2.0e-20;
 static float CosmologySimulationInitialFractionH2II  = 3.0e-14;
 static float CosmologySimulationInitialFractionMetal = 1.0e-10;
+static float CosmologySimulationInitialFractionMetalIa = 1.0e-12;
 static int   CosmologySimulationUseMetallicityField  = FALSE;
  
 static int CosmologySimulationManuallySetParticleMassRatio = FALSE;
   char *DIIName   = "DII_Density";
   char *HDIName   = "HDI_Density";
   char *MetalName = "Metal_Density";
+  char *MetalIaName = "MetalSNIa_Density";
   char *GPotName  = "Grav_Potential";
   char *ForbidName  = "ForbiddenRefinement";
   char *MachName   = "Mach";
 		  &CosmologySimulationInitialFractionH2II);
     ret += sscanf(line, "CosmologySimulationInitialFractionMetal = %"FSYM,
 		  &CosmologySimulationInitialFractionMetal);
+    ret += sscanf(line, "CosmologySimulationInitialFractionMetalIa = %"FSYM,
+		  &CosmologySimulationInitialFractionMetalIa);
     ret += sscanf(line, "CosmologySimulationUseMetallicityField = %"ISYM,
 		  &CosmologySimulationUseMetallicityField);
  
 			     CosmologySimulationInitialFractionH2I,
 			     CosmologySimulationInitialFractionH2II,
 			     CosmologySimulationInitialFractionMetal,
+			     CosmologySimulationInitialFractionMetalIa,
 #ifdef TRANSFER
 			     RadHydroInitialRadiationEnergy,
 #endif
   }
   if (CosmologySimulationUseMetallicityField) {
     DataLabel[i++] = MetalName;
+    if (StarMakerTypeIaSNe)
+      DataLabel[i++] = MetalIaName;
     if(MultiMetals){
       DataLabel[i++] = ExtraNames[0];
       DataLabel[i++] = ExtraNames[1];
 	    CosmologySimulationInitialFractionH2II);
     fprintf(Outfptr, "CosmologySimulationInitialFractionMetal = %"GSYM"\n",
 	    CosmologySimulationInitialFractionMetal);
+    fprintf(Outfptr, "CosmologySimulationInitialFractionMetalIa = %"GSYM"\n",
+	    CosmologySimulationInitialFractionMetalIa);
     fprintf(Outfptr, "CosmologySimulationUseMetallicityField  = %"ISYM"\n\n",
 	    CosmologySimulationUseMetallicityField);
 
 			     CosmologySimulationInitialFractionH2I,
 			     CosmologySimulationInitialFractionH2II,
 			     CosmologySimulationInitialFractionMetal,
+			     CosmologySimulationInitialFractionMetalIa,
 #ifdef TRANSFER
 			     RadHydroInitialRadiationEnergy,
 #endif

src/enzo/DetermineSubgridSizeExtrema.C

 #include "LevelHierarchy.h"
 #include "CommunicationUtilities.h"
 
-#define MINIMUM_EDGE 4
-#define MINIMUM_SIZE 2000
+#define MINIMUM_EDGE 2
+#define MINIMUM_SIZE 64
 
 int DetermineSubgridSizeExtrema(long_int NumberOfCells, int level, int MaximumStaticSubgridLevel)
 {
   /* Now determine subgrid size parameters */
 
   int grids_per_proc = (level > MaximumStaticSubgridLevel) ?
-    OptimalSubgridsPerProcessor : 12;
+    OptimalSubgridsPerProcessor : 2;
 
   MaximumSubgridSize = NumberOfCells / 
     (NumberOfProcessors * grids_per_proc);

src/enzo/EvolveHierarchy.C

 		   Restart);
 
     /* Call inline analysis. */
+
 #ifdef USE_PYTHON
+    LCAPERF_START("CallPython");
     CallPython(LevelArray, &MetaData, 0, 1);
+    LCAPERF_STOP("CallPython");
 #endif
 
     /* Check for resubmission */

src/enzo/EvolveLevel.C

 			  , ImplicitSolver
 #endif
 			  );
+#ifdef USE_PYTHON
+    LCAPERF_START("CallPython");
     CallPython(LevelArray, MetaData, level, 0);
+    LCAPERF_STOP("CallPython");
+#endif
 
     /* Update SubcycleNumber and the timestep counter for the
        streaming data if this is the bottom of the hierarchy -- Note
Add a comment to this file

src/enzo/EvolvePhotons.C

File contents unchanged.

Add a comment to this file

src/enzo/ExternalBoundary.h

File contents unchanged.

Add a comment to this file

src/enzo/ExternalBoundary_IdentifyPhysicalQuantities.C

File contents unchanged.

Add a comment to this file

src/enzo/FluxFix_Grid_CorrectForRefinedFluxes.C

File contents unchanged.

src/enzo/GalaxySimulationInitialize.C

-/***********************************************************************
-/
-/  INITIALIZE A GALAXY SIMULATION
-/
-/  written by: Greg Bryan
-/  date:       May, 1998
-/  modified1:  Elizabeth Tasker, March 2004
-/
-/  PURPOSE:
-/
-/    Set up a number of spherical objects
-/
-/  RETURNS: SUCCESS or FAIL
-/
-************************************************************************/
-
-// This routine intializes a new simulation based on the parameter file.
-//
-
-#ifdef USE_MPI
-#include "mpi.h"
-#endif /* USE_MPI */
-
-#include <string.h>
-#include <stdio.h>
-#include <math.h>
-
-#include "ErrorExceptions.h"
-#include "macros_and_parameters.h"
-#include "typedefs.h"
-#include "global_data.h"
-#include "Fluxes.h"
-#include "GridList.h"
-#include "ExternalBoundary.h"
-#include "Grid.h"
-#include "Hierarchy.h"
-#include "LevelHierarchy.h"
-#include "TopGridData.h"
-
-void WriteListOfFloats(FILE *fptr, int N, float floats[]);
-void WriteListOfFloats(FILE *fptr, int N, FLOAT floats[]);
-void AddLevel(LevelHierarchyEntry *Array[], HierarchyEntry *Grid, int level);
-int RebuildHierarchy(TopGridData *MetaData,
-		     LevelHierarchyEntry *LevelArray[], int level);
-
-
-int GalaxySimulationInitialize(FILE *fptr, FILE *Outfptr, 
-			  HierarchyEntry &TopGrid, TopGridData &MetaData)
-{
-  char *DensName = "Density";
-  char *TEName   = "TotalEnergy";
-  char *GEName   = "GasEnergy";
-  char *Vel1Name = "x-velocity";
-  char *Vel2Name = "y-velocity";
-  char *Vel3Name = "z-velocity";
-  char *MetalName = "Metal_Density";
-
-  /* declarations */
-
-  char  line[MAX_LINE_LENGTH];
-  int   dim, ret, level, disk, i;
-
-  /* make sure it is 3D */
-  
-  if (MetaData.TopGridRank != 3) {
-    ENZO_VFAIL("Cannot do GalaxySimulation in %"ISYM" dimension(s)\n", MetaData.TopGridRank)
-  }
-
-  /* set default parameters */
-
-  float GalaxySimulationGasMass,
-    GalaxySimulationGalaxyMass,
-    GalaxySimulationDiskTemperature,
-    GalaxySimulationAngularMomentum[MAX_DIMENSION],
-    GalaxySimulationUniformVelocity[MAX_DIMENSION],
-    GalaxySimulationUniformDensity,
-    GalaxySimulationUniformEnergy;
-
-  FLOAT GalaxySimulationDiskRadius,
-    GalaxySimulationDiskPosition[MAX_DIMENSION],
-    GalaxySimulationDiskScaleHeightz,
-    GalaxySimulationDiskScaleHeightR;
-
-  float GalaxySimulationInitialTemperature,
-    GalaxySimulationDarkMatterConcentrationParameter,
-    GalaxySimulationInflowTime,
-    GalaxySimulationInflowDensity;
-
-  int   GalaxySimulationRefineAtStart,
-    GalaxySimulationUseMetallicityField;
-  
-  FLOAT LeftEdge[MAX_DIMENSION], RightEdge[MAX_DIMENSION];
-  float ZeroBField[3] = {0.0, 0.0, 0.0};
-
-  /* Default Values */
-
-  GalaxySimulationRefineAtStart      = TRUE;
-  GalaxySimulationUseMetallicityField  = FALSE;
-  GalaxySimulationInitialTemperature = 1000.0;
-  GalaxySimulationDiskRadius         = 0.2;      // [Mpc]
-  GalaxySimulationDiskTemperature    = 1.e4;     // [K]
-  GalaxySimulationDiskScaleHeightz   = 325e-6;
-  GalaxySimulationDiskScaleHeightR   = 3500e-6;
-  GalaxySimulationDarkMatterConcentrationParameter = 12;
-  GalaxySimulationGasMass            = 4.0e10;
-  GalaxySimulationGalaxyMass         = 1.0e12;
-  GalaxySimulationDiskTemperature    = 1000.0;
-  GalaxySimulationInflowTime         = -1;
-  GalaxySimulationInflowDensity      = 0;
-  for (dim = 0; dim < MAX_DIMENSION; dim++) {
-    GalaxySimulationDiskPosition[dim] = 0.5*(DomainLeftEdge[dim] +
-					     DomainRightEdge[dim]);
-    GalaxySimulationAngularMomentum[dim] = 0;
-    GalaxySimulationUniformVelocity[dim] = 0;
-  }
-  GalaxySimulationUniformDensity = 1.0;
-  GalaxySimulationUniformEnergy = 1.0;
-
-  /* read input from file */
-
-  while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) {
-    
-    ret = 0;
-   
-    ret += sscanf(line, "GalaxySimulationRefineAtStart = %"ISYM,
-		  &GalaxySimulationRefineAtStart);
-    ret += sscanf(line, "GalaxySimulationUseMetallicityField = %"ISYM,
-		  &GalaxySimulationUseMetallicityField);
-    ret += sscanf(line, "GalaxySimulationInitialTemperature = %"FSYM,
-		  &GalaxySimulationInitialTemperature);
-    ret += sscanf(line, "GalaxySimulationUniformVelocity = %"FSYM" %"FSYM" %"FSYM,
-                  &GalaxySimulationUniformVelocity[0], &GalaxySimulationUniformVelocity[1],
-                  &GalaxySimulationUniformVelocity[2]);
-    ret += sscanf(line, "GalaxySimulationDiskRadius = %"PSYM,
-		  &GalaxySimulationDiskRadius);
-    ret += sscanf(line, "GalaxySimulationGalaxyMass = %"FSYM,
-		  &GalaxySimulationGalaxyMass);
-    ret += sscanf(line, "GalaxySimulationGasMass = %"FSYM,
-		  &GalaxySimulationGasMass);
-    ret += sscanf(line, "GalaxySimulationDiskPosition = %"PSYM" %"PSYM" %"PSYM, 
-		  &GalaxySimulationDiskPosition[0],
-		  &GalaxySimulationDiskPosition[1],
-		  &GalaxySimulationDiskPosition[2]);
-    ret += sscanf(line, "GalaxySimulationDiskScaleHeightz = %"PSYM,
-		  &GalaxySimulationDiskScaleHeightz);
-    ret += sscanf(line, "GalaxySimulationDiskScaleHeightR = %"PSYM,
-		  &GalaxySimulationDiskScaleHeightR);
-    ret += sscanf(line, "GalaxySimulationDarkMatterConcentrationParameter = %"FSYM,
-		  &GalaxySimulationDarkMatterConcentrationParameter);
-    ret += sscanf(line, "GalaxySimulationDiskTemperature = %"FSYM,
-		  &GalaxySimulationDiskTemperature);
-    ret += sscanf(line, "GalaxySimulationInflowTime = %"FSYM,
-		  &GalaxySimulationInflowTime);
-    ret += sscanf(line, "GalaxySimulationInflowDensity = %"FSYM,
-		  &GalaxySimulationInflowDensity);
-    ret += sscanf(line, "GalaxySimulationAngularMomentum = %"FSYM" %"FSYM" %"FSYM,
-		  &GalaxySimulationAngularMomentum[0],
-		  &GalaxySimulationAngularMomentum[1],
-		  &GalaxySimulationAngularMomentum[2]);
-    
-    /* if the line is suspicious, issue a warning */
-    
-    if (ret == 0 && strstr(line, "=") && strstr(line, "GalaxySimulation") 
-	&& line[0] != '#')
-      fprintf(stderr, "warning: the following parameter line was not interpreted:\n%s\n", line);
-
-  } // end input from parameter file
-
-  /* set up grid */
-
-  if (TopGrid.GridData->GalaxySimulationInitializeGrid(GalaxySimulationDiskRadius,
-						       GalaxySimulationGalaxyMass, 
-						       GalaxySimulationGasMass,
-						       GalaxySimulationDiskPosition, 
-						       GalaxySimulationDiskScaleHeightz,
-						       GalaxySimulationDiskScaleHeightR, 
-						       GalaxySimulationDarkMatterConcentrationParameter,
-						       GalaxySimulationDiskTemperature, 
-						       GalaxySimulationInitialTemperature,
-						       GalaxySimulationAngularMomentum,
-						       GalaxySimulationUniformVelocity,
-						       GalaxySimulationUseMetallicityField,
-						       GalaxySimulationInflowTime,
-						       GalaxySimulationInflowDensity,0)
-	      == FAIL) {
-      ENZO_FAIL("Error in GalaxySimulationInitialize[Sub]Grid.");
-  }// end subgrid if
-
-  /* Convert minimum initial overdensity for refinement to mass
-     (unless MinimumMass itself was actually set). */
-
-  if (MinimumMassForRefinement[0] == FLOAT_UNDEFINED) {
-    MinimumMassForRefinement[0] = MinimumOverDensityForRefinement[0];
-    for (int dim = 0; dim < MetaData.TopGridRank; dim++)
-      MinimumMassForRefinement[0] *=(DomainRightEdge[dim]-DomainLeftEdge[dim])/
-	float(MetaData.TopGridDims[dim]);
-  }
-
-  /* If requested, refine the grid to the desired level. */
-
-  if (GalaxySimulationRefineAtStart) {
-
-    /* Declare, initialize and fill out the LevelArray. */
-
-    LevelHierarchyEntry *LevelArray[MAX_DEPTH_OF_HIERARCHY];
-    for (level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++)
-      LevelArray[level] = NULL;
-    AddLevel(LevelArray, &TopGrid, 0);
-
-    /* Add levels to the maximum depth or until no new levels are created,
-       and re-initialize the level after it is created. */
-
-    for (level = 0; level < MaximumRefinementLevel; level++) {
-      if (RebuildHierarchy(&MetaData, LevelArray, level) == FAIL) {
-	fprintf(stderr, "Error in RebuildHierarchy.\n");
-	return FAIL;
-      }
-      if (LevelArray[level+1] == NULL)
-	break;
-      LevelHierarchyEntry *Temp = LevelArray[level+1];
-      while (Temp != NULL) {
-
-	if (Temp->GridData->GalaxySimulationInitializeGrid(GalaxySimulationDiskRadius,
-						       GalaxySimulationGalaxyMass, 
-						       GalaxySimulationGasMass,
-						       GalaxySimulationDiskPosition, 
-						       GalaxySimulationDiskScaleHeightz,
-						       GalaxySimulationDiskScaleHeightR, 
-						       GalaxySimulationDarkMatterConcentrationParameter,
-						       GalaxySimulationDiskTemperature, 
-						       GalaxySimulationInitialTemperature,
-						       GalaxySimulationAngularMomentum,
-						       GalaxySimulationUniformVelocity,
-						       GalaxySimulationUseMetallicityField,
-						       GalaxySimulationInflowTime,
-						       GalaxySimulationInflowDensity,0)
-	      == FAIL) {
-	    ENZO_FAIL("Error in GalaxySimulationInitialize[Sub]Grid.");
-	}// end subgrid if
-
-	Temp = Temp->NextGridThisLevel;
-      }
-    } // end: loop over levels
-
-    /* Loop back from the bottom, restoring the consistency among levels. */
-
-    for (level = MaximumRefinementLevel; level > 0; level--) {
-      LevelHierarchyEntry *Temp = LevelArray[level];
-      while (Temp != NULL) {
-	if (Temp->GridData->ProjectSolutionToParentGrid(
-				   *LevelArray[level-1]->GridData) == FAIL) {
-	  fprintf(stderr, "Error in grid->ProjectSolutionToParentGrid.\n");
-	  return FAIL;
-	}
-	Temp = Temp->NextGridThisLevel;
-      }
-    }
-
-  } // end: if (GalaxySimulationRefineAtStart)
-
- /* set up field names and units */
-
- int count = 0;
- DataLabel[count++] = DensName;
- DataLabel[count++] = TEName;
- if (DualEnergyFormalism)
-   DataLabel[count++] = GEName;
- DataLabel[count++] = Vel1Name;
- if(MetaData.TopGridRank > 1)
-   DataLabel[count++] = Vel2Name;
- if(MetaData.TopGridRank > 2)
-   DataLabel[count++] = Vel3Name;
- if (GalaxySimulationUseMetallicityField)
-   DataLabel[count++] = MetalName;
-
- for (i = 0; i < count; i++)
-   DataUnits[i] = NULL;
-
- /* Write parameters to parameter output file */
-
- if (MyProcessorNumber == ROOT_PROCESSOR) {
-
-   fprintf(Outfptr, "GalaxySimulationRefineAtStart      = %"ISYM"\n",
-	   GalaxySimulationRefineAtStart);
-   fprintf(Outfptr, "GalaxySimulationUseMetallicityField          = %"ISYM"\n",
-	   GalaxySimulationUseMetallicityField);
-   fprintf(Outfptr, "GalaxySimulationInitialTemperature = %"GOUTSYM"\n",
-	   GalaxySimulationInitialTemperature);
-   fprintf(Outfptr, "GalaxySimulationUniformVelocity    = %"GOUTSYM" %"GOUTSYM" %"GOUTSYM"\n",
-	   GalaxySimulationUniformVelocity[0], GalaxySimulationUniformVelocity[1],
-	   GalaxySimulationUniformVelocity[2]);
-   fprintf(Outfptr, "GalaxySimulationDiskRadius = %"GOUTSYM"\n",
-	   GalaxySimulationDiskRadius);
-   fprintf(Outfptr, "GalaxySimulationGalaxyMass = %"GOUTSYM"\n",
-	   GalaxySimulationGalaxyMass);
-   fprintf(Outfptr, "GalaxySimulationGasMass = %"GOUTSYM"\n",
-	   GalaxySimulationGasMass);
-   fprintf(Outfptr, "GalaxySimulationDiskScaleHeightz = %"GOUTSYM"\n",
-	   GalaxySimulationDiskScaleHeightz);
-   fprintf(Outfptr, "GalaxySimulationDiskScaleHeightR = %"GOUTSYM"\n",
-	   GalaxySimulationDiskScaleHeightR);
-   fprintf(Outfptr, "GalaxySimulationDarkMatterConcentrationParameter = %"GOUTSYM"\n",
-	   GalaxySimulationDarkMatterConcentrationParameter);
-   fprintf(Outfptr, "GalaxySimulationDiskTemperature = %"GOUTSYM"\n",
-	   GalaxySimulationDiskTemperature);
-   fprintf(Outfptr, "GalaxySimulationInflowTime = %"GOUTSYM"\n",
-	   GalaxySimulationInflowTime);
-   fprintf(Outfptr, "GalaxySimulationInflowDensity = %"GOUTSYM"\n",
-	   GalaxySimulationInflowDensity);
-   fprintf(Outfptr, "GalaxySimulationDiskPosition = ");
-   WriteListOfFloats(Outfptr, MetaData.TopGridRank, GalaxySimulationDiskPosition);
-   fprintf(Outfptr, "GalaxySimulationAngularMomentum = ");
-   WriteListOfFloats(Outfptr, MetaData.TopGridRank, GalaxySimulationAngularMomentum);
- }
-
-#ifdef USE_MPI
-
- // BWO: this forces the synchronization of the various point source gravity
- // parameters between processors.  If this is not done, things go to pieces!
-
- MPI_Barrier(MPI_COMM_WORLD);
- MPI_Datatype DataType = (sizeof(float) == 4) ? MPI_FLOAT : MPI_DOUBLE;
- MPI_Bcast(&PointSourceGravityConstant,1,DataType,ROOT_PROCESSOR, MPI_COMM_WORLD);
- MPI_Bcast(&PointSourceGravityCoreRadius,1,DataType,ROOT_PROCESSOR, MPI_COMM_WORLD);
-
-#endif
-
- return SUCCESS;
-
-}
+/***********************************************************************
+/
+/  INITIALIZE A GALAXY SIMULATION
+/
+/  written by: Greg Bryan
+/  date:       May, 1998
+/  modified1:  Elizabeth Tasker, March 2004
+/
+/  PURPOSE:
+/
+/    Set up a number of spherical objects
+/
+/  RETURNS: SUCCESS or FAIL
+/
+************************************************************************/
+
+// This routine intializes a new simulation based on the parameter file.
+//
+
+#ifdef USE_MPI
+#include "mpi.h"
+#endif /* USE_MPI */
+
+#include <string.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "Grid.h"
+#include "Hierarchy.h"
+#include "LevelHierarchy.h"
+#include "TopGridData.h"
+
+void WriteListOfFloats(FILE *fptr, int N, float floats[]);
+void WriteListOfFloats(FILE *fptr, int N, FLOAT floats[]);
+void AddLevel(LevelHierarchyEntry *Array[], HierarchyEntry *Grid, int level);
+int RebuildHierarchy(TopGridData *MetaData,
+		     LevelHierarchyEntry *LevelArray[], int level);
+
+
+int GalaxySimulationInitialize(FILE *fptr, FILE *Outfptr, 
+			  HierarchyEntry &TopGrid, TopGridData &MetaData)
+{
+  char *DensName = "Density";
+  char *TEName   = "TotalEnergy";
+  char *GEName   = "GasEnergy";
+  char *Vel1Name = "x-velocity";
+  char *Vel2Name = "y-velocity";
+  char *Vel3Name = "z-velocity";
+  char *MetalName = "Metal_Density";
+  char *MetalIaName = "MetalSNIa_Density";
+
+  /* declarations */
+
+  char  line[MAX_LINE_LENGTH];
+  int   dim, ret, level, disk, i;
+
+  /* make sure it is 3D */
+  
+  if (MetaData.TopGridRank != 3) {
+    ENZO_VFAIL("Cannot do GalaxySimulation in %"ISYM" dimension(s)\n", MetaData.TopGridRank)
+  }
+
+  /* set default parameters */
+
+  float GalaxySimulationGasMass,
+    GalaxySimulationGalaxyMass,
+    GalaxySimulationDiskTemperature,
+    GalaxySimulationAngularMomentum[MAX_DIMENSION],
+    GalaxySimulationUniformVelocity[MAX_DIMENSION],
+    GalaxySimulationUniformDensity,
+    GalaxySimulationUniformEnergy;
+
+  FLOAT GalaxySimulationDiskRadius,
+    GalaxySimulationDiskPosition[MAX_DIMENSION],
+    GalaxySimulationDiskScaleHeightz,
+    GalaxySimulationDiskScaleHeightR;
+
+  float GalaxySimulationInitialTemperature,
+    GalaxySimulationDarkMatterConcentrationParameter,
+    GalaxySimulationInflowTime,
+    GalaxySimulationInflowDensity;
+
+  int   GalaxySimulationRefineAtStart,
+    GalaxySimulationUseMetallicityField;
+  
+  FLOAT LeftEdge[MAX_DIMENSION], RightEdge[MAX_DIMENSION];
+  float ZeroBField[3] = {0.0, 0.0, 0.0};
+
+  /* Default Values */
+
+  GalaxySimulationRefineAtStart      = TRUE;
+  GalaxySimulationUseMetallicityField  = FALSE;
+  GalaxySimulationInitialTemperature = 1000.0;
+  GalaxySimulationDiskRadius         = 0.2;      // [Mpc]
+  GalaxySimulationDiskTemperature    = 1.e4;     // [K]
+  GalaxySimulationDiskScaleHeightz   = 325e-6;
+  GalaxySimulationDiskScaleHeightR   = 3500e-6;
+  GalaxySimulationDarkMatterConcentrationParameter = 12;
+  GalaxySimulationGasMass            = 4.0e10;
+  GalaxySimulationGalaxyMass         = 1.0e12;
+  GalaxySimulationDiskTemperature    = 1000.0;
+  GalaxySimulationInflowTime         = -1;
+  GalaxySimulationInflowDensity      = 0;
+  for (dim = 0; dim < MAX_DIMENSION; dim++) {
+    GalaxySimulationDiskPosition[dim] = 0.5*(DomainLeftEdge[dim] +
+					     DomainRightEdge[dim]);
+    GalaxySimulationAngularMomentum[dim] = 0;
+    GalaxySimulationUniformVelocity[dim] = 0;
+  }
+  GalaxySimulationUniformDensity = 1.0;
+  GalaxySimulationUniformEnergy = 1.0;
+
+  /* read input from file */
+
+  while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) {
+    
+    ret = 0;
+   
+    ret += sscanf(line, "GalaxySimulationRefineAtStart = %"ISYM,
+		  &GalaxySimulationRefineAtStart);
+    ret += sscanf(line, "GalaxySimulationUseMetallicityField = %"ISYM,
+		  &GalaxySimulationUseMetallicityField);
+    ret += sscanf(line, "GalaxySimulationInitialTemperature = %"FSYM,
+		  &GalaxySimulationInitialTemperature);
+    ret += sscanf(line, "GalaxySimulationUniformVelocity = %"FSYM" %"FSYM" %"FSYM,
+                  &GalaxySimulationUniformVelocity[0], &GalaxySimulationUniformVelocity[1],
+                  &GalaxySimulationUniformVelocity[2]);
+    ret += sscanf(line, "GalaxySimulationDiskRadius = %"PSYM,
+		  &GalaxySimulationDiskRadius);
+    ret += sscanf(line, "GalaxySimulationGalaxyMass = %"FSYM,
+		  &GalaxySimulationGalaxyMass);
+    ret += sscanf(line, "GalaxySimulationGasMass = %"FSYM,
+		  &GalaxySimulationGasMass);
+    ret += sscanf(line, "GalaxySimulationDiskPosition = %"PSYM" %"PSYM" %"PSYM, 
+		  &GalaxySimulationDiskPosition[0],
+		  &GalaxySimulationDiskPosition[1],
+		  &GalaxySimulationDiskPosition[2]);
+    ret += sscanf(line, "GalaxySimulationDiskScaleHeightz = %"PSYM,
+		  &GalaxySimulationDiskScaleHeightz);
+    ret += sscanf(line, "GalaxySimulationDiskScaleHeightR = %"PSYM,
+		  &GalaxySimulationDiskScaleHeightR);
+    ret += sscanf(line, "GalaxySimulationDarkMatterConcentrationParameter = %"FSYM,
+		  &GalaxySimulationDarkMatterConcentrationParameter);
+    ret += sscanf(line, "GalaxySimulationDiskTemperature = %"FSYM,
+		  &GalaxySimulationDiskTemperature);
+    ret += sscanf(line, "GalaxySimulationInflowTime = %"FSYM,
+		  &GalaxySimulationInflowTime);
+    ret += sscanf(line, "GalaxySimulationInflowDensity = %"FSYM,
+		  &GalaxySimulationInflowDensity);
+    ret += sscanf(line, "GalaxySimulationAngularMomentum = %"FSYM" %"FSYM" %"FSYM,
+		  &GalaxySimulationAngularMomentum[0],
+		  &GalaxySimulationAngularMomentum[1],
+		  &GalaxySimulationAngularMomentum[2]);
+    
+    /* if the line is suspicious, issue a warning */
+    
+    if (ret == 0 && strstr(line, "=") && strstr(line, "GalaxySimulation") 
+	&& line[0] != '#')
+      fprintf(stderr, "warning: the following parameter line was not interpreted:\n%s\n", line);
+
+  } // end input from parameter file
+
+  /* set up grid */
+
+  if (TopGrid.GridData->GalaxySimulationInitializeGrid(GalaxySimulationDiskRadius,
+						       GalaxySimulationGalaxyMass, 
+						       GalaxySimulationGasMass,
+						       GalaxySimulationDiskPosition, 
+						       GalaxySimulationDiskScaleHeightz,
+						       GalaxySimulationDiskScaleHeightR, 
+						       GalaxySimulationDarkMatterConcentrationParameter,
+						       GalaxySimulationDiskTemperature, 
+						       GalaxySimulationInitialTemperature,
+						       GalaxySimulationAngularMomentum,
+						       GalaxySimulationUniformVelocity,
+						       GalaxySimulationUseMetallicityField,
+						       GalaxySimulationInflowTime,
+						       GalaxySimulationInflowDensity,0)
+	      == FAIL) {
+      ENZO_FAIL("Error in GalaxySimulationInitialize[Sub]Grid.");
+  }// end subgrid if
+
+  /* Convert minimum initial overdensity for refinement to mass
+     (unless MinimumMass itself was actually set). */
+
+  if (MinimumMassForRefinement[0] == FLOAT_UNDEFINED) {
+    MinimumMassForRefinement[0] = MinimumOverDensityForRefinement[0];
+    for (int dim = 0; dim < MetaData.TopGridRank; dim++)
+      MinimumMassForRefinement[0] *=(DomainRightEdge[dim]-DomainLeftEdge[dim])/
+	float(MetaData.TopGridDims[dim]);
+  }
+
+  /* If requested, refine the grid to the desired level. */
+
+  if (GalaxySimulationRefineAtStart) {
+
+    /* Declare, initialize and fill out the LevelArray. */
+
+    LevelHierarchyEntry *LevelArray[MAX_DEPTH_OF_HIERARCHY];
+    for (level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++)
+      LevelArray[level] = NULL;
+    AddLevel(LevelArray, &TopGrid, 0);
+
+    /* Add levels to the maximum depth or until no new levels are created,
+       and re-initialize the level after it is created. */
+
+    for (level = 0; level < MaximumRefinementLevel; level++) {
+      if (RebuildHierarchy(&MetaData, LevelArray, level) == FAIL) {
+	fprintf(stderr, "Error in RebuildHierarchy.\n");
+	return FAIL;
+      }
+      if (LevelArray[level+1] == NULL)
+	break;
+      LevelHierarchyEntry *Temp = LevelArray[level+1];
+      while (Temp != NULL) {
+
+	if (Temp->GridData->GalaxySimulationInitializeGrid(GalaxySimulationDiskRadius,
+						       GalaxySimulationGalaxyMass, 
+						       GalaxySimulationGasMass,
+						       GalaxySimulationDiskPosition, 
+						       GalaxySimulationDiskScaleHeightz,
+						       GalaxySimulationDiskScaleHeightR, 
+						       GalaxySimulationDarkMatterConcentrationParameter,
+						       GalaxySimulationDiskTemperature, 
+						       GalaxySimulationInitialTemperature,
+						       GalaxySimulationAngularMomentum,
+						       GalaxySimulationUniformVelocity,
+						       GalaxySimulationUseMetallicityField,
+						       GalaxySimulationInflowTime,
+						       GalaxySimulationInflowDensity,0)
+	      == FAIL) {
+	    ENZO_FAIL("Error in GalaxySimulationInitialize[Sub]Grid.");
+	}// end subgrid if
+
+	Temp = Temp->NextGridThisLevel;
+      }
+    } // end: loop over levels
+
+    /* Loop back from the bottom, restoring the consistency among levels. */
+
+    for (level = MaximumRefinementLevel; level > 0; level--) {
+      LevelHierarchyEntry *Temp = LevelArray[level];
+      while (Temp != NULL) {
+	if (Temp->GridData->ProjectSolutionToParentGrid(
+				   *LevelArray[level-1]->GridData) == FAIL) {
+	  fprintf(stderr, "Error in grid->ProjectSolutionToParentGrid.\n");
+	  return FAIL;
+	}
+	Temp = Temp->NextGridThisLevel;
+      }
+    }
+
+  } // end: if (GalaxySimulationRefineAtStart)
+
+ /* set up field names and units */
+
+ int count = 0;
+ DataLabel[count++] = DensName;
+ DataLabel[count++] = TEName;
+ if (DualEnergyFormalism)
+   DataLabel[count++] = GEName;
+ DataLabel[count++] = Vel1Name;
+ if(MetaData.TopGridRank > 1)
+   DataLabel[count++] = Vel2Name;
+ if(MetaData.TopGridRank > 2)
+   DataLabel[count++] = Vel3Name;
+ if (GalaxySimulationUseMetallicityField)
+   DataLabel[count++] = MetalName;
+ if (StarMakerTypeIaSNe)
+   DataLabel[count++] = MetalIaName;
+
+ for (i = 0; i < count; i++)
+   DataUnits[i] = NULL;
+
+ /* Write parameters to parameter output file */
+
+ if (MyProcessorNumber == ROOT_PROCESSOR) {
+
+   fprintf(Outfptr, "GalaxySimulationRefineAtStart      = %"ISYM"\n",
+	   GalaxySimulationRefineAtStart);
+   fprintf(Outfptr, "GalaxySimulationUseMetallicityField          = %"ISYM"\n",
+	   GalaxySimulationUseMetallicityField);
+   fprintf(Outfptr, "GalaxySimulationInitialTemperature = %"GOUTSYM"\n",
+	   GalaxySimulationInitialTemperature);
+   fprintf(Outfptr, "GalaxySimulationUniformVelocity    = %"GOUTSYM" %"GOUTSYM" %"GOUTSYM"\n",
+	   GalaxySimulationUniformVelocity[0], GalaxySimulationUniformVelocity[1],
+	   GalaxySimulationUniformVelocity[2]);
+   fprintf(Outfptr, "GalaxySimulationDiskRadius = %"GOUTSYM"\n",
+	   GalaxySimulationDiskRadius);
+   fprintf(Outfptr, "GalaxySimulationGalaxyMass = %"GOUTSYM"\n",
+	   GalaxySimulationGalaxyMass);
+   fprintf(Outfptr, "GalaxySimulationGasMass = %"GOUTSYM"\n",
+	   GalaxySimulationGasMass);
+   fprintf(Outfptr, "GalaxySimulationDiskScaleHeightz = %"GOUTSYM"\n",
+	   GalaxySimulationDiskScaleHeightz);
+   fprintf(Outfptr, "GalaxySimulationDiskScaleHeightR = %"GOUTSYM"\n",
+	   GalaxySimulationDiskScaleHeightR);
+   fprintf(Outfptr, "GalaxySimulationDarkMatterConcentrationParameter = %"GOUTSYM"\n",
+	   GalaxySimulationDarkMatterConcentrationParameter);
+   fprintf(Outfptr, "GalaxySimulationDiskTemperature = %"GOUTSYM"\n",
+	   GalaxySimulationDiskTemperature);
+   fprintf(Outfptr, "GalaxySimulationInflowTime = %"GOUTSYM"\n",
+	   GalaxySimulationInflowTime);
+   fprintf(Outfptr, "GalaxySimulationInflowDensity = %"GOUTSYM"\n",
+	   GalaxySimulationInflowDensity);
+   fprintf(Outfptr, "GalaxySimulationDiskPosition = ");
+   WriteListOfFloats(Outfptr, MetaData.TopGridRank, GalaxySimulationDiskPosition);
+   fprintf(Outfptr, "GalaxySimulationAngularMomentum = ");
+   WriteListOfFloats(Outfptr, MetaData.TopGridRank, GalaxySimulationAngularMomentum);
+ }
+
+#ifdef USE_MPI
+
+ // BWO: this forces the synchronization of the various point source gravity
+ // parameters between processors.  If this is not done, things go to pieces!
+
+ MPI_Barrier(MPI_COMM_WORLD);
+ MPI_Datatype DataType = (sizeof(float) == 4) ? MPI_FLOAT : MPI_DOUBLE;
+ MPI_Bcast(&PointSourceGravityConstant,1,DataType,ROOT_PROCESSOR, MPI_COMM_WORLD);
+ MPI_Bcast(&PointSourceGravityCoreRadius,1,DataType,ROOT_PROCESSOR, MPI_COMM_WORLD);
+
+#endif
+
+ return SUCCESS;
+
+}
 /* Return time, timestep */
 
    FLOAT ReturnTime() {return Time;};
+   FLOAT ReturnOldTime() {return OldTime;};
    float ReturnTimeStep() {return dtFixed;};
 
   /* Return, set grid ID */
 
    int GadgetComputeTemperatureDEF(FLOAT time, float *temperature);
 
+/* Baryons: compute the dust temperature. */
+
+   int ComputeDustTemperatureField(float *temperature, float *dust_temperature);
+
 /* Baryons: compute X-ray emissivity in specified band. */
 
    int ComputeXrayEmissivity(float *temperature,
 
   /* Identify colour field */
 
-  int IdentifyColourFields(int &SNColourNum, int &MetalNum, int &MBHColourNum,
+  int IdentifyColourFields(int &SNColourNum, int &MetalNum, 
+			   int &MetalIaNum, int &MBHColourNum,
 			   int &Galaxy1ColourNum, int &Galaxy2ColourNum);
 
   /* Identify Multi-species fields. */
 			  float CosmologySimulationInitialFractionH2I,
 			  float CosmologySimulationInitialFractionH2II,
 			  float CosmologySimulationInitialFractionMetal,
+			  float CosmologySimulationInitialFractionMetalIa,
 #ifdef TRANSFER
 			  float RadHydroInitialRadiationEnergy,
 #endif
 			  float CosmologySimulationInitialFractionH2I,
 			  float CosmologySimulationInitialFractionH2II,
 			  float CosmologySimulationInitialFractionMetal,
+			  float CosmologySimulationInitialFractionMetalIa,
 			  int   CosmologySimulationUseMetallicityField,
 			  PINT &CurrentNumberOfParticles,
 			  int CosmologySimulationManuallySetParticleMassRatio,
Add a comment to this file

src/enzo/Grid_AddExternalPotentialField.C

File contents unchanged.

src/enzo/Grid_AddFeedbackSphere.C

   /* Find Metallicity or SNColour field and set flag. */
 
   int SNColourNum, MetalNum, Metal2Num, MBHColourNum, Galaxy1ColourNum, 
-    Galaxy2ColourNum;
+    Galaxy2ColourNum, MetalIaNum;
   int MetallicityField = FALSE;
 
-  if (this->IdentifyColourFields(SNColourNum, Metal2Num, MBHColourNum, 
-				 Galaxy1ColourNum, Galaxy2ColourNum) == FAIL) {
+  if (this->IdentifyColourFields(SNColourNum, Metal2Num, MetalIaNum, 
+				 MBHColourNum, Galaxy1ColourNum, 
+				 Galaxy2ColourNum) == FAIL)
     ENZO_FAIL("Error in grid->IdentifyColourFields.\n");
-  }
 
   MetalNum = max(Metal2Num, SNColourNum);
   MetallicityField = (MetalNum > 0) ? TRUE : FALSE;
   if (MetalNum > 0 && SNColourNum > 0 && cstar->type == PopIII)
     MetalNum = SNColourNum;
 
+  float BubbleVolume = (4.0 * M_PI / 3.0) * radius * radius * radius;
+
   /***********************************************************************
                                 SUPERNOVAE
   ************************************************************************/
   float maxGE, MetalRadius2, PrimordialDensity, metallicity, fhz, fhez;
   float outerRadius2;
 
+  if (cstar->FeedbackFlag == SUPERNOVA || 
+      cstar->FeedbackFlag == CONT_SUPERNOVA) {
+
   // Correct for exaggerated influence radius for pair-instability supernovae
-  if (cstar->FeedbackFlag == SUPERNOVA)
-    radius /= 1.0;
+    if (cstar->FeedbackFlag == SUPERNOVA)
+      radius /= 1.0;
 
   // Correct if the volume with 27 cells is larger than the energy bubble volume
+#ifdef UNUSED
   float BoxVolume = 27 * CellWidth[0][0] * CellWidth[0][0] * CellWidth[0][0];
   float BubbleVolume = (4.0 * M_PI / 3.0) * radius * radius * radius;
   //printf("BoxVolume = %lg, BubbleVolume = %lg\n", BoxVolume, BubbleVolume);
     EjectaMetalDensity *= BubbleVolume / BoxVolume;
     EjectaThermalEnergy *= BubbleVolume / BoxVolume;
   }
+#endif
 //  if (cstar->level > level) {
 //    printf("Reducing ejecta density and energy by 10%% on "
 //	   "level %"ISYM" to avoid crashing.\n", level);
   MetalRadius2 = radius * radius * MetalRadius * MetalRadius;
   outerRadius2 = 1.2 * 1.2 * radius * radius;
 
-  if (cstar->FeedbackFlag == SUPERNOVA || 
-      cstar->FeedbackFlag == CONT_SUPERNOVA) {
-
     /* Remove mass from the star that will now be added to grids. 
        Also, because EjectaDensity will be added with zero net momentum, 
        increase the particle's velocity accordingly. - Ji-hoon Kim, Sep.2009 */

src/enzo/Grid_AddFields.C

 
       }
       BaryonField[n] = new float[size];
-      value = (TypesToAdd[i] == SNColour || TypesToAdd[i] == Metallicity) ?
-	tiny_number : 0.0;
+      value = (TypesToAdd[i] == SNColour || TypesToAdd[i] == Metallicity ||
+	       TypesToAdd[i] == MetalSNIaDensity) ? tiny_number : 0.0;
       for (j = 0; j < size; j++)
 	BaryonField[n][j] = value;
     } // ENDIF this processor
Add a comment to this file

src/enzo/Grid_AddH2Dissociation.C

File contents unchanged.

src/enzo/Grid_CollapseTestInitializeGrid.C

 
   /* Set various units. */
 
-  const float Zsolar = 0.02041;
   const double Mpc = 3.0856e24, SolarMass = 1.989e33, GravConst = 6.67e-8,
     pi = 3.14159, mh = 1.67e-24, kboltz = 1.381e-16, LightSpeed = 2.9979e10;
   float DensityUnits, LengthUnits, TemperatureUnits, TimeUnits, 
 	  /* If there are metals, set it. */
 
 	  if (SphereUseMetals)
-	    BaryonField[MetalNum][n] = metallicity * Zsolar * BaryonField[0][n];
+	    BaryonField[MetalNum][n] = metallicity * CoolData.SolarMetalFractionByMass * 
+	      BaryonField[0][n];
 
 	  /* If there is a colour field, set it. */
 
Add a comment to this file

src/enzo/Grid_ComovingExpansionTerms.C

File contents unchanged.

src/enzo/Grid_ComputeAccelerationFieldExternal.C

       ParticleAcceleration[2][i] += -g*zpos/r;
     }
     
-
   } // end if (ExternalGravity == 1)
 
   /* -----------------------------------------------------------------

src/enzo/Grid_ComputeAccelerationsFromExternalPotential.C

 ************************************************************************/
 
 #include <stdio.h>
+#include "ErrorExceptions.h"
 #include "macros_and_parameters.h"
 #include "typedefs.h"
 #include "global_data.h"
 #include "GridList.h"
 #include "ExternalBoundary.h"
 #include "Grid.h"
-#include "Hierarchy.h"
-#include "LevelHierarchy.h"
-#include "TopGridData.h"
 
 /* function prototypes */
 
 
 
  
-    
   if (DifferenceType != PARTICLES) {
     
     for (dim = 0; dim < GridRank; dim++)
     delete [] Acceleration[dim];
     Acceleration[dim] = NULL;
   }
-   
 
   return SUCCESS;
 }

src/enzo/Grid_ComputeConductionTimeStep.C

   float dt_est;
   double all_units;
 
+  float SpitzerFraction;
+  if (IsotropicConduction && AnisotropicConduction) {
+    SpitzerFraction = max(IsotropicConductionSpitzerFraction, 
+			  AnisotropicConductionSpitzerFraction);
+  }
+  else if (IsotropicConduction) {
+    SpitzerFraction = IsotropicConductionSpitzerFraction;
+  }
+  else if (AnisotropicConduction) {
+    SpitzerFraction = AnisotropicConductionSpitzerFraction;
+  }
+  else {
+    return SUCCESS;
+  }
+
   int size = 1, grid_index, right_side_index; 
   for (int dim = 0; dim < GridRank; dim++) {
     size *= GridDimension[dim];
   // units, scaled correctly. Note that this does NOT contain a factor
   // of 1/mu, since we don't necessarily know anything about the gas in question.
   all_units = POW(dx,2.0)*POW(LengthUnits,2.0)*DensityUnits*kboltz /
-    ( 6.0e-7 * max(IsotropicConductionSpitzerFraction, AnisotropicConductionSpitzerFraction) * mh * TimeUnits );
+    ( 6.0e-7 * SpitzerFraction * mh * TimeUnits );
   
   dt *= float(all_units);
 

src/enzo/Grid_ComputeCoolingTime.C

 
 extern "C" void FORTRAN_NAME(cool_multi_time)(
 	float *d, float *e, float *ge, float *u, float *v, float *w, float *de,
-	   float *HI, float *HII, float *HeI, float *HeII, float *HeIII,
-           float *cooltime,
+	float *HI, float *HII, float *HeI, float *HeII, float *HeIII,
+	float *cooltime,
 	int *in, int *jn, int *kn, int *nratec, int *iexpand,
-           hydro_method *imethod,
-        int *idual, int *ispecies, int *imetal, int *imcool, int *idim,
+	hydro_method *imethod,
+        int *idual, int *ispecies, int *imetal, int *imcool, int *idust, int *idim,
 	int *is, int *js, int *ks, int *ie, int *je, int *ke, int *ih2co,
 	int *ipiht, int *igammah,
 	float *dt, float *aye, float *temstart, float *temend,
 	float *utem, float *uxyz, float *uaye, float *urho, float *utim,
-	float *eta1, float *eta2, float *gamma,
-	float *ceHIa, float *ceHeIa, float *ceHeIIa, float *ciHIa,
-	   float *ciHeIa,
+	float *eta1, float *eta2, float *gamma, float *z_solar,
+	float *ceHIa, float *ceHeIa, float *ceHeIIa, float *ciHIa, float *ciHeIa,
 	float *ciHeISa, float *ciHeIIa, float *reHIIa, float *reHeII1a,
 	float *reHeII2a, float *reHeIIIa, float *brema, float *compa, float *gammaha,
-	float *comp_xraya, float *comp_temp,
-           float *piHI, float *piHeI, float *piHeII,
-	float *HM, float *H2I, float *H2II, float *DI, float *DII, float *HDI,
-           float *metal,
-	float *hyd01ka, float *h2k01a, float *vibha, float *rotha,
-	   float *rotla,
+	float *comp_xraya, float *comp_temp, float *piHI, float *piHeI, float *piHeII,
+	float *HM, float *H2I, float *H2II, float *DI, float *DII, float *HDI, float *metal,
+	float *hyd01ka, float *h2k01a, float *vibha, float *rotha, float *rotla,
 	float *gpldl, float *gphdl, float *HDltea, float *HDlowa,
 	float *gaHIa, float *gaH2a, float *gaHea, float *gaHpa, float *gaela,
-	float *metala, int *n_xe, float *xe_start, float *xe_end,
+	float *gasgra, float *metala, int *n_xe, float *xe_start, float *xe_end,
 	float *inutot, int *iradfield, int *nfreq, int *imetalregen,
 	int *iradshield, float *avgsighp, float *avgsighep, float *avgsighe2p,
 	int *iradtrans, float *photogamma,
-    int *ih2optical, int *iciecool, float *ciecoa,
- 	int *icmbTfloor, int *iClHeat, int *iClMMW,
- 	float *clMetNorm, float *clEleFra, int *clGridRank, int *clGridDim,
+	int *ih2optical, int *iciecool, float *ciecoa,
+ 	int *icmbTfloor, int *iClHeat,
+ 	float *clEleFra, int *clGridRank, int *clGridDim,
  	float *clPar1, float *clPar2, float *clPar3, float *clPar4, float *clPar5,
- 	int *clDataSize, float *clCooling, float *clHeating, float *clMMW);
- 
+ 	int *clDataSize, float *clCooling, float *clHeating);
+
 extern "C" void FORTRAN_NAME(cool_time)(
 	float *d, float *e, float *ge, float *u, float *v, float *w,
            float *cooltime,
        BaryonField[HeINum], BaryonField[HeIINum], BaryonField[HeIIINum],
        cooling_time,
        GridDimension, GridDimension+1, GridDimension+2,
-          &CoolData.NumberOfTemperatureBins, &ComovingCoordinates,
-          &HydroMethod,
+       &CoolData.NumberOfTemperatureBins, &ComovingCoordinates,
+       &HydroMethod,
        &DualEnergyFormalism, &MultiSpecies, &MetalFieldPresent, &MetalCooling, 
+       &H2FormationOnDust,
        &GridRank, GridStartIndex, GridStartIndex+1, GridStartIndex+2,
-          GridEndIndex, GridEndIndex+1, GridEndIndex+2,
+       GridEndIndex, GridEndIndex+1, GridEndIndex+2,
        &CoolData.ih2co, &CoolData.ipiht, &PhotoelectricHeating,
        &dtFixed, &afloat, &CoolData.TemperatureStart,
-          &CoolData.TemperatureEnd,
+       &CoolData.TemperatureEnd,
        &TemperatureUnits, &LengthUnits, &aUnits, &DensityUnits, &TimeUnits,
        &DualEnergyFormalismEta1, &DualEnergyFormalismEta2, &Gamma,
+       &CoolData.SolarMetalFractionByMass,
        CoolData.ceHI, CoolData.ceHeI, CoolData.ceHeII, CoolData.ciHI,
-          CoolData.ciHeI,
+       CoolData.ciHeI,
        CoolData.ciHeIS, CoolData.ciHeII, CoolData.reHII,
-          CoolData.reHeII1,
+       CoolData.reHeII1,
        CoolData.reHeII2, CoolData.reHeIII, CoolData.brem, &CoolData.comp, &CoolData.gammah,
        &CoolData.comp_xray, &CoolData.temp_xray,
-          &CoolData.piHI, &CoolData.piHeI, &CoolData.piHeII,
+       &CoolData.piHI, &CoolData.piHeI, &CoolData.piHeII,
        BaryonField[HMNum], BaryonField[H2INum], BaryonField[H2IINum],
        BaryonField[DINum], BaryonField[DIINum], BaryonField[HDINum],
-          MetalPointer,
+       MetalPointer,
        CoolData.hyd01k, CoolData.h2k01, CoolData.vibh,
-          CoolData.roth, CoolData.rotl,
+       CoolData.roth, CoolData.rotl,
        CoolData.GP99LowDensityLimit, CoolData.GP99HighDensityLimit,
-          CoolData.HDlte, CoolData.HDlow,
+       CoolData.HDlte, CoolData.HDlow,
        CoolData.GAHI, CoolData.GAH2, CoolData.GAHe, CoolData.GAHp,
-       CoolData.GAel,
-          CoolData.metals, &CoolData.NumberOfElectronFracBins, 
-          &CoolData.ElectronFracStart, &CoolData.ElectronFracEnd,
+       CoolData.GAel, CoolData.gas_grain,
+       CoolData.metals, &CoolData.NumberOfElectronFracBins, 
+       &CoolData.ElectronFracStart, &CoolData.ElectronFracEnd,
        RadiationData.Spectrum[0], &RadiationFieldType,
-          &RadiationData.NumberOfFrequencyBins,
-          &RadiationFieldRecomputeMetalRates,
+       &RadiationData.NumberOfFrequencyBins,
+       &RadiationFieldRecomputeMetalRates,
        &RadiationData.RadiationShield, &HIShieldFactor, &HeIShieldFactor, &HeIIShieldFactor,
        &RadiativeTransfer, BaryonField[gammaNum], 
        &H2OpticalDepthApproximation, &CIECooling, CoolData.cieco,
        &CloudyCoolingData.CMBTemperatureFloor,
-       &CloudyCoolingData.IncludeCloudyHeating, &CloudyCoolingData.IncludeCloudyMMW,
-       &CloudyCoolingData.CloudyMetallicityNormalization,
+       &CloudyCoolingData.IncludeCloudyHeating,
        &CloudyCoolingData.CloudyElectronFractionFactor,
        &CloudyCoolingData.CloudyCoolingGridRank,
        CloudyCoolingData.CloudyCoolingGridDimension,
        CloudyCoolingData.CloudyCoolingGridParameters[3],
        CloudyCoolingData.CloudyCoolingGridParameters[4],
        &CloudyCoolingData.CloudyDataSize,
-       CloudyCoolingData.CloudyCooling, CloudyCoolingData.CloudyHeating,
-       CloudyCoolingData.CloudyMeanMolecularWeight);
+       CloudyCoolingData.CloudyCooling, CloudyCoolingData.CloudyHeating);
   } else if (GadgetEquilibriumCooling==1) {  
     int result = GadgetCoolingTime
       (

src/enzo/Grid_ComputeDustTemperatureField.C

+/***********************************************************************
+/
+/  GRID CLASS (COMPUTE THE TEMPERATURE FIELD)
+/
+/  written by: Greg Bryan
+/  date:       April, 1995
+/  modified1:
+/
+/  PURPOSE:
+/
+/  RETURNS:
+/
+************************************************************************/
+ 
+// Compute the pressure at the requested time.  The pressure here is
+//   just the ideal-gas equation-of-state.
+ 
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "fortran.def"
+#include "Grid.h"
+#include "CosmologyParameters.h"
+ 
+/* function prototypes */
+
+int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt);
+int FindField(int f, int farray[], int n);
+int GetUnits(float *DensityUnits, float *LengthUnits,
+	     float *TemperatureUnits, float *TimeUnits,
+	     float *VelocityUnits, FLOAT Time);
+
+extern "C" void FORTRAN_NAME(calc_tdust_3d)(
+	float *d, float *de, float *HI, float *HII, 
+	float *HeI, float *HeII, float *HeIII,
+	float *HM, float *H2I, float *H2II, 
+	int *in, int *jn, int *kn, 
+	int *nratec, int *iexpand,
+	int *ispecies, int *idim,
+	int *is, int *js, int *ks, 
+	int *ie, int *je, int *ke, 
+	float *aye, float *temstart, float *temend,
+	float *gasgra,
+	float *utem, float *uxyz, float *uaye,
+	float *urho, float *utim,
+	float *gas_temp, float *dust_temp);
+
+int grid::ComputeDustTemperatureField(float *temperature, float *dust_temperature)
+{
+  /* Return if this doesn't concern us. */
+ 
+  if (ProcessorNumber != MyProcessorNumber)
+    return SUCCESS;
+
+  /* Return if not using MultiSpecies chemistry. */
+  if (!MultiSpecies) {
+    ENZO_FAIL("Dust temperature calculation requires MultiSpecies > 0.\n");
+  }
+
+  int DensNum;
+  int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
+      DINum, DIINum, HDINum;
+  
+  /* Compute the size of the fields. */
+ 
+  int i, size = 1;
+  for (int dim = 0; dim < GridRank; dim++)
+    size *= GridDimension[dim];
+ 
+  /* Find Density, if possible. */
+ 
+  if ((DensNum = FindField(Density, FieldType, NumberOfBaryonFields)) < 0)
+    ENZO_FAIL("Cannot find density.");
+
+  FLOAT a = 1.0, dadt;
+  float TemperatureUnits = 1, DensityUnits = 1, L