Commits

John Wise committed 97aee6c

Backing out the changeset (594c397dc107) that reversed the premature
merge of the cen and week-of-code branch. When I tried to merge them
together now, all of the changes in the cen branch before this merge
were lost because they were removed in this changeset.

Comments (0)

Files changed (36)

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;
 		  &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
 	    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/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;
+
+}
 
   /* 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,

src/enzo/Grid_AddExternalPotentialField.C

 
 int GetUnits(float *DensityUnits, float *LengthUnits,
        	      float *TemperatureUnits, float *TimeUnits,
-       	      float *VelocityUnits, double *MassUnits, FLOAT Time);
+       	      float *VelocityUnits, float *MassUnits, FLOAT Time);
 
 int grid::AddExternalPotentialField(float *potential)
 {
   /* Get unit conversions */
 
   float DensityUnits = 1 , LengthUnits = 1, TemperatureUnits, 
-    TimeUnits = 1, VelocityUnits = 1, PotentialUnits = 1;
-  double MassUnits = 1;
+    TimeUnits = 1, VelocityUnits = 1, MassUnits = 1, PotentialUnits = 1;
 
    GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
 			&TimeUnits, &VelocityUnits, &MassUnits, Time);

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;

src/enzo/Grid_ComputeAccelerationsFromExternalPotential.C

 
 
  
-    
   if (DifferenceType != PARTICLES) {
     
     for (dim = 0; dim < GridRank; dim++)
     delete [] Acceleration[dim];
     Acceleration[dim] = NULL;
   }
-   
 
   return SUCCESS;
 }

src/enzo/Grid_CosmologySimulationInitializeGrid.C

 			  float CosmologySimulationInitialFractionH2I,
 			  float CosmologySimulationInitialFractionH2II,
 			  float CosmologySimulationInitialFractionMetal,
+			  float CosmologySimulationInitialFractionMetalIa,
 #ifdef TRANSFER
 			  float RadHydroRadiation,
 #endif
  
   int idim, dim, i, j, vel, OneComponentPerFile, ndim, level;
   int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
-      DINum, DIINum, HDINum, MetalNum;
+    DINum, DIINum, HDINum, MetalNum, MetalIaNum;
 #ifdef TRANSFER
   int EgNum;
 #endif
     }
     if (UseMetallicityField) {
       FieldType[MetalNum = NumberOfBaryonFields++] = Metallicity;
+      if (StarMakerTypeIaSNe)
+	FieldType[MetalIaNum = NumberOfBaryonFields++] = MetalSNIaDensity;
       if(MultiMetals){
 	FieldType[ExtraField[0] = NumberOfBaryonFields++] = ExtraType0;
 	FieldType[ExtraField[1] = NumberOfBaryonFields++] = ExtraType1;
   
   // If using metallicity, set the field
  
-  if (UseMetallicityField && ReadData)
-    for (i = 0; i < size; i++) {
+  if (UseMetallicityField && ReadData) {
+    for (i = 0; i < size; i++)
       BaryonField[MetalNum][i] = CosmologySimulationInitialFractionMetal
-	* BaryonField[0][i];  
-      if(MultiMetals){
+	* BaryonField[0][i];
+
+    if (StarMakerTypeIaSNe)
+      for (i = 0; i < size; i++)
+	BaryonField[MetalIaNum][i] = CosmologySimulationInitialFractionMetalIa
+	  * BaryonField[0][i];
+
+    if (MultiMetals) {
+      for (i = 0; i < size; i++) {
 	BaryonField[ExtraField[0]][i] = CosmologySimulationInitialFractionMetal
 	  * BaryonField[0][i];
 	BaryonField[ExtraField[1]][i] = CosmologySimulationInitialFractionMetal
 	  * BaryonField[0][i];
       }
     }
-    if(STARMAKE_METHOD(COLORED_POP3_STAR) && ReadData){
+
+    if (STARMAKE_METHOD(COLORED_POP3_STAR) && ReadData) {
       for (i = 0; i < size; i++)
         BaryonField[ForbidNum][i] = 0.0;
     }
+  } // ENDIF UseMetallicityField
+  
 
 #ifdef EMISSIVITY
     // If using an emissivity field, initialize to zero

src/enzo/Grid_GalaxySimulationInitializeGrid.C

 {
  /* declarations */
 
- int dim, i, j, k, m, field, disk, size, MetalNum, vel;
+  int dim, i, j, k, m, field, disk, size, MetalNum, MetalIaNum, vel;
  int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
    DINum, DIINum, HDINum, B1Num, B2Num, B3Num, PhiNum;
  float DiskDensity, DiskVelocityMag;
 
   if (UseMetallicityField)
     FieldType[MetalNum = NumberOfBaryonFields++] = Metallicity; /* fake it with metals */
+  if (StarMakerTypeIaSNe)
+    FieldType[MetalIaNum = NumberOfBaryonFields++] = MetalSNIaDensity;
 
  /* Return if this doesn't concern us. */
 
  if (ProcessorNumber != MyProcessorNumber) 
    return SUCCESS;
 
-
  /* Set various units. */
 
  float CriticalDensity = 1, BoxLength = 1, mu = 0.6;
 	if (UseMetallicityField)
 	  for (i = 0; i < size; i++)
 	    BaryonField[MetalNum][i] = 1.0e-10;
+	if (StarMakerTypeIaSNe)
+	  for (i = 0; i < size; i++)
+	    BaryonField[MetalIaNum][i] = 1.0e-10;
 
 	/* Set Velocities. */
 

src/enzo/Grid_IdentifyColourFields.C

  
 int FindField(int f, int farray[], int n);
 
-int grid::IdentifyColourFields(int &SNColourNum, int &MetalNum, int &MBHColourNum,
+int grid::IdentifyColourFields(int &SNColourNum, int &MetalNum, 
+			       int &MetalIaNum, int &MBHColourNum,
 			       int &Galaxy1ColourNum, int &Galaxy2ColourNum)
 {
  
-  SNColourNum = MetalNum = MBHColourNum = Galaxy1ColourNum = Galaxy2ColourNum = 0;
+  SNColourNum = MetalNum = MetalIaNum = MBHColourNum = Galaxy1ColourNum = 
+    Galaxy2ColourNum = 0;
  
   SNColourNum = FindField(SNColour, FieldType, NumberOfBaryonFields);
   MetalNum = FindField(Metallicity, FieldType, NumberOfBaryonFields);
+  MetalIaNum = FindField(MetalSNIaDensity, FieldType, NumberOfBaryonFields);
   MBHColourNum = FindField(MBHColour, FieldType, NumberOfBaryonFields);
   Galaxy1ColourNum = FindField(Galaxy1Colour, FieldType, NumberOfBaryonFields);
   Galaxy2ColourNum = FindField(Galaxy2Colour, FieldType, NumberOfBaryonFields);

src/enzo/Grid_InitializeUniformGrid.C

   int dim, i, size, field, GCM;
 
   int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
-    DINum, DIINum, HDINum, MetalNum, B1Num, B2Num, B3Num, PhiNum;
+    DINum, DIINum, HDINum, MetalNum, MetalIaNum, B1Num, B2Num, B3Num, PhiNum;
 
   int CINum, CIINum, OINum, OIINum, SiINum, SiIINum, SiIIINum, CHINum, CH2INum, 
     CH3IINum, C2INum, COINum, HCOIINum, OHINum, H2OINum, O2INum;
   if (TestProblemData.UseMetallicityField) {
     FieldType[MetalNum = NumberOfBaryonFields++] = Metallicity;
 
+    if (StarMakerTypeIaSNe)
+      FieldType[MetalIaNum = NumberOfBaryonFields++] = MetalSNIaDensity;
+
     if(TestProblemData.MultiMetals){
       FieldType[ExtraField[0] = NumberOfBaryonFields++] = ExtraType0;
       FieldType[ExtraField[1] = NumberOfBaryonFields++] = ExtraType1;
     if(TestProblemData.UseMetallicityField){
       BaryonField[MetalNum][i] = TestProblemData.MetallicityField_Fraction* UniformDensity;
 
+      if (StarMakerTypeIaSNe)
+	BaryonField[MetalIaNum][i] = TestProblemData.MetallicitySNIaField_Fraction*
+	  UniformDensity;
+
       if(TestProblemData.MultiMetals){
       BaryonField[ExtraField[0]][i] = TestProblemData.MultiMetalsField1_Fraction* UniformDensity;
       BaryonField[ExtraField[1]][i] = TestProblemData.MultiMetalsField2_Fraction* UniformDensity;

src/enzo/Grid_NestedCosmologySimulationInitializeGrid.C

                           float CosmologySimulationInitialFractionH2I,
                           float CosmologySimulationInitialFractionH2II,
 			  float CosmologySimulationInitialFractionMetal,
+			  float CosmologySimulationInitialFractionMetalIa,
                           int   UseMetallicityField,
                           PINT &CurrentParticleNumber,
                           int CosmologySimulationManuallySetParticleMassRatio,
  
   int idim, ndim, dim, i, j, vel, OneComponentPerFile, level;
   int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
-    DINum, DIINum, HDINum, MetalNum;
+    DINum, DIINum, HDINum, MetalNum, MetalIaNum;
  
   int iTE = ietot;
   int ExtraField[2];
     }
     if (UseMetallicityField) {
       FieldType[MetalNum = NumberOfBaryonFields++] = Metallicity;
+      if (StarMakerTypeIaSNe)
+	FieldType[MetalIaNum = NumberOfBaryonFields++] = MetalSNIaDensity;
       if(MultiMetals){
 	FieldType[ExtraField[0] = NumberOfBaryonFields++] = ExtraType0;
 	FieldType[ExtraField[1] = NumberOfBaryonFields++] = ExtraType1;
  
       // If using metallicity, set the field
       
-      if (UseMetallicityField && ReadData)
-	for (i = 0; i < size; i++) {
+      if (UseMetallicityField && ReadData) {
+	for (i = 0; i < size; i++)
 	  BaryonField[MetalNum][i] = CosmologySimulationInitialFractionMetal
 	    * BaryonField[0][i];
-	  if(MultiMetals){
+
+	if (StarMakerTypeIaSNe)
+	  for (i = 0; i < size; i++)
+	    BaryonField[MetalIaNum][i] = CosmologySimulationInitialFractionMetalIa
+	      * BaryonField[0][i];
+	
+	if (MultiMetals) {
+	  for (i = 0; i < size; i++) {
 	    BaryonField[ExtraField[0]][i] = CosmologySimulationInitialFractionMetal
 	      * BaryonField[0][i];
 	    BaryonField[ExtraField[1]][i] = CosmologySimulationInitialFractionMetal
 	  }
 	}
 
-      if(STARMAKE_METHOD(COLORED_POP3_STAR) && ReadData){
-	for (i = 0; i < size; i++)
-	  BaryonField[ForbidNum][i] = 0.0;
-      }
+	if (STARMAKE_METHOD(COLORED_POP3_STAR) && ReadData) {
+	  for (i = 0; i < size; i++)
+	    BaryonField[ForbidNum][i] = 0.0;
+	}
+      } // ENDIF UseMetallicityField
  
       // If they were not read in above, set the total & gas energy fields now
  

src/enzo/Grid_ParticleSplitter.C

  
   /* Find metallicity field and set flag. */
  
-  int SNColourNum, MetalNum, MBHColourNum, Galaxy1ColourNum, Galaxy2ColourNum; 
+  int SNColourNum, MetalNum, MBHColourNum, Galaxy1ColourNum, Galaxy2ColourNum,
+    MetalIaNum;
   int MetallicityField = FALSE;
 
-  if (this->IdentifyColourFields(SNColourNum, MetalNum, MBHColourNum, 
-				 Galaxy1ColourNum, Galaxy2ColourNum) == FAIL) {
+  if (this->IdentifyColourFields(SNColourNum, MetalNum, MetalIaNum, MBHColourNum,
+				 Galaxy1ColourNum, Galaxy2ColourNum) == FAIL)
     ENZO_FAIL("Error in grid->IdentifyColourFields.\n");
-  }
 
   MetalNum = max(MetalNum, SNColourNum);
   MetallicityField = (MetalNum > 0) ? TRUE : FALSE;

src/enzo/Grid_SolveHydroEquations.C

 
     /* Add "real" colour fields (metallicity, etc.) as colour variables. */
 
-    int SNColourNum, MetalNum, MBHColourNum, Galaxy1ColourNum, Galaxy2ColourNum; 
+    int SNColourNum, MetalNum, MBHColourNum, Galaxy1ColourNum, Galaxy2ColourNum,
+      MetalIaNum; 
 
-    if (this->IdentifyColourFields(SNColourNum, MetalNum, MBHColourNum, 
-				   Galaxy1ColourNum, Galaxy2ColourNum) == FAIL) {
+    if (this->IdentifyColourFields(SNColourNum, MetalNum, MetalIaNum, MBHColourNum, 
+				   Galaxy1ColourNum, Galaxy2ColourNum) == FAIL)
       ENZO_FAIL("Error in grid->IdentifyColourFields.\n");
-    }
 
     if (MetalNum != -1) {
       colnum[NumberOfColours++] = MetalNum;
       }
     }
 
+    if (MetalIaNum       != -1) colnum[NumberOfColours++] = MetalIaNum;
     if (SNColourNum      != -1) colnum[NumberOfColours++] = SNColourNum;
     if (MBHColourNum     != -1) colnum[NumberOfColours++] = MBHColourNum;
     if (Galaxy1ColourNum != -1) colnum[NumberOfColours++] = Galaxy1ColourNum;

src/enzo/Grid_StarParticleHandler.C

              float *odthresh, float *massff, float *smthrest, int *level,
 		 int *np, 
              FLOAT *xp, FLOAT *yp, FLOAT *zp, float *up, float *vp, float *wp,
-             float *mp, float *tdp, float *tcp, float *metalf);
+		float *mp, float *tdp, float *tcp, float *metalf,
+	     int *imetalSNIa, float *metalSNIa, float *metalfSNIa);
  
 extern "C" void FORTRAN_NAME(star_maker3)(int *nx, int *ny, int *nz,
              float *d, float *dm, float *temp, float *u, float *v, float *w,
              float *odthresh, float *massff, float *smthrest, int *level,
 		 int *np, 
              FLOAT *xp, FLOAT *yp, FLOAT *zp, float *up, float *vp, float *wp,
-             float *mp, float *tdp, float *tcp, float *metalf);
- 
+		 float *mp, float *tdp, float *tcp, float *metalf,
+ 	     int *imetalSNIa, float *metalSNIa, float *metalfSNIa);
+
 extern "C" void FORTRAN_NAME(star_maker4)(int *nx, int *ny, int *nz,
              float *d, float *dm, float *temp, float *u, float *v, float *w,
                 float *cooltime,
 	     float *mp, float *tdp, float *tcp, float *metalf, int *type,
 			float *justburn, int *ctype, float *mbhradius);
 
+extern "C" void FORTRAN_NAME(star_feedback_pn_snia)
+  (int *nx, int *ny, int *nz,
+   float *d, float *dm, float *te, float *ge, float *u, float *v,
+       float *w, float *metal,
+   int *idual, int *imetal, hydro_method *imethod, float *dt,
+       float *r, float *dx, FLOAT *t, float *z,
+   float *d1, float *x1, float *v1, float *t1,
+       float *sn_param, float *m_eject, float *yield,
+   int *distrad, int *diststep, int *distcells,
+       int *nmax, FLOAT *xstart, FLOAT *ystart, FLOAT *zstart, int *ibuff,
+   FLOAT *xp, FLOAT *yp, FLOAT *zp, float *up, float *vp, float *wp,
+       float *mp, float *tdp, float *tcp, float *metalf, int *type,
+   float *justburn, int *iPN, int *imetalSNIa, float *metalSNIa);
+
 extern "C" void FORTRAN_NAME(pop3_maker)
   (int *nx, int *ny, int *nz, 
    float *d, float *dm, float *h2d, float *temp, 
 
   /* Find metallicity field and set flag. */
  
-  int SNColourNum, MetalNum, MBHColourNum, Galaxy1ColourNum, Galaxy2ColourNum; 
+  int SNColourNum, MetalNum, MBHColourNum, Galaxy1ColourNum, Galaxy2ColourNum,
+    MetalIaNum;
 
-  if (this->IdentifyColourFields(SNColourNum, MetalNum, MBHColourNum, 
-				 Galaxy1ColourNum, Galaxy2ColourNum) == FAIL) {
+  if (this->IdentifyColourFields(SNColourNum, MetalNum, MetalIaNum, MBHColourNum, 
+				 Galaxy1ColourNum, Galaxy2ColourNum) == FAIL)
     ENZO_FAIL("Error in grid->IdentifyColourFields.\n");
-  }
 
   /* Set variables to type defines to pass to FORTRAN routines */
 
        tg->ParticleVelocity[0], tg->ParticleVelocity[1],
           tg->ParticleVelocity[2],
        tg->ParticleMass, tg->ParticleAttribute[1], tg->ParticleAttribute[0],
-          tg->ParticleAttribute[2]);
+       tg->ParticleAttribute[2],
+       &StarMakerTypeIaSNe, BaryonField[MetalIaNum], tg->ParticleAttribute[3]);
 
       for (i = NumberOfNewParticlesSoFar; i < NumberOfNewParticles; i++)
           tg->ParticleType[i] = NormalStarType;
        tg->ParticleVelocity[0], tg->ParticleVelocity[1],
           tg->ParticleVelocity[2],
        tg->ParticleMass, tg->ParticleAttribute[1], tg->ParticleAttribute[0],
-          tg->ParticleAttribute[2]);
+       tg->ParticleAttribute[2],
+       &StarMakerTypeIaSNe, BaryonField[MetalIaNum], tg->ParticleAttribute[3]);
 
       for (i = NumberOfNewParticlesSoFar; i < NumberOfNewParticles; i++)
           tg->ParticleType[i] = NormalStarType;
        ParticleVelocity[0], ParticleVelocity[1],
           ParticleVelocity[2],
        ParticleMass, ParticleAttribute[1], ParticleAttribute[0],
-          ParticleAttribute[2], ParticleType, &RadiationData.IntegratedStarFormation);
+       ParticleAttribute[2], ParticleType, &RadiationData.IntegratedStarFormation);
  
   } // end: if NORMAL_STAR
  
        ParticleVelocity[0], ParticleVelocity[1],
           ParticleVelocity[2],
        ParticleMass, ParticleAttribute[1], ParticleAttribute[0],
-          ParticleAttribute[2], ParticleType, &RadiationData.IntegratedStarFormation);
+       ParticleAttribute[2], ParticleType, &RadiationData.IntegratedStarFormation);
  
   } // end: if UNIGRID_STAR
  
 
   } // end: if (StarParticleFeedback == 5)
 
+  if (StarMakerTypeIaSNe == 1 || StarMakerPlanetaryNebulae == 1) {
+
+      FORTRAN_NAME(star_feedback_pn_snia)(
+       GridDimension, GridDimension+1, GridDimension+2,
+          BaryonField[DensNum], dmfield,
+          BaryonField[TENum], BaryonField[GENum], BaryonField[Vel1Num],
+          BaryonField[Vel2Num], BaryonField[Vel3Num], MetalPointer,
+       &DualEnergyFormalism, &MetallicityField, &HydroMethod,
+       &dtFixed, BaryonField[NumberOfBaryonFields], &CellWidthTemp,
+          &Time, &zred,
+       &DensityUnits, &LengthUnits, &VelocityUnits, &TimeUnits,
+          &StarEnergyToThermalFeedback, &StarMassEjectionFraction,
+          &StarMetalYield, &StarFeedbackDistRadius, &StarFeedbackDistCellStep, 
+       &StarFeedbackDistTotalCells,
+       &NumberOfParticles,
+          CellLeftEdge[0], CellLeftEdge[1], CellLeftEdge[2], &GhostZones,
+       ParticlePosition[0], ParticlePosition[1],
+          ParticlePosition[2],
+       ParticleVelocity[0], ParticleVelocity[1],
+          ParticleVelocity[2],
+       ParticleMass, ParticleAttribute[1], ParticleAttribute[0],
+       ParticleAttribute[2], ParticleType, &RadiationData.IntegratedStarFormation,
+       &StarMakerPlanetaryNebulae, &StarMakerTypeIaSNe, BaryonField[MetalIaNum]);
+
+  }
+
   /* Convert the species back from fractional densities to real densities. */
  
   for (field = 0; field < NumberOfBaryonFields; field++) {

src/enzo/Grid_SubtractAccretedMassFromSphere.C

 
   /* Find Metallicity or SNColour field and set flag. */
 
-  int SNColourNum, MetalNum, MBHColourNum, Galaxy1ColourNum, Galaxy2ColourNum; 
+  int SNColourNum, MetalNum, MBHColourNum, Galaxy1ColourNum, Galaxy2ColourNum,
+    MetalIaNum;
   int MetallicityField = FALSE;
 
-  if (this->IdentifyColourFields(SNColourNum, MetalNum, MBHColourNum, 
-				 Galaxy1ColourNum, Galaxy2ColourNum) == FAIL) {
+  if (this->IdentifyColourFields(SNColourNum, MetalNum, MetalIaNum, MBHColourNum, 
+				 Galaxy1ColourNum, Galaxy2ColourNum) == FAIL)
     ENZO_FAIL("Error in grid->IdentifyColourFields.\n");
-  }
 
   MetalNum = max(MetalNum, SNColourNum);
   MetallicityField = (MetalNum > 0) ? TRUE : FALSE;

src/enzo/InitializeNew.C

   // Set the number of particle attributes, if left unset
  
   if (NumberOfParticleAttributes == INT_UNDEFINED)
-    if (StarParticleCreation || StarParticleFeedback)
+    if (StarParticleCreation || StarParticleFeedback) {
       NumberOfParticleAttributes = 3;
-    else
+      if (StarMakerTypeIaSNe) NumberOfParticleAttributes++;
+    } else {
       NumberOfParticleAttributes = 0;
+    }
  
   // Give unset parameters their default values
  

src/enzo/Make.config.objects

         FastSiblingLocatorFinalize.o \
         FastSiblingLocatorInitialize.o \
         FastSiblingLocatorInitializeStaticChainingMesh.o \
-        fft66.o \
+        feedback_formulae.o \
+	fft66.o \
         fft90.o \
         ffte4X.o \
         ffte_st1.o \
         Star_SetFeedbackFlag.o \
 	Star_SphereContained.o \
 	Star_SubtractAccretedMassFromCell.o \
-        star_maker1.o \
+	star_feedback_pn_snia.o \
+	star_maker1.o \
         star_maker2.o \
         star_maker3.o \
         star_maker4.o \

src/enzo/NestedCosmologySimulationInitialize.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;
 		  &CosmologySimulationInitialFractionH2II);
     ret += sscanf(line, "CosmologySimulationInitialFractionMetal = %"FSYM,
 		  &CosmologySimulationInitialFractionMetal);
+    ret += sscanf(line, "CosmologySimulationInitialFractionMetalIa = %"FSYM,
+		  &CosmologySimulationInitialFractionMetalIa);
     ret += sscanf(line, "CosmologySimulationUseMetallicityField = %"ISYM,
 		  &CosmologySimulationUseMetallicityField);
  
 			     CosmologySimulationInitialFractionH2I,
 			     CosmologySimulationInitialFractionH2II,
 			     CosmologySimulationInitialFractionMetal,
+			     CosmologySimulationInitialFractionMetalIa,
 			     CosmologySimulationUseMetallicityField,
 			     MetaData.NumberOfParticles,
 			     CosmologySimulationManuallySetParticleMassRatio,
 	    CosmologySimulationInitialFractionH2II);
     fprintf(Outfptr, "CosmologySimulationInitialFractionMetal = %"GSYM"\n",
 	    CosmologySimulationInitialFractionMetal);
+    fprintf(Outfptr, "CosmologySimulationInitialFractionMetalIa = %"GSYM"\n",
+	    CosmologySimulationInitialFractionMetalIa);
     fprintf(Outfptr, "CosmologySimulationUseMetallicityField  = %"ISYM"\n\n",
 	    CosmologySimulationUseMetallicityField);
 
 	   CosmologySimulationInitialFractionH2I,
 	   CosmologySimulationInitialFractionH2II,
 	   CosmologySimulationInitialFractionMetal,
+	   CosmologySimulationInitialFractionMetalIa,
 	   CosmologySimulationUseMetallicityField,
 	   ParticleTempCount,
 	   CosmologySimulationManuallySetParticleMassRatio,

src/enzo/RadiationFieldCalculateRates.C

 
   float exp_arg = -1.0 * POW(Redshift-2.3, 2);
 
+  /* Above z = 2.3, increase the width of the Gaussian by a factor of
+     3 for HI and HeI. */
+
+  float beta2 = (AdjustUVBackgroundHighRedshift && Redshift > 2.3) ? 9.0 : 1.0;
+
   /* ------------------------------------------------------------------ */
   /* 1) For the Haardt and Madau (1996) quasar spectrum (alpha_q = 1.5) */
 
   if (RadiationFieldType == 1) {
 
-    RateData.k24 = 6.7e-13 * POW(1.0+Redshift, 0.43) * exp(exp_arg/1.95)
+    RateData.k24 = 6.7e-13 * POW(1.0+Redshift, 0.43) * exp(exp_arg/1.95/beta2)
                      * TimeUnits * Ramp;
     RateData.k25 = 6.3e-15 * POW(1.0+Redshift, 0.51) * exp(exp_arg/2.35) 
                      * TimeUnits * Ramp;
-    RateData.k26 = 3.2e-13 * POW(1.0+Redshift, 0.50) * exp(exp_arg/2.00) 
+    RateData.k26 = 3.2e-13 * POW(1.0+Redshift, 0.50) * exp(exp_arg/2.00/beta2) 
                      * TimeUnits * Ramp;
-    CoolData.piHI   = 4.7e-24 * POW(1.0+Redshift, 0.43) * exp(exp_arg/1.95) 
+    CoolData.piHI   = 4.7e-24 * POW(1.0+Redshift, 0.43) * exp(exp_arg/1.95/beta2) 
                      / CoolingUnits * Ramp;
-    CoolData.piHeI  = 8.2e-24 * POW(1.0+Redshift, 0.50) * exp(exp_arg/2.00) 
+    CoolData.piHeI  = 8.2e-24 * POW(1.0+Redshift, 0.50) * exp(exp_arg/2.00/beta2) 
                      / CoolingUnits * Ramp;
     CoolData.piHeII = 1.6e-25 * POW(1.0+Redshift, 0.51) * exp(exp_arg/2.35) 
                      / CoolingUnits * Ramp;
   /* 2) For the Haardt and Madau (1996) quasar spectrum (alpha_q = 1.8) */
 
   if (RadiationFieldType == 2) {
-    RateData.k24 = 5.6e-13 * POW(1.0+Redshift, 0.43) * exp(exp_arg/1.95)
+    RateData.k24 = 5.6e-13 * POW(1.0+Redshift, 0.43) * exp(exp_arg/1.95/beta2)
                  * TimeUnits * Ramp;
     RateData.k25 = 3.2e-15 * POW(1.0+Redshift, 0.30) * exp(exp_arg/2.60)
                  * TimeUnits * Ramp;
-    RateData.k26 = 4.8e-13 * POW(1.0+Redshift, 0.43) * exp(exp_arg/1.95)
+    RateData.k26 = 4.8e-13 * POW(1.0+Redshift, 0.43) * exp(exp_arg/1.95/beta2)
                  * TimeUnits * Ramp;
-    CoolData.piHI   = 3.9e-24 * POW(1.0+Redshift, 0.43) * exp(exp_arg/1.95)
+    CoolData.piHI   = 3.9e-24 * POW(1.0+Redshift, 0.43) * exp(exp_arg/1.95/beta2)
                  / CoolingUnits * Ramp;
-    CoolData.piHeI  = 6.4e-24 * POW(1.0+Redshift, 0.43) * exp(exp_arg/2.10)
+    CoolData.piHeI  = 6.4e-24 * POW(1.0+Redshift, 0.43) * exp(exp_arg/2.10/beta2)
                  / CoolingUnits * Ramp;
     CoolData.piHeII = 8.7e-26 * POW(1.0+Redshift, 0.30) * exp(exp_arg/2.70)
                  / CoolingUnits * Ramp;

src/enzo/ReadParameterFile.C

     ret += sscanf(line, "RadiationFieldRedshift = %"FSYM, &RadiationFieldRedshift);
     ret += sscanf(line, "TabulatedLWBackground = %"ISYM, &TabulatedLWBackground);
     ret += sscanf(line, "AdjustUVBackground = %"ISYM, &AdjustUVBackground);
+    ret += sscanf(line, "AdjustUVBackgroundHighRedshift = %"ISYM, &AdjustUVBackgroundHighRedshift);
     ret += sscanf(line, "SetUVBAmplitude = %"FSYM, &SetUVBAmplitude);
     ret += sscanf(line, "SetHeIIHeatingScale = %"FSYM, &SetHeIIHeatingScale);
     ret += sscanf(line, "RadiationFieldLevelRecompute = %"ISYM, &RadiationFieldLevelRecompute);    

src/enzo/SetDefaultGlobalValues.C

   RadiationFieldLevelRecompute = 0;
   RadiationData.RadiationShield = 0;
   AdjustUVBackground          = 1;
+  AdjustUVBackgroundHighRedshift = 0;
   SetUVBAmplitude             = 1.0;
   SetHeIIHeatingScale         = 1.8;
   PhotoelectricHeating	      = 0;
 
   TestProblemData.UseMetallicityField = 0;
   TestProblemData.MetallicityField_Fraction = tiny_number;
+  TestProblemData.MetallicitySNIaField_Fraction = tiny_number;
 
   TestProblemData.UseMassInjection = 0;
   TestProblemData.InitialHydrogenMass = tiny_number;

src/enzo/StarParticleData.h

 /* Star particle parameters. */
 
 SPEXTERN int StarFeedbackType;
+SPEXTERN int StarMakerTypeIaSNe;
+SPEXTERN int StarMakerPlanetaryNebulae;
 SPEXTERN float StarMakerOverDensityThreshold;
 SPEXTERN float StarMakerSHDensityThreshold;
 SPEXTERN float StarMakerMassEfficiency;

src/enzo/Star_SubtractAccretedMassFromCell.C

 
   /* Find Metallicity or SNColour field and set flag. */
 
-  int SNColourNum, MetalNum, MBHColourNum, Galaxy1ColourNum, Galaxy2ColourNum; 
+  int SNColourNum, MetalNum, MBHColourNum, Galaxy1ColourNum, Galaxy2ColourNum,
+    MetalIaNum;
   int MetallicityField = FALSE;
 
-  if (CurrentGrid->IdentifyColourFields(SNColourNum, MetalNum, MBHColourNum, 
-					Galaxy1ColourNum, Galaxy2ColourNum) == FAIL) {
+  if (CurrentGrid->IdentifyColourFields(SNColourNum, MetalNum, MetalIaNum, MBHColourNum, 
+					Galaxy1ColourNum, Galaxy2ColourNum) == FAIL)
     ENZO_FAIL("Error in grid->IdentifyColourFields.\n");
-  }
 
   MetalNum = max(MetalNum, SNColourNum);
   MetallicityField = (MetalNum > 0) ? TRUE : FALSE;

src/enzo/TestProblemData.h

   /*  metallicity fields */
   int UseMetallicityField;
   float MetallicityField_Fraction;
+  float MetallicitySNIaField_Fraction;
+  float MetallicityNormalization;
 
   float InitialMetalMass;
 

src/enzo/WriteParameterFile.C

   fprintf(fptr, "RadiationFieldType             = %"ISYM"\n", RadiationFieldType);
   fprintf(fptr, "TabulatedLWBackground          = %"ISYM"\n", TabulatedLWBackground);
   fprintf(fptr, "AdjustUVBackground             = %"ISYM"\n", AdjustUVBackground);
+  fprintf(fptr, "AdjustUVBackgroundHighRedshift = %"ISYM"\n", AdjustUVBackgroundHighRedshift);
   fprintf(fptr, "SetUVBAmplitude                = %"GSYM"\n", SetUVBAmplitude);
   fprintf(fptr, "SetHeIIHeatingScale            = %"GSYM"\n", SetHeIIHeatingScale);
   fprintf(fptr, "RadiationFieldLevelRecompute   = %"ISYM"\n", RadiationFieldLevelRecompute);

src/enzo/cool1d_multi.F

 !  Parameters
 
       double precision mh
-      real mu_metal
+      real    ZSOLAR, mu_metal, boost_factor
+! MKRJ 12/21/07 boost gammaha by this factor (JHW reduced by 1.0 because
+! it has been included in v2.0
+      parameter (boost_factor = 3.5) 
       parameter (mh = mass_h)      !DPC
       parameter (mu_metal = 16.d0)    ! approx. mean molecular weight of metals
 
+
 !  Locals
 
       integer i, j1, iter, iradfield
          ndi = min(max(1, ndi), NID)
          nti = (log10(0.5d0*(tgas(i)+tgasold(i)))-TEMMIN)/DELT + 1
          nti = min(max(1, nti), NIT)
-         xi = min(3.d0, metal(i,j,k)/(d(i,j,k)*z_solar))
-         metal_cool = cbovcool(nti,ndi) * xi / dom *
-     &                (HI(i,j,k)+HII(i,j,k))
-         metal_heat = cbovheat(nti,ndi) * xi / dom *
-     &                (HI(i,j,k)+HII(i,j,k))
-         if (tgas(i) .le. 1e4) then
-            metal_cool = metal_cool*(tgas(i)/1.d4)**3
-            metal_heat = 0.d0
-         endif
+         xi = min(3.d0, metal(i,j,k)/(d(i,j,k)*ZSOLAR))
+
+c
+c        MKRJ 12/18/06
+c        Use the previous metal cooling/heating rate if T >= 10^4 K
+         if (tgas(i) .gt. 1.e4) then
+            metal_cool = cbovcool(nti,ndi) * xi * dom_inv *
+     &           (HI(i,j,k)+HII(i,j,k))
+            metal_heat = cbovheat(nti,ndi) * xi * dom_inv *
+     &           (HI(i,j,k)+HII(i,j,k))
+         else
+
+c           MKRJ: For low T, adopt the x=10^-2 cooling curve from
+c           Dalgarno & McCray 1972
+	    if (tgas(i) .ge. 1.e3) then
+               if (tgas(i) .ge. 6.31e3) then
+                  metal_cool = 3.13e-27*(tgas(i)**0.396)
+               else
+                  if (tgas(i) .ge. 3.16e3) then
+                     metal_cool = 2.64e-26*(tgas(i)**0.152)
+                  else
+                     metal_cool = 5.28e-27*(tgas(i)**0.352)
+                  endif
+               endif
+            else
+               if (tgas(i) .ge. 3.98e1) then
+                  if (tgas(i) .ge. 2.00e2) then
+                     metal_cool = 3.06e-27*(tgas(i)**0.431)
+                  else
+                     metal_cool = 1.52e-28*(tgas(i)**0.997)
+                  endif
+               else
+                  if (tgas(i) .ge. 2.51e1) then
+                     metal_cool = 2.39e-29*(tgas(i)**1.50)
+                  else
+                     metal_cool = 1.095e-32*(tgas(i)**3.885)
+                  endif
+               endif
+            endif
+
+c	    MKRJ 1/3/06 -- convert metal_cool (n^2 Lambda) from cgs to sim units 
+            metal_cool = metal_cool/coolunit * xi * 
+     &           ((HI(i,j,k)+HII(i,j,k)) * dom_inv)**2
+
+c	    MKRJ 12/18/06 -- Photoelectric heating rate (Gamma_pe)
+c            8.5e-26 came from Gamma_pe, assuming epsilon=0.05, G_0=1.7
+c            (see Gerritsen & Icke '97, or Joung & Mac Low '06, Sec. 2.2)
+c
+	    metal_heat = float(igammah) * boost_factor * 
+     &           gammaha * xi * (HI(i,j,k)+HII(i,j,k)) * dom_inv
+c
+c           MKRJ 10/15/08 -- Gamma (in erg/s) should weakly depend on 
+c            the gas density. The value 8.5e-26 is really for the solar 
+c            neighborhood, where n \approx 0.5 cm^-3.
+c            So assume:  emissivity \propto rho^1.4 &  
+c                        F_uv = emissivity*l_mfp \propto rho^0.4
+c            (n_H = 0.76*d(i,j,k)*dom)
+c
+            if (0.76*d(i,j,k)*dom .gt. 5.0e-1) then
+	       metal_heat = metal_heat * 
+     &                     (0.76*d(i,j,k)*dom/5.0e-1)**0.4
+            endif
+         endif ! tgas(i) .gt. 1.e4
 
          edot(i) = edot(i) + metal_heat - metal_cool
 

src/enzo/feedback_formulae.F

+c=======================================================================
+c//////////////////////////////////////\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
+c
+      subroutine SNII_feedback(mtot, age, Edot, Mdot)
+c
+c  Fitting formula of Type II SNe feedback
+c 
+c  written by: Shikui Tang
+c  date:  Febuary 9 2009
+c
+c  dE/dt = 1.e52 * (Tot_mass / 1.e11 Msun) * A(t),
+c   where A(t) equals:
+c     (age[i]/3.16e6 yr) * exp(1-age[i]/3.16e6 yr) , for 0 < t < 3.16e6 yr
+c     pow(age[i]/3.16e6 yr, -1.1) , for 3.16e6 yr <= t < 3.16e7 yr
+c     0 , for t >= 3.16e7 yr. The integration of A(t/3.16e6) is
+c  \intA(x) = e*(1-exp(-x)*(x+1)) (x=t/3.16E6<=1)
+c           = e*(1-2/e) + 10(1-x^(-0.1)) (x>1)
+c
+c  dM/dt = B * dE/dt ,
+c  where B is chosen such that 25% of the original star particle mass
+c  gets ejected after 3.16e7 yr. Revised it according to A(t):
+c  dM/dt = B' * Tot_mass * A(x) = B' * A(t/3.16E6)
+c  Given \intA(x<=10)=2.775, B' = 2.85095E-8 /yr/Msun
+c         ??? seems 15% is more reasonable, which gives B'=1.71E-8
+c
+c  Thus dE/dt = 1D41/B' * dM/dt = 3.5075E48 * dM/dt  (25% ejected at x=10)
+c                               = 5.8458E48 * dM/dt  (15% ejected at x=10)
+c  dM_ZII/dt = YII * dM/dt ,
+c  where YII is the metal yield for SN II.
+c  But metal is not handled in this subroutine.
+c
+c  INPUTS:
+c
+c    mtot  - total mass of the star particle, in the unit of solar mass
+c          - Should it be the initial mass of the gas particle ???
+c    age   - the age of the star particle, in the unit of year
+c
+c  OUTPUTS:
+c
+c    Edot - energy input rate of the star particle, ergs/yr
+c    Mdot - mass input rate of the star particle, Msun/yr
+c
+c  REMARKS:
+c    We may use the integral form to make result independant on time step
+c-----------------------------------------------------------------------
+       implicit none
+c-----------------------------------------------------------------------
+c
+c  Arguments:
+c
+       real mtot, age, Edot, Mdot
+c
+c
+c  Locals:
+c
+       real tt, At
+c
+       tt = age/3.16E6       
+c       if ( age .le. 0 ) stop !!Hope this never happens.
+       if (age .lt. 3.16E6) then
+          At = tt * exp(1-tt)
+       else if (age .lt. 3.16E7) then
+          At = tt**(-1.1)
+       else
+          At = 0.
+       endif
+c
+       Edot = 1.D41 * mtot * At        !ergs/yr
+       !!Mdot = 2.851E-8 * mtot * At     !Msun/yr,  25% return at tt=10
+       Mdot = 1.71D-8 * mtot * At      !Msun/yr, 15% return at tt=10
+       return
+       end
+c       
+c
+c=======================================================================
+c//////////////////////////////////////\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
+c
+c      subroutine SNIa_feedback(mtot, age, Edot, Mdot, tau)
+      subroutine SNIa_feedback(age, Mdot, tau) 
+c 
+c  Fitting formula of Type Ia SNe feedback
+c
+c  written by: Shikui Tang
+c  date:  Febuary 7 2009
+c  modified1: March 7 2009
+c      Keep the form of Mdot, Edot=Mdot*eSNIa*c^2
+c
+c  dE/dt = 3.e41 * (Tot_mass/1.e12 Msun) * pow(age[i]/1.5e10 yr,-1.1)
+c         * F(t) * Year_in_sec,
+c  where F(t) = pow(age[i]/tau, 1.5) / (1 + pow(age[i]/tau, 1.5)), and
+c  age[i] is the age in yrs from t=0, i.e., the birth epoch of a stellar
+c  population with a stellar mass "Tot_mass".  'tau' is what's called tdyn
+c  or tdp in the Fortran source code.
+c  dM/dt  = ??
+c   1 solar mass of metals per 10^51 erg ??? Iron
+c  dM_ZIa = YIa * dM/dt
+c
+c  The Arguments have the same meaning as those in SNII_feedback
+c
+c-----------------------------------------------------------------------
+       implicit none
+c-----------------------------------------------------------------------
+c
+c  Arguments:
+c
+c       real mtot, age, Edot, Mdot, tau
+       real age, Mdot, tau
+c
+c
+c  Locals:
+c
+       real tt
+       if (age .le. 3.16E7) then
+         !!Edot = 0.0
+         Mdot = 0.0
+       else
+         tt = (age/tau)**1.5
+         !!Edot = mtot * 9.467E36 * (age/1.5E10)**(-1.1) * tt/(1+tt) !!ergs/yr
+         !!Edot = Edot *2 !!In order to match Ken's plot
+         !!9.467E36 = 3.e41/1.e12 * Year_in_sec
+         !!Mdot = 1.4* Edot/1D51  !!Msun/yr, 1.4Msun per 1D51 ergs
+         !!To match Ken's plot of SNIa mass loss, in unit of Msun/yr
+         Mdot = 2.65076E-14*(age/1.5E10)**(-1.1)*tt/(1+tt)
+       endif
+c
+       return
+       end
+c
+c
+c=======================================================================
+c//////////////////////////////////////\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
+c
+      subroutine PN_feedback(mtot, age, Edot, Mdot)
+c
+c  Fitting formula of PN  feedback
+c
+c  written by: Shikui Tang
+c  date:  Febuary 7 2009
+c  dE/dt = 0
+c  dM/dt = (Tot_mass / 1.e11 Msun) * P(t), where P(t) is given by:
+c    pow(10, -0.9 * (x-9) + 0.6) , for 7.5 < x < 9.2
+c    pow(10, -1.4 * (x-9) + 0.7) , for x > 9.2
+c  where x=log(t [yr]), and t=0 is the formation time of a given star
+c  particle in our simulations.
+c
+c  dM_Z/dt = 0 (metal output negligible)
+c
+c  The Arguments have the same meaning as those in SNII_feedback
+c
+c  REMARKS:
+c    We may use the integral form to make result independant on time step
+c
+c-----------------------------------------------------------------------
+       implicit none
+c-----------------------------------------------------------------------
+c
+c  Arguments:
+c
+       real mtot, age, Edot, Mdot
+c
+c  Locals:
+c
+      real tt
+c      
+      Edot = 0.0
+      tt = log10(age)
+      if ( tt .le. 7.5) then
+         Mdot = 0.0
+      else if (tt .le. 9.2) then
+         Mdot = mtot/1.E11 * 10**(-0.9*(tt-9)+0.6)
+      else
+         Mdot = mtot/1.E11 * 10**(-1.4*(tt-9)+0.7)
+      endif
+c
+      return
+      end
+c
+c=======================================================================
+c//////////////////////////////////////\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
+c      function rmf(age, tau)
+      subroutine cal_rmf(age, rmf, mlossSNII, mlossPN)
+c
+c  The residual mass fraction of a star particle due to mass ejection
+c  by SNe II and PN, using the differential dM/dt listed above,
+c  normalized to a one Msun stellar population
+c
+c  INPUTS:
+c     age - the age of a stellar population
+c  OUTPUT:
+c     rmf - return the residual mass of a one Msun stellar population
+c     mlossSNII - the cumulative mass ejection due to Type II SNe
+c     mlossPN - the cumulative mass ejection due to PN.
+c
+c-----------------------------------------------------------------------
+       implicit none
+c-----------------------------------------------------------------------
+c
+c  Arguments:
+c
+       real rmf, age
+c
+c  Locals:
+c
+       real x, mlossSNII, mlossPN
+       real e
+c
+c      cumulative mass loss due to SNe II, 15% mass loss when x=10 
+c
+       x = age/3.16E6
+       e = exp(1.0)
+       if ( x .le. 1.0) then
+          mlossSNII = 0.054 * e * (1 - (x+1)/exp(x))
+       else if (x .le. 10.0) then
+          mlossSNII = 0.054 * (e*(1-2/e) + 10*(1-x**(-0.1)))
+       else
+          mlossSNII = 0.1498
+       endif
+c
+c      cumulative mass loss due to PN
+c
+       x = alog10(age)
+       if ( x .le. 7.5) then
+          mlossPN = 0
+       else if ( x .le. 9.2) then
+          mlossPN = 10**(-1.3) * (age**0.1 - 10**0.75)
+       else
+          mlossPN = 0.239 - 498.82 * age**(-0.4)
+       endif
+c
+       rmf = 1.0 - mlossSNII - mlossPN
+       return
+       end       

src/enzo/global_data.h

    10 - homogenous internal radiation field (a la Renyue's work) */
 
 EXTERN int RadiationFieldType;
-EXTERN int AdjustUVBackground; 
+EXTERN int AdjustUVBackground;
+EXTERN int AdjustUVBackgroundHighRedshift;
 EXTERN float SetUVBAmplitude;
 EXTERN float SetHeIIHeatingScale;
 EXTERN RadiationFieldDataType RadiationData;

src/enzo/hydro_rk/Grid_ReturnHydroRKPointers.C

 
   /* Add the colours (NColor is determined in EvolveLevel_RK2) */  
 
-  int SNColourNum, MetalNum, MBHColourNum, Galaxy1ColourNum, Galaxy2ColourNum; 
+  int SNColourNum, MetalNum, MetalIaNum, MBHColourNum, Galaxy1ColourNum, 
+    Galaxy2ColourNum; 
 
-  if (this->IdentifyColourFields(SNColourNum, MetalNum, MBHColourNum, 
+  if (this->IdentifyColourFields(SNColourNum, MetalNum, MetalIaNum, MBHColourNum, 
 				 Galaxy1ColourNum, Galaxy2ColourNum) == FAIL) {
     fprintf(stderr, "Error in grid->IdentifyColourFields.\n");
     return FAIL;
   
   if (MetalNum != -1) {
     Prim[nfield++] = BaryonField[MetalNum];
+    if (StarMakerTypeIaSNe)
+      Prim[nfield++] = OldBaryonField[MetalIaNum];
     if (MultiMetals || TestProblemData.MultiMetals) {
       Prim[nfield++] = BaryonField[MetalNum+1];
       Prim[nfield++] = BaryonField[MetalNum+2];

src/enzo/hydro_rk/Grid_ReturnOldHydroRKPointers.C

 
   /* Add the colours (treat them as species) */
 
-  int SNColourNum, MetalNum, MBHColourNum, Galaxy1ColourNum, Galaxy2ColourNum; 
+  int SNColourNum, MetalNum, MetalIaNum, MBHColourNum, Galaxy1ColourNum, 
+    Galaxy2ColourNum; 
 
-  if (this->IdentifyColourFields(SNColourNum, MetalNum, MBHColourNum, 
+  if (this->IdentifyColourFields(SNColourNum, MetalNum, MetalIaNum, MBHColourNum, 
 				 Galaxy1ColourNum, Galaxy2ColourNum) == FAIL) {
     fprintf(stderr, "Error in grid->IdentifyColourFields.\n");
     return FAIL;
   
   if (MetalNum != -1) {
     Prim[nfield++] = OldBaryonField[MetalNum];
+    if (StarMakerTypeIaSNe)
+      Prim[nfield++] = OldBaryonField[MetalIaNum];
     if (MultiMetals || TestProblemData.MultiMetals) {
       Prim[nfield++] = OldBaryonField[MetalNum+1];
       Prim[nfield++] = OldBaryonField[MetalNum+2];

src/enzo/hydro_rk/Grid_SetNumberOfColours.C

 
   /* Count colours */  
 
-  int SNColourNum, MetalNum, MBHColourNum, Galaxy1ColourNum, Galaxy2ColourNum; 
+  int SNColourNum, MetalNum, MetalIaNum, MBHColourNum, Galaxy1ColourNum, 
+    Galaxy2ColourNum; 
 
-  if (this->IdentifyColourFields(SNColourNum, MetalNum, MBHColourNum, 
+  if (this->IdentifyColourFields(SNColourNum, MetalNum, MetalIaNum, MBHColourNum, 
 				 Galaxy1ColourNum, Galaxy2ColourNum) == FAIL) {
     fprintf(stderr, "Error in grid->IdentifyColourFields.\n");
     return FAIL;
 
   if (MetalNum != -1) {
     _nc++;
+    if (StarMakerTypeIaSNe) _nc++;
     if (MultiMetals || TestProblemData.MultiMetals) {
       _nc += 2;
     }

src/enzo/star_feedback_pn_snia.F

+#include "fortran.def"
+#include "phys_const.def"
+c
+c=======================================================================
+c/////////////////////  SUBROUTINE STAR_FEEDBACK \\\\\\\\\\\\\\\\\\\\\\\
+c
+      subroutine star_feedback_pn_snia(nx, ny, nz,
+     &                      d, dm, te, ge, u, v, w, metal,
+     &                      idual, imetal, imethod, dt, r, dx, t, z,
+     &                      d1, x1, v1, t1, sn_param, m_eject, yield,
+     &                      distrad, diststep, distcells,
+     &                      npart, xstart, ystart, zstart, ibuff,
+     &                      xp, yp, zp, up, vp, wp,
+     &                      mp, tdp, tcp, metalf, type, justburn,
+     &                      iPN, imetalSNIa, metalSNIa)
+
+c
+c  RELEASES "STAR" PARTICLE ENERGY, MASS AND METALS
+c
+c  written by: Chris Loken & Greg Bryan
+c  date:       3 March 1997
+c  modified1:  BWO
+c              13 Nov 2002
+c              Many changes have been made between the date of
+c              initial creation and today - all of them unlogged,
+c              unfortunately.  Bugs were fixed in calculation of 
+c              total energy and conservation of metal and gas density.
+c              The code has been cleaned up in general to enhance
+c              readability.  This is the stable version - star_maker1
+c              and star_maker3 should be used for experimentation.
+c  modified2: 20 Feb 2008 by M. Ryan Joung
+c    added metalSNIa & metalfSNIa; included feedback from SN Ia/PN
+c  modified3: 3 Apr 2008 by M. Ryan Joung
+c    added smthresh to use as proxy for initial SP mass
+c  modified4: 10 Feb 2009 by Shikui Tang
+c    use the KENTARO NAGAMINE's stellar feedback fitting formula
+c    http://www.physics.unlv.edu/~kn/SNIa_2/
+c  modified5: 21 Aug 2011 by John Wise
+c    Removed Type II SNe, so this can be used in conjunction with any
+c    other feedback routine