From bd6c84a8fe9c92c15937df054f95f19481f368bd Mon Sep 17 00:00:00 2001
From: tabel <none@none>
Date: Tue, 29 Sep 2009 09:09:29 -0700
Subject: [PATCH 01/12] Since we can read in Mu, the mean molecular weight, we
 use that now in Grid_ComputeTemperatureField to be consistent between all
 solvers. The default remains 0.6 as set in SetDefaulGlobalValues.

--HG--
branch : week-of-code
---
 src/enzo/Grid_ComputeTemperatureField.C | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/enzo/Grid_ComputeTemperatureField.C b/src/enzo/Grid_ComputeTemperatureField.C
index e7de3a4ef..e3658c9e0 100644
--- a/src/enzo/Grid_ComputeTemperatureField.C
+++ b/src/enzo/Grid_ComputeTemperatureField.C
@@ -30,7 +30,7 @@
  
 /* Set the mean molecular mass. */
  
-#define DEFAULT_MU 0.6
+//#define DEFAULT_MU 0.6
 #define MU_METAL 16.0
  
 /* This is minimum returned temperature. (K) */
@@ -122,7 +122,9 @@ int grid::ComputeTemperatureField(float *temperature)
 
   /* For Sedov Explosion compute temperature without floor */
 
-  float mol_weight = DEFAULT_MU, min_temperature = 1.0;
+  float mol_weight = Mu, min_temperature = 1.0;   // Mu is now read from parameter file
+
+  // this : would not be needed as we can specify it in parameter file:
   if (ProblemType == 7) {//AK for Sedov explosion test
     mol_weight = 1.0;
     min_temperature = tiny_number;
@@ -131,7 +133,7 @@ int grid::ComputeTemperatureField(float *temperature)
   if (MultiSpecies == FALSE)
  
     /* If the multi-species flag is not set,
-       Compute temperature T = p/d and assume mu = DEFAULT_MU. */
+       Compute temperature T = p/d and assume mu = Mu (global data). */
  
     for (i = 0; i < size; i++)
       temperature[i] = max((TemperatureUnits*temperature[i]*mol_weight

From 214f02f0e3f714c69fb3e72e44f340e612dd1caa Mon Sep 17 00:00:00 2001
From: tabel <none@none>
Date: Tue, 6 Oct 2009 16:54:15 -0700
Subject: [PATCH 02/12] Added a PutSink option to TurbulenceInitialize problem
 type parameter file.

--HG--
branch : week-of-code
---
 src/enzo/Grid.h                               |  4 +--
 src/enzo/Grid_DepositBaryons.C                |  2 ++
 src/enzo/Grid_Group_WriteGrid.C               |  2 +-
 src/enzo/Grid_PreparePotentialField.C         |  2 +-
 src/enzo/hydro_rk/Collapse3DInitialize.C      | 32 +++++++++----------
 .../hydro_rk/Grid_Collapse3DInitializeGrid.C  |  2 +-
 .../Grid_MHDTurbulenceInitializeGrid.C        |  2 +-
 .../hydro_rk/Grid_TurbulenceInitializeGrid.C  | 32 ++++++++++++-------
 src/enzo/hydro_rk/MHDTurbulenceInitialize.C   |  3 +-
 src/enzo/hydro_rk/TurbulenceInitialize.C      |  6 ++--
 src/enzo/macros_and_parameters.h              |  2 +-
 11 files changed, 51 insertions(+), 38 deletions(-)

diff --git a/src/enzo/Grid.h b/src/enzo/Grid.h
index 997f59529..fbc0fb121 100644
--- a/src/enzo/Grid.h
+++ b/src/enzo/Grid.h
@@ -1948,8 +1948,8 @@ int CollapseTestInitializeGrid(int NumberOfSpheres,
 	      float fluxcoef, int fallback);
   int TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FLOAT CloudRadius, 
 			       float CloudMachNumber, float CloudAngularVelocity, float InitialBField,
-			       int SetTurbulence, int CloudType, int TurbulenceSeed, int level, 
-			       int SetBaryonFields);
+			       int SetTurbulence, int CloudType, int TurbulenceSeed, int PutSink, 
+			       int level, int SetBaryonFields);
   int Collapse3DInitializeGrid(int n_sphere,
 			       FLOAT r_sphere[MAX_SPHERES],
 			       FLOAT rc_sphere[MAX_SPHERES],
diff --git a/src/enzo/Grid_DepositBaryons.C b/src/enzo/Grid_DepositBaryons.C
index 059325a04..3fd07f4b8 100644
--- a/src/enzo/Grid_DepositBaryons.C
+++ b/src/enzo/Grid_DepositBaryons.C
@@ -145,6 +145,8 @@ int grid::DepositBaryons(grid *TargetGrid, FLOAT DepositTime)
  
     if (RegionDim[dim] < 2) {
       fprintf(stderr, "RegionDim[%"ISYM"] = %"ISYM" < 2\n", dim, RegionDim[dim]);
+      fprintf(stderr, "GridOffsetEnd[%"ISYM"] = %"ISYM" < 2\n", dim, GridOffsetEnd[dim]);
+      fprintf(stderr, "GridOffset[%"ISYM"] = %"ISYM" < 2\n", dim, GridOffset[dim]);
       ENZO_FAIL("");
     }
  
diff --git a/src/enzo/Grid_Group_WriteGrid.C b/src/enzo/Grid_Group_WriteGrid.C
index 3c554dde1..147aa9132 100644
--- a/src/enzo/Grid_Group_WriteGrid.C
+++ b/src/enzo/Grid_Group_WriteGrid.C
@@ -82,7 +82,7 @@ int grid::Group_WriteGrid(FILE *fptr, char *base_name, int grid_id, HDF5_hid_t f
   char *ParticleVelocityLabel[] =
      {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"};
   char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time",
-				    "metallicity_fraction", "alpha_fraction"};
+				    "metallicity_fraction", "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"};
   char *SmoothedDMLabel[] = {"Dark_Matter_Density", "Velocity_Dispersion",
 			     "Particle_x-velocity", "Particle_y-velocity",
 			     "Particle_z-velocity"};
diff --git a/src/enzo/Grid_PreparePotentialField.C b/src/enzo/Grid_PreparePotentialField.C
index 65a70198a..55539f35c 100644
--- a/src/enzo/Grid_PreparePotentialField.C
+++ b/src/enzo/Grid_PreparePotentialField.C
@@ -75,7 +75,7 @@ int grid::PreparePotentialField(grid *ParentGrid)
   // Only done in COMMUNICATION_SEND because
   // CommunicationMethodShouldExit() will exit in other modes when the
   // grids are on the same processor.
-  if (ProcessorNumber == ProcessorNumber &&
+  if (MyProcessorNumber == ProcessorNumber &&
       CommunicationDirection != COMMUNICATION_POST_RECEIVE) {
     if (PotentialField != NULL)
       delete [] PotentialField;
diff --git a/src/enzo/hydro_rk/Collapse3DInitialize.C b/src/enzo/hydro_rk/Collapse3DInitialize.C
index 84d6ad5d5..2a171eaff 100644
--- a/src/enzo/hydro_rk/Collapse3DInitialize.C
+++ b/src/enzo/hydro_rk/Collapse3DInitialize.C
@@ -132,10 +132,10 @@ int Collapse3DInitialize(FILE *fptr, FILE *Outfptr,
       ret += sscanf(line, "SphereType[%d] = %d", &sphere,
 		    &SphereType[sphere]);
     if (sscanf(line, "SphereRadius[%d]", &sphere) > 0)
-      ret += sscanf(line, "SphereRadius[%d] = %"FSYM, &sphere,
+      ret += sscanf(line, "SphereRadius[%d] = %"PSYM, &sphere,
 		    &SphereRadius[sphere]);
     if (sscanf(line, "SphereCoreRadius[%d]", &sphere) > 0)
-      ret += sscanf(line, "SphereCoreRadius[%d] = %"FSYM, &sphere,
+      ret += sscanf(line, "SphereCoreRadius[%d] = %"PSYM, &sphere,
 		    &SphereCoreRadius[sphere]);
     if (sscanf(line, "SphereDensity[%d]", &sphere) > 0)
       ret += sscanf(line, "SphereDensity[%d] = %f", &sphere,
@@ -147,7 +147,7 @@ int Collapse3DInitialize(FILE *fptr, FILE *Outfptr,
       ret += sscanf(line, "SphereSoundVelocity[%d] = %f", &sphere,
 		    &SphereSoundVelocity[sphere]);
     if (sscanf(line, "SpherePosition[%d]", &sphere) > 0)
-      ret += sscanf(line, "SpherePosition[%d] = %"FSYM" %"FSYM" %"FSYM, 
+      ret += sscanf(line, "SpherePosition[%d] = %"PSYM" %"PSYM" %"PSYM, 
 		    &sphere, &SpherePosition[sphere][0],
 		    &SpherePosition[sphere][1],
 		    &SpherePosition[sphere][2]);
@@ -354,25 +354,25 @@ int Collapse3DInitialize(FILE *fptr, FILE *Outfptr,
 
   /* Write parameters to parameter output file */
 
-  /*if (MyProcessorNumber == ROOT_PROCESSOR) {
+  if (MyProcessorNumber == ROOT_PROCESSOR) {
     fprintf(Outfptr, "NumberOfSpheres    = %d\n",
 	    n_sphere);
     fprintf(Outfptr, "RefineAtStart      = %d\n",
 	    RefineAtStart);
     fprintf(Outfptr, "UseParticles       = %d\n",
 	    UseParticles);
-    fprintf(Outfptr, "UseColour          = %d\n",
-	    UseColour);
-    fprintf(Outfptr, "InitialTemperature = %f\n",
-	    InitialTemperature);
+    //    fprintf(Outfptr, "UseColour          = %d\n",
+    //	    UseColour);
+    //    fprintf(Outfptr, "InitialTemperature = %f\n",
+    //	    InitialTemperature);
     fprintf(Outfptr, "UniformVelocity    = %f %f %f\n",
 	    UniformVelocity[0], UniformVelocity[1],
 	    UniformVelocity[2]);
     fprintf(Outfptr, "LengthUnit = %f\n",
-            LengthUnit);
+            lenu);
     fprintf(Outfptr, "DensityUnit = %f\n",
-            DensityUnit);
-    for (sphere = 0; sphere < NumberOfSpheres; sphere++) {
+            rhou);
+    for (sphere = 0; sphere < n_sphere; sphere++) {
       fprintf(Outfptr, "SphereType[%d] = %d\n", sphere,
 	      SphereType[sphere]);
       fprintf(Outfptr, "SphereRadius[%d] = %"GOUTSYM"\n", sphere,
@@ -381,16 +381,16 @@ int Collapse3DInitialize(FILE *fptr, FILE *Outfptr,
 	      SphereCoreRadius[sphere]);
       fprintf(Outfptr, "SphereDensity[%d] = %f\n", sphere,
 	      SphereDensity[sphere]);
-      fprintf(Outfptr, "SphereTemperature[%d] = %f\n", sphere,
-	      SphereTemperature[sphere]);
+      fprintf(Outfptr, "SphereSoundVelocity[%d] = %f\n", sphere,
+	      SphereSoundVelocity[sphere]);
       fprintf(Outfptr, "SpherePosition[%d] = ", sphere);
       WriteListOfFloats(Outfptr, MetaData.TopGridRank,
 			SpherePosition[sphere]);
       fprintf(Outfptr, "SphereVelocity[%d] = ", sphere);
       WriteListOfFloats(Outfptr, MetaData.TopGridRank,
 			SphereVelocity[sphere]);
-      fprintf(Outfptr, "FracKeplarianRot[%d] = %"GOUTSYM"\n", sphere,
-              FracKeplarianRot[sphere]);
+      //      fprintf(Outfptr, "FracKeplarianRot[%d] = %"GOUTSYM"\n", sphere,
+      //              FracKeplarianRot[sphere]);
       fprintf(Outfptr, "SphereTurbulence[%d] = %"GOUTSYM"\n", sphere,
               SphereTurbulence[sphere]);
       fprintf(Outfptr, "SphereCutOff[%d] = %"GOUTSYM"\n", sphere,
@@ -402,7 +402,7 @@ int Collapse3DInitialize(FILE *fptr, FILE *Outfptr,
       fprintf(Outfptr, "SphereNumShells[%d] = %d\n", sphere,
               SphereNumShells[sphere]);
     }
-    }*/
+    }
 
   return SUCCESS;
 
diff --git a/src/enzo/hydro_rk/Grid_Collapse3DInitializeGrid.C b/src/enzo/hydro_rk/Grid_Collapse3DInitializeGrid.C
index d674e59fc..9adb96b29 100644
--- a/src/enzo/hydro_rk/Grid_Collapse3DInitializeGrid.C
+++ b/src/enzo/hydro_rk/Grid_Collapse3DInitializeGrid.C
@@ -214,7 +214,7 @@ int grid::Collapse3DInitializeGrid(int n_sphere,
 	      eint = pow(cs_sphere[sphere], 2)/(Gamma-1.0);
 	      vel[0] = -omega_sphere[sphere]*ypos;
 	      vel[1] = omega_sphere[sphere]*xpos;
-	      //printf("cs=%g, vel[0]=%g ", cs, vel[0]);
+	      printf("cs=%g, vel[0]=%g ", cs, vel[0]);
 	      double mach_turb = 0.0;
 	      vel[0] += mach_turb*Gaussian(cs);
 	      vel[1] += mach_turb*Gaussian(cs);
diff --git a/src/enzo/hydro_rk/Grid_MHDTurbulenceInitializeGrid.C b/src/enzo/hydro_rk/Grid_MHDTurbulenceInitializeGrid.C
index 1bb206533..a213e2e04 100644
--- a/src/enzo/hydro_rk/Grid_MHDTurbulenceInitializeGrid.C
+++ b/src/enzo/hydro_rk/Grid_MHDTurbulenceInitializeGrid.C
@@ -143,7 +143,7 @@ int grid::MHDTurbulenceInitializeGrid(float rho_medium, float cs_medium, float m
   float eint, h, dpdrho, dpde, cs;
   eint = cs_medium*cs_medium/(Gamma-1.0);
   FLOAT xc = 0.5, yc = 0.5, zc = 0.5, x, y, z, r;
-  FLOAT rs = 0.3;
+  FLOAT rs = 1.; // 0.3;
   n=0;
   for (int k = 0; k < GridDimension[2]; k++) {
     for (int j = 0; j < GridDimension[1]; j++) {
diff --git a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C
index 53a0c504c..a6d83a4cd 100644
--- a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C
+++ b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C
@@ -32,7 +32,7 @@ void Turbulence_Generator(float **vel, int dim0, int dim1, int dim2, int ind,
 
 int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FLOAT CloudRadius, 
 				   float CloudMachNumber, float CloudAngularVelocity, float InitialBField,
-				   int SetTurbulence, int CloudType, int TurbulenceSeed, int level,
+				   int SetTurbulence, int CloudType, int TurbulenceSeed, int PutSink, int level,
 				   int SetBaryonFields)
 {
 
@@ -100,16 +100,17 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL
         FieldType[kdissH2INum    = NumberOfBaryonFields++] = kdissH2I;
       }
     }
+    
+    if (RadiationPressure) {
+      FieldType[RPresNum1 = NumberOfBaryonFields++] = RadPressure0;
+      FieldType[RPresNum2 = NumberOfBaryonFields++] = RadPressure1;
+      FieldType[RPresNum3 = NumberOfBaryonFields++] = RadPressure2;
+    }
+    
+    NumberOfPhotonPackages = 0;
+    PhotonPackages-> NextPackage= NULL;
+    
   }
-
-  if (RadiationPressure && RadiativeTransfer) {
-    FieldType[RPresNum1 = NumberOfBaryonFields++] = RadPressure0;
-    FieldType[RPresNum2 = NumberOfBaryonFields++] = RadPressure1;
-    FieldType[RPresNum3 = NumberOfBaryonFields++] = RadPressure2;
-  }
-
-  NumberOfPhotonPackages = 0;
-  PhotonPackages-> NextPackage= NULL;
 #endif
 
   int idrivex, idrivey, idrivez;
@@ -551,9 +552,8 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL
 
   /* Put a sink particle if we are studying massive star formation */
 
-  int PutSink = 0;
   if (PutSink == 1 && level == MaximumRefinementLevel) {
-
+      NumberOfParticleAttributes = 6;
     double mass_p = 20.0*1.989e33;
     mass_p /= MassUnits;
     double dx = CellWidth[0][0];
@@ -563,6 +563,8 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL
 
     double dxm = dx / pow(2.0, MaximumRefinementLevel);
 
+    printf("Adding a Sink Particle. \n");
+
     NumberOfParticles = 1;
     NumberOfStars = 1;
     //    MaximumParticleNumber = 1;
@@ -580,6 +582,12 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL
     ParticleAttribute[0][0] = 0.0; // creation time             
     ParticleAttribute[1][0] = 0;
     ParticleAttribute[2][0] = mass_p;
+
+    if (StellarWindFeedback) {
+      ParticleAttribute[3][0] = 1.0;  
+      ParticleAttribute[4][0] = 0.0;
+      ParticleAttribute[5][0] = 0.0;
+    }
   }
 
   /*  printf("XXX PutSinkParticle = %d\n", PutSinkParticle);
diff --git a/src/enzo/hydro_rk/MHDTurbulenceInitialize.C b/src/enzo/hydro_rk/MHDTurbulenceInitialize.C
index be92b5506..17becd5d1 100644
--- a/src/enzo/hydro_rk/MHDTurbulenceInitialize.C
+++ b/src/enzo/hydro_rk/MHDTurbulenceInitialize.C
@@ -61,7 +61,7 @@ int MHDTurbulenceInitialize(FILE *fptr, FILE *Outfptr,
 
   int RefineAtStart   = TRUE;
   int RandomSeed = 1;
-  FLOAT rho_medium=1.0, cs=1.0, mach=1.0, Bnaught=0.0;
+  float rho_medium=1.0, cs=1.0, mach=1.0, Bnaught=0.0;
 
   /* read input from file */
   rewind(fptr);
@@ -139,6 +139,7 @@ int MHDTurbulenceInitialize(FILE *fptr, FILE *Outfptr,
     
     v_rms = sqrt(v_rms/Volume); // actuall v_rms
     fac = cs*mach/v_rms;
+
     CurrentGrid = &TopGrid;
     while (CurrentGrid != NULL) {
       if (CurrentGrid->GridData->NormalizeVelocities(fac) == FAIL) {
diff --git a/src/enzo/hydro_rk/TurbulenceInitialize.C b/src/enzo/hydro_rk/TurbulenceInitialize.C
index 9b0eccf47..b257d432f 100644
--- a/src/enzo/hydro_rk/TurbulenceInitialize.C
+++ b/src/enzo/hydro_rk/TurbulenceInitialize.C
@@ -89,6 +89,7 @@ int TurbulenceInitialize(FILE *fptr, FILE *Outfptr,
   /* set default parameters */
 
   int RefineAtStart   = TRUE;
+  int PutSink         = FALSE;
   int SetTurbulence = TRUE;
   int RandomSeed = 52761;
   float CloudDensity=1.0, CloudSoundSpeed=1.0, CloudMachNumber=1.0, CloudAngularVelocity = 0.0, InitialBField = 0.0;
@@ -102,6 +103,7 @@ int TurbulenceInitialize(FILE *fptr, FILE *Outfptr,
     ret = 0;
 
     ret += sscanf(line, "RefineAtStart = %"ISYM, &RefineAtStart);
+    ret += sscanf(line, "PutSink = %"ISYM, &PutSink);
     ret += sscanf(line, "Density = %"FSYM, &CloudDensity);
     ret += sscanf(line, "SoundVelocity = %"FSYM, &CloudSoundSpeed);
     ret += sscanf(line, "MachNumber = %"FSYM, &CloudMachNumber);
@@ -143,7 +145,7 @@ printf("Plasma beta=%g\n", CloudDensity*CloudSoundSpeed*CloudSoundSpeed/(Initial
   while (CurrentGrid != NULL) {
     if (CurrentGrid->GridData->TurbulenceInitializeGrid(
                 CloudDensity, CloudSoundSpeed, CloudRadius, CloudMachNumber, CloudAngularVelocity, InitialBField, 
-		SetTurbulence, CloudType, RandomSeed, 0, SetBaryonFields) == FAIL) {
+		SetTurbulence, CloudType, RandomSeed, PutSink, 0, SetBaryonFields) == FAIL) {
       fprintf(stderr, "Error in TurbulenceInitializeGrid.\n");
       return FAIL;
     }
@@ -223,7 +225,7 @@ printf("Plasma beta=%g\n", CloudDensity*CloudSoundSpeed*CloudSoundSpeed/(Initial
       while (Temp != NULL) {
 	if (Temp->GridData->TurbulenceInitializeGrid(
 		  CloudDensity, CloudSoundSpeed, CloudRadius, fac, CloudAngularVelocity, InitialBField,
-		  SetTurbulence, CloudType, RandomSeed, level+1, SetBaryonFields) == FAIL) {
+		  SetTurbulence, CloudType, RandomSeed, PutSink, level+1, SetBaryonFields) == FAIL) {
 	  fprintf(stderr, "Error in TurbulenceInitializeGrid.\n");
 	  return FAIL;
 	}
diff --git a/src/enzo/macros_and_parameters.h b/src/enzo/macros_and_parameters.h
index b712182f9..e4f69d4f3 100644
--- a/src/enzo/macros_and_parameters.h
+++ b/src/enzo/macros_and_parameters.h
@@ -71,7 +71,7 @@
 
 #define MAX_STATIC_REGIONS               1000
 
-#ifdef WINDS
+#ifdef WINDS 
 #define MAX_NUMBER_OF_PARTICLE_ATTRIBUTES  6
 #else
 #define MAX_NUMBER_OF_PARTICLE_ATTRIBUTES  3

From 4f9e20d9cffe471d9b58b4cd2b61ef83ace1fe49 Mon Sep 17 00:00:00 2001
From: tabel <none@none>
Date: Tue, 6 Oct 2009 18:22:05 -0700
Subject: [PATCH 03/12] cleaned up a bit the use of UsePhysicalUnit so that
 also the EOSSoundSpeed is specified in cgs and then just converted once and
 not everytime in EOS.h

--HG--
branch : week-of-code
---
 src/enzo/Grid_ComputeTimeStep.C               |  2 +-
 src/enzo/ReadParameterFile.C                  | 51 ++++++++++---------
 src/enzo/hydro_rk/EOS.h                       |  4 +-
 .../hydro_rk/Grid_TurbulenceInitializeGrid.C  | 18 +++++--
 4 files changed, 44 insertions(+), 31 deletions(-)

diff --git a/src/enzo/Grid_ComputeTimeStep.C b/src/enzo/Grid_ComputeTimeStep.C
index 1d6f018c9..1a1dcc817 100644
--- a/src/enzo/Grid_ComputeTimeStep.C
+++ b/src/enzo/Grid_ComputeTimeStep.C
@@ -420,7 +420,7 @@ float grid::ComputeTimeStep()
     if (HydroMethod != MHD_RK && NumberOfBaryonFields > 0)
       printf("Bar = %"FSYM" ", dtBaryons);
     if (HydroMethod == MHD_RK)
-      printf("dtMHD = %"FSYM" ", dtMHD);
+      printf("dtMHD = %e ", dtMHD);
     if (HydroMethod == Zeus_Hydro)
       printf("Vis = %"FSYM" ", dtViscous);
     if (ComovingCoordinates)
diff --git a/src/enzo/ReadParameterFile.C b/src/enzo/ReadParameterFile.C
index 46d37f4cd..5b3222a87 100644
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -823,32 +823,35 @@ int ReadParameterFile(FILE *fptr, TopGridData &MetaData, float *Initialdt)
     GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, &TimeUnits, &VelocityUnits, 
 	     &MassUnits, MetaData.Time);
     PressureUnits = DensityUnits*pow(VelocityUnits,2);
-  }
 
-  /* Change input physical parameters into code units */
-
-  double mh = 1.6726e-24;
-  double uheat = pow(VelocityUnits,2)*2.0*mh/TimeUnits;
-  PhotoelectricHeating /= uheat;
-  StarMakerOverDensityThreshold /= DensityUnits;
-  //  StarEnergyFeedbackRate = StarEnergyFeedbackRate/pow(LengthUnits,2)*pow(TimeUnits,3);
-
-  if (SinkMergeDistance > 1.0)
-    SinkMergeDistance /= LengthUnits;
-  //printf(" \n SinkMergeDistance = %"FSYM"\n \n", SinkMergeDistance);
-  SmallRho /= DensityUnits;
-  SmallP /= PressureUnits;
-  SmallT /= TemperatureUnits;
-  MaximumAlvenSpeed /= VelocityUnits;
-  float h, cs, dpdrho, dpde;
-  EOS(SmallP, SmallRho, SmallEint, h, cs, dpdrho, dpde, EOSType, 1);
-  if (debug && (HydroMethod == HD_RK || HydroMethod == MHD_RK))
-    printf("smallrho=%g, smallp=%g, smalleint=%g, PressureUnits=%g, MaximumAlvenSpeed=%g\n",
-	   SmallRho, SmallP, SmallEint, PressureUnits, MaximumAlvenSpeed);
-  for (int i = 0; i < MAX_FLAGGING_METHODS; i++) {
-    if (MinimumMassForRefinement[i] != FLOAT_UNDEFINED) {
-      MinimumMassForRefinement[i] /= MassUnits;
+
+    /* Change input physical parameters into code units */
+    
+    double mh = 1.6726e-24;
+    double uheat = pow(VelocityUnits,2)*2.0*mh/TimeUnits;
+    PhotoelectricHeating /= uheat;
+    StarMakerOverDensityThreshold /= DensityUnits;
+    //  StarEnergyFeedbackRate = StarEnergyFeedbackRate/pow(LengthUnits,2)*pow(TimeUnits,3);
+    
+    if (SinkMergeDistance > 1.0)
+      SinkMergeDistance /= LengthUnits;
+    //printf(" \n SinkMergeDistance = %"FSYM"\n \n", SinkMergeDistance);
+    SmallRho /= DensityUnits;
+    SmallP /= PressureUnits;
+    SmallT /= TemperatureUnits;
+    MaximumAlvenSpeed /= VelocityUnits;
+    EOSSoundSpeed /=  VelocityUnits;
+    float h, cs, dpdrho, dpde;
+    EOS(SmallP, SmallRho, SmallEint, h, cs, dpdrho, dpde, EOSType, 1);
+    if (debug && (HydroMethod == HD_RK || HydroMethod == MHD_RK))
+      printf("smallrho=%g, smallp=%g, smalleint=%g, PressureUnits=%g, MaximumAlvenSpeed=%g\n",
+	     SmallRho, SmallP, SmallEint, PressureUnits, MaximumAlvenSpeed);
+    for (int i = 0; i < MAX_FLAGGING_METHODS; i++) {
+      if (MinimumMassForRefinement[i] != FLOAT_UNDEFINED) {
+	MinimumMassForRefinement[i] /= MassUnits;
+      }
     }
+    
   }
 
   if (!ComovingCoordinates && UsePhysicalUnit) {
diff --git a/src/enzo/hydro_rk/EOS.h b/src/enzo/hydro_rk/EOS.h
index ef4e8c4fc..287258236 100644
--- a/src/enzo/hydro_rk/EOS.h
+++ b/src/enzo/hydro_rk/EOS.h
@@ -42,7 +42,7 @@ inline void EOS(float &p, float &rho, float &e, float &h, float &cs, float &dpdr
     GetUnits(&denu, &lenu, &tempu, &tu, &velu, 1);
     double c_s = EOSSoundSpeed;
     double rho_cr = EOSCriticalDensity;
-    c_s /= velu;
+    //    c_s /= velu;
     rho_cr /= denu;
 
     cs = c_s*sqrt(1.0 + EOSGamma*pow(rho/rho_cr, EOSGamma-1.0));
@@ -59,7 +59,7 @@ inline void EOS(float &p, float &rho, float &e, float &h, float &cs, float &dpdr
     GetUnits(&denu, &lenu, &tempu, &tu, &velu, 1);
     double c_s = EOSSoundSpeed;
     double rho_cr = EOSCriticalDensity;
-    c_s /= velu;
+    //    c_s /= velu;
     rho_cr /= denu;
     
     if (rho <= rho_cr) {
diff --git a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C
index a6d83a4cd..c7644846b 100644
--- a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C
+++ b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C
@@ -193,7 +193,9 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL
   /* assume isothermal initially */
 
   CloudPressure = CloudDensity * pow(CloudSoundSpeed,2);
-  CloudInternalEnergy = pow(CloudSoundSpeed,2) / (Gamma - 1.0);
+  //  CloudInternalEnergy = pow(CloudSoundSpeed,2) / (Gamma - 1.0);
+  float h, cs, dpdrho, dpde;
+  EOS(CloudPressure, CloudDensity, CloudInternalEnergy, h, cs, dpdrho, dpde, EOSType, 1);
 
   float InitialFractionHII = 1.2e-5;
   float InitialFractionHeII = 1.0e-14;
@@ -293,7 +295,7 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL
 
 	  if (CloudType == 0) {
 	    Density = CloudDensity/10.0;
-	    eint = CloudInternalEnergy;
+	    eint = CloudInternalEnergy*10.;
 	  }
 
 	  if (CloudType == 1) {
@@ -553,7 +555,7 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL
   /* Put a sink particle if we are studying massive star formation */
 
   if (PutSink == 1 && level == MaximumRefinementLevel) {
-      NumberOfParticleAttributes = 6;
+
     double mass_p = 20.0*1.989e33;
     mass_p /= MassUnits;
     double dx = CellWidth[0][0];
@@ -561,14 +563,16 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL
     double t_dyn = sqrt(3*M_PI/(6.672e-8*den_p*DensityUnits));
     t_dyn /= TimeUnits;
 
-    double dxm = dx / pow(2.0, MaximumRefinementLevel);
+    double dxm = dx / pow(RefineBy, MaximumRefinementLevel);
 
     printf("Adding a Sink Particle. \n");
 
     NumberOfParticles = 1;
     NumberOfStars = 1;
     //    MaximumParticleNumber = 1;
+    if (StellarWindFeedback) NumberOfParticleAttributes = 6;
     this->AllocateNewParticles(NumberOfParticles);
+
     ParticleMass[0] = den_p;
     ParticleNumber[0] = 0;
     ParticleType[0] = PARTICLE_TYPE_MUST_REFINE;
@@ -588,6 +592,12 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL
       ParticleAttribute[4][0] = 0.0;
       ParticleAttribute[5][0] = 0.0;
     }
+
+    for (i = 0; i< MAX_DIMENSION+1; i++){
+      ParticleAcceleration[i] = NULL;
+    }
+    this->ClearParticleAccelerations();
+
   }
 
   /*  printf("XXX PutSinkParticle = %d\n", PutSinkParticle);

From 8a5ae83847cf38c54c5430b0c66752eb1d9b6f55 Mon Sep 17 00:00:00 2001
From: tabel <none@none>
Date: Tue, 6 Oct 2009 18:46:39 -0700
Subject: [PATCH 04/12] Changed condition that determines when the jet
 direction for stellar wind feedback gets set to allow one to specify the jet
 direction in test problems.

--HG--
branch : week-of-code
---
 src/enzo/star_maker8.C | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/enzo/star_maker8.C b/src/enzo/star_maker8.C
index bb6b2089c..01616e6ba 100644
--- a/src/enzo/star_maker8.C
+++ b/src/enzo/star_maker8.C
@@ -323,8 +323,8 @@ int star_maker8(int *nx, int *ny, int *nz, int *size, float *d, float *te, float
       if (mpold[bb]*pow(*dx,3)*umass < StellarWindTurnOnMass && (*t - tcpold[bb])*(*t1) < 1e5*3.1557e7) continue;
 
       int first = 0;
-      if (dmold[bb] > 0.99*mpold[bb]*pow(*dx,3)) first = 1;
-      
+      //      if (dmold[bb] > 0.99*mpold[bb]*pow(*dx,3)) first = 1;
+      if (nx_jet[bb]+ny_jet[bb]+nz_jet[bb] < 0.1) first = 1;
 
       /* Decide whether the current grid contains the whole supercell */
 

From 90fd0bd7691c8ff4622e00c8ca6af2ed58cf6bc6 Mon Sep 17 00:00:00 2001
From: tabel <none@none>
Date: Sat, 10 Oct 2009 14:46:40 -0700
Subject: [PATCH 05/12] Some cleanup in EvolveLevel_RK2

--HG--
branch : week-of-code
---
 input/metal_cool.dat                          |   2 +-
 input/metal_cool_ratios.dat                   |  18 ++--
 src/enzo/ComputePotentialFieldLevelZero.C     |   2 +-
 src/enzo/DepositBaryons.C                     |   5 +-
 src/enzo/EvolveLevel.C                        |   5 +-
 src/enzo/Grid_CopyPotentialToBaryonField.C    |   2 +-
 src/enzo/Grid_ZeroSolutionUnderSubgrid.C      |   2 +-
 src/enzo/PrepareGravitatingMassField.C        |   2 +-
 src/enzo/hydro_rk/EvolveLevel_RK2.C           | 102 +++++++++---------
 .../hydro_rk/Grid_TurbulenceInitializeGrid.C  |  21 ++--
 10 files changed, 84 insertions(+), 77 deletions(-)

diff --git a/input/metal_cool.dat b/input/metal_cool.dat
index dd4be6d22..fc7446103 100644
--- a/input/metal_cool.dat
+++ b/input/metal_cool.dat
@@ -168,7 +168,7 @@
   6.9499E-29  6.9503E-29  6.9509E-29  6.9516E-29  6.9525E-29  6.9536E-29  6.9550E-29  6.9567E-29  6.9590E-29  6.9618E-29  6.9654E-29  6.9699E-29  6.9756E-29  6.9828E-29  6.9918E-29  7.0033E-29  7.0179E-29  7.0362E-29  7.0594E-29  7.0887E-29  7.1257E-29  7.1725E-29  7.2317E-29  7.3064E-29  7.4009E-29  7.5202E-29  7.6711E-29  7.8618E-29  8.1028E-29  8.4073E-29  8.7922E-29  9.2786E-29  9.8934E-29  1.0670E-28  1.1652E-28  1.2893E-28  1.4461E-28  1.6444E-28  1.8949E-28  2.2115E-28  2.6116E-28  3.1172E-28  3.7563E-28  4.5639E-28  5.5846E-28  6.8745E-28  8.5045E-28  1.0564E-27  1.3167E-27  1.6457E-27  2.0613E-27  2.5864E-27  3.2498E-27  4.0879E-27  5.1465E-27  6.4834E-27  8.1716E-27  1.0303E-26  1.2992E-26  1.6384E-26
   7.3342E-29  7.3347E-29  7.3353E-29  7.3360E-29  7.3370E-29  7.3382E-29  7.3397E-29  7.3416E-29  7.3441E-29  7.3471E-29  7.3510E-29  7.3559E-29  7.3621E-29  7.3699E-29  7.3798E-29  7.3923E-29  7.4081E-29  7.4280E-29  7.4532E-29  7.4851E-29  7.5253E-29  7.5762E-29  7.6405E-29  7.7217E-29  7.8244E-29  7.9541E-29  8.1181E-29  8.3253E-29  8.5873E-29  8.9183E-29  9.3366E-29  9.8653E-29  1.0534E-28  1.1378E-28  1.2445E-28  1.3794E-28  1.5499E-28  1.7653E-28  2.0376E-28  2.3817E-28  2.8166E-28  3.3662E-28  4.0608E-28  4.9387E-28  6.0481E-28  7.4501E-28  9.2218E-28  1.1461E-27  1.4290E-27  1.7865E-27  2.2383E-27  2.8090E-27  3.5302E-27  4.4411E-27  5.5918E-27  7.0449E-27  8.8799E-27  1.1196E-26  1.4119E-26  1.7807E-26
   7.7359E-29  7.7364E-29  7.7371E-29  7.7379E-29  7.7389E-29  7.7402E-29  7.7419E-29  7.7439E-29  7.7466E-29  7.7499E-29  7.7541E-29  7.7594E-29  7.7661E-29  7.7746E-29  7.7853E-29  7.7989E-29  7.8159E-29  7.8376E-29  7.8649E-29  7.8994E-29  7.9430E-29  7.9981E-29  8.0678E-29  8.1558E-29  8.2671E-29  8.4077E-29  8.5854E-29  8.8100E-29  9.0939E-29  9.4526E-29  9.9060E-29  1.0479E-28  1.1203E-28  1.2118E-28  1.3275E-28  1.4737E-28  1.6584E-28  1.8919E-28  2.1869E-28  2.5599E-28  3.0312E-28  3.6268E-28  4.3796E-28  5.3309E-28  6.5332E-28  8.0526E-28  9.9726E-28  1.2399E-27  1.5465E-27  1.9340E-27  2.4235E-27  3.0421E-27  3.8236E-27  4.8108E-27  6.0578E-27  7.6327E-27  9.6212E-27  1.2132E-26  1.5300E-26  1.9296E-26
-  8.1555E-29  8.1560E-29  8.1567E-29  8.1576E-29  8.1588E-29  8.1602E-29  8.1619E-29  8.1642E-29  8.1670E-29  8.1706E-29  8.1752E-29  8.1809E-29  8.1882E-29  8.1973E-29  8.2089E-29  8.2235E-29  8.2420E-29  8.2654E-29  8.2949E-29  8.3322E-29  8.3793E-29  8.4388E-29  8.5141E-29  8.6092E-29  8.7295E-29  8.8814E-29  9.0735E-29  9.3161E-29  9.6229E-29  1.0010E-28  1.0500E-28  1.1119E-28  1.1902E-28  1.2891E-28  1.4141E-28  1.5720E-28  1.7716E-28  2.0239E-28  2.3428E-28  2.7458E-28  3.2550E-28  3.8987E-28  4.7121E-28  5.7401E-28  7.0393E-28  8.6811E-28  1.0756E-27  1.3378E-27  1.6691E-27  2.0878E-27  2.6168E-27  3.2852E-27  4.1297E-27  5.1965E-27  6.5440E-27  8.2458E-27  1.0395E-26  1.3107E-26  1.6531E-26  2.0850E-26
+  8.1555E-29  8.1560E-29  8.1568E-29  8.1576E-29  8.1588E-29  8.1602E-29  8.1619E-29  8.1642E-29  8.1670E-29  8.1706E-29  8.1752E-29  8.1809E-29  8.1882E-29  8.1973E-29  8.2089E-29  8.2235E-29  8.2420E-29  8.2654E-29  8.2949E-29  8.3322E-29  8.3793E-29  8.4388E-29  8.5141E-29  8.6092E-29  8.7295E-29  8.8814E-29  9.0735E-29  9.3161E-29  9.6229E-29  1.0010E-28  1.0500E-28  1.1119E-28  1.1902E-28  1.2891E-28  1.4141E-28  1.5720E-28  1.7716E-28  2.0239E-28  2.3428E-28  2.7458E-28  3.2550E-28  3.8987E-28  4.7121E-28  5.7401E-28  7.0393E-28  8.6811E-28  1.0756E-27  1.3378E-27  1.6691E-27  2.0878E-27  2.6168E-27  3.2852E-27  4.1297E-27  5.1965E-27  6.5440E-27  8.2458E-27  1.0395E-26  1.3107E-26  1.6531E-26  2.0850E-26
   8.5935E-29  8.5941E-29  8.5949E-29  8.5958E-29  8.5970E-29  8.5986E-29  8.6005E-29  8.6029E-29  8.6060E-29  8.6098E-29  8.6147E-29  8.6209E-29  8.6287E-29  8.6386E-29  8.6511E-29  8.6668E-29  8.6867E-29  8.7119E-29  8.7437E-29  8.7839E-29  8.8347E-29  8.8988E-29  8.9800E-29  9.0825E-29  9.2121E-29  9.3757E-29  9.5827E-29  9.8442E-29  1.0175E-28  1.0592E-28  1.1120E-28  1.1787E-28  1.2631E-28  1.3696E-28  1.5043E-28  1.6745E-28  1.8896E-28  2.1615E-28  2.5051E-28  2.9393E-28  3.4881E-28  4.1817E-28  5.0582E-28  6.1659E-28  7.5658E-28  9.3351E-28  1.1571E-27  1.4396E-27  1.7966E-27  2.2478E-27  2.8179E-27  3.5381E-27  4.4481E-27  5.5977E-27  7.0497E-27  8.8836E-27  1.1199E-26  1.4122E-26  1.7811E-26  2.2465E-26
   9.0505E-29  9.0511E-29  9.0519E-29  9.0529E-29  9.0542E-29  9.0559E-29  9.0579E-29  9.0605E-29  9.0638E-29  9.0680E-29  9.0732E-29  9.0799E-29  9.0883E-29  9.0989E-29  9.1123E-29  9.1292E-29  9.1506E-29  9.1777E-29  9.2118E-29  9.2550E-29  9.3096E-29  9.3786E-29  9.4657E-29  9.5759E-29  9.7151E-29  9.8910E-29  1.0113E-28  1.0394E-28  1.0750E-28  1.1198E-28  1.1766E-28  1.2483E-28  1.3389E-28  1.4534E-28  1.5981E-28  1.7810E-28  2.0122E-28  2.3044E-28  2.6736E-28  3.1402E-28  3.7300E-28  4.4753E-28  5.4172E-28  6.6077E-28  8.1121E-28  1.0013E-27  1.2416E-27  1.5452E-27  1.9289E-27  2.4137E-27  3.0263E-27  3.8004E-27  4.7783E-27  6.0136E-27  7.5741E-27  9.5448E-27  1.2033E-26  1.5175E-26  1.9139E-26  2.4141E-26
   9.5268E-29  9.5275E-29  9.5283E-29  9.5294E-29  9.5308E-29  9.5326E-29  9.5348E-29  9.5376E-29  9.5411E-29  9.5456E-29  9.5512E-29  9.5583E-29  9.5673E-29  9.5787E-29  9.5930E-29  9.6112E-29  9.6341E-29  9.6631E-29  9.6997E-29  9.7460E-29  9.8045E-29  9.8784E-29  9.9719E-29  1.0090E-28  1.0239E-28  1.0428E-28  1.0666E-28  1.0967E-28  1.1348E-28  1.1829E-28  1.2437E-28  1.3205E-28  1.4177E-28  1.5404E-28  1.6955E-28  1.8916E-28  2.1393E-28  2.4524E-28  2.8482E-28  3.3483E-28  3.9804E-28  4.7793E-28  5.7888E-28  7.0647E-28  8.6771E-28  1.0715E-27  1.3290E-27  1.6544E-27  2.0657E-27  2.5853E-27  3.2419E-27  4.0715E-27  5.1196E-27  6.4437E-27  8.1162E-27  1.0228E-26  1.2896E-26  1.6263E-26  2.0512E-26  2.5873E-26
diff --git a/input/metal_cool_ratios.dat b/input/metal_cool_ratios.dat
index 91739c750..43f800093 100644
--- a/input/metal_cool_ratios.dat
+++ b/input/metal_cool_ratios.dat
@@ -15229,7 +15229,7 @@
       269.60  0.00026478     0.26249    0.056946  3.3441E-08     0.14765     0.49034  1.1984E-05    0.023418    0.013208  6.2553E-08   0.0056664
       278.02  0.00026417     0.26498    0.055478  3.2602E-08     0.14403     0.49358  1.1526E-05    0.023874    0.012383  5.9447E-08   0.0054052
       286.70  0.00026364     0.26727    0.054067  3.1795E-08     0.14054     0.49676  1.1094E-05    0.024328    0.011609  5.6483E-08   0.0051542
-      295.65  0.00026317     0.26936    0.052711  3.1020E-08     0.13719     0.49988  1.0686E-05    0.024780    0.010884  5.3655E-08   0.0049133
+      295.65  0.00026317     0.26936    0.052711  3.1020E-08     0.13719     0.49988  1.0686E-05    0.024780    0.010884  5.3655E-08   0.0049132
       304.89  0.00026277     0.27127    0.051408  3.0274E-08     0.13397     0.50297  1.0301E-05    0.025230    0.010205  5.0959E-08   0.0046821
       314.41  0.00026245     0.27299    0.050155  2.9557E-08     0.13087     0.50601  9.9368E-06    0.025679   0.0095678  4.8389E-08   0.0044605
       324.23  0.00026219     0.27453    0.048951  2.8867E-08     0.12788     0.50902  9.5927E-06    0.026127   0.0089710  4.5941E-08   0.0042483
@@ -22456,7 +22456,7 @@
       427.61  0.00011962     0.67268    0.017967  1.0661E-08    0.047460     0.24324  3.2624E-06    0.013828   0.0027875  1.9419E-08   0.0019153
       440.97  0.00012088     0.66985    0.017754  1.0541E-08    0.046954     0.24669  3.1993E-06    0.014150   0.0026424  1.8585E-08   0.0018359
       454.74  0.00012221     0.66688    0.017553  1.0429E-08    0.046477     0.25022  3.1401E-06    0.014479   0.0025056  1.7787E-08   0.0017599
-      468.94  0.00012361     0.66375    0.017362  1.0322E-08    0.046027     0.25385  3.0844E-06    0.014816   0.0023765  1.7025E-08   0.0016871
+      468.94  0.00012361     0.66375    0.017363  1.0322E-08    0.046027     0.25385  3.0844E-06    0.014816   0.0023765  1.7025E-08   0.0016871
       483.59  0.00012507     0.66048    0.017182  1.0222E-08    0.045603     0.25758  3.0320E-06    0.015160   0.0022545  1.6296E-08   0.0016172
       498.69  0.00012660     0.65706    0.017012  1.0127E-08    0.045203     0.26139  2.9827E-06    0.015513   0.0021393  1.5599E-08   0.0015503
       514.26  0.00012820     0.65351    0.016850  1.0037E-08    0.044827     0.26529  2.9363E-06    0.015873   0.0020304  1.4932E-08   0.0014861
@@ -24214,7 +24214,7 @@
       107.16  0.00010882     0.64946    0.049681  2.8485E-08     0.12330     0.12351  2.0322E-05   0.0036492    0.037913  1.4578E-07    0.012358
       110.51  0.00010536     0.66364    0.046751  2.6830E-08     0.11623     0.12276  1.8486E-05   0.0037098    0.034990  1.3792E-07    0.011797
       113.96  0.00010209     0.67686    0.044033  2.5293E-08     0.10966     0.12200  1.6852E-05   0.0037696    0.032301  1.3048E-07    0.011255
-      117.52  9.8990E-05     0.68917    0.041515  2.3869E-08     0.10356     0.12125  1.5397E-05   0.0038291    0.029828  1.2343E-07    0.010733
+      117.52  9.8990E-05     0.68917    0.041515  2.3869E-08     0.10356     0.12125  1.5397E-05   0.0038290    0.029828  1.2343E-07    0.010733
       121.19  9.6072E-05     0.70060    0.039184  2.2549E-08    0.097906     0.12053  1.4098E-05   0.0038884    0.027558  1.1676E-07    0.010232
       124.98  9.3329E-05     0.71119    0.037028  2.1327E-08    0.092669     0.11983  1.2938E-05   0.0039481    0.025474  1.1047E-07   0.0097521
       128.88  9.0760E-05     0.72100    0.035034  2.0196E-08    0.087821     0.11917  1.1900E-05   0.0040084    0.023563  1.0454E-07   0.0092934
@@ -24890,7 +24890,7 @@
       1075.8  0.00011369     0.71539   0.0088105  5.3297E-09    0.024090     0.23408  1.4432E-06    0.016526  0.00048825  4.8598E-09  0.00050420
       1109.4  0.00011631     0.71000   0.0088266  5.3428E-09    0.024160     0.23898  1.4454E-06    0.016961  0.00046715  4.6699E-09  0.00048474
       1144.0  0.00011900     0.70451   0.0088435  5.3562E-09    0.024233     0.24398  1.4481E-06    0.017405  0.00044693  4.4869E-09  0.00046596
-      1179.7  0.00012176     0.69892   0.0088608  5.3701E-09    0.024307     0.24906  1.4511E-06    0.017859  0.00042757  4.3104E-09  0.00044785
+      1179.7  0.00012176     0.69892   0.0088608  5.3701E-09    0.024307     0.24906  1.4511E-06    0.017859  0.00042757  4.3105E-09  0.00044785
       1216.6  0.00012459     0.69323   0.0088787  5.3842E-09    0.024382     0.25422  1.4545E-06    0.018322  0.00040902  4.1404E-09  0.00043039
       1254.6  0.00012748     0.68744   0.0088969  5.3985E-09    0.024459     0.25947  1.4582E-06    0.018794  0.00039124  3.9765E-09  0.00041354
       1293.8  0.00013044     0.68156   0.0089153  5.4130E-09    0.024536     0.26481  1.4623E-06    0.019275  0.00037421  3.8185E-09  0.00039729
@@ -26089,7 +26089,7 @@
       980.96  7.6485E-05     0.80639   0.0062468  3.7715E-09    0.017021     0.15814  1.0359E-06    0.011069  0.00049967  5.2892E-09  0.00055252
       1011.6  7.8440E-05     0.80227   0.0062708  3.7889E-09    0.017110     0.16187  1.0393E-06    0.011390  0.00047957  5.0964E-09  0.00053260
       1043.2  8.0456E-05     0.79803   0.0062982  3.8078E-09    0.017203     0.16569  1.0431E-06    0.011720  0.00046030  4.9103E-09  0.00051337
-      1075.8  8.2534E-05     0.79368   0.0063268  3.8274E-09    0.017300     0.16961  1.0472E-06    0.012058  0.00044180  4.7307E-09  0.00049479
+      1075.8  8.2534E-05     0.79369   0.0063268  3.8274E-09    0.017300     0.16961  1.0472E-06    0.012058  0.00044180  4.7307E-09  0.00049479
       1109.4  8.4675E-05     0.78923   0.0063563  3.8476E-09    0.017400     0.17362  1.0518E-06    0.012406  0.00042405  4.5574E-09  0.00047684
       1144.0  8.6879E-05     0.78468   0.0063867  3.8685E-09    0.017502     0.17772  1.0567E-06    0.012762  0.00040701  4.3900E-09  0.00045951
       1179.7  8.9149E-05     0.78001   0.0064180  3.8898E-09    0.017607     0.18191  1.0619E-06    0.013128  0.00039065  4.2285E-09  0.00044276
@@ -29031,7 +29031,7 @@
       141.34  1.7226E-05     0.93573   0.0055510  3.2086E-09    0.013984    0.028297  1.8709E-06   0.0012453   0.0088311  6.3863E-08   0.0063383
       145.75  1.6717E-05     0.93815   0.0052483  3.0362E-09    0.013241    0.027893  1.7262E-06   0.0012471   0.0081889  6.0409E-08   0.0060152
       150.30  1.6253E-05     0.94035   0.0049724  2.8789E-09    0.012564    0.027530  1.5973E-06   0.0012500   0.0076046  5.7191E-08   0.0057123
-      154.99  1.5832E-05     0.94236   0.0047207  2.7354E-09    0.011946    0.027206  1.4822E-06   0.0012540   0.0070722  5.4190E-08   0.0054284
+      154.99  1.5832E-05     0.94236   0.0047207  2.7354E-09    0.011946    0.027206  1.4822E-06   0.0012540   0.0070722  5.4189E-08   0.0054284
       159.84  1.5451E-05     0.94418   0.0044907  2.6043E-09    0.011381    0.026920  1.3791E-06   0.0012590   0.0065863  5.1387E-08   0.0051619
       164.83  1.5106E-05     0.94585   0.0042803  2.4843E-09    0.010863    0.026670  1.2867E-06   0.0012651   0.0061424  4.8768E-08   0.0049117
       169.97  1.4794E-05     0.94737   0.0040877  2.3743E-09    0.010390    0.026454  1.2036E-06   0.0012723   0.0057360  4.6318E-08   0.0046766
@@ -29679,7 +29679,7 @@
       599.74  1.5098E-05     0.95760   0.0015931  9.5223E-10   0.0042641    0.032648  3.0429E-07   0.0022661  0.00071453  8.5010E-09  0.00089856
       618.47  1.5460E-05     0.95685   0.0015967  9.5502E-10   0.0042787    0.033373  3.0356E-07   0.0023274  0.00068759  8.2121E-09  0.00086834
       637.79  1.5840E-05     0.95607   0.0016013  9.5834E-10   0.0042958    0.034127  3.0307E-07   0.0023910  0.00066185  7.9343E-09  0.00083926
-      657.70  1.6236E-05     0.95525   0.0016067  9.6219E-10   0.0043152    0.034909  3.0281E-07   0.0024570  0.00063724  7.6670E-09  0.00081126
+      657.70  1.6236E-05     0.95525   0.0016067  9.6219E-10   0.0043152    0.034909  3.0281E-07   0.0024571  0.00063724  7.6670E-09  0.00081126
       678.24  1.6651E-05     0.95439   0.0016129  9.6653E-10   0.0043368    0.035721  3.0277E-07   0.0025255  0.00061369  7.4097E-09  0.00078429
       699.43  1.7083E-05     0.95349   0.0016199  9.7137E-10   0.0043607    0.036564  3.0295E-07   0.0025966  0.00059116  7.1620E-09  0.00075832
       721.27  1.7535E-05     0.95256   0.0016278  9.7670E-10   0.0043868    0.037438  3.0333E-07   0.0026703  0.00056959  6.9235E-09  0.00073328
@@ -31435,7 +31435,7 @@
       141.34  7.0775E-06     0.96363   0.0019639  1.1353E-09   0.0049487    0.015167  7.6498E-07  0.00081142   0.0073771  6.0230E-08   0.0060910
       145.75  6.8604E-06     0.96516   0.0018551  1.0734E-09   0.0046818    0.014859  7.0560E-07  0.00080706   0.0068488  5.6977E-08   0.0057774
       150.30  6.6632E-06     0.96656   0.0017562  1.0170E-09   0.0044389    0.014577  6.5279E-07  0.00080345   0.0063683  5.3947E-08   0.0054840
-      154.99  6.4845E-06     0.96785   0.0016661  9.6557E-10   0.0042174    0.014321  6.0569E-07  0.00080056   0.0059304  5.1124E-08   0.0052094
+      154.99  6.4845E-06     0.96785   0.0016661  9.6557E-10   0.0042174    0.014321  6.0569E-07  0.00080056   0.0059303  5.1124E-08   0.0052094
       159.84  6.3227E-06     0.96902   0.0015839  9.1869E-10   0.0040154    0.014088  5.6359E-07  0.00079836   0.0055307  4.8489E-08   0.0049521
       164.83  6.1766E-06     0.97010   0.0015089  8.7587E-10   0.0038308    0.013877  5.2588E-07  0.00079686   0.0051653  4.6028E-08   0.0047109
       169.97  6.0451E-06     0.97109   0.0014403  8.3672E-10   0.0036619    0.013687  4.9201E-07  0.00079604   0.0048307  4.3726E-08   0.0044844
@@ -35631,7 +35631,7 @@
       100.77  2.4882E-06     0.95884  2.8558E-06  1.9735E-12  1.0182E-05    0.012437  4.7718E-07  0.00070138    0.016630  1.1540E-07    0.011373
       103.92  2.3486E-06     0.96162  2.6597E-06  1.8421E-12  9.4969E-06    0.011861  4.2399E-07  0.00068617    0.015143  1.0785E-07    0.010677
       107.16  2.2221E-06     0.96413  2.4838E-06  1.7238E-12  8.8813E-06    0.011330  3.7832E-07  0.00067183    0.013817  1.0092E-07    0.010033
-      110.51  2.1074E-06     0.96642  2.3258E-06  1.6173E-12  8.3272E-06    0.010838  3.3895E-07  0.00065832    0.012633  9.4550E-08   0.0094364
+      110.51  2.1074E-06     0.96642  2.3258E-06  1.6173E-12  8.3272E-06    0.010838  3.3895E-07  0.00065832    0.012633  9.4551E-08   0.0094364
       113.96  2.0034E-06     0.96850  2.1836E-06  1.5212E-12  7.8278E-06    0.010384  3.0489E-07  0.00064562    0.011572  8.8688E-08   0.0088843
       117.52  1.9090E-06     0.97040  2.0557E-06  1.4344E-12  7.3771E-06   0.0099639  2.7532E-07  0.00063368    0.010622  8.3287E-08   0.0083724
       121.19  1.8232E-06     0.97213  1.9403E-06  1.3559E-12  6.9699E-06   0.0095752  2.4957E-07  0.00062248   0.0097676  7.8304E-08   0.0078975
diff --git a/src/enzo/ComputePotentialFieldLevelZero.C b/src/enzo/ComputePotentialFieldLevelZero.C
index 652c68a9a..9bc217ae4 100644
--- a/src/enzo/ComputePotentialFieldLevelZero.C
+++ b/src/enzo/ComputePotentialFieldLevelZero.C
@@ -277,7 +277,7 @@ int ComputePotentialFieldLevelZeroPer(TopGridData *MetaData,
   float coef = GravitationalConstant/a;
 
   //  for (int dim = 0; dim < MetaData->TopGridRank; dim++)
-  //    coef *= (DomainRightEdge[dim] - DomainLeftEdge[dim])/float(DomainDim[dim]);
+  //     coef *= (DomainRightEdge[dim] - DomainLeftEdge[dim])/float(DomainDim[dim]);
 			
   /* Multiply density by Green's function to get potential. */
  
diff --git a/src/enzo/DepositBaryons.C b/src/enzo/DepositBaryons.C
index 93c643cf6..5044612ff 100644
--- a/src/enzo/DepositBaryons.C
+++ b/src/enzo/DepositBaryons.C
@@ -36,8 +36,7 @@ int DepositBaryonsChildren(HierarchyEntry *DepositGrid,
  
 int DepositBaryons(HierarchyEntry *Grid, FLOAT When)
 {
- 
-  /* Get the time and dt for this grid.  Compute time+1/2 dt. */
+   /* Get the time and dt for this grid.  Compute time+1/2 dt. */
  
   FLOAT TimeMidStep =     Grid->GridData->ReturnTime() 
     + When*Grid->GridData->ReturnTimeStep();
@@ -47,7 +46,7 @@ int DepositBaryons(HierarchyEntry *Grid, FLOAT When)
  
     /* Set the under_subgrid field (indicating if a cell is refined or not)
        on this grid. */
- 
+    //  printf("DepositBaryons:\n");
     HierarchyEntry *Temp = Grid->NextGridNextLevel;
     Grid->GridData->ZeroSolutionUnderSubgrid(NULL, ZERO_UNDER_SUBGRID_FIELD);
     while (Temp != NULL) {
diff --git a/src/enzo/EvolveLevel.C b/src/enzo/EvolveLevel.C
index 2bd647a6c..60d9060fb 100644
--- a/src/enzo/EvolveLevel.C
+++ b/src/enzo/EvolveLevel.C
@@ -461,7 +461,7 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 
     //dcc cut second potential cut: Duplicate?
  
-    if (ComovingCoordinates && SelfGravity && WritePotential) {
+    if (SelfGravity && WritePotential) {
       CopyGravPotential = TRUE;
       When = 0.0;
  
@@ -471,7 +471,6 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
       PrepareDensityField(LevelArray, level, MetaData, When);
 #endif  // end FAST_SIB
  
-      CopyGravPotential = FALSE;
  
       for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
         if (level <= MaximumGravityRefinementLevel) {
@@ -483,6 +482,8 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
           Grids[grid1]->GridData->CopyPotentialToBaryonField();
         }
       } //  end loop over grids
+      CopyGravPotential = FALSE;
+
     } // if WritePotential
  
 
diff --git a/src/enzo/Grid_CopyPotentialToBaryonField.C b/src/enzo/Grid_CopyPotentialToBaryonField.C
index c63c6c769..57e2ea132 100644
--- a/src/enzo/Grid_CopyPotentialToBaryonField.C
+++ b/src/enzo/Grid_CopyPotentialToBaryonField.C
@@ -101,7 +101,7 @@ int grid::CopyPotentialToBaryonField()
  
       for (i = 0; i < GridDimension[0]; i++, index++)
       {
-	//	BaryonField[field+1][jj] = GravitatingMassField[index]; // use this for debugging 
+		BaryonField[field+1][jj] = GravitatingMassField[index]; // use this for debugging 
 	BaryonField[field][jj++] = PotentialField[index];
 	// debuggin:
 	maxPot = max(maxPot,PotentialField[index]);
diff --git a/src/enzo/Grid_ZeroSolutionUnderSubgrid.C b/src/enzo/Grid_ZeroSolutionUnderSubgrid.C
index e2f575463..5bc4f685a 100644
--- a/src/enzo/Grid_ZeroSolutionUnderSubgrid.C
+++ b/src/enzo/Grid_ZeroSolutionUnderSubgrid.C
@@ -30,7 +30,7 @@ int grid::ZeroSolutionUnderSubgrid(grid *Subgrid, int FieldsToZero,
 {
  
   /* Return if this grid is not on this processor. */
- 
+  //  printf("Valeu:%f\n", Value); 
   if (AllProcessors == FALSE)
     if (MyProcessorNumber != ProcessorNumber)
       return SUCCESS;
diff --git a/src/enzo/PrepareGravitatingMassField.C b/src/enzo/PrepareGravitatingMassField.C
index b310fa1b1..538ebfab0 100644
--- a/src/enzo/PrepareGravitatingMassField.C
+++ b/src/enzo/PrepareGravitatingMassField.C
@@ -101,7 +101,7 @@ int PrepareGravitatingMassField2(HierarchyEntry *Grid, TopGridData *MetaData,
   // if (CurrentGrid->AddBaryonsToGravitatingMassField() == FAIL) {
       //fprintf(stderr, "  PGMF - DepositBaryons\n");
   if (DepositBaryons(Grid, When) == FAIL) {
-    fprintf(stderr, "Error in grid->AddBaryonsToGravitatingMassField\n");
+    fprintf(stderr, "Error in DepositBaryons\n");
     ENZO_FAIL("");
   }
  
diff --git a/src/enzo/hydro_rk/EvolveLevel_RK2.C b/src/enzo/hydro_rk/EvolveLevel_RK2.C
index 141eeaa0e..959b9aef5 100644
--- a/src/enzo/hydro_rk/EvolveLevel_RK2.C
+++ b/src/enzo/hydro_rk/EvolveLevel_RK2.C
@@ -368,85 +368,63 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 //     }
 
     When = 0.5;
-
+    if (SelfGravity) {
 #ifdef FAST_SIB
-    if (SelfGravity)
-      if (PrepareDensityField(LevelArray, SiblingList,
-			      level, MetaData, When) == FAIL) {
-	fprintf(stderr, "Error in PrepareDensityField.\n");
-	ENZO_FAIL("");
-      }
+      PrepareDensityField(LevelArray, SiblingList, level, MetaData, When);
 #else   // !FAST_SIB
-    if (SelfGravity)
-      if (PrepareDensityField(LevelArray, level, MetaData, When) == FAIL) {
-        fprintf(stderr, "Error in PrepareDensityField.\n");
-        ENZO_FAIL("");
-      }
+      PrepareDensityField(LevelArray, level, MetaData, When);
 #endif  // end FAST_SIB
+    }
 
     /* Solve the radiative transfer */
 
 #ifdef TRANSFER
-    //Grids[grid1]->GridData->SetTimePreviousTimestep();
     FLOAT GridTime = Grids[0]->GridData->ReturnTime();
     EvolvePhotons(MetaData, LevelArray, AllStars, GridTime, level);
-    //Grids[grid1]->GridData->SetTimeNextTimestep();
 #endif /* TRANSFER */
 
     /* Compute particle-particle acceleration */
 
-    //if (NBodyDirectSummation == TRUE) 
-    //for (grid = 0; grid < NumberOfGrids; grid++)
-    //Grids[grid1]->GridData->ComputeParticleParticleAcceleration(level);
+    if (NBodyDirectSummation == TRUE) 
+      for (grid1 = 0; grid1 < NumberOfGrids; grid1++)
+	Grids[grid1]->GridData->ComputeParticleParticleAcceleration(level);
 
     /* ------------------------------------------------------- */
-    /* Evolve all grids by timestep dtThisLevel. */
+    /* Evolve all grids by timestep dtThisLevel. Predictor Step*/
 
     for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
+      
+      CallProblemSpecificRoutines(MetaData, Grids[grid1], grid1, &norm, 
+				  TopGridTimeStep, level, LevelCycleCount);
 
       /* Gravity: compute acceleration field for grid and particles. */
-
       if (SelfGravity) {
 	int Dummy;
 	if (level <= MaximumGravityRefinementLevel) {
 	  if (level > 0) 
 	    Grids[grid1]->GridData->SolveForPotential(level) ;
-
 	  Grids[grid1]->GridData->ComputeAccelerations(level) ;
 	}
-	// otherwise, interpolate potential from coarser grid, which is
-	//   now done in PrepareDensity.
+	// otherwise, use interpolated potential from coarser grid, which is
+	// done in PrepareDensity.
       } // end: if (SelfGravity)
 
       Grids[grid1]->GridData->ComputeAccelerationFieldExternal() ;
 
-
 #ifdef TRANSFER
-
       /* Radiation Pressure: add to acceleration field */
-
-      if (RadiativeTransfer && RadiationPressure)
-	if (Grids[grid1]->GridData->AddRadiationPressureAcceleration() == FAIL) {
-	  fprintf(stderr,"Error in grid->AddRadiationPressureAcceleration.\n");
-	  ENZO_FAIL("");
-	}
-
+      Grids[grid1]->GridData->AddRadiationPressureAcceleration();
 #endif /* TRANSFER */
 
-
       Grids[grid1]->GridData->CopyBaryonFieldToOldBaryonField();
 
       if (UseHydro) {
-
 	if (HydroMethod == HD_RK)
 	  Grids[grid1]->GridData->RungeKutta2_1stStep
-	    (SubgridFluxesEstimate[grid1], 
-	     NumberOfSubgrids[grid1], level, Exterior);
-
+	    (SubgridFluxesEstimate[grid1], NumberOfSubgrids[grid1], level, Exterior);
 	else if (HydroMethod == MHD_RK) {
 	  Grids[grid1]->GridData->MHDRK2_1stStep
-	    (SubgridFluxesEstimate[grid1], 
-	     NumberOfSubgrids[grid1], level, Exterior);
+	    (SubgridFluxesEstimate[grid1], NumberOfSubgrids[grid1], level, Exterior);
 	}
       } // ENDIF UseHydro
 	
@@ -461,34 +439,30 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
       ENZO_FAIL("");
       }*/
 #ifdef FAST_SIB
-  SetBoundaryConditions(Grids, NumberOfGrids, SiblingList,
-			level, MetaData, Exterior, LevelArray[level]);
+    SetBoundaryConditions(Grids, NumberOfGrids, SiblingList, level, MetaData, Exterior, LevelArray[level]);
 #else
-  SetBoundaryConditions(Grids, NumberOfGrids, level, MetaData,
-			Exterior, LevelArray[level]);
+    SetBoundaryConditions(Grids, NumberOfGrids, level, MetaData, Exterior, LevelArray[level]);
 #endif
 
+
+
     for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
 
       if (UseHydro) {
 	if (HydroMethod == HD_RK)
 	  Grids[grid1]->GridData->RungeKutta2_2ndStep
-	    (SubgridFluxesEstimate[grid1], 
-	     NumberOfSubgrids[grid1], level, Exterior);
+	    (SubgridFluxesEstimate[grid1], NumberOfSubgrids[grid1], level, Exterior);
 
 	else if (HydroMethod == MHD_RK) {
 
 	  Grids[grid1]->GridData->MHDRK2_2ndStep
-	    (SubgridFluxesEstimate[grid1], 
-	     NumberOfSubgrids[grid1], level, Exterior);
+	    (SubgridFluxesEstimate[grid1], NumberOfSubgrids[grid1], level, Exterior);
 	  
-	  if (UseAmbipolarDiffusion) {
+	  if (UseAmbipolarDiffusion) 
 	    Grids[grid1]->GridData->AddAmbipolarDiffusion();
-	  }
-
-	  if (UseResistivity) {
+	  
+	  if (UseResistivity) 
 	    Grids[grid1]->GridData->AddResistivity();
-	  }
 	
 	 
 	  time1 = ReturnWallTime();
@@ -552,6 +526,32 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
     StarParticleFinalize(Grids, MetaData, NumberOfGrids, LevelArray,
 			 level, AllStars);
 
+    /* Reompute potential if output is requested */
+    if (SelfGravity && WritePotential) {
+      CopyGravPotential = TRUE;
+      When = 0.0;
+ 
+#ifdef FAST_SIB
+      PrepareDensityField(LevelArray, SiblingList, level, MetaData, When);
+#else   // !FAST_SIB
+      PrepareDensityField(LevelArray, level, MetaData, When);
+#endif  // end FAST_SIB
+ 
+      //      CopyGravPotential = FALSE;
+ 
+      for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
+        if (level <= MaximumGravityRefinementLevel) {
+ 
+          /* Compute the potential. */
+ 
+          if (level > 0)
+            Grids[grid1]->GridData->SolveForPotential(level);
+          Grids[grid1]->GridData->CopyPotentialToBaryonField();
+        }
+      } //  end loop over grids
+    } // if WritePotential
+
+
     OutputFromEvolveLevel(LevelArray,MetaData,level,Exterior);
     CallPython(LevelArray, MetaData, level);
 
diff --git a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C
index c7644846b..24b6bbb97 100644
--- a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C
+++ b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C
@@ -231,6 +231,11 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL
 
 	Velx = Vely = Velz = 0.0;
 
+	// default values:
+	Density = CloudDensity;
+	eint = CloudInternalEnergy;
+
+
 	if (r < CloudRadius) {
 
 	  /* Type 0: uniform cloud */
@@ -294,8 +299,8 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL
 	} else {
 
 	  if (CloudType == 0) {
-	    Density = CloudDensity/10.0;
-	    eint = CloudInternalEnergy*10.;
+	    Density = CloudDensity/100.0;
+	    eint = CloudInternalEnergy*100.;
 	  }
 
 	  if (CloudType == 1) {
@@ -554,9 +559,11 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL
 
   /* Put a sink particle if we are studying massive star formation */
 
-  if (PutSink == 1 && level == MaximumRefinementLevel) {
+  //  if (PutSink == 1 && level == MaximumRefinementLevel) {
+  if (PutSink == 1 && level == 0) {  // set it up on level zero and make it mustrefine
 
-    double mass_p = 20.0*1.989e33;
+    //    double mass_p = 20.0*1.989e33;
+    double mass_p = 0.01*1.989e33;
     mass_p /= MassUnits;
     double dx = CellWidth[0][0];
     double den_p = mass_p / pow(dx,3);
@@ -576,9 +583,9 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL
     ParticleMass[0] = den_p;
     ParticleNumber[0] = 0;
     ParticleType[0] = PARTICLE_TYPE_MUST_REFINE;
-    ParticlePosition[0][0] = 0.5+0.5*dx;
-    ParticlePosition[1][0] = 0.5+0.5*dx;
-    ParticlePosition[2][0] = 0.5+0.5*dx;
+    ParticlePosition[0][0] = 0.5;//+0.5*dx;
+    ParticlePosition[1][0] = 0.5;//+0.5*dx;
+    ParticlePosition[2][0] = 0.5;//+0.5*dx;
 
     ParticleVelocity[0][0] = 0.0;
     ParticleVelocity[1][0] = 0.0;

From 86567d82755f00b5a35afb5b922c0e5219a70346 Mon Sep 17 00:00:00 2001
From: tabel <none@none>
Date: Sat, 10 Oct 2009 15:17:21 -0700
Subject: [PATCH 06/12] Put SetBoundaryConditions in a separate file.

--HG--
branch : week-of-code
---
 src/enzo/EvolveLevelRoutinesOptimized.C | 198 ----------------------
 src/enzo/Make.config.objects            |   1 +
 src/enzo/PrepareDensityField.C          |  28 ++-
 src/enzo/SetBoundaryConditions.C        | 215 ++++++++++++++++++++++++
 4 files changed, 239 insertions(+), 203 deletions(-)
 create mode 100644 src/enzo/SetBoundaryConditions.C

diff --git a/src/enzo/EvolveLevelRoutinesOptimized.C b/src/enzo/EvolveLevelRoutinesOptimized.C
index b8a95bc72..b29c5a7e5 100644
--- a/src/enzo/EvolveLevelRoutinesOptimized.C
+++ b/src/enzo/EvolveLevelRoutinesOptimized.C
@@ -38,211 +38,13 @@
 
 /* function prototypes */
  
-int DepositParticleMassField(HierarchyEntry *Grid, FLOAT Time = -1.0);
-
-int CommunicationBufferPurge(void);
 int CommunicationReceiveHandler(fluxes **SubgridFluxesEstimate[] = NULL,
 				int NumberOfSubgrids[] = NULL,
 				int FluxFlag = FALSE,
 				TopGridData* MetaData = NULL);
- 
-int PrepareGravitatingMassField1(HierarchyEntry *Grid);
-#ifdef FAST_SIB
-int PrepareGravitatingMassField2(HierarchyEntry *Grid, int grid1,
-				 SiblingGridList SiblingList[],
-				 TopGridData *MetaData, int level,
-				 FLOAT When);
-#else
-int PrepareGravitatingMassField2(HierarchyEntry *Grid, TopGridData *MetaData,
-				 LevelHierarchyEntry *LevelArray[], int level,
-				 FLOAT When);
-#endif
- 
-#ifdef FAST_SIB
-int ComputePotentialFieldLevelZero(TopGridData *MetaData,
-				   SiblingGridList SiblingList[],
-				   HierarchyEntry *Grids[], int NumberOfGrids);
-#else
-int ComputePotentialFieldLevelZero(TopGridData *MetaData,
-				   HierarchyEntry *Grids[], int NumberOfGrids);
-#endif
-
-int GenerateGridArray(LevelHierarchyEntry *LevelArray[], int level,
-		      HierarchyEntry **Grids[]);
- 
- 
- 
 
 #define GRIDS_PER_LOOP 20000
  
-/* ======================================================================= */
-/* This routine sets all the boundary conditions for Grids by either
-   interpolating from their parents or copying from sibling grids. */
- 
-
-
-#ifdef FAST_SIB
-int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids,
-			  SiblingGridList SiblingList[],
-			  int level, TopGridData *MetaData,
-			  ExternalBoundary *Exterior, LevelHierarchyEntry *Level)
-#else
-int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids,
-			  int level, TopGridData *MetaData,
-			  ExternalBoundary *Exterior, LevelHierarchyEntry *Level)
-#endif
-{
- 
-  int loopEnd = (ShearingBoundaryDirection != 1) ? 2 : 1;
-  
-  int grid1, grid2, StartGrid, EndGrid, loop;
-  
-  JBPERF_START("SetBoundaryConditions");
-    
-  for (loop = 0; loop < loopEnd; loop++){
-    
-#ifdef FORCE_MSG_PROGRESS 
-    CommunicationBarrier();
-#endif
-  
-    TIME_MSG("Interpolating boundaries from parent");
-      
-    for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) {
-	
-      if (traceMPI) fprintf(tracePtr, "SBC loop\n");
-	
-      EndGrid = min(StartGrid + GRIDS_PER_LOOP, NumberOfGrids);
-	
-      /* -------------- FIRST PASS ----------------- */
-      /* Here, we just generate the calls to generate the receive buffers,
-	 without actually doing anything. */
-      
-      CommunicationDirection = COMMUNICATION_POST_RECEIVE;
-      CommunicationReceiveIndex = 0;
-      
-      for (grid1 = StartGrid; grid1 < EndGrid; grid1++) {
-	
-	/* a) Interpolate boundaries from the parent grid or set external
-	   boundary conditions. */
-	
-	CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE;
-	
-	if (level == 0) {
-	  if (loop == 0) 
-	    Grids[grid1]->GridData->SetExternalBoundaryValues(Exterior);
-	} else {
-	  Grids[grid1]->GridData->InterpolateBoundaryFromParent
-	    (Grids[grid1]->ParentGrid->GridData);
-	}
-
-      } // ENDFOR grids
-
-	/* -------------- SECOND PASS ----------------- */
-	/* Now we generate all the sends, and do all the computation
-	   for grids which are on the same processor as well. */
-
-      CommunicationDirection = COMMUNICATION_SEND;
-      for (grid1 = StartGrid; grid1 < EndGrid; grid1++) {
-
-	/* a) Interpolate boundaries from the parent grid or set
-	   external boundary conditions. */
-
-	if (level > 0)
-	  Grids[grid1]->GridData->InterpolateBoundaryFromParent
-	    (Grids[grid1]->ParentGrid->GridData);
-     
-      
-      } // ENDFOR grids
-
-	//Grids[StartGrid]->GridData->PrintToScreenBoundaries(0);
-
-	/* -------------- THIRD PASS ----------------- */
-	/* In this final step, we get the messages as they come in and
-	   then match them to the methods which generate the receive
-	   handle. */
-
-      if (CommunicationReceiveHandler() == FAIL)
-	ENZO_FAIL("");
-
-    } // ENDFOR grid batches
-
-    TIME_MSG("Copying zones in SetBoundaryConditions");
-    for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) {
-      EndGrid = min(StartGrid + GRIDS_PER_LOOP, NumberOfGrids);
-
-      /* -------------- FIRST PASS ----------------- */
-      /* b) Copy any overlapping zones for sibling grids.  */
-
-      CommunicationDirection = COMMUNICATION_POST_RECEIVE;
-      CommunicationReceiveIndex = 0;
- 
-#ifdef FAST_SIB
-      for (grid1 = StartGrid; grid1 < EndGrid; grid1++)
-	for (grid2 = 0; grid2 < SiblingList[grid1].NumberOfSiblings; grid2++)
-	  Grids[grid1]->GridData->
-	    CheckForOverlap(SiblingList[grid1].GridList[grid2],
-			    MetaData->LeftFaceBoundaryCondition,
-			    MetaData->RightFaceBoundaryCondition,
-			    &grid::CopyZonesFromGrid);
-#else
-      for (grid1 = StartGrid; grid1 < EndGrid; grid1++)
-	for (grid2 = 0; grid2 < NumberOfGrids; grid2++)
-	  Grids[grid1]->GridData->
-	    CheckForOverlap(Grids[grid2]->GridData,
-			    MetaData->LeftFaceBoundaryCondition,
-			    MetaData->RightFaceBoundaryCondition,
-			    &grid::CopyZonesFromGrid);
-#endif
-
-      /* -------------- SECOND PASS ----------------- */
-      /* b) Copy any overlapping zones for sibling grids.  */
-
-      CommunicationDirection = COMMUNICATION_SEND;
-#ifdef FAST_SIB
-      for (grid1 = StartGrid; grid1 < EndGrid; grid1++)
-	for (grid2 = 0; grid2 < SiblingList[grid1].NumberOfSiblings; grid2++)
-	  Grids[grid1]->GridData->
-	    CheckForOverlap(SiblingList[grid1].GridList[grid2],
-			    MetaData->LeftFaceBoundaryCondition,
-			    MetaData->RightFaceBoundaryCondition,
-			    &grid::CopyZonesFromGrid);
-#else
-      for (grid1 = StartGrid; grid1 < EndGrid; grid1++)
-	for (grid2 = 0; grid2 < NumberOfGrids; grid2++)
-	  Grids[grid1]->GridData->
-	    CheckForOverlap(Grids[grid2]->GridData,
-			    MetaData->LeftFaceBoundaryCondition,
-			    MetaData->RightFaceBoundaryCondition,
-			    &grid::CopyZonesFromGrid);
-#endif
-
-      /* -------------- THIRD PASS ----------------- */
-
-      if (CommunicationReceiveHandler() == FAIL)
-	ENZO_FAIL("");
-
-    } // end loop over batchs of grids
-
-  } // ENDFOR loop (for ShearinBox)
- 
-    /* c) Apply external reflecting boundary conditions, if needed.  */
-
-  for (grid1 = 0; grid1 < NumberOfGrids; grid1++)
-    Grids[grid1]->GridData->CheckForExternalReflections
-      (MetaData->LeftFaceBoundaryCondition,
-       MetaData->RightFaceBoundaryCondition);
-  
-#ifdef FORCE_MSG_PROGRESS 
-  CommunicationBarrier();
-#endif
- 
-  CommunicationDirection = COMMUNICATION_SEND_RECEIVE;
- 
-  JBPERF_STOP("SetBoundaryConditions");
-
-  return SUCCESS;
-  
-}
  
   /* ======================================================================= */
   /* This routines does the flux correction and project for all grids on this
diff --git a/src/enzo/Make.config.objects b/src/enzo/Make.config.objects
index bbf37bdd7..fc06d42da 100644
--- a/src/enzo/Make.config.objects
+++ b/src/enzo/Make.config.objects
@@ -544,6 +544,7 @@ OBJS_CONFIG_LIB = \
         s90_st1.o \
         SedovBlastInitialize.o \
         select_fft.o \
+	SetBoundaryConditions.o \
         SetDefaultGlobalValues.o \
         SetLevelTimeStep.o \
         sgi_st1_fft64.o \
diff --git a/src/enzo/PrepareDensityField.C b/src/enzo/PrepareDensityField.C
index 3bea2b751..69c0e2258 100644
--- a/src/enzo/PrepareDensityField.C
+++ b/src/enzo/PrepareDensityField.C
@@ -1,3 +1,26 @@
+/***********************************************************************
+/
+/  PREPARE DENSITY FIELD (CALLED BY EVOLVE LEVEL)
+/
+/  written by: Greg Bryan
+/  date:       June, 1999
+/  modifiedN:  Robert Harkness
+/  date:       February, 2008
+/
+/ ======================================================================= 
+/ This routine prepares the density field for all the grids on this level,
+/ both particle and baryonic densities.  It also calculates the potential
+/ field if this is level 0 (since this involves communication). 
+/
+/   This is part of a collection of routines called by EvolveLevel.
+/   These have been optimized for enhanced message passing
+/   performance by performing two passes -- one which generates
+/   sends and the second which receives them.
+/
+/  modified: Robert Harkness, December 2007
+/
+************************************************************************/
+
 #ifdef USE_MPI
 #include <mpi.h>
 #endif /* USE_MPI */
@@ -59,11 +82,6 @@ extern int CopyPotentialFieldAverage;
 #define GRIDS_PER_LOOP 20000
 
  
-  /* ======================================================================= */
-  /* This routine prepares the density field for all the grids on this level,
-     both particle and baryonic densities.  It also calculates the potential
-     field if this is level 0 (since this involves communication). */
- 
 #ifdef FAST_SIB
 int PrepareDensityField(LevelHierarchyEntry *LevelArray[],
 			SiblingGridList SiblingList[],
diff --git a/src/enzo/SetBoundaryConditions.C b/src/enzo/SetBoundaryConditions.C
new file mode 100644
index 000000000..f34f04e25
--- /dev/null
+++ b/src/enzo/SetBoundaryConditions.C
@@ -0,0 +1,215 @@
+/***********************************************************************
+/
+/  SET BOUNDARY CONDITIONS (CALLED BY EVOLVE LEVEL)
+/
+/  written by: Greg Bryan
+/  date:       June, 1999
+/  modifiedN:  Robert Harkness
+/  date:       February, 2008
+/
+/  PURPOSE:   ======================================================================= 
+/  This routine sets all the boundary conditions for Grids by either
+/   interpolating from their parents or copying from sibling grids. 
+/
+/   This is part of a collection of routines called by EvolveLevel.
+/   These have been optimized for enhanced message passing
+/   performance by performing two passes -- one which generates
+/   sends and the second which receives them.
+/
+/  modified: Robert Harkness, December 2007
+/
+************************************************************************/
+ 
+#ifdef USE_MPI
+#include <mpi.h>
+#endif /* USE_MPI */
+ 
+#include <stdio.h>
+#include "ErrorExceptions.h"
+#include "performance.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 "TopGridData.h"
+#include "LevelHierarchy.h"
+#include "communication.h"
+#include "CommunicationUtilities.h"
+
+/* function prototypes */
+ 
+int CommunicationReceiveHandler(fluxes **SubgridFluxesEstimate[] = NULL,
+				int NumberOfSubgrids[] = NULL,
+				int FluxFlag = FALSE,
+				TopGridData* MetaData = NULL);
+
+#define GRIDS_PER_LOOP 20000
+ 
+
+
+#ifdef FAST_SIB
+int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids,
+			  SiblingGridList SiblingList[],
+			  int level, TopGridData *MetaData,
+			  ExternalBoundary *Exterior, LevelHierarchyEntry *Level)
+#else
+int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids,
+			  int level, TopGridData *MetaData,
+			  ExternalBoundary *Exterior, LevelHierarchyEntry *Level)
+#endif
+{
+ 
+  int loopEnd = (ShearingBoundaryDirection != 1) ? 2 : 1;
+  
+  int grid1, grid2, StartGrid, EndGrid, loop;
+  
+  JBPERF_START("SetBoundaryConditions");
+    
+  for (loop = 0; loop < loopEnd; loop++){
+    
+#ifdef FORCE_MSG_PROGRESS 
+    CommunicationBarrier();
+#endif
+  
+    TIME_MSG("Interpolating boundaries from parent");
+      
+    for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) {
+	
+      if (traceMPI) fprintf(tracePtr, "SBC loop\n");
+	
+      EndGrid = min(StartGrid + GRIDS_PER_LOOP, NumberOfGrids);
+	
+      /* -------------- FIRST PASS ----------------- */
+      /* Here, we just generate the calls to generate the receive buffers,
+	 without actually doing anything. */
+      
+      CommunicationDirection = COMMUNICATION_POST_RECEIVE;
+      CommunicationReceiveIndex = 0;
+      
+      for (grid1 = StartGrid; grid1 < EndGrid; grid1++) {
+	
+	/* a) Interpolate boundaries from the parent grid or set external
+	   boundary conditions. */
+	
+	CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE;
+	
+	if (level == 0) {
+	  if (loop == 0) 
+	    Grids[grid1]->GridData->SetExternalBoundaryValues(Exterior);
+	} else {
+	  Grids[grid1]->GridData->InterpolateBoundaryFromParent
+	    (Grids[grid1]->ParentGrid->GridData);
+	}
+
+      } // ENDFOR grids
+
+	/* -------------- SECOND PASS ----------------- */
+	/* Now we generate all the sends, and do all the computation
+	   for grids which are on the same processor as well. */
+
+      CommunicationDirection = COMMUNICATION_SEND;
+      for (grid1 = StartGrid; grid1 < EndGrid; grid1++) {
+
+	/* a) Interpolate boundaries from the parent grid or set
+	   external boundary conditions. */
+
+	if (level > 0)
+	  Grids[grid1]->GridData->InterpolateBoundaryFromParent
+	    (Grids[grid1]->ParentGrid->GridData);
+     
+      
+      } // ENDFOR grids
+
+	//Grids[StartGrid]->GridData->PrintToScreenBoundaries(0);
+
+	/* -------------- THIRD PASS ----------------- */
+	/* In this final step, we get the messages as they come in and
+	   then match them to the methods which generate the receive
+	   handle. */
+
+      if (CommunicationReceiveHandler() == FAIL)
+	ENZO_FAIL("");
+
+    } // ENDFOR grid batches
+
+    TIME_MSG("Copying zones in SetBoundaryConditions");
+    for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) {
+      EndGrid = min(StartGrid + GRIDS_PER_LOOP, NumberOfGrids);
+
+      /* -------------- FIRST PASS ----------------- */
+      /* b) Copy any overlapping zones for sibling grids.  */
+
+      CommunicationDirection = COMMUNICATION_POST_RECEIVE;
+      CommunicationReceiveIndex = 0;
+ 
+#ifdef FAST_SIB
+      for (grid1 = StartGrid; grid1 < EndGrid; grid1++)
+	for (grid2 = 0; grid2 < SiblingList[grid1].NumberOfSiblings; grid2++)
+	  Grids[grid1]->GridData->
+	    CheckForOverlap(SiblingList[grid1].GridList[grid2],
+			    MetaData->LeftFaceBoundaryCondition,
+			    MetaData->RightFaceBoundaryCondition,
+			    &grid::CopyZonesFromGrid);
+#else
+      for (grid1 = StartGrid; grid1 < EndGrid; grid1++)
+	for (grid2 = 0; grid2 < NumberOfGrids; grid2++)
+	  Grids[grid1]->GridData->
+	    CheckForOverlap(Grids[grid2]->GridData,
+			    MetaData->LeftFaceBoundaryCondition,
+			    MetaData->RightFaceBoundaryCondition,
+			    &grid::CopyZonesFromGrid);
+#endif
+
+      /* -------------- SECOND PASS ----------------- */
+      /* b) Copy any overlapping zones for sibling grids.  */
+
+      CommunicationDirection = COMMUNICATION_SEND;
+#ifdef FAST_SIB
+      for (grid1 = StartGrid; grid1 < EndGrid; grid1++)
+	for (grid2 = 0; grid2 < SiblingList[grid1].NumberOfSiblings; grid2++)
+	  Grids[grid1]->GridData->
+	    CheckForOverlap(SiblingList[grid1].GridList[grid2],
+			    MetaData->LeftFaceBoundaryCondition,
+			    MetaData->RightFaceBoundaryCondition,
+			    &grid::CopyZonesFromGrid);
+#else
+      for (grid1 = StartGrid; grid1 < EndGrid; grid1++)
+	for (grid2 = 0; grid2 < NumberOfGrids; grid2++)
+	  Grids[grid1]->GridData->
+	    CheckForOverlap(Grids[grid2]->GridData,
+			    MetaData->LeftFaceBoundaryCondition,
+			    MetaData->RightFaceBoundaryCondition,
+			    &grid::CopyZonesFromGrid);
+#endif
+
+      /* -------------- THIRD PASS ----------------- */
+
+      if (CommunicationReceiveHandler() == FAIL)
+	ENZO_FAIL("");
+
+    } // end loop over batchs of grids
+
+  } // ENDFOR loop (for ShearinBox)
+ 
+    /* c) Apply external reflecting boundary conditions, if needed.  */
+
+  for (grid1 = 0; grid1 < NumberOfGrids; grid1++)
+    Grids[grid1]->GridData->CheckForExternalReflections
+      (MetaData->LeftFaceBoundaryCondition,
+       MetaData->RightFaceBoundaryCondition);
+  
+#ifdef FORCE_MSG_PROGRESS 
+  CommunicationBarrier();
+#endif
+ 
+  CommunicationDirection = COMMUNICATION_SEND_RECEIVE;
+ 
+  JBPERF_STOP("SetBoundaryConditions");
+
+  return SUCCESS;
+  
+}

From ba57ee6786c2ba005ea136714bd68f66cf663922 Mon Sep 17 00:00:00 2001
From: tabel <none@none>
Date: Sat, 10 Oct 2009 15:29:31 -0700
Subject: [PATCH 07/12] Split also GenerateGridArray and UpdateFromFinerGrids
 into separate file. EvolveLevelRoutinesOptimized is now obsolete and has been
 removed after all its contents were split into separate files.

--HG--
branch : week-of-code
---
 src/enzo/GenerateGridArray.C                  | 49 +++++++++++++++++
 src/enzo/Make.config.objects                  |  2 +
 ...inesOptimized.C => UpdateFromFinerGrids.C} | 52 ++++---------------
 3 files changed, 61 insertions(+), 42 deletions(-)
 create mode 100644 src/enzo/GenerateGridArray.C
 rename src/enzo/{EvolveLevelRoutinesOptimized.C => UpdateFromFinerGrids.C} (87%)

diff --git a/src/enzo/GenerateGridArray.C b/src/enzo/GenerateGridArray.C
new file mode 100644
index 000000000..12574e123
--- /dev/null
+++ b/src/enzo/GenerateGridArray.C
@@ -0,0 +1,49 @@
+/***********************************************************************
+/
+/  EVOLVE LEVEL ROUTINES (CALLED BY EVOLVE LEVEL)
+/
+/  written by: Greg Bryan
+/  date:       June, 1999
+/
+/  ======================================================================= 
+/  This routine simply converts a linked list of grids into an array of
+/   pointers. 
+/
+************************************************************************/ 
+
+#include <stdio.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"
+
+int GenerateGridArray(LevelHierarchyEntry *LevelArray[], int level,
+		      HierarchyEntry **Grids[])
+{
+  /* Count the number of grids on this level. */
+ 
+  int NumberOfGrids = 0, counter = 0;
+  LevelHierarchyEntry *Temp = LevelArray[level];
+  while (Temp != NULL) {
+    NumberOfGrids++;
+    Temp             = Temp->NextGridThisLevel;
+  }
+  /* Create a list of pointers and number of subgrids (and fill it out). */
+ 
+  typedef HierarchyEntry* HierarchyEntryPointer;
+  *Grids = new HierarchyEntryPointer[NumberOfGrids];
+  Temp = LevelArray[level];
+  while (Temp != NULL) {
+    (*Grids)[counter++] = Temp->GridHierarchyEntry;
+    Temp              = Temp->NextGridThisLevel;
+  }
+ 
+  return NumberOfGrids;
+}
+ 
+  
diff --git a/src/enzo/Make.config.objects b/src/enzo/Make.config.objects
index fc06d42da..802e5c8c9 100644
--- a/src/enzo/Make.config.objects
+++ b/src/enzo/Make.config.objects
@@ -182,6 +182,7 @@ OBJS_CONFIG_LIB = \
 	GalaxySimulationInitialize.o \
         F_ColumnFormat.o \
         f_message.o \
+	GenerateGridArray.o \
         GetNodeFreeMemory.o \
         GetUnits.o \
         GravityBdryExchange.o \
@@ -606,6 +607,7 @@ OBJS_CONFIG_LIB = \
         tscint3d.o \
         TurbulenceSimulationInitialize.o \
         twoshock.o \
+	UpdateFromFinerGrids.o \
         UpdateParticlePositions.o \
         utilities.o \
         wall_clock.o \
diff --git a/src/enzo/EvolveLevelRoutinesOptimized.C b/src/enzo/UpdateFromFinerGrids.C
similarity index 87%
rename from src/enzo/EvolveLevelRoutinesOptimized.C
rename to src/enzo/UpdateFromFinerGrids.C
index b29c5a7e5..714b0bffa 100644
--- a/src/enzo/EvolveLevelRoutinesOptimized.C
+++ b/src/enzo/UpdateFromFinerGrids.C
@@ -1,16 +1,21 @@
 /***********************************************************************
 /
-/  EVOLVE LEVEL ROUTINES (CALLED BY EVOLVE LEVEL)
+/  UPDATE FROM FINER GRIDS (CALLED BY EVOLVE LEVEL)
 /
 /  written by: Greg Bryan
 /  date:       June, 1999
 /  modifiedN:  Robert Harkness
 /  date:       February, 2008
 /
-/  PURPOSE:  This is a collection of routines called by EvolveLevel.
-/            These have been optimized for enhanced message passing
-/            performance by performing two passes -- one which generates
-/            sends and the second which receives them.
+/  PURPOSE:  
+/  ======================================================================= 
+/  This routines does the flux correction and project for all grids on this
+/  level from the list of subgrids. 
+/
+/  This is part of a collection of routines called by EvolveLevel.
+/  These have been optimized for enhanced message passing
+/  performance by performing two passes -- one which generates
+/  sends and the second which receives them.
 /
 /  modified: Robert Harkness, December 2007
 /
@@ -46,10 +51,6 @@ int CommunicationReceiveHandler(fluxes **SubgridFluxesEstimate[] = NULL,
 #define GRIDS_PER_LOOP 20000
  
  
-  /* ======================================================================= */
-  /* This routines does the flux correction and project for all grids on this
-     level from the list of subgrids. */
- 
 #ifdef FLUX_FIX
 int UpdateFromFinerGrids(int level, HierarchyEntry *Grids[], int NumberOfGrids,
 			 int NumberOfSubgrids[],
@@ -325,36 +326,3 @@ int UpdateFromFinerGrids(int level, HierarchyEntry *Grids[], int NumberOfGrids,
   return SUCCESS;
 }
  
- 
- 
-/* ======================================================================= */
-/* This routine simply converts a linked list of grids into an array of
-   pointers. */
- 
-int GenerateGridArray(LevelHierarchyEntry *LevelArray[], int level,
-		      HierarchyEntry **Grids[])
-{
- 
-  /* Count the number of grids on this level. */
- 
-  int NumberOfGrids = 0, counter = 0;
-  LevelHierarchyEntry *Temp = LevelArray[level];
-  while (Temp != NULL) {
-    NumberOfGrids++;
-    Temp             = Temp->NextGridThisLevel;
-  }
- 
-  /* Create a list of pointers and number of subgrids (and fill it out). */
- 
-  typedef HierarchyEntry* HierarchyEntryPointer;
-  *Grids = new HierarchyEntryPointer[NumberOfGrids];
-  Temp = LevelArray[level];
-  while (Temp != NULL) {
-    (*Grids)[counter++] = Temp->GridHierarchyEntry;
-    Temp              = Temp->NextGridThisLevel;
-  }
- 
-  return NumberOfGrids;
-}
- 
-  

From 34eb43ad19c46b10d2e27c2bfbed6376e8970ac7 Mon Sep 17 00:00:00 2001
From: tabel <none@none>
Date: Sat, 10 Oct 2009 19:10:11 -0700
Subject: [PATCH 08/12] In EvolveLevel_RK2 there is now an option to use the
 time averaged density field for the baryon gravity. This should be more
 accurate but also will cost to calculate the potential and acceleration twice
 and will run consequently slower. In simple collapse problems this does not
 seem worth it. Differences in collapse time are approximately of order a few
 percent of the dynamical time. However, I'll leave this in for particular
 cases when one might worry about even better accuracy. To turn it on one has
 to change (currently line 457: RK2SecondStepBaryonDeposit = 0 to
 RK2SecondStepBaryonDeposit = 1

--HG--
branch : week-of-code
---
 src/enzo/DepositBaryons.C                     |  2 -
 src/enzo/EvolveLevel.C                        | 10 +--
 src/enzo/Grid_AccelerationBoundaryRoutines.C  | 17 ++--
 .../Grid_CopyBaryonFieldToOldBaryonField.C    |  5 --
 src/enzo/Grid_DepositBaryons.C                | 60 ++++++++++----
 src/enzo/Make.config.objects                  |  1 -
 src/enzo/hydro_rk/EvolveLevel_RK2.C           | 79 ++++++++++---------
 src/enzo/hydro_rk/Grid_MHDRK2_2ndStep.C       |  7 --
 8 files changed, 98 insertions(+), 83 deletions(-)

diff --git a/src/enzo/DepositBaryons.C b/src/enzo/DepositBaryons.C
index 5044612ff..d44fcee01 100644
--- a/src/enzo/DepositBaryons.C
+++ b/src/enzo/DepositBaryons.C
@@ -77,8 +77,6 @@ int DepositBaryons(HierarchyEntry *Grid, FLOAT When)
 }
  
  
-/* this doesn't quite work yet (subgrid zeroing needs to be worked out). */
- 
 int DepositBaryonsChildren(HierarchyEntry *DepositGrid,
 			   HierarchyEntry *Grid, FLOAT DepositTime)
 {
diff --git a/src/enzo/EvolveLevel.C b/src/enzo/EvolveLevel.C
index 4691345a3..9f0999e54 100644
--- a/src/enzo/EvolveLevel.C
+++ b/src/enzo/EvolveLevel.C
@@ -219,6 +219,8 @@ static int StaticLevelZero = 1;
 static int StaticLevelZero = 0;
 #endif
 
+extern int RK2SecondStepBaryonDeposit;
+
 int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 		int level, float dtLevelAbove, ExternalBoundary *Exterior)
 {
@@ -253,7 +255,6 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
   LevelHierarchyEntry **SUBlingList;
 #endif
 
-
   /* Initialize the chaining mesh used in the FastSiblingLocator. */
 
   if (dbx) fprintf(stderr, "EL: Initialize FSL \n"); 
@@ -493,13 +494,6 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
       Grids[grid1]->GridData->DeleteGravitatingMassFieldParticles();
 
 
-    /* Run the Divergence Cleaing                */
-
-    for (grid1 = 0; grid1 < NumberOfGrids; grid1++)
-      Grids[grid1]->GridData->PoissonSolver(level);
-    
-
- 
     /* ----------------------------------------- */
     /* Evolve the next level down (recursively). */
  
diff --git a/src/enzo/Grid_AccelerationBoundaryRoutines.C b/src/enzo/Grid_AccelerationBoundaryRoutines.C
index 2bd337e27..434a8fa17 100644
--- a/src/enzo/Grid_AccelerationBoundaryRoutines.C
+++ b/src/enzo/Grid_AccelerationBoundaryRoutines.C
@@ -122,29 +122,28 @@ int SetAccelerationBoundary(HierarchyEntry *Grids[], int NumberOfGrids,
 			    int CycleNumber)
 {
 
-  if ( ! (SelfGravity || UniformGravity || PointSourceGravity) && level > 0 )
+  if ( ! ( (SelfGravity || UniformGravity || PointSourceGravity) && level > 0 ))
     return SUCCESS;
 
   //Set the boundary on the Acceleration field.  Reuse SetBoundaryConditions.  
   //Juggle pointers around.
 
-  int grid, ConservativeTruth;
+  int grid1, ConservativeTruth;
   char basename[30];  
 
   //We don't want conservative interpolation actually being done for the acceleration field.
   ConservativeTruth = ConservativeInterpolation;
   ConservativeInterpolation = FALSE;
 
-  for (grid = 0; grid < NumberOfGrids; grid++) {
-    if( Grids[grid]->GridData->AttachAcceleration() == FAIL ) {
+  for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
+    if( Grids[grid1]->GridData->AttachAcceleration() == FAIL ) {
       fprintf(stderr,"Error in AttachAcceleration \n");
       ENZO_FAIL("");
     }
-    if( Grids[grid]->ParentGrid->GridData->AttachAcceleration() ==FAIL ){
+    if( Grids[grid1]->ParentGrid->GridData->AttachAcceleration() ==FAIL ){
       fprintf(stderr,"Error in AttachAcceleration, Parent \n");
       ENZO_FAIL("");
     }
-
   }
 
 #ifdef FAST_SIB
@@ -159,13 +158,13 @@ int SetAccelerationBoundary(HierarchyEntry *Grids[], int NumberOfGrids,
   
 
   
-  for (grid = 0; grid < NumberOfGrids; grid++) {
+  for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
 
-    if( Grids[grid]->GridData->DetachAcceleration() == FAIL ) {
+    if( Grids[grid1]->GridData->DetachAcceleration() == FAIL ) {
       fprintf(stderr,"Error in DetachAcceleration\n");
       ENZO_FAIL("");
     }
-    if( Grids[grid]->ParentGrid->GridData->DetachAcceleration() == FAIL ) {
+    if( Grids[grid1]->ParentGrid->GridData->DetachAcceleration() == FAIL ) {
       fprintf(stderr,"Error in DetachAcceleration, parent\n");
       ENZO_FAIL("");
     }
diff --git a/src/enzo/Grid_CopyBaryonFieldToOldBaryonField.C b/src/enzo/Grid_CopyBaryonFieldToOldBaryonField.C
index 4d167a9c1..5d71ac161 100644
--- a/src/enzo/Grid_CopyBaryonFieldToOldBaryonField.C
+++ b/src/enzo/Grid_CopyBaryonFieldToOldBaryonField.C
@@ -48,11 +48,6 @@ int grid::CopyBaryonFieldToOldBaryonField()
  
   int size = 1;
 
-/* BUG reported 4th June 2006
-  for (int dim = 0; dim < GridDimension[dim]; dim++)
-    size *= GridDimension[dim];
-*/
-
   for (int dim = 0; dim < GridRank; dim++) {
     size *= GridDimension[dim];
   }
diff --git a/src/enzo/Grid_DepositBaryons.C b/src/enzo/Grid_DepositBaryons.C
index 3fd07f4b8..ac227184d 100644
--- a/src/enzo/Grid_DepositBaryons.C
+++ b/src/enzo/Grid_DepositBaryons.C
@@ -51,6 +51,8 @@ extern "C" void FORTRAN_NAME(dep_grid_cic)(
 			       int *ddim1, int *ddim2, int *ddim3,
 			       int *refine1, int *refine2, int *refine3);
  
+int RK2SecondStepBaryonDeposit = 0;
+
 /* InterpolateBoundaryFromParent function */
  
 int grid::DepositBaryons(grid *TargetGrid, FLOAT DepositTime)
@@ -174,8 +176,9 @@ int grid::DepositBaryons(grid *TargetGrid, FLOAT DepositTime)
       ENZO_FAIL("");
     }
     float dt = (DepositTime - Time)/a;
+    //    if (HydroMethod < 3) dt = 0;
     dt = 0;
- 
+
     /* Set up a float version of cell size to pass to fortran. */
  
     float dxfloat[MAX_DIMENSION] = {0,0,0};
@@ -190,21 +193,50 @@ int grid::DepositBaryons(grid *TargetGrid, FLOAT DepositTime)
        velocity field. */
  
 //  fprintf(stderr, "Grid_DepositBaryons - call dep_grid_cic\n");
- 
-    FORTRAN_NAME(dep_grid_cic)(BaryonField[DensNum], dens_field, vel_field,
-			 BaryonField[Vel1Num], BaryonField[Vel2Num],
-			     BaryonField[Vel3Num], &dt,
-			 BaryonField[NumberOfBaryonFields], &GridRank,
-                             &HydroMethod,
-			 dxfloat, dxfloat+1, dxfloat+2,
-			 GridDimension, GridDimension+1, GridDimension+2,
-			 GridStartIndex, GridStartIndex+1, GridStartIndex+2,
-			 GridEndIndex, GridEndIndex+1, GridEndIndex+2,
-			 GridStart, GridStart+1, GridStart+2,
-			 RegionDim, RegionDim+1, RegionDim+2,
-			 Refinement, Refinement+1, Refinement+2);
+
+    float *input_density = BaryonField[DensNum];
+    float *input_velx    = BaryonField[Vel1Num];
+    float *input_vely    = BaryonField[Vel2Num];
+    float *input_velz    = BaryonField[Vel3Num];
+  
+    float *av_dens = NULL;
+    // supply zero velocity field and the time averaged density to the deposit routine for second step
+    if (RK2SecondStepBaryonDeposit == 1) { 
+      int current_size = 1;
+      dt = 0;
+      for (dim = 0; dim < GridRank; dim++)
+	current_size *= GridDimension[dim];
+
+      if (OldBaryonField[DensNum] != NULL) {
+	av_dens = new float[current_size];
+	for (int i=0; i<current_size; i++) 
+	  av_dens[i] = (BaryonField[DensNum][i] + OldBaryonField[DensNum][i]) / 2.;
+      } else 	
+	av_dens = BaryonField[DensNum];
+
+      input_density = av_dens;
+    }
+
+    printf("DepositBaryons, %i\n", RK2SecondStepBaryonDeposit);
+    
+    FORTRAN_NAME(dep_grid_cic)(input_density, dens_field, vel_field,
+			       input_velx, input_vely, input_velz,
+			       &dt,
+			       BaryonField[NumberOfBaryonFields], &GridRank,
+			       &HydroMethod,
+			       dxfloat, dxfloat+1, dxfloat+2,
+			       GridDimension, GridDimension+1, GridDimension+2,
+			       GridStartIndex, GridStartIndex+1, GridStartIndex+2,
+			       GridEndIndex, GridEndIndex+1, GridEndIndex+2,
+			       GridStart, GridStart+1, GridStart+2,
+			       RegionDim, RegionDim+1, RegionDim+2,
+			       Refinement, Refinement+1, Refinement+2);
  
     delete [] vel_field;
+    if ( RK2SecondStepBaryonDeposit ) 
+      if (OldBaryonField[DensNum] != NULL) 
+	delete [] av_dens;
+    
  
   } // end: if (ProcessorNumber == MyProcessorNumber)
  
diff --git a/src/enzo/Make.config.objects b/src/enzo/Make.config.objects
index 802e5c8c9..f24393146 100644
--- a/src/enzo/Make.config.objects
+++ b/src/enzo/Make.config.objects
@@ -135,7 +135,6 @@ OBJS_CONFIG_LIB = \
         FastFourierTransformPrepareComplex.o \
         FastFourierTransformSGIMATH.o \
         EvolveLevel.o \
-        EvolveLevelRoutinesOptimized.o \
         FastSiblingLocatorFinalize.o \
         FastSiblingLocatorInitialize.o \
         FastSiblingLocatorInitializeStaticChainingMesh.o \
diff --git a/src/enzo/hydro_rk/EvolveLevel_RK2.C b/src/enzo/hydro_rk/EvolveLevel_RK2.C
index aa7f19461..3b3477b8c 100644
--- a/src/enzo/hydro_rk/EvolveLevel_RK2.C
+++ b/src/enzo/hydro_rk/EvolveLevel_RK2.C
@@ -2,7 +2,7 @@
 /
 /  EVOLVE LEVEL USING RUNGE-KUTTA2 METHOD
 /
-/  written by: Peng Wang
+/  written by: Peng Wang & Tom Abel 
 /  date:       May, 2007
 /  modified1:
 /
@@ -77,9 +77,20 @@ int CallProblemSpecificRoutines(TopGridData * MetaData, HierarchyEntry *ThisGrid
 int PrepareDensityField(LevelHierarchyEntry *LevelArray[],
 			SiblingGridList SiblingList[],
 			int level, TopGridData *MetaData, FLOAT When);
+int SetAccelerationBoundary(HierarchyEntry *Grids[], int NumberOfGrids,
+			    SiblingGridList SiblingList[],
+			    int level, TopGridData *MetaData,
+			    ExternalBoundary *Exterior,
+			    LevelHierarchyEntry * Level,
+			    int CycleNumber);
 #else  // !FAST_SIB
 int PrepareDensityField(LevelHierarchyEntry *LevelArray[],
                         int level, TopGridData *MetaData, FLOAT When);
+int SetAccelerationBoundary(HierarchyEntry *Grids[], int NumberOfGrids,
+			    int level, TopGridData *MetaData, 
+			    ExternalBoundary *Exterior,
+			    LevelHierarchyEntry * Level,
+			    int CycleNumber);
 #endif  // end FAST_SIB
 
 int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids,
@@ -182,6 +193,8 @@ static int StaticLevelZero = 1;
 static int StaticLevelZero = 0;
 #endif
 
+extern int RK2SecondStepBaryonDeposit;
+
 /* EvolveGrid function */
 
 int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
@@ -365,6 +378,7 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
     }
 
     When = 0.5;
+    RK2SecondStepBaryonDeposit = 0;
     if (SelfGravity) {
 #ifdef FAST_SIB
       PrepareDensityField(LevelArray, SiblingList, level, MetaData, When);
@@ -425,33 +439,52 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 	}
       } // ENDIF UseHydro
 	
-      /* Do this here so that we can get the correct
-	 time interpolated boundary condition */
+      /* Do this here so that we can get the correct time interpolated boundary condition */
       Grids[grid1]->GridData->SetTimeNextTimestep();
       
     }  // end loop over grids
       
-    /*if (SetBoundaryConditions(Grids, NumberOfGrids,SiblingList,level, MetaData,
-      Exterior) == FAIL) {
-      ENZO_FAIL("");
-      }*/
+
 #ifdef FAST_SIB
     SetBoundaryConditions(Grids, NumberOfGrids, SiblingList, level, MetaData, Exterior, LevelArray[level]);
 #else
     SetBoundaryConditions(Grids, NumberOfGrids, level, MetaData, Exterior, LevelArray[level]);
 #endif
 
+    // Recompute potential and accelerations with time centered baryon Field
+    // this also does the particles again at the moment so could be made more efficient.
 
+    RK2SecondStepBaryonDeposit = 0; // set this to (0/1) to (not use/use) this extra step
+    printf("SECOND STEP\n");
+    if (RK2SecondStepBaryonDeposit & SelfGravity && UseHydro) {
+      When = 0.5;
+#ifdef FAST_SIB
+      PrepareDensityField(LevelArray, SiblingList, level, MetaData, When);
+#else  
+      PrepareDensityField(LevelArray, level, MetaData, When);
+#endif  // end FAST_SIB
+    }
 
     for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
 
+      /* Gravity: compute acceleration field for grid and particles. */
+      if (RK2SecondStepBaryonDeposit) {
+	int Dummy;
+	if (level <= MaximumGravityRefinementLevel) {
+	  if (level > 0) 
+	    Grids[grid1]->GridData->SolveForPotential(level) ;
+	  Grids[grid1]->GridData->ComputeAccelerations(level) ;
+	}
+	Grids[grid1]->GridData->ComputeAccelerationFieldExternal() ;
+      } // end: if (SelfGravity)
+
       if (UseHydro) {
 	if (HydroMethod == HD_RK)
 	  Grids[grid1]->GridData->RungeKutta2_2ndStep
 	    (SubgridFluxesEstimate[grid1], NumberOfSubgrids[grid1], level, Exterior);
 
 	else if (HydroMethod == MHD_RK) {
-
+	  
 	  Grids[grid1]->GridData->MHDRK2_2ndStep
 	    (SubgridFluxesEstimate[grid1], NumberOfSubgrids[grid1], level, Exterior);
 	  
@@ -461,7 +494,6 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 	  if (UseResistivity) 
 	    Grids[grid1]->GridData->AddResistivity();
 	
-	 
 	  time1 = ReturnWallTime();
 	  
 	  Grids[grid1]->GridData->PoissonSolver(level);
@@ -471,9 +503,8 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 
       /* Add viscosity */
 
-      if (UseViscosity) {
+      if (UseViscosity) 
 	Grids[grid1]->GridData->AddViscosity();
-      }
 
       /* Solve the cooling and species rate equations. */
  
@@ -523,32 +554,6 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
     StarParticleFinalize(Grids, MetaData, NumberOfGrids, LevelArray,
 			 level, AllStars);
 
-    /* Reompute potential if output is requested */
-    if (SelfGravity && WritePotential) {
-      CopyGravPotential = TRUE;
-      When = 0.0;
- 
-#ifdef FAST_SIB
-      PrepareDensityField(LevelArray, SiblingList, level, MetaData, When);
-#else   // !FAST_SIB
-      PrepareDensityField(LevelArray, level, MetaData, When);
-#endif  // end FAST_SIB
- 
-      //      CopyGravPotential = FALSE;
- 
-      for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
-        if (level <= MaximumGravityRefinementLevel) {
- 
-          /* Compute the potential. */
- 
-          if (level > 0)
-            Grids[grid1]->GridData->SolveForPotential(level);
-          Grids[grid1]->GridData->CopyPotentialToBaryonField();
-        }
-      } //  end loop over grids
-    } // if WritePotential
-
-
     OutputFromEvolveLevel(LevelArray,MetaData,level,Exterior);
     CallPython(LevelArray, MetaData, level);
 
diff --git a/src/enzo/hydro_rk/Grid_MHDRK2_2ndStep.C b/src/enzo/hydro_rk/Grid_MHDRK2_2ndStep.C
index a8911ad7d..739bc0cec 100644
--- a/src/enzo/hydro_rk/Grid_MHDRK2_2ndStep.C
+++ b/src/enzo/hydro_rk/Grid_MHDRK2_2ndStep.C
@@ -58,13 +58,6 @@ int grid::MHDRK2_2ndStep(fluxes *SubgridFluxes[],
   this->ReturnHydroRKPointers(Prim, false);
   this->ReturnOldHydroRKPointers(OldPrim, false);
 
-  // ok. this line should absolutely not be needed.
-  // however there is one piece of hardware i've compiled on where it will not work 
-  // without. There is absoltuely no reasons I can see how this can happen ... argh.
-  // this same line is in ReturnHydroRKPointers but for some reason n the kolob cluster
-  // in Heidlberg that is not enough. On my laptop no problem.
-  Prim[0] = BaryonField[0];
-
 #ifdef ECUDADEBUG
   printf("in Grid_MHDRK_2ndStep.C.\n");
   for (int j=30; j < 33; j++) 

From d9de205112da90784637fb6d26a19b0ef5fa4d6f Mon Sep 17 00:00:00 2001
From: tabel <none@none>
Date: Sat, 10 Oct 2009 19:14:31 -0700
Subject: [PATCH 09/12] More cleaning of EvolveLevel routines and adding the
 StreamingFormat write calls to EvolveLevel_RK2. The StreamingFormat for
 EvolveLevel_RK2 is not yet tested.

--HG--
branch : week-of-code
---
 src/enzo/EvolveLevel.C              | 15 +--------------
 src/enzo/hydro_rk/EvolveLevel_RK2.C | 30 +++++++++++++++--------------
 2 files changed, 17 insertions(+), 28 deletions(-)

diff --git a/src/enzo/EvolveLevel.C b/src/enzo/EvolveLevel.C
index 9f0999e54..72b21adf6 100644
--- a/src/enzo/EvolveLevel.C
+++ b/src/enzo/EvolveLevel.C
@@ -559,12 +559,7 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 
     /* Recompute radiation field, if requested. */
  
-    if (RadiationFieldType >= 10 && RadiationFieldType <= 11 &&
-	level <= RadiationFieldLevelRecompute)
-      if (RadiationFieldUpdate(LevelArray, level, MetaData) == FAIL) {
-	fprintf(stderr, "Error in RecomputeRadiationField.\n");
-	ENZO_FAIL("");
-      }
+    RadiationFieldUpdate(LevelArray, level, MetaData);
  
     /* Rebuild the Grids on the next level down.
        Don't bother on the last cycle, as we'll rebuild this grid soon. */
@@ -584,8 +579,6 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 	LevelZoneCycleCountPerProc[level] += NumberOfCells;
     }
  
-    
-
     cycle++;
     LevelCycleCount[level]++;
  
@@ -605,12 +598,6 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 
   /* Clean up. */
  
-#ifdef UNUSED
-  if (level > MaximumGravityRefinementLevel &&
-      level == MaximumRefinementLevel)
-    ZEUSQuadraticArtificialViscosity /= 1;
-#endif
- 
   delete [] NumberOfSubgrids;
   delete [] Grids;
   delete [] SubgridFluxesEstimate;
diff --git a/src/enzo/hydro_rk/EvolveLevel_RK2.C b/src/enzo/hydro_rk/EvolveLevel_RK2.C
index 3b3477b8c..79cd9c37e 100644
--- a/src/enzo/hydro_rk/EvolveLevel_RK2.C
+++ b/src/enzo/hydro_rk/EvolveLevel_RK2.C
@@ -177,9 +177,9 @@ int CallPython(LevelHierarchyEntry *LevelArray[], TopGridData *MetaData,
 
 static int LevelCycleCount[MAX_DEPTH_OF_HIERARCHY];
 static int MovieCycleCount[MAX_DEPTH_OF_HIERARCHY];
-//double LevelWallTime[MAX_DEPTH_OF_HIERARCHY];
-//double LevelZoneCycleCount[MAX_DEPTH_OF_HIERARCHY];
-//double LevelZoneCycleCountPerProc[MAX_DEPTH_OF_HIERARCHY];
+static double LevelWallTime[MAX_DEPTH_OF_HIERARCHY];
+static double LevelZoneCycleCount[MAX_DEPTH_OF_HIERARCHY];
+static double LevelZoneCycleCountPerProc[MAX_DEPTH_OF_HIERARCHY];
  
 static float norm = 0.0;            //AK
 static float TopGridTimeStep = 0.0; //AK
@@ -263,6 +263,11 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
   for (grid1 = 0; grid1 < NumberOfGrids; grid1++)
     Grids[grid1]->GridData->ClearBoundaryFluxes();
 
+  /* After we calculate the ghost zones, we can initialize streaming
+     data files (only on level 0) */
+
+  if (MetaData->FirstTimestepAfterRestart == TRUE && level == 0)
+    WriteStreamData(LevelArray, level, MetaData, MovieCycleCount);
 
   /* ================================================================== */
   /* Loop over grid timesteps until the elapsed time equals the timestep 
@@ -273,6 +278,11 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
     SetLevelTimeStep(Grids, NumberOfGrids, level, 
         &dtThisLevelSoFar, &dtThisLevel, dtLevelAbove);
 
+    /* Streaming movie output (write after all parent grids are
+       updated) */
+
+    WriteStreamData(LevelArray, level, MetaData, MovieCycleCount);
+
    /* Initialize the star particles */
 
     Star *AllStars = NULL;
@@ -529,7 +539,6 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
  
       if (UseFloor) 
 	Grids[grid1]->GridData->SetFloor();
-      
 
       /* If using comoving co-ordinates, do the expansion terms now. */
 
@@ -580,7 +589,6 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 
 #ifdef USE_JBPERF
     // Update lcaperf "level" attribute
-
     jbPerf.attribute ("level",&jb_level,JB_INT);
 #endif
 
@@ -651,9 +659,9 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
     for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
       Grids[grid1]->GridData->CollectGridInformation
         (GridMemory, GridVolume, NumberOfCells, AxialRatio, CellsTotal, Particles);
-      //      LevelZoneCycleCount[level] += NumberOfCells;
-      //if (MyProcessorNumber == Grids[grid1]->GridData->ReturnProcessorNumber())
-      //	LevelZoneCycleCountPerProc[level] += NumberOfCells;
+      LevelZoneCycleCount[level] += NumberOfCells;
+      if (MyProcessorNumber == Grids[grid1]->GridData->ReturnProcessorNumber())
+      	LevelZoneCycleCountPerProc[level] += NumberOfCells;
     }
 #endif
 
@@ -673,12 +681,6 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 
   /* Clean up. */
 
-#ifdef UNUSED
-  if (level > MaximumGravityRefinementLevel &&
-      level == MaximumRefinementLevel)
-    ZEUSQuadraticArtificialViscosity /= 1;
-#endif
-
   delete [] NumberOfSubgrids;
   delete [] Grids;
   delete [] SubgridFluxesEstimate;

From 6463db1badd32fb95184b9552eb68a84dfe23e02 Mon Sep 17 00:00:00 2001
From: tabel <none@none>
Date: Sat, 10 Oct 2009 19:27:09 -0700
Subject: [PATCH 10/12] Added AutoAdjustRefineRegion also to EvolveLevel_RK2
 and did a litte more cleaning up.

--HG--
branch : week-of-code
---
 src/enzo/EvolveLevel.C              |  2 --
 src/enzo/hydro_rk/EvolveLevel_RK2.C | 34 +++++++++++++++++++----------
 2 files changed, 23 insertions(+), 13 deletions(-)

diff --git a/src/enzo/EvolveLevel.C b/src/enzo/EvolveLevel.C
index 72b21adf6..f3c9c113e 100644
--- a/src/enzo/EvolveLevel.C
+++ b/src/enzo/EvolveLevel.C
@@ -281,8 +281,6 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
     ENZO_FAIL("");
 #endif
  
-
-
   /* Clear the boundary fluxes for all Grids (this will be accumulated over
      the subcycles below (i.e. during one current grid step) and used to by the
      current grid to correct the zones surrounding this subgrid (step #18). */
diff --git a/src/enzo/hydro_rk/EvolveLevel_RK2.C b/src/enzo/hydro_rk/EvolveLevel_RK2.C
index 79cd9c37e..4c25d0941 100644
--- a/src/enzo/hydro_rk/EvolveLevel_RK2.C
+++ b/src/enzo/hydro_rk/EvolveLevel_RK2.C
@@ -40,6 +40,8 @@ int StarParticleInitialize(LevelHierarchyEntry *LevelArray[], int ThisLevel,
 int StarParticleFinalize(HierarchyEntry *Grids[], TopGridData *MetaData,
 			 int NumberOfGrids, LevelHierarchyEntry *LevelArray[], 
 			 int level, Star *&AllStars);
+int AdjustRefineRegion(LevelHierarchyEntry *LevelArray[], 
+		       TopGridData *MetaData, int EL_level);
 
 #ifdef TRANSFER
 int EvolvePhotons(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
@@ -220,6 +222,7 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 #endif
 
   FLOAT When;
+
   float DensityUnits = 1.0, LengthUnits = 1.0, TemperatureUnits = 1, TimeUnits, 
     VelocityUnits, CriticalDensity = 1, BoxLength = 1, MagneticUnits;
   double MassUnits;
@@ -235,14 +238,23 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
   int *NumberOfSubgrids = new int[NumberOfGrids];
   fluxes ***SubgridFluxesEstimate = new fluxes **[NumberOfGrids];
 
+  /* Initialize the chaining mesh used in the FastSiblingLocator. */
+
+
   SiblingGridList *SiblingList = new SiblingGridList[NumberOfGrids];
   CreateSiblingList(Grids, NumberOfGrids, SiblingList, StaticLevelZero, 
 		    MetaData, level);
 
-  /*if (SetBoundaryConditions(Grids, NumberOfGrids,SiblingList,level, MetaData,
-			      Exterior) == FAIL) {
-    ENZO_FAIL("");
-    }*/
+  /* On the top grid, adjust the refine region so that only the finest
+     particles are included.  We don't want the more massive particles
+     to contaminate the high-resolution region. */
+
+  AdjustRefineRegion(LevelArray, MetaData, level);
+
+  /* ================================================================== */
+  /* For each grid: a) interpolate boundaries from its parent.
+                    b) copy any overlapping zones.  */
+ 
 #ifdef FAST_SIB
   if (SetBoundaryConditions(Grids, NumberOfGrids, SiblingList,
 			    level, MetaData, Exterior, LevelArray[level]) == FAIL)
@@ -252,14 +264,16 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
                             Exterior, LevelArray[level]) == FAIL)
     ENZO_FAIL("");
 #endif
-
-
-  
+ 
   /* Count the number of colours in the first grid (to define NColor) */
 
   Grids[0]->GridData->SetNumberOfColours();
   //fprintf(stdout, "EvolveLevel_RK2: NColor =%d, NSpecies = %d\n", NColor, NSpecies);
 
+  /* Clear the boundary fluxes for all Grids (this will be accumulated over
+     the subcycles below (i.e. during one current grid step) and used to by the
+     current grid to correct the zones surrounding this subgrid (step #18). */
+
   for (grid1 = 0; grid1 < NumberOfGrids; grid1++)
     Grids[grid1]->GridData->ClearBoundaryFluxes();
 
@@ -439,15 +453,13 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 
       Grids[grid1]->GridData->CopyBaryonFieldToOldBaryonField();
 
-      if (UseHydro) {
+      if (UseHydro) 
 	if (HydroMethod == HD_RK)
 	  Grids[grid1]->GridData->RungeKutta2_1stStep
 	    (SubgridFluxesEstimate[grid1], NumberOfSubgrids[grid1], level, Exterior);
-	else if (HydroMethod == MHD_RK) {
+	else if (HydroMethod == MHD_RK) 
 	  Grids[grid1]->GridData->MHDRK2_1stStep
 	    (SubgridFluxesEstimate[grid1], NumberOfSubgrids[grid1], level, Exterior);
-	}
-      } // ENDIF UseHydro
 	
       /* Do this here so that we can get the correct time interpolated boundary condition */
       Grids[grid1]->GridData->SetTimeNextTimestep();

From 5f10ffea5a6596b74f46a3a95208b47dce7813e4 Mon Sep 17 00:00:00 2001
From: tabel <none@none>
Date: Sat, 10 Oct 2009 20:00:59 -0700
Subject: [PATCH 11/12] Split computation of Dedner Wave speeds into a separate
 file.

--HG--
branch : week-of-code
---
 src/enzo/Make.config.objects                |  1 +
 src/enzo/hydro_rk/ComputeDednerWaveSpeeds.C | 77 +++++++++++++++++++++
 src/enzo/hydro_rk/EvolveLevel_RK2.C         | 55 +++------------
 3 files changed, 89 insertions(+), 44 deletions(-)
 create mode 100644 src/enzo/hydro_rk/ComputeDednerWaveSpeeds.C

diff --git a/src/enzo/Make.config.objects b/src/enzo/Make.config.objects
index f24393146..f253ade8d 100644
--- a/src/enzo/Make.config.objects
+++ b/src/enzo/Make.config.objects
@@ -701,6 +701,7 @@ OBJS_HYDRO_RK = \
       	hydro_rk/AGNDiskInitialize.o \
         hydro_rk/Collapse1DInitialize.o \
         hydro_rk/Collapse3DInitialize.o \
+        hydro_rk/ComputeDednerWaveSpeeds.o \
         hydro_rk/EvolveLevel_RK2.o \
         hydro_rk/GalaxyDiskInitialize.o \
         hydro_rk/Grid_AddSelfGravity.o \
diff --git a/src/enzo/hydro_rk/ComputeDednerWaveSpeeds.C b/src/enzo/hydro_rk/ComputeDednerWaveSpeeds.C
new file mode 100644
index 000000000..c66dcecf0
--- /dev/null
+++ b/src/enzo/hydro_rk/ComputeDednerWaveSpeeds.C
@@ -0,0 +1,77 @@
+/***********************************************************************
+/
+/  COMPUTE DEDNER WAVE SPEEDS
+/
+/  written by: Tom Abel
+/  date:       October 2009
+/
+/  ======================================================================= 
+/  This routine computes the wave speeds used for the Dedner formalism.
+/
+/       Reference: e.g.  Matsumoto, PASJ, 2007, 59, 905 
+/
+/  Output: C_h and C_p global variables defined in global_data.h
+/
+************************************************************************/ 
+
+#include <stdio.h>
+#include <math.h>
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "TopGridData.h"
+#include "Grid.h"
+#include "Hierarchy.h"
+#include "LevelHierarchy.h"
+#include "../hydro_rk/tools.h"
+
+int GetUnits(float *DensityUnits, float *LengthUnits,
+	     float *TemperatureUnits, float *TimeUnits,
+	     float *VelocityUnits, FLOAT Time);
+
+int ComputeDednerWaveSpeeds(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], 
+			    int level, FLOAT dt0)
+{
+  /* Count the number of grids on this level. */
+ 
+  if (HydroMethod != MHD_RK) 
+    return SUCCESS;
+
+  float DensityUnits = 1.0, LengthUnits = 1.0, TemperatureUnits = 1, TimeUnits, 
+    VelocityUnits, CriticalDensity = 1, BoxLength = 1, MagneticUnits;
+  double MassUnits;
+  GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, &TimeUnits, &VelocityUnits, 1.0);
+  MassUnits = DensityUnits*pow(LengthUnits,3);
+
+  int lmax;
+  LevelHierarchyEntry *Temp;
+  for (lmax = MAX_DEPTH_OF_HIERARCHY-1; lmax >= 0; lmax--) {
+    Temp = LevelArray[lmax];
+    if (Temp != NULL) 
+      break;
+  }
+
+  //      lmax = 0; // <- Pengs version had lmax = 6
+  FLOAT dx0, dy0, dz0, h_min, DivBDampingLength = 1.0;
+  
+  dx0 = (DomainRightEdge[0] - DomainLeftEdge[0]) / MetaData->TopGridDims[0];
+  dy0 = (MetaData->TopGridRank > 1) ? 
+    (DomainRightEdge[1] - DomainLeftEdge[1]) / MetaData->TopGridDims[1] : 1e8;
+  dz0 = (MetaData->TopGridRank > 2) ? 
+    (DomainRightEdge[2] - DomainLeftEdge[2]) / MetaData->TopGridDims[2] : 1e8;
+  h_min = my_MIN(dx0, dy0, dz0);
+  h_min /= pow(RefineBy, lmax);
+  C_h = 0.5*MetaData->CourantSafetyNumber*(h_min/dt0);
+  C_h = min( C_h, 1e6/VelocityUnits); // never faster than __ cm/s (for very small dt0 a problems)
+  if (EOSType == 3)  // for isothermal runs just use the constant sound speed
+    C_h = EOSSoundSpeed;
+  
+  C_p = sqrt(0.18*DivBDampingLength*C_h);
+
+  return SUCCESS;
+}
+ 
+  
diff --git a/src/enzo/hydro_rk/EvolveLevel_RK2.C b/src/enzo/hydro_rk/EvolveLevel_RK2.C
index 4c25d0941..0224be8cc 100644
--- a/src/enzo/hydro_rk/EvolveLevel_RK2.C
+++ b/src/enzo/hydro_rk/EvolveLevel_RK2.C
@@ -42,7 +42,8 @@ int StarParticleFinalize(HierarchyEntry *Grids[], TopGridData *MetaData,
 			 int level, Star *&AllStars);
 int AdjustRefineRegion(LevelHierarchyEntry *LevelArray[], 
 		       TopGridData *MetaData, int EL_level);
-
+int ComputeDednerWaveSpeeds(TopGridData *MetaData,LevelHierarchyEntry *LevelArray[], 
+			    int level, FLOAT dt0);
 #ifdef TRANSFER
 int EvolvePhotons(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 		  Star *AllStars, FLOAT GridTime, int level, int LoopTime = TRUE);
@@ -161,9 +162,6 @@ int FastSiblingLocatorFinalize(ChainingMeshStructure *Mesh);
 int FastSiblingLocatorInitializeStaticChainingMesh(ChainingMeshStructure *Mesh, int Rank,
 						   int TopGridDims[]); 
 
-int GetUnits(float *DensityUnits, float *LengthUnits,
-	     float *TemperatureUnits, float *TimeUnits,
-	     float *VelocityUnits, FLOAT Time);
 double ReturnWallTime();
 int CallPython(LevelHierarchyEntry *LevelArray[], TopGridData *MetaData,
                int level);
@@ -223,12 +221,6 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 
   FLOAT When;
 
-  float DensityUnits = 1.0, LengthUnits = 1.0, TemperatureUnits = 1, TimeUnits, 
-    VelocityUnits, CriticalDensity = 1, BoxLength = 1, MagneticUnits;
-  double MassUnits;
-  GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
-	   &TimeUnits, &VelocityUnits, 1.0);
-  MassUnits = DensityUnits*pow(LengthUnits,3);
   
   /* Create an array (Grids) of all the grids. */
 
@@ -366,40 +358,15 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 
     } // end loop over grids (create Subgrid list)
 
-    /* compute wave speed
-       Reference: Matsumoto, PASJ, 2007, 59, 905 */
-
-    if (HydroMethod == MHD_RK) {
-      int lmax;
-      LevelHierarchyEntry *Temp;
-      for (lmax = MAX_DEPTH_OF_HIERARCHY-1; lmax >= 0; lmax--) {
-	Temp = LevelArray[lmax];
-	if (Temp != NULL) {
-	  break;
-	} 
-      }
-      //      lmax = 0; // <- Pengs version had lmax = 6
-      //      lmax = 1;
-      FLOAT dx0, dy0, dz0, h_min, DivBDampingLength = 1.0;
-
-      dx0 = (DomainRightEdge[0] - DomainLeftEdge[0]) / MetaData->TopGridDims[0];
-      dy0 = (MetaData->TopGridRank > 1) ? 
-	(DomainRightEdge[1] - DomainLeftEdge[1]) / MetaData->TopGridDims[1] : 1e8;
-      dz0 = (MetaData->TopGridRank > 2) ? 
-	(DomainRightEdge[2] - DomainLeftEdge[2]) / MetaData->TopGridDims[2] : 1e8;
-      h_min = my_MIN(dx0, dy0, dz0);
-      h_min /= pow(RefineBy, lmax);
-      C_h = 0.5*MetaData->CourantSafetyNumber*(h_min/dt0);
-      C_h = min( C_h, 1e6/VelocityUnits); // never faster than __ cm/s (for very small dt0 a problems)
-      if (EOSType == 3)  // for isothermal runs just use the constant sound speed
-	C_h = EOSSoundSpeed;
-
-      C_p = sqrt(0.18*DivBDampingLength*C_h);
-      //      C_p = sqrt(0.18*DivBDampingLength)*C_h;
-      if (debug) 
-	fprintf(stderr, "lengthscale %g timestep: %g  C_h: %g  C_p: %g\n ", 
-		h_min, dt0, C_h, C_p);
-    }
+    /* compute wave speed Reference: Matsumoto, PASJ, 2007, 59, 905 */
+
+    if (HydroMethod == MHD_RK) 
+      ComputeDednerWaveSpeeds(MetaData, LevelArray, level, dt0);
+	
+    if (debug && HydroMethod == MHD_RK) 
+      fprintf(stderr, "wave speeds: timestep: %g  C_h: %g  C_p: %g\n ", 
+	       dt0, C_h, C_p);
+
 
     When = 0.5;
     RK2SecondStepBaryonDeposit = 0;

From 64be9f1ddc0d6df2051c4f220dd7415d40f27f17 Mon Sep 17 00:00:00 2001
From: tabel <none@none>
Date: Sat, 10 Oct 2009 20:08:33 -0700
Subject: [PATCH 12/12] remove some debugging printf's

--HG--
branch : week-of-code
---
 src/enzo/Grid_DepositBaryons.C      | 2 +-
 src/enzo/hydro_rk/EvolveLevel_RK2.C | 3 +--
 2 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/enzo/Grid_DepositBaryons.C b/src/enzo/Grid_DepositBaryons.C
index ac227184d..07afc1c0c 100644
--- a/src/enzo/Grid_DepositBaryons.C
+++ b/src/enzo/Grid_DepositBaryons.C
@@ -217,7 +217,7 @@ int grid::DepositBaryons(grid *TargetGrid, FLOAT DepositTime)
       input_density = av_dens;
     }
 
-    printf("DepositBaryons, %i\n", RK2SecondStepBaryonDeposit);
+    //    printf("DepositBaryons, %i\n", RK2SecondStepBaryonDeposit);
     
     FORTRAN_NAME(dep_grid_cic)(input_density, dens_field, vel_field,
 			       input_velx, input_vely, input_velz,
diff --git a/src/enzo/hydro_rk/EvolveLevel_RK2.C b/src/enzo/hydro_rk/EvolveLevel_RK2.C
index 0224be8cc..1e91c4f64 100644
--- a/src/enzo/hydro_rk/EvolveLevel_RK2.C
+++ b/src/enzo/hydro_rk/EvolveLevel_RK2.C
@@ -299,7 +299,6 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
     RadiativeTransferPrepare(LevelArray, level, MetaData, AllStars, dtLevelAbove);
 #endif /* TRANSFER */
 
-
     /* For each grid, compute the number of it's subgrids. */
 
     for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
@@ -444,7 +443,7 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
     // this also does the particles again at the moment so could be made more efficient.
 
     RK2SecondStepBaryonDeposit = 0; // set this to (0/1) to (not use/use) this extra step
-    printf("SECOND STEP\n");
+    //    printf("SECOND STEP\n");
     if (RK2SecondStepBaryonDeposit & SelfGravity && UseHydro) {
       When = 0.5;
 #ifdef FAST_SIB