Commits

John Wise  committed b86d52b Merge

Merging heads.

  • Participants
  • Parent commits 7844600, 4a52471

Comments (0)

Files changed (194)

    * Greg Bryan             gbryan@astro.columbia.edu
    * Renyue Cen             cen@astro.princeton.edu
    * Dave Collins           dcollins@physics.ucsd.edu
+   * Brian Crosby           crosby.bd@gmail.com
    * Nathan Goldbaum        goldbaum@ucolick.org
+   * Oliver Hahn            hahn@phys.ethz.ch
    * Robert Harkness        harkness@sdsc.edu
    * Elizabeth Harper-Clark h-clark@astro.utoronto.ca
    * Cameron Hummels        chummels@gmail.com
+   * Ji-hoon Kim            mornkr@ucolick.org
+   * Daegene Koh            dkoh30@gatech.edu
    * Alexei Kritsuk         akritsuk@ucsd.edu
    * Michael Kuhlen         mqk@astro.berkeley.edu
    * Eve Lee                elee@cita.utoronto.ca
    * JS Oishi               jsoishi@gmail.com
    * Brian O'Shea           oshea@msu.edu
    * Pascal Paschos         ppaschos@minbari.ucsd.edu
+   * Carolyn Peruta         perutaca@msu.edu
    * Alex Razoumov          razoumov@gmail.com
    * Dan Reynolds           reynolds@smu.edu
+   * Munier Salem           msalem@astro.columbia.edu
    * Christine Simpson      csimpson@astro.columbia.edu
    * Samuel Skillman        samskillman@gmail.com 
    * Stephen Skory          s@skory.us
    * Britton Smith          brittonsmith@gmail.com
+   * Geoffrey So            gsiisg@gmail.com
    * Elizabeth Tasker       tasker@astro1.sci.hokudai.ac.jp
    * Matthew Turk           matthewturk@gmail.com
    * Rick Wagner            rwagner@physics.ucsd.edu
    * Peng Wang              pengw@slac.stanford.edu
    * John Wise              jwise@physics.gatech.edu
+   * Hao Xu                 haxu@ucsd.edu
    * Fen Zhao               fenzhao@stanford.edu
 
 

File doc/manual/source/developer_guide/FineGrainedOutput.rst

+Fine Grained Output
+===================
+
+When making significant changes to Enzo that have non-local impact, such as adding a new accretion
+mechanism for sink particles or face centered magnetic fields, there are many
+places to introduce errors.  In order to examine the effect of changes at
+specific points in the code, one can use ``ExtraOutputs``.  This run time
+parameter makes a call to ``WriteAllData`` at various points in ``EvolveLevel``.
+For instance, putting::
+
+    ExtraOutput = 1
+    StopCycle = 2
+    MaximumRefinementLevel = 0
+
+will cause::
+
+    ED01_0000
+    ED01_0001
+
+to be written, along with your regular outputs.  With one level of refinement,
+six outputs will be written.  The relation between output number and position is
+below.  Up to 10 output points can be specified.  
+
+
+Unraveling what output gets written when can be a challenge.  One technique is
+to run with ``-d``, and use the following command::
+
+    egrep "^Level||ExtraOutput" output.log
+
+on the output log, which will show what output gets called on which level, and a string indicating
+at which point in ``EvolveLevel`` it was called.
+
+It should be noted that ``ExtraOutputs`` is not written into parameter files on
+data dumps, though it can be added to restart parameter files.  This is to prevent absurd amounts of data being written to disk.
+By design, this technique outputs many data dumps for each root grid timestep,
+following the W cycle.
+This has the added disadvantage of making the code slower, as disk access is
+rarely the fastest part of any machine.   
+
+In the code, overhead is minimized by wrapping the full function signature in a macro.  New calls can be added with::
+
+    EXTRA_OUTPUT_MACRO(42, "After My Special Purpose")
+
+where, of course, 42 is replaced by an integer not used by another output, and the string represents the location in the code.  It is often instructive to include this output mechanism in ``EvolveHierarchy`` as well, though this has not been done in the current checkin.
+
+
+
+Here's a table of output number vs. position in ``EvolveLevel``.  Please refer to the :ref:`FlowChart` 
+to understand each entry. The non-continuity represents some outputs that will be introduced when 
+MHDCT is merged, but not relevant for pure hydro.
+
+===== ===================================
+Index Position in ``EvolveLevel``
+----- -----------------------------------
+1     Before time loop
+2     After SolveHydroEquations grid loop
+25     After SetBoundaryConditions
+3     Before UpdateFromFinerGrids
+4     After UpdateFromFinerGrids
+6     After the time loop
+===== ===================================

File doc/manual/source/developer_guide/FloatIsDouble.rst

 .. code-block:: c
 
     c     Declarations
-          R_PREC     third
+          R_PREC     third, tenth
           INTG_PREC  one
           P_PREC     fifth
           CMPLX_PREC two_i
 
     c     Calculations
           third = 1._RKIND / 3._RKIND
+	  tenth = 1.e-1_RKIND
           one   = 1_IKIND
           fifth = real(1, PKIND) / 5._PKIND
           two_i = (0._RKIND, 2._RKIND)

File doc/manual/source/developer_guide/index.rst

    ProgrammingGuide.rst
    FilenameConventions.rst
    DebuggingWithGDB.rst
+   FineGrainedOutput.rst
    AddingNewParameters.rst
    HowToAddNewBaryonField.rst
    FloatIsDouble.rst

File doc/manual/source/parameters/problemtypes.rst

 102          :ref:`1dcollapse_param`
 106          :ref:`mhdhydro_param`
 107          :ref:`putsink_param`
-108          Cluster Cooling Flow 
+108          :ref:`clustercoolingflow_param` 
 200          :ref:`mhd1d_param`
 201          :ref:`mhd2d_param`
 202          :ref:`mhd3d_param`
 ``PutSinkRestartName`` (external)
      Filename to restart from. 
 
+
+.. _clustercoolingflow_param:
+
+Cluster Cooling Flow (108)
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+``ClusterSMBHFeedback`` (external)
+    Boolean flag. Default: FALSE
+``ClusterSMBHJetMdot`` (external)
+    Mdot of one Jet. Units: Solar mass per year. Default: 3.0
+``ClusterSMBHJetVelocity`` (external)
+    Units:km/s. Default: 10000.0
+``ClusterSMBHJetRadius`` (external)
+    The radius of the jet launching region. Units: cell width. Default: 6.0
+``ClusterSMBHJetLaunchOffset`` (external)
+    The distance of the jet launching plane to the center of the cluster. Units: cell width. Default: 10.0
+``ClusterSMBHStartTime`` (external)
+    The time to start feedback in code unit. Default: 1.0
+``ClusterSMBHTramp`` (external)
+    The ramp time in Myr. Default: 0.1
+``ClusterSMBHJetOpenAngleRadius`` (external)
+    Default: 0.0
+``ClusterSMBHFastJetRadius`` (external)
+    Default: 0.1
+``ClusterSMBHFastJetVelocity`` (external)
+    Unit: km/s. Default: 10000.0
+``ClusterSMBHJetEdot`` (external)
+    Unit: 10^44 ergs/s. Default: 1.0
+``ClusterSMBHKineticFraction`` (external)
+    The fraction of kinetic energy feedback; the rest is thermal feedback. Default: 1.0
+``ClusterSMBHJetAngleTheta`` (external)
+    The angle of the jet direction with respect to z-axis. Default: 0.0 (along the axis)
+``ClusterSMBHJetAnglePhi`` (external)
+    Default: 0.0
+``ClusterSMBHJetPrecessionPeriod`` (external)
+    Unit: Myr. Default: 0.0 (not precessing)
+``ClusterSMBHCalculateGasMass`` (external)
+    Type: integer. 1--Calculate the amount of cold gas around the SMBH and remove it at the rate of 2*Mdot; 2--Calculate Mdot based on the amount of cold gas around the SMBH; 0--off (do not remove cold gas). Default: 1.
+``ClusterSMBHFeedbackSwitch`` (external)
+    Boolean flag. When ClusterSMBHCalculateGasMass=1, ClusterSMBHFeedbackSwitch is turned on when there is enough cold gas (ClusterSMBHEnoughColdGas) around the SMBH. Default: FALSE
+``ClusterSMBHEnoughColdGas`` (external)
+    Unit: Solar mass. Default: 1.0e7
+``ClusterSMBHAccretionTime`` (external)
+    When ClusterSMBHCalculateGasMass = 2, Mdot = Mcold/ClusterSMBHAccretionTime. Default: 5.0 (Myr)
+``ClusterSMBHJetDim`` (external)
+    0--x; 1--y; 2--z. Default: 2
+``ClusterSMBHAccretionEpsilon`` (external)
+    Jet Edot = ClusterSMBHAccretionEpsilon * Mdot * c^2. Default: 0.001
+
+
 .. _mhd1d_param:
 
 1D MHD Test (200)

File doc/manual/source/user_guide/CosmologicalInitialConditions.rst

      RefineRegionLeftEdge          = 0.15523 0.14551 0.30074
      RefineRegionRightEdge         = 0.38523 0.37551 0.53074
      NewCenterFloat                = 0.270230055 0.260508984 0.415739357
+     AutomaticSubgridBuffer        = 4
 
 MaximumInitialRefinementLevel indicates how many extra levels you want
 to generate (in this case two additional levels, or 3 in total,
 including the root grid).  The next two parameters
 (RefineRegionLeftEdge and RefineRegionRightEdge) describe the region
 to be refined.  The fourth (optional) parameter re-centers the grid on
-the halo to be resimulated.
+the halo to be resimulated.  The fifth parameter (AutomaticSubgridBuffer)
+indicates how many course cells should be added around each refined
+region.
 
 Once you have added these parameters, run inits once on the new
-parameter file.  It will give you a progress report as it runs (note
+parameter file in the standard way:  
+
+::
+
+     inits [-d] MultiGridParameterFile
+
+It will give you a progress report as it runs (note
 that if MaximumInitialRefinementLevel is large, this can take a long
 time), and generate all of the necessary files (e.g.  GridDensity.0,
 GridDensity.1, etc.).
 **NewCenterFloat**
     Indicates that the final grid should be recenter so that this point
     is the new center (0.5 0.5 0.5) of the grid.
+**AutomaticSubgridBuffer**
+    For multi-grid (nested) initial code generation (with the above
+    parameters).  This parameter controls how many coarse cells are
+    added around each refined region as buffer zones.  The value
+    of 1 is probably ok, but larger values (4?) are probably safer.
+    Default: 1
 **MaxDims**
     All dimensions are specified as one to three numbers deliminated by
     spaces (and for those familiar with the KRONOS or ZEUS method of

File doc/manual/source/user_guide/FlowChart.rst

+.. _FlowChart:
+
 Enzo Flow Chart, Source Browser
 ===============================
 

File run/Cooling/OneZoneFreefallTest/OneZoneFreefallTest.enzo

 #
 TopGridRank = 2
 
-TopGridDimensions = 3 5
+TopGridDimensions = 3 3
 
 OneZoneFreefallTestInitialDensity = 1.0
-OneZoneFreefallTestMinimumEnergy = 1.0
-OneZoneFreefallTestMaximumEnergy = 10.0
+OneZoneFreefallTestMinimumEnergy = 2.0
+OneZoneFreefallTestMaximumEnergy = 2.0
 OneZoneFreefallTestMinimumMetallicity = 1.0e-7
-OneZoneFreefallTestMaximumMetallicity = 1.0e-3
+OneZoneFreefallTestMaximumMetallicity = 1.0e-5
 
 # Set timestep as a fraction of free-fall time
 OneZoneFreefallTimestepFraction = 1e-2
+OneZoneFreefallUseEffectiveGamma = 1
+Gamma = 1.333
 
 #
 #  set I/O and stop/start parameters
 #
-StopTime                  = 50
+StopTime                  = 100
 StopCycle                 = 3000
 CycleSkipDataDump         = 150
 DataDumpDir               = DD

File run/Cooling/OneZoneFreefallTest/plot.py

 
 for pf in es:
     T.append(pf.h.grids[0]['Temperature'])
-    n.append(pf.h.grids[0]['NumberDensity'][0,0,0])
+    n.append(pf.h.grids[0]['NumberDensity'])
     Z.append(pf.h.grids[0]['Metallicity'])
     if do_fH2: fH2.append(pf.h.grids[0]['H2I_Fraction'])
     if do_t_cool:
         t_cool.append(pf.h.grids[0]['Cooling_Time'])
-        t_dyn.append(pf.h.grids[0]['DynamicalTime'][0,0,0])
+        t_dyn.append(pf.h.grids[0]['DynamicalTime'])
     if dust: Tdust.append(pf.h.grids[0]['Dust_Temperature'])
     del pf
 
 
 met = na.round(na.log10(Z[0,0,:,0]))
 for i in range(T.shape[2]):
-    pylab.loglog(n, T[:, 0, i, 0], 
+    pylab.loglog(n[:, 0, i, 0], T[:, 0, i, 0], 
                  label='log (Z/Z$_{\odot}$) = %d' % met[i],
                  color=colors[i], linestyle='-')
     if dust:
-        pylab.loglog(n, Tdust[:, 0, i, 0], 
+        pylab.loglog(n[:, 0, i, 0], Tdust[:, 0, i, 0], 
                      color=colors[i], linestyle='--')
 pylab.xlim(xmin=1.0)
 pylab.ylim(1e0, 1e4)
 pylab.xlabel('n [cm$^{-3}$]')
 pylab.ylabel('T [K]')
-pylab.legend(labelspacing=0.0, loc='upper left')
+pylab.legend(labelspacing=0.0, loc='lower right')
 pylab.savefig('n-T_%s.png' % keyword)
 pylab.clf()
 
 if do_fH2:
     for i in range(T.shape[2]):
-        pylab.loglog(n, fH2[:, 0, i, 0], 
+        pylab.loglog(n[:, 0, i, 0], fH2[:, 0, i, 0], 
                      label='log (Z/Z$_{\odot}$) = %d' % met[i],
                      color=colors[i])
     pylab.xlim(xmin=1.0)
 
 if do_t_cool:
     for i in range(T.shape[2]):
-        pylab.loglog(n, (t_cool[:, 0, i, 0]/t_dyn), 
+        pylab.loglog(n[:, 0, i, 0], (t_cool[:, 0, i, 0]/t_dyn[:, 0, i, 0]), 
                      label='log (Z/Z$_{\odot}$) = %d' % met[i],
                      color=colors[i])
     pylab.xlim(xmin=1.0)

File run/DrivenTurbulence3D/DrivenTurbulence3D.enzo

File contents unchanged.

File run/GravitySolver/GravityTest/GravityTest__test_gravity.py

 from yt.frontends.enzo.answer_testing_support import \
      requires_outputlog
 
+tolerance = ytcfg.getint("yt", "answer_testing_tolerance")
+     
 _pf_name = os.path.basename(os.path.dirname(__file__)) + ".enzo"
 _dir_name = os.path.dirname(__file__)
 
         return rmsError
 
     def compare(self, new_result, old_result):
-        assert_allclose(new_result, old_result, rtol=1e-1, atol=0)
+        assert_allclose(new_result, old_result, rtol=10.**-tolerance, atol=0)
 
 @requires_outputlog(_dir_name, _pf_name)
 def test_gravity_test():

File run/RadiationTransportFLD/RHIonization2_sp/iliev2_check.py

     tval, vol, gamma, dUnit, tUnit, lUnit = get_params(pfile)
     f = h5py.File(hfile,'r')
     Eg = f.get('/Grid00000001/Grey_Radiation_Energy')
-    energy = f.get('/Grid00000001/Total_Energy')
+    energy = f.get('/Grid00000001/TotalEnergy')
     HI = f.get('/Grid00000001/HI_Density')
     HII = f.get('/Grid00000001/HII_Density')
     rho = f.get('/Grid00000001/Density')

File run/RadiationTransportFLD/RHIonization2_sp/iliev2_makeplots.py

     tval, vol, gamma, dUnit, tUnit, lUnit = get_params(pfile)
     f = h5py.File(hfile,'r')
     Eg = f.get('/Grid00000001/Grey_Radiation_Energy')
-    energy = f.get('/Grid00000001/Total_Energy')
+    energy = f.get('/Grid00000001/TotalEnergy')
     HI = f.get('/Grid00000001/HI_Density')
     HII = f.get('/Grid00000001/HII_Density')
     rho = f.get('/Grid00000001/Density')

File run/test_runner.py

         self.test_dir = test_dir
         self.test_data = test_data
         self.exe_path = exe_path
+        self.finished = False
         self.results = {}
         if self.exe_path is None:
             self.local_exe = None

File src/enzo/CUDAUtil.cu

+/***********************************************************************
+/
+/  GPU ROUTINES
+/
+/  written by: Peng Wang 
+/  date:       June, 2012
+/  modified1:
+/
+/  PURPOSE:
+/         Init GPU.
+/
+************************************************************************/
+
+#include <stdio.h>
+#include "macros_and_parameters.h"
+#include "CUDAUtil.h"
+
+extern "C"
+int InitGPU(int ProcessRank)
+{
+  int NumGPU;
+  CUDA_SAFE_CALL( cudaGetDeviceCount(&NumGPU) );
+  if (NumGPU < 1) {
+    printf("No NVIDIA GPU found. Cannot use CUDA version.\n");
+    exit(1);
+  }
+  CUDA_SAFE_CALL( cudaSetDevice( ProcessRank % NumGPU ) );
+  return SUCCESS;
+}

File src/enzo/CUDAUtil.h

+#ifndef GPU_UTIL_H_
+#define GPU_UTIL_H_
+
+#include "cuda_runtime.h"
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#define CUDA_SAFE_CALL(call) {                                    \
+    cudaError err = call;                                                    \
+    if( cudaSuccess != err) {                                                \
+      fprintf(stderr, "Cuda error in file '%s' in line %i : %s.\n",        \
+              __FILE__, __LINE__, cudaGetErrorString( err) );              \
+      fflush(stderr); \
+      exit(1);                                                  \
+    } }
+
+template <class T>
+class gpu_array
+{
+public:
+  gpu_array() {
+    size_ = 0;
+    ptr_ = NULL;
+  }
+  gpu_array(size_t size) {
+    size_ = size;
+    CUDA_SAFE_CALL( cudaMalloc((void**)&ptr_, size*sizeof(T)) );
+  }
+  gpu_array(T* data, size_t size) {
+    size_ = size;
+    CUDA_SAFE_CALL( cudaMalloc((void**)&ptr_, size*sizeof(T)) );
+    CUDA_SAFE_CALL( cudaMemcpy(ptr_, data, size*sizeof(T), cudaMemcpyHostToDevice) );
+  }
+  ~gpu_array() {
+    if (ptr_ != NULL) {
+      CUDA_SAFE_CALL( cudaFree(ptr_) );
+    }
+  }
+  void malloc(size_t size) {
+    CUDA_SAFE_CALL( cudaMalloc((void**)&ptr_, size*sizeof(T)) );
+  }
+  void free() {
+    if (ptr_ != NULL) {
+      CUDA_SAFE_CALL( cudaFree(ptr_) );
+      ptr_ = NULL;
+    }
+  }
+  void SetPtr(T *ptr) {
+    ptr_ = ptr;
+  }
+  void D2H(T* data) {
+    assert( ptr_ != NULL);
+    CUDA_SAFE_CALL( cudaMemcpy(data, ptr_, size_*sizeof(T), cudaMemcpyDeviceToHost) );
+  }
+  void D2H(T* data, size_t size) {
+    assert( ptr_ != NULL);
+    CUDA_SAFE_CALL( cudaMemcpy(data, ptr_, size*sizeof(T), cudaMemcpyDeviceToHost) );
+  }
+  void H2D(T *data) {
+    assert( ptr_ != NULL);
+    CUDA_SAFE_CALL( cudaMemcpy(ptr_, data, size_*sizeof(T), cudaMemcpyHostToDevice) );
+  }
+  void H2D(T *data, size_t size) {
+    assert( ptr_ != NULL);
+    CUDA_SAFE_CALL( cudaMemcpy(ptr_, data, size*sizeof(T), cudaMemcpyHostToDevice) );
+  }
+  void D2D(T *data) {
+    assert( ptr_ != NULL);
+    CUDA_SAFE_CALL( cudaMemcpy(ptr_, data, size_*sizeof(T), cudaMemcpyDeviceToDevice) );
+  }
+  void D2D(T *data, size_t size) {
+    assert( ptr_ != NULL);
+    CUDA_SAFE_CALL( cudaMemcpy(ptr_, data, size*sizeof(T), cudaMemcpyDeviceToDevice) );
+  }
+  T* ptr() { return ptr_; }
+  size_t size() { return size_; }
+private:
+  T *ptr_;
+  size_t size_;
+};
+
+inline void cumalloc(void **p, size_t size)
+{
+  assert(size > 0);
+  CUDA_SAFE_CALL( cudaMalloc(p, size) );
+}
+
+inline void cufree(void *p)
+{
+  if (p != NULL)
+    CUDA_SAFE_CALL( cudaFree(p) );
+}
+
+inline void cuh2df(float *dst, float *src, size_t size)
+{
+  assert(dst != NULL);
+  assert(src != NULL);
+  assert(size > 0);
+  CUDA_SAFE_CALL( cudaMemcpy(dst, src, size, cudaMemcpyHostToDevice) );
+}
+
+inline void cud2hf(float *dst, float *src, size_t size)
+{
+  assert(dst != NULL);
+  assert(src != NULL);
+  assert(size > 0);
+  CUDA_SAFE_CALL( cudaMemcpy(dst, src, size, cudaMemcpyDeviceToHost) );
+}
+
+class CUDATimer
+{
+protected:
+  cudaEvent_t start, stop;
+public:
+  CUDATimer() {
+    cudaEventCreate(&start);
+    cudaEventCreate(&stop);
+  }
+  ~CUDATimer() {
+    cudaEventDestroy(start);
+    cudaEventDestroy(stop);
+  }
+  void Start() {
+    cudaEventRecord(start, 0);
+  }
+  double GetET() {
+    float gpuTime;
+    cudaEventRecord(stop, 0);
+    cudaEventSynchronize(stop);
+    cudaEventElapsedTime(&gpuTime, start, stop);
+    return (double)1e-3*gpuTime;
+  }
+};
+
+inline void gpu_print_sum(float *a, size_t size)
+{
+  float *tmp = (float*)malloc(size*sizeof(float));
+  cudaMemcpy(tmp, a, size*sizeof(float), cudaMemcpyDeviceToHost);
+  float sum = 0;
+  for (int i = 0; i < size; i++)
+    sum += tmp[i];
+  printf("%f ", sum);
+  free(tmp);
+}
+
+extern "C"
+int InitGPU(int ProcessRank);
+
+#endif

File src/enzo/CallProblemSpecificRoutines.C

     ThisGrid->GridData->SolveOneZoneFreefall();
   }
 
+  /* Add radio-mode jet feedback */
+  if (ClusterSMBHFeedback == TRUE)
+   ThisGrid->GridData->ClusterSMBHFeedback(level);
+
   return SUCCESS;
 }

File src/enzo/ClusterInitialize.C

   int ClusterUseParticles    = FALSE;
   int ClusterUseColour       = FALSE;
   float ClusterInitialTemperature = 1000;
+  float ClusterInitialSpinParameter = 0.05;
   int   ClusterSphereType[MAX_SPHERES];
   float ClusterSphereDensity[MAX_SPHERES],
         ClusterSphereTemperature[MAX_SPHERES],
                   &ClusterUseColour);
     ret += sscanf(line, "ClusterInitialTemperature = %"FSYM, 
                   &ClusterInitialTemperature);
+    ret += sscanf(line, "ClusterInitialSpinParameter = %"FSYM,
+		  &ClusterInitialSpinParameter);
     ret += sscanf(line, "ClusterUniformVelocity = %"FSYM" %"FSYM" %"FSYM, 
                   ClusterUniformVelocity, ClusterUniformVelocity+1,
                   ClusterUniformVelocity+2);
   if (TopGrid.GridData->ClusterInitializeGrid(
              ClusterNumberOfSpheres, ClusterSphereRadius,
              ClusterSphereCoreRadius, ClusterSphereDensity,
-             ClusterSphereTemperature,
+             ClusterSphereTemperature, 
              ClusterSpherePosition, ClusterSphereVelocity,
              ClusterSphereType, ClusterUseParticles,
              ClusterUniformVelocity, ClusterUseColour,
-             ClusterInitialTemperature, 0) == FAIL) {
+             ClusterInitialTemperature, ClusterInitialSpinParameter, 0) == FAIL) {
     fprintf(stderr, "Error in ClusterInitializeGrid.\n");
     return FAIL;
   }
         if (Temp->GridData->ClusterInitializeGrid(
              ClusterNumberOfSpheres, ClusterSphereRadius,
              ClusterSphereCoreRadius, ClusterSphereDensity,
-             ClusterSphereTemperature,
+             ClusterSphereTemperature, 
              ClusterSpherePosition, ClusterSphereVelocity,
              ClusterSphereType, ClusterUseParticles,
              ClusterUniformVelocity, ClusterUseColour,
-             ClusterInitialTemperature, level+1) == FAIL) {
+             ClusterInitialTemperature, ClusterInitialSpinParameter, level+1) == FAIL) {
           fprintf(stderr, "Error in ClusterInitializeGrid.\n");
           return FAIL;
         }
             ClusterUseColour);
     fprintf(Outfptr, "ClusterInitialTemperature = %"GOUTSYM"\n",
             ClusterInitialTemperature);
+    fprintf(fptr, "ClusterInitialSpinParameter   = %"FSYM"\n",
+	    ClusterInitialSpinParameter);
     fprintf(Outfptr, "ClusterUniformVelocity    = %"GOUTSYM" %"GOUTSYM" %"GOUTSYM"\n",
             ClusterUniformVelocity[0], ClusterUniformVelocity[1],
             ClusterUniformVelocity[2]);

File src/enzo/ClusterSMBHSumGasMass.C

+/***********************************************************************
+/
+/  LOOPS OVER GRID AND SUMS COLD GAS COMPONENT
+/
+/  written by: Yuan Li
+/  date:       May 2012
+/  modified1: 
+/
+/  NOTES:
+/
+************************************************************************/
+
+#ifdef USE_MPI
+#include "mpi.h"
+#endif /* USE_MPI */
+ 
+#include <stdlib.h>
+#include <stdio.h>
+#include <unistd.h>
+#include <math.h>
+
+#include "performance.h"
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "Grid.h"
+#include "Hierarchy.h"
+#include "TopGridData.h"
+#include "LevelHierarchy.h"
+#include "CommunicationUtilities.h"
+#include "CosmologyParameters.h"
+#include "phys_constants.h"
+
+int GetUnits(float *DensityUnits, float *LengthUnits,
+             float *TemperatureUnits, float *TimeUnits,
+             float *VelocityUnits, FLOAT Time);
+
+float ClusterSMBHColdGasMass;
+
+int ClusterSMBHSumGasMass(HierarchyEntry *Grids[], int NumberOfGrids, int level)
+{
+
+  /* Return if not running Cluster SMBH feedback */
+  if (ProblemType != 108)
+    return SUCCESS;
+
+  /* Return if we do not want to calculate the cold gas mass */
+  if (ClusterSMBHCalculateGasMass == 0)
+    return SUCCESS;
+
+  /* Return if not on most-refined level. */
+  if (level != MaximumRefinementLevel)
+    return SUCCESS;
+
+  int grid;
+
+  /* Sum over all grids on this processor. */
+
+  ClusterSMBHColdGasMass = 0;
+  for (grid = 0; grid < NumberOfGrids; grid++) {
+//    if (level == MaximumRefinementLevel)
+      Grids[grid]->GridData->ClusterSMBHEachGridGasMass(level);
+  }
+
+  /* Sum over all processors. */
+
+  CommunicationAllSumValues(&ClusterSMBHColdGasMass, 1);
+  FLOAT Time = Grids[0]->GridData->ReturnTime();
+  float DensityUnits = 1.0, LengthUnits = 1.0, TemperatureUnits = 1,
+    TimeUnits = 1.0, VelocityUnits = 1.0;
+  double MassUnits = 1.0;
+
+  if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
+               &TimeUnits, &VelocityUnits, Time) == FAIL) {
+    fprintf(stderr, "Error in GetUnits.\n");
+    return FAIL;
+  }
+  MassUnits = DensityUnits*pow(LengthUnits,3);
+
+  float ColdGasMassMsun=ClusterSMBHColdGasMass*MassUnits/SolarMass;
+  if (ClusterSMBHCalculateGasMass == 2){
+    if (ColdGasMassMsun  > 0.001) {
+       ClusterSMBHFeedbackSwitch = TRUE;
+       ClusterSMBHJetMdot = (ColdGasMassMsun/(ClusterSMBHAccretionTime*1e6))/2.0;  // AccretionTime from Myr to yr; reset Mdot, still in Msun/yr. Devide it by 2 because Mdot is for only one jet.
+       float epsilon=0.001;
+       ClusterSMBHJetEdot = (epsilon*ClusterSMBHJetMdot * SolarMass/3.1557e7) * pow(clight,2)/1.0e44;   //for one jet
+       }
+    else
+       ClusterSMBHFeedbackSwitch = FALSE;    // if there is not enough ColdGas, then do not turn jet on.
+  } // end if ClusterSMBHCalculateGasMass == 2
+  if (ClusterSMBHCalculateGasMass == 1) {
+    int LastClusterSMBHFeedbackSwitch = ClusterSMBHFeedbackSwitch;
+    if (ColdGasMassMsun < 1.0e5)
+       ClusterSMBHFeedbackSwitch = FALSE;
+    else
+       ClusterSMBHFeedbackSwitch = (ColdGasMassMsun <= ClusterSMBHEnoughColdGas && LastClusterSMBHFeedbackSwitch == FALSE) ? FALSE : TRUE;
+
+    if (LastClusterSMBHFeedbackSwitch == FALSE && ClusterSMBHFeedbackSwitch == TRUE) {
+       ClusterSMBHStartTime = Time + ClusterSMBHTramp*1.0e6*3.1557e7/TimeUnits;
+       if (ClusterSMBHJetPrecessionPeriod < 0.00001 & ClusterSMBHJetPrecessionPeriod > -0.00001)  //if precession off (set to 0), change the angle of the jets
+       ClusterSMBHJetAnglePhi += 0.5;
+    if (ClusterSMBHJetPrecessionPeriod < -0.00001)  //if precession negative, change the jet dimension
+       ClusterSMBHJetDim += 1;
+    }
+  } // end if ClusterSMBHCalculateGasMass == 1
+
+  if (MyProcessorNumber == ROOT_PROCESSOR) {
+    FILE *fptr=fopen("MT.out","a");
+    fprintf(fptr,"Time, ClusterSMBHStartTime, Switch, and Total ClusterSMBHColdGasMass in Msun = %g %g %d %g \n", Time, ClusterSMBHStartTime, ClusterSMBHFeedbackSwitch, ColdGasMassMsun);
+    fclose(fptr);
+  }
+  return SUCCESS;
+}

File src/enzo/CommunicationBufferedSend.C

File contents unchanged.

File src/enzo/CommunicationReceiveHandler.C

 	  igrid = CommunicationReceiveArgumentInt[0][index];
 	  isubgrid = CommunicationReceiveArgumentInt[1][index];
 	  SUBling = CommunicationReceiveArgumentInt[2][index];
-#ifdef FLUX_FIX
 	  if ((errcode = grid_two->CorrectForRefinedFluxes
 	      (SubgridFluxesEstimate[igrid][isubgrid], &SubgridFluxesRefined, 
 	       SubgridFluxesEstimate[igrid][NumberOfSubgrids[igrid] - 1],
 	       SUBling, MetaData)) == FAIL) {
 	    ENZO_FAIL("Error in grid->CorrectForRefinedFluxes.\n");
 	  }
-#else
-	  if ((errcode = grid_two->CorrectForRefinedFluxes
-	      (SubgridFluxesEstimate[igrid][isubgrid], &SubgridFluxesRefined, 
-	       SubgridFluxesEstimate[igrid][NumberOfSubgrids[igrid] - 1]     ))
-	      == FAIL) {
-	    ENZO_FAIL("Error in grid->CorrectForRefinedFluxes.\n");
-	  }
-#endif
 	  break;
 
 	case 12:

File src/enzo/CreateSUBlingList.C

 #include "LevelHierarchy.h"
 #include "CommunicationUtilities.h"
 
-#ifdef FLUX_FIX
 int GenerateGridArray(LevelHierarchyEntry *LevelArray[], int level,
 		      HierarchyEntry **Grids[]);
 
   return SUCCESS;
  
 }
-#endif

File src/enzo/DeleteSUBlingList.C

 #include "TopGridData.h"
 #include "LevelHierarchy.h"
  
-#ifdef FLUX_FIX
  
 int DeleteSUBlingList( int NumberOfGrids,
 		      LevelHierarchyEntry **SUBlingList)
   return SUCCESS;
  
 }
-#endif

File src/enzo/DetermineSubgridSizeExtrema.C

File contents unchanged.

File src/enzo/EvolveLevel.C

  
 /* function prototypes */
  
+#ifdef TRANSFER
+#define IMPLICIT_MACRO , ImplicitSolver
+#else
+#define IMPLICIT_MACRO 
+#endif
+
+#define EXTRA_OUTPUT_MACRO(A,B) ExtraOutput(A,LevelArray,MetaData,level,Exterior IMPLICIT_MACRO,B);
+int ExtraOutput(int output_flag, LevelHierarchyEntry *LevelArray[],TopGridData *MetaData, int level, ExternalBoundary *Exterior
+#ifdef TRANSFER
+			  , ImplicitProblemABC *ImplicitSolver
+#endif
+        , char * output_string);
 
 int  RebuildHierarchy(TopGridData *MetaData,
 		      LevelHierarchyEntry *LevelArray[], int level);
 #endif
 #endif
 
-#ifdef FLUX_FIX
 int UpdateFromFinerGrids(int level, HierarchyEntry *Grids[], int NumberOfGrids,
 			 int NumberOfSubgrids[],
 			 fluxes **SubgridFluxesEstimate[],
 			 LevelHierarchyEntry *SUBlingList[],
 			 TopGridData *MetaData);
-#else
-int UpdateFromFinerGrids(int level, HierarchyEntry *Grids[], int NumberOfGrids,
-			 int NumberOfSubgrids[],
-			 fluxes **SubgridFluxesEstimate[]);
-#endif
+
 int CreateFluxes(HierarchyEntry *Grids[],fluxes **SubgridFluxesEstimate[],
 		 int NumberOfGrids,int NumberOfSubgrids[]);		 
 int FinalizeFluxes(HierarchyEntry *Grids[],fluxes **SubgridFluxesEstimate[],
 int ComputeRandomForcingNormalization(LevelHierarchyEntry *LevelArray[],
                                       int level, TopGridData *MetaData,
                                       float * norm, float * pTopGridTimeStep);
+int ClusterSMBHSumGasMass(HierarchyEntry *Grids[], int NumberOfGrids, int level);
 int CreateSiblingList(HierarchyEntry ** Grids, int NumberOfGrids, SiblingGridList *SiblingList, 
 		      int StaticLevelZero,TopGridData * MetaData,int level);
 
-#ifdef FLUX_FIX
 #ifdef FAST_SIB 
 int CreateSUBlingList(TopGridData *MetaData,
 		      LevelHierarchyEntry *LevelArray[], int level,
 #endif /* FAST_SIB */
 int DeleteSUBlingList(int NumberOfGrids,
 		      LevelHierarchyEntry **SUBlingList);
-#endif
 
 int StarParticleInitialize(HierarchyEntry *Grids[], TopGridData *MetaData,
 			   int NumberOfGrids, LevelHierarchyEntry *LevelArray[], 
   int *TotalStarParticleCountPrevious = new int[NumberOfGrids];
   RunEventHooks("EvolveLevelTop", Grids, *MetaData);
 
-#ifdef FLUX_FIX
   /* Create a SUBling list of the subgrids */
   LevelHierarchyEntry **SUBlingList;
-#endif
 
   /* Initialize the chaining mesh used in the FastSiblingLocator. */
 
   /* Loop over grid timesteps until the elapsed time equals the timestep
      from the level above (or loop once for the top level). */
  
+  EXTRA_OUTPUT_MACRO(1, "Before Time Loop")
+
   while ((CheckpointRestart == TRUE)
         || (dtThisLevelSoFar[level] < dtLevelAbove)) {
     if(CheckpointRestart == FALSE) {
     StarParticleInitialize(Grids, MetaData, NumberOfGrids, LevelArray,
 			   level, AllStars, TotalStarParticleCountPrevious);
 
+    /* Calculate ClusterSMBHColdGasMass */
+
+    ClusterSMBHSumGasMass(Grids, NumberOfGrids, level);
+
 #ifdef TRANSFER
     /* Initialize the radiative transfer */
 
     /* For each grid: a) interpolate boundaries from the parent grid.
                       b) copy any overlapping zones from siblings. */
  
+    EXTRA_OUTPUT_MACRO(2,"After SolveHydroEquations grid loop")
 #ifdef FAST_SIB
     SetBoundaryConditions(Grids, NumberOfGrids, SiblingList,
 			  level, MetaData, Exterior, LevelArray[level]);
     SetBoundaryConditions(Grids, NumberOfGrids, level, MetaData,
 			  Exterior, LevelArray[level]);
 #endif
+    EXTRA_OUTPUT_MACRO(25,"After SBC")
 
     /* If cosmology, then compute grav. potential for output if needed. */
 
      *     subgrid's fluxes. (step #19)
      */
  
-#ifdef FLUX_FIX
     SUBlingList = new LevelHierarchyEntry*[NumberOfGrids];
 #ifdef FAST_SIB
     CreateSUBlingList(MetaData, LevelArray, level, SiblingList,
 #else
     CreateSUBlingList(MetaData, LevelArray, level, &SUBlingList);
 #endif /* FAST_SIB */
-#endif /* FLUX_FIX */
 
-#ifdef FLUX_FIX
+
+    EXTRA_OUTPUT_MACRO(3,"Before UFG")
+
     UpdateFromFinerGrids(level, Grids, NumberOfGrids, NumberOfSubgrids,
 			     SubgridFluxesEstimate,SUBlingList,MetaData);
-#else
-    UpdateFromFinerGrids(level, Grids, NumberOfGrids, NumberOfSubgrids,
-			 SubgridFluxesEstimate);
-#endif
 
-#ifdef FLUX_FIX
     DeleteSUBlingList( NumberOfGrids, SUBlingList );
-#endif
 
+    EXTRA_OUTPUT_MACRO(4,"After UFG")
   /* ------------------------------------------------------- */
   /* Add the saved fluxes (in the last subsubgrid entry) to the exterior
      fluxes for this subgrid .
  
   } // end of loop over subcycles
  
+    EXTRA_OUTPUT_MACRO(6, "After Subcycle Loop")
   if (debug)
     fprintf(stdout, "EvolveLevel[%"ISYM"]: NumberOfSubCycles = %"ISYM" (%"ISYM" total)\n", level,
            cycle, LevelCycleCount[level]);

File src/enzo/ExtraOutput.C

+
+//
+// ExtraOutput
+// 
+// Written by: David Collins 
+// date      : May 25, 2011.  3:36 pm.  Cold in my office.
+// 
+// Purpose   : To provide a horrific amount of data.  Calls to ExtraOutput(), thus WriteAllData, 
+//             are sprinkled throughout the code,which helps debug events throughout the evolution.
+//             Output positions are indicated by an integer, and outputs are of the form 
+//             EDXX_YYYY/ExtraOutputXX_YYYY, where XX is the identifier of the dump and YYYY is a counter
+//             for that dump location.  
+//
+//             Output is controlled by the parameter ExtraOutput, a list of integers.  Each integer
+//             corresponds to a location in the code.  Each dump is performed every time the location is reached,
+//             which will result in qute a lot of data if one is not careful.
+//
+//             ExtraOuput is, by design, not re-written into the parameter file when data dumps are made,
+//             to prevent writing too much data.  
+//
+#include "preincludes.h"
+ 
+#ifdef USE_MPI
+#include "mpi.h"
+#endif /* USE_MPI */
+ 
+#include <stdlib.h>
+#include <stdio.h>
+#include <unistd.h>
+#include <math.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 "CommunicationUtilities.h"
+#ifdef TRANSFER
+#include "ImplicitProblemABC.h"
+#endif
+int Group_WriteAllData(char *basename, int filenumber, HierarchyEntry *TopGrid,
+		       TopGridData &MetaData, ExternalBoundary *Exterior,
+#ifdef TRANSFER
+		       ImplicitProblemABC *ImplicitSolver,
+#endif
+		       FLOAT WriteTime = -1, int CheckpointDump = FALSE);
+ int *output_number = NULL;
+int ExtraOutput(int output_flag, LevelHierarchyEntry *LevelArray[],TopGridData *MetaData, int level, ExternalBoundary *Exterior
+#ifdef TRANSFER
+			  , ImplicitProblemABC *ImplicitSolver
+#endif
+       ,char * message ){
+    //initialize output_number array.
+    int n_outputs=100;
+    if( output_number == NULL){
+        output_number = new int[n_outputs];
+        for( int nout=0;nout<n_outputs;nout++){
+            output_number[nout]=0;
+        }
+    }
+
+    //Check for ouptut
+    int WriteOut = FALSE;
+    for( int i=0; i<MAX_EXTRA_OUTPUTS; i++){
+        if( output_flag == ExtraOutputs[i]){
+            WriteOut = TRUE;
+            break;
+        }
+    }
+
+    if( WriteOut ){
+        fflush(stdout);
+        fprintf(stderr,"Extra Output ED%02"ISYM"_%04d %s\n",output_flag,output_number[output_flag],message);
+        LevelHierarchyEntry *Temp2 = LevelArray[0];
+        while (Temp2->NextGridThisLevel != NULL)
+          Temp2 = Temp2->NextGridThisLevel; /* ugh: find last in linked list */
+        //#ifdef USE_HDF5_GROUPS
+        sprintf(MetaData->ExtraDumpName,"Extra%02"ISYM"_",output_flag);
+        sprintf(MetaData->ExtraDumpDir,"ED%02"ISYM"_",output_flag);
+        if (Group_WriteAllData(MetaData->ExtraDumpName, output_number[output_flag]++,
+                   Temp2->GridHierarchyEntry, *MetaData, Exterior,
+#ifdef TRANSFER
+                   ImplicitSolver,
+#endif
+                   -1, FALSE) == FAIL) {
+                ENZO_FAIL("Error in Group_WriteAllData.");
+        }
+        fflush(stdout);
+    }
+  return SUCCESS;
+}

File src/enzo/FluxFix_Grid_CorrectForRefinedFluxes.C

-#ifdef FLUX_FIX
-/***********************************************************************
-/
-/  GRID CLASS (CORRECT SOLUTION GIVEN ORIGINAL AND REFINED FLUXES)
-/
-/  written by: Greg Bryan
-/  date:       November, 1994
-/  modified1:  Robert Harkness
-/  date:       January, 2003
-/	       Include extra fields beyond Metallicity!
-/  modified2:  David Collins & Rick Wagner
-/  date:       May, 2005
-/	       Include flux correction for outside grids.
-/              Re-instated CorrectBoundaryFlux code.
-/  modified2: David Collins, 2005
-/              Updated algebra so Cosmological Expansion is also
-/              conservative.  This fix also came with fixes to euler.src and
-/              Grid_GetProjectedBoundaryFluxes.C, so make sure you get those.
-/
-/  PURPOSE:    Ensures conservation of stuff.
-/
-/  RETURNS: SUCCESS or FAIL 
-/
-************************************************************************/
- 
-// Given both the original and refined fluxes for a subgrid in the current
-//   grid, correct the solution in the cells just outside the solution to
-//   reflect the new fluxes (which are determined from the solution of
-//   the subgrid).
-//   Note that if the subgrid is on the boundary of the current grid, we
-//     do not correct the values but instead replace the boundary fluxes
-//     for the current time step (BoundaryFluxesThisTimeStep).
-//   Also note that subgrids sharing a face with This Grid, but not proper subgrids,
-//   also need to be taken into account.
- 
-#include <stdio.h>
-#include <math.h>
-#include "ErrorExceptions.h"
-#include "macros_and_parameters.h"
-#include "typedefs.h"
-#include "global_data.h"
-#include "Fluxes.h"
-#include "GridList.h"
-#include "ExternalBoundary.h"
-#include "Grid.h"
- 
-/* function prototypes */
- 
-int FindField(int f, int farray[], int n);
-int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt);
-int MakeFieldConservative(field_type field); 
- 
-int grid::CorrectForRefinedFluxes(fluxes *InitialFluxes,
-				  fluxes *RefinedFluxes,
-				  fluxes *BoundaryFluxesThisTimeStep,
-				  int SUBlingGrid,
-				  TopGridData *MetaData)
-{
- 
-  // Return if this doesn't concern us.
- 
-  if (ProcessorNumber != MyProcessorNumber || !UseHydro)
-    return SUCCESS;
- 
-  // declarations
- 
-  int i1, i2, i, j, k, dim, field, ffield, index;
-  int FieldIndex, FluxIndex, GridFluxIndex, Offset, RefinedFluxIndex;
-  int End[MAX_DIMENSION], Start[MAX_DIMENSION];
- 
-  //Dimensions of Initial and Refined fluxes.  ("Dim" should be "InitialDim")
-  int Dim[MAX_DIMENSION],RefinedDim[MAX_DIMENSION] = {0,0,0};
- 
-  //Dimension and Start Index for the ParentGrid (this grid) boundary flux.
-  int GridFluxDim[MAX_DIMENSION], GridFluxStartIndex[MAX_DIMENSION];
- 
-  // Used for calculating position in the RefinedFlux and InitialFlux data structure.
-  int RefinedOffset[MAX_DIMENSION] ={0,0,0}, InitialOffset[MAX_DIMENSION] = {0,0,0};
- 
-  //Internal flags: Correct the BaryonField
-  int CorrectLeftBaryonField, CorrectRightBaryonField;
-  //For correction of the Parent Grid (this grid) boundary flux.
-  int CorrectLeftBoundaryFlux, CorrectRightBoundaryFlux;
-  float CorrectionAmountLeft, CorrectionAmountRight; 
- 
-  long_int GlobalDim;
- 
-  /* If there are no fields, don't do anything. */
- 
-  if (NumberOfBaryonFields > 0) {
- 
-    /* Find fields: density, total energy, velocity1-3. */
- 
-    int DensNum, GENum, Vel1Num, Vel2Num, Vel3Num, TENum, B1Num, B2Num, B3Num;
-    if (this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, 
-					 Vel3Num, TENum, B1Num, B2Num, B3Num) == FAIL) {
-            ENZO_FAIL("Error in grid->IdentifyPhysicalQuantities.");
-    }
-
-    //dcc kludge:  Just remove a(t)? 09/06/05 
-    /* If using comoving coordinates, compute a(t) because we'll need it
-       to multiply against the CellWidth. */
-
-    //  DC revision September 16th 2005 
-    //    FLOAT a = 1, dadt;
-    //    if (ComovingCoordinates)
-    //      if (CosmologyComputeExpansionFactor(Time, &a, &dadt) == FAIL) {
-    //        ENZO_FAIL("Error in CosmologyComputeExpansionFactors.\n");
-    //      }
- 
- 
-    /* Main loop over all faces. */
- 
-    for (dim = 0; dim < GridRank; dim++) {
-      if (GridDimension[dim] > 1) {
-	
-	/* Check that the dims of InitialFluxes & RefinedFluxes are the same */
- 
-	/* don't do this for SUBlings */
-	if( SUBlingGrid == FALSE ){
-	  for (j = 0; j < GridRank; j++)
-	    if ((InitialFluxes->LeftFluxStartGlobalIndex[dim][j] !=
-		 RefinedFluxes->LeftFluxStartGlobalIndex[dim][j])  ||
-		(InitialFluxes->LeftFluxEndGlobalIndex[dim][j] !=
-		 RefinedFluxes->LeftFluxEndGlobalIndex[dim][j])) {
-//	      printf("dim=%d / j=%d //// %d == %d :: %d == %d\n", 
-//		     dim, j, 
-//		     InitialFluxes->LeftFluxStartGlobalIndex[dim][j],
-//		     RefinedFluxes->LeftFluxStartGlobalIndex[dim][j],
-//		     InitialFluxes->LeftFluxEndGlobalIndex[dim][j],
-//		     RefinedFluxes->LeftFluxEndGlobalIndex[dim][j]);
-	      ENZO_FAIL("InitialFluxes & RefinedFluxes are different.\n");
-	    }
-	  /* Error check Fluxes to make sure they all exist. */
-	  for (field = 0; field < NumberOfBaryonFields; field++)
-	    if ((InitialFluxes->LeftFluxes[field][dim] == NULL) ||
-		(RefinedFluxes->LeftFluxes[field][dim] == NULL) ||
-		(InitialFluxes->RightFluxes[field][dim] == NULL) ||
-		(RefinedFluxes->RightFluxes[field][dim] == NULL)) {
-	      fprintf(stderr,"Some Flux data is not present.\n");
-	    return FAIL;
-	  }
-	}
- 
-	//by default, we want to correct the flux.
-	CorrectLeftBaryonField = CorrectRightBaryonField = TRUE;
-	if( SUBlingGrid == TRUE ){
-	
-	  /* calculate Global dimensions on this level */
-	  GlobalDim = 	nlongint(( DomainRightEdge[dim] - DomainLeftEdge[dim])
-				 / CellWidth[dim][0]);
-	
-	  /* get the dims of the refined fluxes to calculate
-	     array indices */
-	
-	  for (i = 0; i < MAX_DIMENSION; i++){
-	    RefinedDim[i] = RefinedFluxes->LeftFluxEndGlobalIndex[dim][i] -
-	      RefinedFluxes->LeftFluxStartGlobalIndex[dim][i] + 1;
-	  }
-	
-	  /* check if SUBling left or right edge lies on
-	     Domain boundary, and if so, modulo the indices,
-	     otherwise, bump the indices to match the initial flux's.
-	  */
-	
-	  if( RefinedFluxes->LeftFluxStartGlobalIndex[dim][dim] == 0 &&
-	      MetaData->LeftFaceBoundaryCondition[dim] == periodic ){
-	    RefinedFluxes->LeftFluxStartGlobalIndex[dim][dim] = GlobalDim - 1;
-	    RefinedFluxes->LeftFluxEndGlobalIndex[dim][dim] = GlobalDim - 1;
-	  }else{
-	    RefinedFluxes->LeftFluxStartGlobalIndex[dim][dim]--;
-	    RefinedFluxes->LeftFluxEndGlobalIndex[dim][dim]--;
-	  }
-	
-	
-	  if( RefinedFluxes->RightFluxStartGlobalIndex[dim][dim] == GlobalDim - 1 &&
-	      MetaData->RightFaceBoundaryCondition[dim] == periodic){
-	    RefinedFluxes->RightFluxStartGlobalIndex[dim][dim] = 0;
-	    RefinedFluxes->RightFluxEndGlobalIndex[dim][dim] = 0;
-	  }else{
-	    RefinedFluxes->RightFluxStartGlobalIndex[dim][dim]++;
-	    RefinedFluxes->RightFluxEndGlobalIndex[dim][dim]++;
-	  }
-	
-	  /* check to see if we're doing this dimension at all.
-	     only the dimension of contact needs to be checked,
-	     since SUBling grids can only have contact along a
-	     single axis. any corrections to this statement
-	     earns a beer */
-	
-	
-	  //note these are integers, so comparing them directly is ok.
-	  //Also note that here we don't do the complete check for SUBling-ness.
-	  //More logic is necessary for any two arbitrary grids, but it's been done
-	  //already in populating the SUBling list.  Here, all we need
-	  //to do is determine which face needs the correction, as we already know one exists.
-	
-	  if( InitialFluxes->RightFluxStartGlobalIndex[dim][dim] ==
-	      RefinedFluxes->LeftFluxStartGlobalIndex[dim][dim] ){
-	    CorrectRightBaryonField = TRUE;
-	  }else{
-	    CorrectRightBaryonField = FALSE;
-	  }
-	
-	  if( InitialFluxes->LeftFluxStartGlobalIndex[dim][dim]==
-	      RefinedFluxes->RightFluxStartGlobalIndex[dim][dim] ){
-	    CorrectLeftBaryonField = TRUE;
-	  }else{
-	    CorrectLeftBaryonField = FALSE;
-	  }
-	
-	  for (i = 0; i < MAX_DIMENSION; i++) {
-	    /* calculate the offset, so the index of the refined fluxes can
-	       be determined from the grid's index */
-	    RefinedOffset[i] = max(InitialFluxes->LeftFluxStartGlobalIndex[dim][i]-
-				   RefinedFluxes->LeftFluxStartGlobalIndex[dim][i],0);
-	  }
- 
-	  RefinedFluxes->LeftFluxStartGlobalIndex[dim][dim]=0;
-	  RefinedOffset[dim]=0;
- 
-	}//Subling == TRUE	
- 
-	if( CorrectLeftBaryonField || CorrectRightBaryonField){
-	
-	  /* Compute Start and end indicies of flux region (with respect to
-	     the current grid's flux region). */
-	
-	  for (i = 0; i < MAX_DIMENSION; i++) {
-	    Start[i] = 0;
-	    End[i] = 0;
-	  }
-	
-	  /* start index = subgrid flux left edge global index -
-	     grid far left edge global index
-	     end index = subgrid flux right edge -
-	     grid far left edge global index. */
-	
-	  /* modified to account for different dimensions to the
-	     initial and refined fluxes */
-	
-	  for (i = 0; i < GridRank; i++) {
-	    Start[i] = max(InitialFluxes->LeftFluxStartGlobalIndex[dim][i],
-			   RefinedFluxes->LeftFluxStartGlobalIndex[dim][i]) -
-	      nlongint((CellLeftEdge[i][0] - DomainLeftEdge[i])/
-		       CellWidth[i][0]);
-	    End[i] = min(InitialFluxes->LeftFluxEndGlobalIndex[dim][i],
-			 RefinedFluxes->LeftFluxEndGlobalIndex[dim][i]) -
-	      nlongint((CellLeftEdge[i][0] - DomainLeftEdge[i])/
-		       CellWidth[i][0]);
-	
-	
-	    if (Start[i] < 0 || End[i] > GridDimension[i]) {
-	      fprintf(stderr, "Start/End[%"ISYM"] = %"ISYM"/%"ISYM"\n",
-		      dim, Start[i], End[i]);
-	      fprintf(stderr, "%"GOUTSYM" %"GOUTSYM" %lld\n",
-		      CellLeftEdge[i][0], CellWidth[i][0],
-		      InitialFluxes->LeftFluxStartGlobalIndex[dim][i]);
-	      ENZO_FAIL("Error in FluxFix_Grid_CorrectForRefinedFluxes!\n");
-	    }
-	  }
-	
-	  /* Correct vector to point at cells just outside the left face.
-	     Start[dim] and End[dim] should be the same because the
-	     layer to be corrected is but one cell thick. */
-	
-	  Start[dim] = max(Start[dim] - 1, 0);
-	  End[dim]   = Start[dim];
- 
-	  /* Compute Dimensions of InitialFluxes */
-	
-	  for (i = 0; i < MAX_DIMENSION; i++)
-	    Dim[i] = End[i] - Start[i] + 1;
-	
-	  /* Compute Offset (in baryon field) for right side of face.
-	     The +2 is there because we want to correct the cells just the
-	     right face.*/
-	
-	  Offset = InitialFluxes->RightFluxStartGlobalIndex[dim][dim] -
-	    InitialFluxes->LeftFluxStartGlobalIndex[dim][dim] + 2;
-	  Offset = min(Offset, GridDimension[dim]-1);  // this isn't needed (?)
- 
- 
-	  //For SUBling grids, alter Offset, Start, and End to reflect that we're
-	  //adjusting the INNER edge of the grid if the SUBgrid is outside of it.
- 
-	  if( SUBlingGrid ){
-	    Offset -= 2;
-	    Start[dim]++;
-	    End[dim]++;
-	
-	    //Also correct Dim, the size of the Initial Flux: it comes from This Grid,
-	    //not the Subgrid.
-	
-	    for(i=0;i<GridRank;i++)
-	      if(i != dim){
-		Dim[i] = GridEndIndex[i]-GridStartIndex[i]+1;
-		InitialOffset[i] = max( RefinedFluxes->LeftFluxStartGlobalIndex[dim][i]-
-					InitialFluxes->LeftFluxStartGlobalIndex[dim][i],
-					0);
-	      }else{
-		Dim[i] = 1;
-		InitialOffset[i] = 0;
-	      }
-	  }
-	
-	
-	  /* Check to see if we should correct BoundaryFluxesThisTimeStep
-	     instead of the fields themselves. */
-	
-	  CorrectLeftBoundaryFlux = FALSE;
-	  CorrectRightBoundaryFlux = FALSE;
-	
-	  if (Start[dim] == GridStartIndex[dim]-1){
-	    CorrectLeftBoundaryFlux = TRUE;
-	    CorrectLeftBaryonField  = FALSE;
-	  }
-	  if (Start[dim] + Offset == GridEndIndex[dim]+1){
-	    CorrectRightBoundaryFlux = TRUE;
-	    CorrectRightBaryonField  = FALSE;
-	  }
- 
-	  /* Set GridFluxStartIndex to the starting position of the flux
-	     plane (i.e. exclude grid boundary zones), except for the direction
-	     of the flux is set such that GridStartIndex[dim] - Start[dim] = 0 */
-	
-	  for (i = 0; i < MAX_DIMENSION; i++) {
-	    GridFluxStartIndex[i] = GridStartIndex[i];
-	    GridFluxDim[i] = GridEndIndex[i] - GridStartIndex[i] + 1;
-	  }
-	
-	  GridFluxStartIndex[dim] = Start[dim];
-	  GridFluxDim[dim] = 1;
-	
-	  /* Turn Offset (in dim direction) to Offset (in field array) */
-	
-	  for (i = 0; i < dim; i++)
-	    Offset *= GridDimension[i];
-	
-	  /* Multiply faces by density to get conserved quantities
-	     (only multiply fields which we are going to correct) */
-	
-	  if (HydroMethod != Zeus_Hydro) {
-
-	    for (field = 0; field < NumberOfBaryonFields; field++) 
-	      if ( MakeFieldConservative(FieldType[field]) ) {
-		for (k = Start[2]; k <= End[2]; k++) {
-		  for (j = Start[1]; j <= End[1]; j++) {
-		    index = (k*GridDimension[1] + j)*GridDimension[0] + Start[0];
-		    for (i = Start[0]; i <= End[0]; i++, index++) {
-		      BaryonField[field][index] *= BaryonField[DensNum][index];
-		      BaryonField[field][index+Offset] *= 
-			BaryonField[DensNum][index+Offset];
-		    }
-		  }
-		}
-	      }
-	  }
-
-	  
-	
-	  /* Divide species by densities so that at the end we can multiply
-	     them by the new density (species are not otherwise modified --
-	     see the next comment).  This ensures that the species are changed
-	     to keep the same fractional density. */
-	
-	  for (field = 0; field < NumberOfBaryonFields; field++)
-	    if (FieldType[field] >= ElectronDensity &&
-		FieldType[field] < Metallicity &&
-		FieldTypeNoInterpolate(FieldType[field]) == FALSE &&
-		FieldTypeIsRadiation(FieldType[field]) == FALSE)
-	      for (k = Start[2]; k <= End[2]; k++)
-		for (j = Start[1]; j <= End[1]; j++) {
-		  index = (k*GridDimension[1] + j)*GridDimension[0] + Start[0];
-		  for (i = Start[0]; i <= End[0]; i++, index++) {
-		    BaryonField[field][index] /= BaryonField[DensNum][index];
-		    BaryonField[field][index+Offset] /=
-		      BaryonField[DensNum][index+Offset];
-		  }
-		}
-	
-	  /* Correct face for difference between refined and initial fluxes.
-	     (Don't do this for energy if radiative cooling is on because it's
-	     no longer conserved.  Similarly, don't do it for the species
-	     because they are not individually conserved either -- in fact,
-	     this could be done for the conserved quantities like charge,
-	     total number density summed over ionization, etc.) */
-	
-	  if (Coordinate == Cartesian) {
-	  for (field = 0; field < NumberOfBaryonFields; field++){
-	    if ((FieldTypeNoInterpolate(FieldType[field]) == FALSE) &&
-		(RadiativeCooling == 0 || (FieldType[field] != TotalEnergy &&
-					   FieldType[field] != InternalEnergy)) &&
-		(FieldType[field] < ElectronDensity) && 
-		FieldType[field] != DrivingField1 &&
-		FieldType[field] != DrivingField2 &&
-		FieldType[field] != DrivingField3 &&
-		FieldType[field] != GravPotential &&
-		FieldType[field] != DebugField) {
-	      for (k = Start[2]; k <= End[2]; k++){
-		for (j = Start[1]; j <= End[1]; j++){
-		  for (i = Start[0]; i <= End[0]; i++) {
-		
-		    /* Compute indexes. */
-		
-		    FieldIndex = (k*GridDimension[1] + j)*GridDimension[0] + i;
-		    FluxIndex  = ((k - Start[2]+InitialOffset[2])*Dim[1] +
-				  (j - Start[1]+InitialOffset[1]))*Dim[0] +
-		      (i - Start[0]+InitialOffset[0]);
-		
-		    if( SUBlingGrid ){
-		      RefinedFluxIndex = ((k - Start[2] + RefinedOffset[2])*RefinedDim[1] +
-					  (j - Start[1] + RefinedOffset[1]))*RefinedDim[0] +
-			(i - Start[0] + RefinedOffset[0]);
-		
-		    }else{
-		      RefinedFluxIndex = FluxIndex;
-		    }
-		
-		    GridFluxIndex =
-		      (i - GridFluxStartIndex[0])
-		      + (j - GridFluxStartIndex[1])*GridFluxDim[0]
-		      + (k - GridFluxStartIndex[2])*GridFluxDim[1]*GridFluxDim[0];
-		
-		
-// 		    if (CorrectLeftBoundaryFlux)
-//		      BoundaryFluxesThisTimeStep->LeftFluxes[field][dim][GridFluxIndex] =
-//			RefinedFluxes->LeftFluxes[field][dim][FluxIndex];
-		
-		    if(CorrectLeftBaryonField){
-		
-		      if( SUBlingGrid == FALSE ){
-			CorrectionAmountLeft = 
-			  InitialFluxes->LeftFluxes[field][dim][FluxIndex] -
-			  RefinedFluxes->LeftFluxes[field][dim][FluxIndex];
-			BaryonField[field][FieldIndex] += CorrectionAmountLeft;
-			
-		      }else{ /* if( SUBlingGrid == False) */
-			
-			CorrectionAmountLeft = 
-			  -(InitialFluxes->LeftFluxes[field][dim][FluxIndex] -
-			    RefinedFluxes->RightFluxes[field][dim][RefinedFluxIndex]);
-			BaryonField[field][FieldIndex] += CorrectionAmountLeft;
-			
-		      }
-		    }
-		
-// 		    if (CorrectRightBoundaryFlux)
-//		      BoundaryFluxesThisTimeStep->RightFluxes[field][dim] [GridFluxIndex] =
-//			RefinedFluxes->RightFluxes[field][dim][FluxIndex];
-		
-		    /* update only if necessary */
-		    if(CorrectRightBaryonField){
-		
-		      if( SUBlingGrid == FALSE ){
-
-			CorrectionAmountRight = 
-			  -(InitialFluxes->RightFluxes[field][dim][FluxIndex] -
-			    RefinedFluxes->RightFluxes[field][dim][FluxIndex]);
-			  
-			BaryonField[field][FieldIndex + Offset] += CorrectionAmountRight;
-			
-		      }else{ /* if( SUBlingGrid == FALSE ){ */
-			  CorrectionAmountRight = 
-			    InitialFluxes->RightFluxes[field][dim][FluxIndex] -
-			    RefinedFluxes->LeftFluxes[field][dim][RefinedFluxIndex];
-			  BaryonField[field][FieldIndex + Offset] += CorrectionAmountRight;
-			
-		      } // else{ /* if( SUBlingGrid == FALSE ){ */
-		    } // if(CorrectRightBaryonField)
-		
-		    if ((FieldTypeIsDensity(FieldType[field]) == TRUE ||
-			 FieldType[field] == TotalEnergy ||
-			 FieldType[field] == InternalEnergy)) {
-
-		      /* If new density & energy is < 0 then undo the
-			 flux correction. */
-
-		      if (CorrectLeftBaryonField &&
-			  BaryonField[field][FieldIndex] <= 0) {
-			BaryonField[field][FieldIndex] -= CorrectionAmountLeft;
-
-			if (SUBlingGrid == FALSE) {
-			  if (debug)
-			    printf("P(%d) -- CFRFl warn: %e %e %e %e %"ISYM
-				   " %"ISYM" %"ISYM" %"ISYM" [%"ISYM"]\n",
-				   MyProcessorNumber, BaryonField[field][FieldIndex],
-				   InitialFluxes->LeftFluxes[field][dim][FluxIndex],
-				   RefinedFluxes->LeftFluxes[field][dim][FluxIndex],
-				   CorrectionAmountLeft,
-				   i, j, k, dim, field);
-			  for (ffield = 0; ffield < NumberOfBaryonFields; ffield++)
-			    RefinedFluxes->LeftFluxes[ffield][dim][FluxIndex] =
-			      InitialFluxes->LeftFluxes[ffield][dim][FluxIndex];
-			} else {
-			  if (debug)
-			    printf("P(%d) -- CFRFlS warn: %e %e %e %e %"ISYM
-				   " %"ISYM" %"ISYM" %"ISYM" [%"ISYM"]\n",
-				   MyProcessorNumber, BaryonField[field][FieldIndex],
-				   InitialFluxes->LeftFluxes[field][dim][FluxIndex],
-				   RefinedFluxes->RightFluxes[field][dim][RefinedFluxIndex],
-				   CorrectionAmountLeft,
-				   i, j, k, dim, field);
-			  for (ffield = 0; ffield < NumberOfBaryonFields; ffield++)
-			    RefinedFluxes->RightFluxes[ffield][dim][RefinedFluxIndex] =
-			      InitialFluxes->LeftFluxes[ffield][dim][FluxIndex];
-			} // ENDELSE (SUBlingGrid == FALSE)
-		      } // ENDIF CorrectLeftBaryonField
-
-		      if (CorrectRightBaryonField &&
-			  BaryonField[field][FieldIndex+Offset] <= 0) {
-			BaryonField[field][FieldIndex+Offset] -= CorrectionAmountRight;
-
-			if (SUBlingGrid == FALSE) {
-			  if (debug)
-			    printf("P(%d) -- CFRFr warn: %e %e %e %e %"ISYM
-				   " %"ISYM" %"ISYM" %"ISYM" [%"ISYM"]\n",
-				   MyProcessorNumber, BaryonField[field][FieldIndex],
-				   InitialFluxes->RightFluxes[field][dim][FluxIndex],
-				   RefinedFluxes->RightFluxes[field][dim][FluxIndex],
-				   CorrectionAmountRight,
-				   i, j, k, dim, field);
-			  for (ffield = 0; ffield < NumberOfBaryonFields; ffield++)
-			    RefinedFluxes->RightFluxes[ffield][dim][FluxIndex] =
-			      InitialFluxes->RightFluxes[ffield][dim][FluxIndex];
-			} else {
-			  if (debug)
-			    printf("P(%d) -- CFRFrS warn: %e %e %e %e %"ISYM
-				   " %"ISYM" %"ISYM" %"ISYM" [%"ISYM"]\n",
-				   MyProcessorNumber, BaryonField[field][FieldIndex],
-				   InitialFluxes->LeftFluxes[field][dim][FluxIndex],
-				   RefinedFluxes->RightFluxes[field][dim][RefinedFluxIndex],
-				   CorrectionAmountRight,
-				   i, j, k, dim, field);
-			  for (ffield = 0; ffield < NumberOfBaryonFields; ffield++)
-			    RefinedFluxes->LeftFluxes[ffield][dim][RefinedFluxIndex] =
-			      InitialFluxes->RightFluxes[ffield][dim][FluxIndex];
-			}
-		      } // ENDIF CorrectRightBaryonField
-		    }
-		  }// for (i = Start[0]; i <= End[0]; i++) {
-		} // for (j = Start[1]; j <= End[1]; j++){
-	      } // for (k = Start[2]; k <= End[2]; k++){
-	    }	// if ((RadiativeCooling == 0 || (FieldType[field] != TotalEnergy && etc
-	  } // for (field = 0; field < NumberOfBaryonFields; field++){
-	  }
-
-
-	if (Coordinate == Cylindrical) {
-	FLOAT xr, xl, xc, geofacr, geofacl;
-	for (field = 0; field < NumberOfBaryonFields; field++) {
-	  /*if ((RadiativeCooling == 0 || (FieldType[field] != TotalEnergy && 
-	    FieldType[field] != InternalEnergy))
-	    && (FieldType[field] < ElectronDensity) ) {*/
-	  if (FieldType[field] < ElectronDensity &&
-	      FieldType[field] != DrivingField1 &&
-	      FieldType[field] != DrivingField2 &&
-	      FieldType[field] != DrivingField3 &&
-	      FieldType[field] != GravPotential &&
-	      FieldType[field] != DebugField
-	      ) {
-	  for (k = Start[2]; k <= End[2]; k++) {
-	    for (j = Start[1]; j <= End[1]; j++) {
-	      for (i = Start[0]; i <= End[0]; i++) {
-		/* Compute indexes. */
-		FieldIndex = (k*GridDimension[1] + j)*GridDimension[0] + i;
-		FluxIndex  = ((k - Start[2])*Dim[1] + (j - Start[1]))*Dim[0] +
-		              (i - Start[0]);
-		GridFluxIndex = 
-                     (i - GridFluxStartIndex[0]) 
-		   + (j - GridFluxStartIndex[1])*GridFluxDim[0]
-		   + (k - GridFluxStartIndex[2])*GridFluxDim[1]*GridFluxDim[0];
-		if (dim == 0) {
-		  /* Left side */
-		  xr = CellLeftEdge[0][i] + CellWidth[0][i];
-                  xc = CellLeftEdge[0][i] + 0.5*CellWidth[0][i];
-                  geofacr = xr/xc;
-		  BaryonField[field][FieldIndex] +=
-		    (InitialFluxes->LeftFluxes[field][dim][FluxIndex] -
-		     RefinedFluxes->LeftFluxes[field][dim][FluxIndex] )*geofacr;
-
-		  /* Right side */
-		  xl = CellLeftEdge[0][i+Offset];
-                  xc = xl + 0.5*CellWidth[0][i+Offset];
-                  geofacl = xl/xc;
-		  BaryonField[field][FieldIndex + Offset] -=
-		    (InitialFluxes->RightFluxes[field][dim][FluxIndex] -
-		     RefinedFluxes->RightFluxes[field][dim][FluxIndex] )*geofacl;
-		}		
-		if (dim == 1) {
-		  /* Left side */
-		  BaryonField[field][FieldIndex] +=
-		    (InitialFluxes->LeftFluxes[field][dim][FluxIndex] -
-		     RefinedFluxes->LeftFluxes[field][dim][FluxIndex] );
-		  /* Right side */
-		  BaryonField[field][FieldIndex + Offset] -=
-		    (InitialFluxes->RightFluxes[field][dim][FluxIndex] -
-		     RefinedFluxes->RightFluxes[field][dim][FluxIndex] );
-		}		
-		if (dim == 2) {
-
-		  /* Left side */
-		  xc = CellLeftEdge[0][i] + 0.5*CellWidth[0][i]; 
-		  BaryonField[field][FieldIndex] +=
-		    (InitialFluxes->LeftFluxes[field][dim][FluxIndex] -
-		     RefinedFluxes->LeftFluxes[field][dim][FluxIndex] )*xc;
-		  /* Right side */
-		  xl = CellLeftEdge[0][i+Offset];
-                  xc = xl + 0.5*CellWidth[0][i+Offset];
-		  BaryonField[field][FieldIndex + Offset] -=
-		    (InitialFluxes->RightFluxes[field][dim][FluxIndex] -
-		     RefinedFluxes->RightFluxes[field][dim][FluxIndex] )*xc;
-
-		}		
-		/* Check for posivtivity and undo flux correction if negative */
-		if ((FieldType[field] == Density || 
-		     FieldType[field] == TotalEnergy ||
-		     FieldType[field] == InternalEnergy) &&
-		    BaryonField[field][FieldIndex] <= 0) {
-		  /*if (debug) {
-		    printf("CFRFl warn: %e %e %e %d %d %d %d [%d]\n",
-			   BaryonField[field][FieldIndex],
-			   InitialFluxes->LeftFluxes[field][dim][FluxIndex],
-			   RefinedFluxes->LeftFluxes[field][dim][FluxIndex],
-			   i, j, k, dim, field);
-			   }*/
-		  /* If new density & energy is < 0 then undo the flux correction. */
-		  BaryonField[field][FieldIndex] -=
-		     (InitialFluxes->LeftFluxes[field][dim][FluxIndex] -
-		      RefinedFluxes->LeftFluxes[field][dim][FluxIndex] );
-		  for (ffield = 0; ffield < NumberOfBaryonFields; ffield++)
-		    RefinedFluxes->LeftFluxes[ffield][dim][FluxIndex] =
-		      InitialFluxes->LeftFluxes[ffield][dim][FluxIndex];
-		}
-		if ((FieldType[field] == Density || 
-		     FieldType[field] == TotalEnergy ||
-		     FieldType[field] == InternalEnergy) &&
-		    BaryonField[field][FieldIndex + Offset] <= 0.0) {
-		  /*if (debug) {
-		    printf("CFRFr warn: %e %e %e %d %d %d %d (%d) [%d]\n",
-			   BaryonField[field][FieldIndex + Offset],
-			   InitialFluxes->RightFluxes[field][dim][FluxIndex],
-			   RefinedFluxes->RightFluxes[field][dim][FluxIndex],
-			   i, j, k, dim, Offset, field);
-			   }*/
-		  /* If new density & energy is < 0 then undo the flux correction. */
-		  BaryonField[field][FieldIndex + Offset] +=