Commits

John Wise committed c598e84

Imported my old C++ euler wrappers. Less overhead and the coding is
more clear than the current wrappers.

Comments (0)

Files changed (8)

 
 /* PPM Direct Euler Solver. */
 
-  int PPMDirectEuler(int CycleNumber, int NumberOfSubgrids, 
-                     fluxes *SubgridFluxes[], float *CellWidthTemp[], 
-                     long_int GridGlobalStart[], int GravityOn,
-		     int NumberOfColours, int colnum[]);
-
-  int euler_sweep(int dim, int iter, int CycleNumber, int NumberOfSubgrids, 
-                  fluxes *SubgridFluxes[], float *CellWidthTemp[],
-                  long_int GridGlobalStart[], int GravityOn,
-		  int NumberOfColours, int colnum[], float *temp, int tempsize);
+int SolvePPM_DE(int CycleNumber, int NumberOfSubgrids, 
+		fluxes *SubgridFluxes[], float *CellWidthTemp[], 
+		Elong_int GridGlobalStart[], int GravityOn, 
+		int NumberOfColours, int colnum[]);
+
+int xEulerSweep(int k, int NumberOfSubgrids, fluxes *SubgridFluxes[], 
+		Elong_int GridGlobalStart[], float *CellWidthTemp[], 
+		int GravityOn, int NumberOfColours, int colnum[]);
+
+int yEulerSweep(int i, int NumberOfSubgrids, fluxes *SubgridFluxes[], 
+		Elong_int GridGlobalStart[], float *CellWidthTemp[], 
+		int GravityOn, int NumberOfColours, int colnum[]);
+
+int zEulerSweep(int j, int NumberOfSubgrids, fluxes *SubgridFluxes[], 
+		Elong_int GridGlobalStart[], float *CellWidthTemp[], 
+		int GravityOn, int NumberOfColours, int colnum[]);
 
 // AccelerationHack
 

src/enzo/Grid_SolveHydroEquations.C

 
     /* initialize */
 
-    int dim, i, j, field, size, subgrid, n, colnum[MAX_COLOR];   // MAX_COLOR is defined in fortran.def
+    // MAX_COLOR is defined in fortran.def
+    int dim, i, j, field, size, subgrid, n, colnum[MAX_COLOR];
     Elong_int GridGlobalStart[MAX_DIMENSION];
     FLOAT a = 1, dadt;
 
     /* note: Start/EndIndex are zero based */
         
     if (HydroMethod == PPM_DirectEuler)
-      this->PPMDirectEuler(CycleNumber, NumberOfSubgrids, 
-			   SubgridFluxes, CellWidthTemp, 
-			   GridGlobalStart, GravityOn,
-			   NumberOfColours, colnum);
+      this->SolvePPM_DE(CycleNumber, NumberOfSubgrids, SubgridFluxes, 
+			CellWidthTemp, GridGlobalStart, GravityOn, 
+			NumberOfColours, colnum);
 
     /* PPM LR has been withdrawn. */
 

src/enzo/Grid_SolvePPM_DE.C

+/***********************************************************************
+/
+/  GRID CLASS (WRAPPER FOR EULER SOLVER)
+/
+/  written by: John Wise
+/  date:       May, 2007
+/  modified1:
+/
+/  PURPOSE:
+/
+/  RETURNS:
+/    SUCCESS or FAIL
+/
+************************************************************************/
+
+// Solve the hydro equations with the solver, saving the subgrid fluxes
+//
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "Grid.h"
+//#include "fortran.def"
+
+int grid::SolvePPM_DE(int CycleNumber, int NumberOfSubgrids, 
+		      fluxes *SubgridFluxes[], float *CellWidthTemp[], 
+		      Elong_int GridGlobalStart[], int GravityOn, 
+		      int NumberOfColours, int colnum[])
+{
+
+  int nxz, nyz, nzz, ixyz;
+
+  nxz = GridEndIndex[0] - GridStartIndex[0] + 1;
+  nyz = GridEndIndex[1] - GridStartIndex[1] + 1;
+  nzz = GridEndIndex[2] - GridStartIndex[2] + 1;
+
+  ixyz = CycleNumber % 3;
+
+  int i,j,k,n;
+  for (n = ixyz; n <= ixyz+2; n++) {
+
+    // Update in x-direction
+    if ((n % 3 == 0) && nxz > 1) {
+      for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) {
+	if (this->xEulerSweep(k, NumberOfSubgrids, SubgridFluxes, 
+			      GridGlobalStart, CellWidthTemp, GravityOn, 
+			      NumberOfColours, colnum) == FAIL) {
+	  fprintf(stderr, "Error in xEulerSweep.  k = %d\n", k);
+	  ENZO_FAIL("");
+	}
+      } // ENDFOR k
+    } // ENDIF x-direction
+
+    // Update in y-direction
+    if ((n % 3 == 1) && nyz > 1) {
+      for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++) {
+	if (this->yEulerSweep(i, NumberOfSubgrids, SubgridFluxes, 
+			      GridGlobalStart, CellWidthTemp, GravityOn, 
+			      NumberOfColours, colnum) == FAIL) {
+	  fprintf(stderr, "Error in yEulerSweep.  i = %d\n", i);
+	  ENZO_FAIL("");
+	}
+      } // ENDFOR i
+    } // ENDIF y-direction
+
+    // Update in z-direction
+    if ((n % 3 == 2) && nzz > 1) {
+      for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) {
+	if (this->zEulerSweep(j, NumberOfSubgrids, SubgridFluxes, 
+			      GridGlobalStart, CellWidthTemp, GravityOn, 
+			      NumberOfColours, colnum) == FAIL) {
+	  fprintf(stderr, "Error in zEulerSweep.  j = %d\n", j);
+	  ENZO_FAIL("");
+	}
+      } // ENDFOR j
+    } // ENDIF z-direction
+
+  } // ENDFOR n
+
+  return SUCCESS;
+
+}

src/enzo/Grid_xEulerSweep.C

+/***********************************************************************
+/
+/  GRID CLASS (WRAPPER FOR EULERIAN PPM SOLVER)
+/
+/  written by: John H. Wise
+/  date:       May 2007
+/  modified1:
+/
+/  PURPOSE:
+/
+/  RETURNS:
+/    SUCCESS or FAIL
+/
+************************************************************************/
+
+// Solve the hydro equations with the solver, saving the subgrid fluxes
+//
+
+#include <stdlib.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 "Grid.h"
+#include "euler_sweep.h"
+//#include "fortran.def"
+
+int grid::xEulerSweep(int k, int NumberOfSubgrids, fluxes *SubgridFluxes[], 
+		      Elong_int GridGlobalStart[], float *CellWidthTemp[], 
+		      int GravityOn, int NumberOfColours, int colnum[])
+{
+
+  int dim = 0, idim = 1, jdim = 2;
+  int ierr = 0;
+
+  /* Find fields: density, total energy, velocity1-3. */
+
+  int DensNum, GENum, Vel1Num, Vel2Num, Vel3Num, TENum;
+  
+  this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, 
+				   Vel3Num, TENum);
+
+  int nxz, nyz, nzz, ixyz;
+
+  nxz = GridEndIndex[0] - GridStartIndex[0] + 1;
+  nyz = GridEndIndex[1] - GridStartIndex[1] + 1;
+  nzz = GridEndIndex[2] - GridStartIndex[2] + 1;
+
+  float MinimumPressure = tiny_number;
+  
+  // Copy from field to slice
+
+  float *dslice, *eslice, *uslice, *vslice, *wslice, *grslice, *geslice, 
+    *colslice, *pslice;
+
+  int size = GridDimension[0] * GridDimension[1];
+  dslice = new float[size];
+  eslice = new float[size];
+  uslice = new float[size];
+  vslice = new float[size];
+  wslice = new float[size];
+  pslice = new float[size];
+  if (GravityOn) {
+    grslice = new float[size];  
+  }
+  if (DualEnergyFormalism) {
+    geslice = new float[size];  
+  }
+  if (NumberOfColours > 0) {
+    colslice = new float[NumberOfColours * size];  
+  }
+
+  int i, j, n, ncolour, index2, index3;
+  for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) {
+
+    index2 = j * GridDimension[0];
+
+    for (i = 0; i < GridDimension[0]; i++) {
+      index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+      dslice[index2+i] = BaryonField[DensNum][index3];
+      eslice[index2+i] = BaryonField[TENum][index3];
+      uslice[index2+i] = BaryonField[Vel1Num][index3];
+    } // ENDFOR i
+
+    // Set velocities to zero if rank < 3 since hydro routines are
+    // hard-coded for 3-d
+
+    if (GridRank > 1) 
+      for (i = 0; i < GridDimension[0]; i++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	vslice[index2+i] = BaryonField[Vel2Num][index3];
+      }
+    else
+      for (i = 0; i < GridDimension[0]; i++)
+	vslice[index2+i] = 0;
+  
+    if (GridRank > 2)
+      for (i = 0; i < GridDimension[0]; i++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	wslice[index2+i] = BaryonField[Vel3Num][index3];
+      }
+    else
+      for (i = 0; i < GridDimension[0]; i++)
+	wslice[index2+i] = 0;
+    
+    if (GravityOn)
+      for (i = 0; i < GridDimension[0]; i++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	grslice[index2+i] = AccelerationField[dim][index3];
+      }
+
+    if (DualEnergyFormalism)
+      for (i = 0; i < GridDimension[0]; i++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	geslice[index2+i] = BaryonField[GENum][index3];
+      }
+
+    for (n = 0; n < NumberOfColours; n++) {
+      index2 = (n*GridDimension[1] + j) * GridDimension[0];
+      for (i = 0; i < GridDimension[0]; i++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	colslice[index2+i] = BaryonField[colnum[n]][index3];
+      }
+    } // ENDFOR colours
+  } // ENDFOR j
+
+  /* Allocate memory for fluxes */
+
+  float *dls, *drs, *flatten, *pbar, *pls, *prs, *ubar, *uls, *urs, *vls, 
+    *vrs, *gels, *gers, *wls, *wrs, *diffcoef, *df, *ef, *uf, *vf, *wf, *gef,
+    *colf, *colls, *colrs;
+
+  dls = new float[size];	
+  drs = new float[size];	
+  flatten = new float[size];	
+  pbar = new float[size];	
+  pls = new float[size];	
+  prs = new float[size];	
+  ubar = new float[size];	
+  uls = new float[size];	
+  urs = new float[size];	
+  vls = new float[size];	
+  vrs = new float[size];	
+  gels = new float[size];	
+  gers = new float[size];	
+  wls = new float[size];	
+  wrs = new float[size];	
+  diffcoef = new float[size];	
+  df = new float[size];		
+  ef = new float[size];		
+  uf = new float[size];		
+  vf = new float[size];		
+  wf = new float[size];		
+  gef = new float[size];	
+  colf = new float[NumberOfColours*size];  
+  colls = new float[NumberOfColours*size];  
+  colrs = new float[NumberOfColours*size];  
+
+  /* Convert start and end indexes into 1-based for FORTRAN */
+
+  int is, ie, js, je, is_m3, ie_p3, ie_p1, k_p1;
+
+  is = GridStartIndex[0] + 1;
+  ie = GridEndIndex[0] + 1;
+  js = GridStartIndex[1] + 1;
+  je = GridEndIndex[1] + 1;
+  is_m3 = is - 3;
+  ie_p1 = ie + 1;
+  ie_p3 = ie + 3;
+  k_p1 = k + 1;
+
+  /* Compute the pressure on a slice */
+
+  if (DualEnergyFormalism)
+    FORTRAN_NAME(pgas2d_dual)(dslice, eslice, geslice, pslice, uslice, vslice, 
+			      wslice, &DualEnergyFormalismEta1, 
+			      &DualEnergyFormalismEta2, &GridDimension[0], 
+			      &GridDimension[1], &is_m3, &ie_p3, &js, &je, 
+			      &Gamma, &MinimumPressure);
+  else
+    FORTRAN_NAME(pgas2d)(dslice, eslice, pslice, uslice, vslice, 
+			 wslice, &GridDimension[0], &GridDimension[1], 
+			 &is_m3, &ie_p3, &js, &je, &Gamma, &MinimumPressure);
+
+  /* If requested, compute diffusion and slope flattening coefficients */
+
+  if (PPMDiffusionParameter != 0 || PPMFlatteningParameter != 0)
+    FORTRAN_NAME(calcdiss)(dslice, eslice, uslice, BaryonField[Vel2Num],
+			   BaryonField[Vel3Num], pslice, CellWidthTemp[0],
+			   CellWidthTemp[1], CellWidthTemp[2], &GridDimension[0],
+			   &GridDimension[1], &GridDimension[2],
+			   &is, &ie, &js, &je, &k_p1,
+			   &nzz, &dim, &GridDimension[0],
+			   &GridDimension[1], &GridDimension[2],
+			   &dtFixed, &Gamma, &PPMDiffusionParameter,
+			   &PPMFlatteningParameter, diffcoef, flatten);
+
+  /* Compute Eulerian left and right states at zone edges via interpolation */
+
+  FORTRAN_NAME(inteuler)(dslice, pslice, &GravityOn, grslice, geslice, uslice,
+			 vslice, wslice, CellWidthTemp[0], flatten,
+			 &GridDimension[0], &GridDimension[1],
+			 &is, &ie, &js, &je, &DualEnergyFormalism, 
+			 &DualEnergyFormalismEta1, &DualEnergyFormalismEta2,
+			 &PPMSteepeningParameter, &PPMFlatteningParameter,
+			 &dtFixed, &Gamma, &PressureFree, 
+			 dls, drs, pls, prs, gels, gers, uls, urs, vls, vrs,
+			 wls, wrs, &NumberOfColours, colslice, colls, colrs);
+
+  /* Compute (Lagrangian part of the) Riemann problem at each zone boundary */
+
+  FORTRAN_NAME(twoshock)(dls, drs, pls, prs, uls, urs,
+			 &GridDimension[0], &GridDimension[1],
+			 &is, &ie_p1, &js, &je,
+			 &dtFixed, &Gamma, &MinimumPressure, &PressureFree,
+			 pbar, ubar, &GravityOn, grslice,
+			 &DualEnergyFormalism, &DualEnergyFormalismEta1);
+
+  /* Compute Eulerian fluxes and update zone-centered quantities */
+
+  FORTRAN_NAME(euler)(dslice, eslice, grslice, geslice, uslice, vslice, wslice,
+		      CellWidthTemp[0], diffcoef, 
+		      &GridDimension[0], &GridDimension[1], 
+		      &is, &ie, &js, &je, &dtFixed, &Gamma, 
+		      &PPMDiffusionParameter, &GravityOn, &DualEnergyFormalism, 
+		      &DualEnergyFormalismEta1, &DualEnergyFormalismEta2,
+		      dls, drs, pls, prs, gels, gers, uls, urs, vls, vrs,
+		      wls, wrs, pbar, ubar, df, ef, uf, vf, wf, gef,
+		      &NumberOfColours, colslice, colls, colrs, colf);
+
+  /* If necessary, recompute the pressure to correctly set ge and e */
+
+  if (DualEnergyFormalism)
+    FORTRAN_NAME(pgas2d_dual)(dslice, eslice, geslice, pslice, uslice, vslice, 
+			      wslice, &DualEnergyFormalismEta1, 
+			      &DualEnergyFormalismEta2, &GridDimension[0], 
+			      &GridDimension[1], &is_m3, &ie_p3, &js, &je, 
+			      &Gamma, &MinimumPressure);
+
+  /* Check this slice against the list of subgrids (all subgrid
+     quantities are zero-based) */
+
+  int jstart, jend, offset, nfi, lface, rface, lindex, rindex, 
+    fistart, fiend, fjstart, fjend, clindex, crindex;
+  
+  for (n = 0; n < NumberOfSubgrids; n++) {
+
+    fistart = SubgridFluxes[n]->RightFluxStartGlobalIndex[dim][idim] - 
+      GridGlobalStart[idim];
+    fiend = SubgridFluxes[n]->RightFluxEndGlobalIndex[dim][idim] -
+      GridGlobalStart[idim];
+    fjstart = SubgridFluxes[n]->RightFluxStartGlobalIndex[dim][jdim] - 
+      GridGlobalStart[jdim];
+    fjend = SubgridFluxes[n]->RightFluxEndGlobalIndex[dim][jdim] -
+      GridGlobalStart[jdim];
+
+    if (k >= fjstart && k <= fjend) {
+
+      nfi = fiend - fistart + 1;
+      for (j = fistart; j <= fiend; j++) {
+
+	offset = (j-fistart) + (k-fjstart)*nfi;
+
+	lface = SubgridFluxes[n]->LeftFluxStartGlobalIndex[dim][dim] -
+	  GridGlobalStart[dim];
+	lindex = j * GridDimension[dim] + lface;
+
+	rface = SubgridFluxes[n]->RightFluxStartGlobalIndex[dim][dim] -
+	  GridGlobalStart[dim] + 1;
+	rindex = j * GridDimension[dim] + rface;	
+
+	SubgridFluxes[n]->LeftFluxes [DensNum][dim][offset] = df[lindex];
+	SubgridFluxes[n]->RightFluxes[DensNum][dim][offset] = df[rindex];
+	SubgridFluxes[n]->LeftFluxes [TENum][dim][offset]   = ef[lindex];
+	SubgridFluxes[n]->RightFluxes[TENum][dim][offset]   = ef[rindex];
+	SubgridFluxes[n]->LeftFluxes [Vel1Num][dim][offset] = uf[lindex];
+	SubgridFluxes[n]->RightFluxes[Vel1Num][dim][offset] = uf[rindex];
+
+	if (nyz > 1) {
+	  SubgridFluxes[n]->LeftFluxes [Vel2Num][dim][offset] = vf[lindex];
+	  SubgridFluxes[n]->RightFluxes[Vel2Num][dim][offset] = vf[rindex];
+	} // ENDIF y-data
+
+	if (nzz > 1) {
+	  SubgridFluxes[n]->LeftFluxes [Vel3Num][dim][offset] = wf[lindex];
+	  SubgridFluxes[n]->RightFluxes[Vel3Num][dim][offset] = wf[rindex];
+	} // ENDIF z-data
+
+	if (DualEnergyFormalism) {
+	  SubgridFluxes[n]->LeftFluxes [GENum][dim][offset] = gef[lindex];
+	  SubgridFluxes[n]->RightFluxes[GENum][dim][offset] = gef[rindex];
+	} // ENDIF DualEnergyFormalism
+
+	for (ncolour = 0; ncolour < NumberOfColours; ncolour++) {
+	  clindex = (j + ncolour * GridDimension[1]) * GridDimension[dim] +
+	    lface;
+	  crindex = (j + ncolour * GridDimension[1]) * GridDimension[dim] +
+	    rface;
+
+	  SubgridFluxes[n]->LeftFluxes [colnum[ncolour]][dim][offset] = 
+	    colf[clindex];
+	  SubgridFluxes[n]->RightFluxes[colnum[ncolour]][dim][offset] = 
+	    colf[crindex];
+	} // ENDFOR ncolour
+
+      } // ENDFOR J
+
+    } // ENDIF k inside
+
+  } // ENDFOR n
+
+  /* Copy from slice to field */
+
+  for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) {
+
+    index2 = j * GridDimension[0];
+
+    for (i = 0; i < GridDimension[0]; i++) {
+      index3 = (k*GridDimension[1] + j)*GridDimension[0] + i;
+      BaryonField[DensNum][index3] = dslice[index2+i];
+      BaryonField[TENum][index3] = eslice[index2+i];
+      BaryonField[Vel1Num][index3] = uslice[index2+i];
+    } // ENDFOR i
+
+    if (GridRank > 1)
+      for (i = 0; i < GridDimension[0]; i++) {
+	index3 = (k*GridDimension[1] + j)*GridDimension[0] + i;
+	BaryonField[Vel2Num][index3] = vslice[index2+i];
+      }
+
+    if (GridRank > 2)
+      for (i = 0; i < GridDimension[0]; i++) {
+	index3 = (k*GridDimension[1] + j)*GridDimension[0] + i;
+	BaryonField[Vel3Num][index3] = wslice[index2+i];
+      }
+
+    if (DualEnergyFormalism)
+      for (i = 0; i < GridDimension[0]; i++) {
+	index3 = (k*GridDimension[1] + j)*GridDimension[0] + i;
+	BaryonField[GENum][index3] = geslice[index2+i];
+      }
+
+    for (n = 0; n < NumberOfColours; n++) {
+      index2 = (n*GridDimension[1] + j) * GridDimension[0];
+      for (i = 0; i < GridDimension[0]; i++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	BaryonField[colnum[n]][index3] = colslice[index2+i];
+      }
+    } // ENDFOR colours
+  } // ENDFOR j
+
+  /* Delete all temporary slices */
+
+  delete [] dslice;
+  delete [] eslice;
+  delete [] uslice;
+  delete [] vslice;
+  delete [] wslice;
+  delete [] pslice;
+  if (GravityOn)
+    delete [] grslice;
+  if (DualEnergyFormalism)
+    delete [] geslice;
+  if (NumberOfColours > 0)
+    delete [] colslice;
+
+  delete [] dls;
+  delete [] drs;
+  delete [] flatten;
+  delete [] pbar;
+  delete [] pls;
+  delete [] prs;
+  delete [] ubar;
+  delete [] uls;
+  delete [] urs;
+  delete [] vls;
+  delete [] vrs;
+  delete [] gels;
+  delete [] gers;
+  delete [] wls;
+  delete [] wrs;
+  delete [] diffcoef;
+  delete [] df;
+  delete [] ef;
+  delete [] uf;
+  delete [] vf;
+  delete [] wf;
+  delete [] gef;
+  delete [] colf;
+  delete [] colls;
+  delete [] colrs;
+
+  return SUCCESS;
+
+}

src/enzo/Grid_yEulerSweep.C

+/***********************************************************************
+/
+/  GRID CLASS (WRAPPER FOR EULERIAN PPM SOLVER)
+/
+/  written by: John H. Wise
+/  date:       May 2007
+/  modified1:
+/
+/  PURPOSE:
+/
+/  RETURNS:
+/    SUCCESS or FAIL
+/
+************************************************************************/
+
+// Solve the hydro equations with the solver, saving the subgrid fluxes
+//
+
+#include <stdlib.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 "Grid.h"
+#include "euler_sweep.h"
+//#include "fortran.def"
+
+int grid::yEulerSweep(int i, int NumberOfSubgrids, fluxes *SubgridFluxes[], 
+		      Elong_int GridGlobalStart[], float *CellWidthTemp[], 
+		      int GravityOn, int NumberOfColours, int colnum[])
+{
+
+  int dim = 1, idim = 0, jdim = 2;
+
+  /* Find fields: density, total energy, velocity1-3. */
+  
+  int DensNum, GENum, Vel1Num, Vel2Num, Vel3Num, TENum;
+
+  this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, 
+				   Vel3Num, TENum);
+
+  int nxz, nyz, nzz, ixyz;
+
+  nxz = GridEndIndex[0] - GridStartIndex[0] + 1;
+  nyz = GridEndIndex[1] - GridStartIndex[1] + 1;
+  nzz = GridEndIndex[2] - GridStartIndex[2] + 1;
+
+  float MinimumPressure = tiny_number;
+  
+  // Copy from field to slice
+
+  float *dslice, *eslice, *uslice, *vslice, *wslice, *grslice, *geslice, 
+    *colslice, *pslice;
+
+  int size = GridDimension[1] * GridDimension[2];
+  dslice = new float[size];
+  eslice = new float[size];
+  uslice = new float[size];
+  vslice = new float[size];
+  wslice = new float[size];
+  pslice = new float[size];
+  if (GravityOn) {
+    grslice = new float[size];  
+  }
+  if (DualEnergyFormalism) {
+    geslice = new float[size];  
+  }
+  if (NumberOfColours > 0) {
+    colslice = new float[NumberOfColours * size];  
+  }
+
+  int j, k, n, ncolour, index2, index3;
+  for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) {
+
+    index2 = k * GridDimension[1];
+
+    for (j = 0; j < GridDimension[1]; j++) {
+      index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+      dslice[index2+j] = BaryonField[DensNum][index3];
+      eslice[index2+j] = BaryonField[TENum][index3];
+      wslice[index2+j] = BaryonField[Vel1Num][index3];
+    } // ENDFOR i
+
+    // Set velocities to zero if rank < 3 since hydro routines are
+    // hard-coded for 3-d
+
+    if (GridRank > 1) 
+      for (j = 0; j < GridDimension[1]; j++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	uslice[index2+j] = BaryonField[Vel2Num][index3];
+      }
+    else
+      for (j = 0; j < GridDimension[1]; j++)
+	uslice[index2+j] = 0;
+  
+    if (GridRank > 2)
+      for (j = 0; j < GridDimension[1]; j++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	vslice[index2+j] = BaryonField[Vel3Num][index3];
+      }
+    else
+      for (j = 0; j < GridDimension[1]; j++)
+	vslice[index2+j] = 0;
+    
+    if (GravityOn)
+      for (j = 0; j < GridDimension[1]; j++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	grslice[index2+j] = AccelerationField[dim][index3];
+      }
+
+    if (DualEnergyFormalism)
+      for (j = 0; j < GridDimension[1]; j++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	geslice[index2+j] = BaryonField[GENum][index3];
+      }
+
+    for (n = 0; n < NumberOfColours; n++) {
+      index2 = (n*GridDimension[2] + k) * GridDimension[1];
+      for (j = 0; j < GridDimension[1]; j++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	colslice[index2+j] = BaryonField[colnum[n]][index3];
+      }
+    } // ENDFOR colours
+
+  } // ENDFOR j
+
+  /* Allocate memory for fluxes */
+
+  float *dls, *drs, *flatten, *pbar, *pls, *prs, *ubar, *uls, *urs, *vls, 
+    *vrs, *gels, *gers, *wls, *wrs, *diffcoef, *df, *ef, *uf, *vf, *wf, *gef,
+    *colf, *colls, *colrs;
+
+  dls = new float[size];	
+  drs = new float[size];	
+  flatten = new float[size];	
+  pbar = new float[size];	
+  pls = new float[size];	
+  prs = new float[size];	
+  ubar = new float[size];	
+  uls = new float[size];	
+  urs = new float[size];	
+  vls = new float[size];	
+  vrs = new float[size];	
+  gels = new float[size];	
+  gers = new float[size];	
+  wls = new float[size];	
+  wrs = new float[size];	
+  diffcoef = new float[size];	
+  df = new float[size];		
+  ef = new float[size];		
+  uf = new float[size];		
+  vf = new float[size];		
+  wf = new float[size];		
+  gef = new float[size];	
+  colf = new float[NumberOfColours*size];  
+  colls = new float[NumberOfColours*size];  
+  colrs = new float[NumberOfColours*size];  
+
+  /* Convert start and end indexes into 1-based for FORTRAN */
+
+  int is, ie, js, je, is_m3, ie_p3, ie_p1, k_p1;
+
+  is = GridStartIndex[1] + 1;
+  ie = GridEndIndex[1] + 1;
+  js = GridStartIndex[2] + 1;
+  je = GridEndIndex[2] + 1;
+  is_m3 = is - 3;
+  ie_p1 = ie + 1;
+  ie_p3 = ie + 3;
+  k_p1 = i + 1;
+
+  /* Compute the pressure on a slice */
+
+  if (DualEnergyFormalism)
+    FORTRAN_NAME(pgas2d_dual)(dslice, eslice, geslice, pslice, uslice, vslice, 
+			      wslice, &DualEnergyFormalismEta1, 
+			      &DualEnergyFormalismEta2, &GridDimension[1], 
+			      &GridDimension[2], &is_m3, &ie_p3, &js, &je, 
+			      &Gamma, &MinimumPressure);
+  else
+    FORTRAN_NAME(pgas2d)(dslice, eslice, pslice, uslice, vslice,
+			 wslice, &GridDimension[1], &GridDimension[2], 
+			 &is_m3, &ie_p3, &js, &je, &Gamma, &MinimumPressure);
+
+  /* If requested, compute diffusion and slope flattening coefficients */
+
+  if (PPMDiffusionParameter != 0 || PPMFlatteningParameter != 0)
+    FORTRAN_NAME(calcdiss)(dslice, eslice, uslice, BaryonField[Vel3Num],
+			   BaryonField[Vel1Num], pslice, CellWidthTemp[1],
+			   CellWidthTemp[2], CellWidthTemp[0], 
+			   &GridDimension[1], &GridDimension[2],
+			   &GridDimension[0], &is, &ie, &js, &je, &k_p1,
+			   &nxz, &dim, &GridDimension[0],
+			   &GridDimension[1], &GridDimension[2],
+			   &dtFixed, &Gamma, &PPMDiffusionParameter,
+			   &PPMFlatteningParameter, diffcoef, flatten);
+
+  /* Compute Eulerian left and right states at zone edges via interpolation */
+
+  FORTRAN_NAME(inteuler)(dslice, pslice, &GravityOn, grslice, geslice, uslice,
+			 vslice, wslice, CellWidthTemp[1], flatten,
+			 &GridDimension[1], &GridDimension[2],
+			 &is, &ie, &js, &je, &DualEnergyFormalism, 
+			 &DualEnergyFormalismEta1, &DualEnergyFormalismEta2,
+			 &PPMSteepeningParameter, &PPMFlatteningParameter,
+			 &dtFixed, &Gamma, &PressureFree, 
+			 dls, drs, pls, prs, gels, gers, uls, urs, vls, vrs,
+			 wls, wrs, &NumberOfColours, colslice, colls, colrs);
+
+  /* Compute (Lagrangian part of the) Riemann problem at each zone boundary */
+
+  FORTRAN_NAME(twoshock)(dls, drs, pls, prs, uls, urs,
+			 &GridDimension[1], &GridDimension[2],
+			 &is, &ie_p1, &js, &je,
+			 &dtFixed, &Gamma, &MinimumPressure, &PressureFree,
+			 pbar, ubar, &GravityOn, grslice,
+			 &DualEnergyFormalism, &DualEnergyFormalismEta1);
+
+  /* Compute Eulerian fluxes and update zone-centered quantities */
+
+  FORTRAN_NAME(euler)(dslice, eslice, grslice, geslice, uslice, vslice, wslice,
+		      CellWidthTemp[1], diffcoef, 
+		      &GridDimension[1], &GridDimension[2], 
+		      &is, &ie, &js, &je, &dtFixed, &Gamma, 
+		      &PPMDiffusionParameter, &GravityOn, &DualEnergyFormalism, 
+		      &DualEnergyFormalismEta1, &DualEnergyFormalismEta2,
+		      dls, drs, pls, prs, gels, gers, uls, urs, vls, vrs,
+		      wls, wrs, pbar, ubar, df, ef, uf, vf, wf, gef,
+		      &NumberOfColours, colslice, colls, colrs, colf);
+
+  /* If necessary, recompute the pressure to correctly set ge and e */
+
+  if (DualEnergyFormalism)
+    FORTRAN_NAME(pgas2d_dual)(dslice, eslice, geslice, pslice, uslice, vslice, 
+			      wslice, &DualEnergyFormalismEta1, 
+			      &DualEnergyFormalismEta2, &GridDimension[1], 
+			      &GridDimension[2], &is_m3, &ie_p3, &js, &je, 
+			      &Gamma, &MinimumPressure);
+
+  /* Check this slice against the list of subgrids (all subgrid
+     quantities are zero-based) */
+
+  int jstart, jend, offset, nfi, lface, rface, lindex, rindex, 
+    fistart, fiend, fjstart, fjend, clindex, crindex;
+  
+  for (n = 0; n < NumberOfSubgrids; n++) {
+
+    fistart = SubgridFluxes[n]->RightFluxStartGlobalIndex[dim][idim] - 
+      GridGlobalStart[idim];
+    fiend = SubgridFluxes[n]->RightFluxEndGlobalIndex[dim][idim] -
+      GridGlobalStart[idim];
+    fjstart = SubgridFluxes[n]->RightFluxStartGlobalIndex[dim][jdim] - 
+      GridGlobalStart[jdim];
+    fjend = SubgridFluxes[n]->RightFluxEndGlobalIndex[dim][jdim] -
+      GridGlobalStart[jdim];
+
+    if (i >= fistart && i <= fiend) {
+
+      nfi = fiend - fistart + 1;
+      for (k = fjstart; k <= fjend; k++) {
+
+	offset = (i-fistart) + (k-fjstart)*nfi;
+
+	lface = SubgridFluxes[n]->LeftFluxStartGlobalIndex[dim][dim] -
+	  GridGlobalStart[dim];
+	lindex = k * GridDimension[dim] + lface;
+
+	rface = SubgridFluxes[n]->RightFluxStartGlobalIndex[dim][dim] -
+	  GridGlobalStart[dim] + 1;
+	rindex = k * GridDimension[dim] + rface;	
+
+	SubgridFluxes[n]->LeftFluxes [DensNum][dim][offset] = df[lindex];
+	SubgridFluxes[n]->RightFluxes[DensNum][dim][offset] = df[rindex];
+	SubgridFluxes[n]->LeftFluxes [TENum][dim][offset]   = ef[lindex];
+	SubgridFluxes[n]->RightFluxes[TENum][dim][offset]   = ef[rindex];
+
+	if (nxz > 1) {
+	  SubgridFluxes[n]->LeftFluxes [Vel1Num][dim][offset] = wf[lindex];
+	  SubgridFluxes[n]->RightFluxes[Vel1Num][dim][offset] = wf[rindex];
+	} // ENDIF x-data
+
+	SubgridFluxes[n]->LeftFluxes [Vel2Num][dim][offset] = uf[lindex];
+	SubgridFluxes[n]->RightFluxes[Vel2Num][dim][offset] = uf[rindex];
+
+	if (nzz > 1) {
+	  SubgridFluxes[n]->LeftFluxes [Vel3Num][dim][offset] = vf[lindex];
+	  SubgridFluxes[n]->RightFluxes[Vel3Num][dim][offset] = vf[rindex];
+	} // ENDIF z-data
+
+	if (DualEnergyFormalism) {
+	  SubgridFluxes[n]->LeftFluxes [GENum][dim][offset] = gef[lindex];
+	  SubgridFluxes[n]->RightFluxes[GENum][dim][offset] = gef[rindex];
+	} // ENDIF DualEnergyFormalism
+
+	for (ncolour = 0; ncolour < NumberOfColours; ncolour++) {
+	  clindex = (k + ncolour * GridDimension[2]) * GridDimension[dim] +
+	    lface;
+	  crindex = (k + ncolour * GridDimension[2]) * GridDimension[dim] +
+	    rface;
+
+	  SubgridFluxes[n]->LeftFluxes [colnum[ncolour]][dim][offset] = 
+	    colf[clindex];
+	  SubgridFluxes[n]->RightFluxes[colnum[ncolour]][dim][offset] = 
+	    colf[crindex];
+	} // ENDFOR ncolour
+
+      } // ENDFOR J
+
+    } // ENDIF k inside
+
+  } // ENDFOR n
+
+  /* Copy from slice to field */
+
+  for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) {
+    index2 = k * GridDimension[1];
+    for (j = 0; j < GridDimension[1]; j++) {
+      index3 = (k*GridDimension[1] + j)*GridDimension[0] + i;
+      BaryonField[DensNum][index3] = dslice[index2+j];
+      BaryonField[TENum][index3] = eslice[index2+j];
+      BaryonField[Vel1Num][index3] = wslice[index2+j];
+    } // ENDFOR i
+
+    if (GridRank > 1)
+      for (j = 0; j < GridDimension[1]; j++) {
+	index3 = (k*GridDimension[1] + j)*GridDimension[0] + i;
+	BaryonField[Vel2Num][index3] = uslice[index2+j];
+      }
+
+    if (GridRank > 2)
+      for (j = 0; j < GridDimension[1]; j++) {
+	index3 = (k*GridDimension[1] + j)*GridDimension[0] + i;
+	BaryonField[Vel3Num][index3] = vslice[index2+j];
+      }
+
+    if (DualEnergyFormalism)
+      for (j = 0; j < GridDimension[1]; j++) {
+	index3 = (k*GridDimension[1] + j)*GridDimension[0] + i;
+	BaryonField[GENum][index3] = geslice[index2+j];
+      }
+
+    for (n = 0; n < NumberOfColours; n++) {
+      index2 = (n*GridDimension[2] + k) * GridDimension[1];
+      for (j = 0; j < GridDimension[1]; j++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	BaryonField[colnum[n]][index3] = colslice[index2+j];
+      }
+    } // ENDFOR colours    
+
+  } // ENDFOR j
+
+  /* Delete all temporary slices */
+
+  delete [] dslice;
+  delete [] eslice;
+  delete [] uslice;
+  delete [] vslice;
+  delete [] wslice;
+  delete [] pslice;
+  if (GravityOn)
+    delete [] grslice;
+  if (DualEnergyFormalism)
+    delete [] geslice;
+  if (NumberOfColours > 0)
+    delete [] colslice;
+
+  delete [] dls;
+  delete [] drs;
+  delete [] flatten;
+  delete [] pbar;
+  delete [] pls;
+  delete [] prs;
+  delete [] ubar;
+  delete [] uls;
+  delete [] urs;
+  delete [] vls;
+  delete [] vrs;
+  delete [] gels;
+  delete [] gers;
+  delete [] wls;
+  delete [] wrs;
+  delete [] diffcoef;
+  delete [] df;
+  delete [] ef;
+  delete [] uf;
+  delete [] vf;
+  delete [] wf;
+  delete [] gef;
+  delete [] colf;
+  delete [] colls;
+  delete [] colrs;
+
+  return SUCCESS;
+
+}

src/enzo/Grid_zEulerSweep.C

+/***********************************************************************
+/
+/  GRID CLASS (WRAPPER FOR EULERIAN PPM SOLVER)
+/
+/  written by: John H. Wise
+/  date:       May 2007
+/  modified1:
+/
+/  PURPOSE:
+/
+/  RETURNS:
+/    SUCCESS or FAIL
+/
+************************************************************************/
+
+// Solve the hydro equations with the solver, saving the subgrid fluxes
+//
+
+#include <stdlib.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 "Grid.h"
+#include "euler_sweep.h"
+//#include "fortran.def"
+
+int grid::zEulerSweep(int j, int NumberOfSubgrids, fluxes *SubgridFluxes[], 
+		      Elong_int GridGlobalStart[], float *CellWidthTemp[], 
+		      int GravityOn, int NumberOfColours, int colnum[])
+{
+
+  int dim = 2, idim = 0, jdim = 1;
+
+  /* Find fields: density, total energy, velocity1-3. */
+  
+  int DensNum, GENum, Vel1Num, Vel2Num, Vel3Num, TENum;
+
+  this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, 
+				   Vel3Num, TENum);
+
+  int nxz, nyz, nzz, ixyz;
+
+  nxz = GridEndIndex[0] - GridStartIndex[0] + 1;
+  nyz = GridEndIndex[1] - GridStartIndex[1] + 1;
+  nzz = GridEndIndex[2] - GridStartIndex[2] + 1;
+
+  float MinimumPressure = tiny_number;
+  
+  // Copy from field to slice
+
+  float *dslice, *eslice, *uslice, *vslice, *wslice, *grslice, *geslice, 
+    *colslice, *pslice;
+
+  int size = GridDimension[2] * GridDimension[0];
+  dslice = new float[size];  
+  eslice = new float[size];  
+  uslice = new float[size];  
+  vslice = new float[size];  
+  wslice = new float[size];  
+  pslice = new float[size];  
+  if (GravityOn) {
+    grslice = new float[size];  
+  }
+  if (DualEnergyFormalism) {
+    geslice = new float[size];  
+  }
+  if (NumberOfColours > 0) {
+    colslice = new float[NumberOfColours * size];  
+  }
+
+  int i, k, n, ncolour, index2, index3;
+
+  for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++) {
+    index2 = i * GridDimension[2];
+    for (k = 0; k < GridDimension[2]; k++) {
+      index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+      dslice[index2+k] = BaryonField[DensNum][index3];
+      eslice[index2+k] = BaryonField[TENum][index3];
+      vslice[index2+k] = BaryonField[Vel1Num][index3];
+    } // ENDFOR i
+
+    // Set velocities to zero if rank < 3 since hydro routines are
+    // hard-coded for 3-d
+
+    if (GridRank > 1)
+      for (k = 0; k < GridDimension[2]; k++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	wslice[index2+k] = BaryonField[Vel2Num][index3];
+      }
+    else
+      for (k = 0; k < GridDimension[2]; k++)
+	wslice[index2+k] = 0;
+  
+    if (GridRank > 2)
+      for (k = 0; k < GridDimension[2]; k++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	uslice[index2+k] = BaryonField[Vel3Num][index3];
+      }
+    else
+      for (k = 0; k < GridDimension[2]; k++)
+	uslice[index2+k] = 0;
+
+    if (GravityOn)
+      for (k = 0; k < GridDimension[2]; k++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	grslice[index2+k] = AccelerationField[dim][index3];
+      }
+
+    if (DualEnergyFormalism)
+      for (k = 0; k < GridDimension[2]; k++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	geslice[index2+k] = BaryonField[GENum][index3];
+      }
+
+    for (n = 0; n < NumberOfColours; n++) {
+      index2 = (n*GridDimension[0] + i) * GridDimension[2];
+      for (k = 0; k < GridDimension[2]; k++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	colslice[index2+k] = BaryonField[colnum[n]][index3];
+      }
+    } // ENDFOR colours
+  } // ENDFOR j
+
+  /* Allocate memory for temporaries used in solver */
+
+  float *dls, *drs, *flatten, *pbar, *pls, *prs, *ubar, *uls, *urs, *vls, 
+    *vrs, *gels, *gers, *wls, *wrs, *diffcoef, *df, *ef, *uf, *vf, *wf, *gef,
+    *colf, *colls, *colrs;
+
+  dls = new float[size];	
+  drs = new float[size];	
+  flatten = new float[size];	
+  pbar = new float[size];	
+  pls = new float[size];	
+  prs = new float[size];	
+  ubar = new float[size];	
+  uls = new float[size];	
+  urs = new float[size];	
+  vls = new float[size];	
+  vrs = new float[size];	
+  gels = new float[size];	
+  gers = new float[size];	
+  wls = new float[size];	
+  wrs = new float[size];	
+  diffcoef = new float[size];	
+  df = new float[size];		
+  ef = new float[size];		
+  uf = new float[size];		
+  vf = new float[size];		
+  wf = new float[size];		
+  gef = new float[size];	
+  colf = new float[NumberOfColours*size];  
+  colls = new float[NumberOfColours*size];  
+  colrs = new float[NumberOfColours*size];  
+
+  /* Convert start and end indexes into 1-based for FORTRAN */
+
+  int is, ie, js, je, is_m3, ie_p3, ie_p1, k_p1;
+
+  is = GridStartIndex[2] + 1;
+  ie = GridEndIndex[2] + 1;
+  js = GridStartIndex[0] + 1;
+  je = GridEndIndex[0] + 1;
+  is_m3 = is - 3;
+  ie_p1 = ie + 1;
+  ie_p3 = ie + 3;
+  k_p1 = j + 1;
+
+  /* Compute the pressure on a slice */
+
+  if (DualEnergyFormalism)
+    FORTRAN_NAME(pgas2d_dual)(dslice, eslice, geslice, pslice, uslice, vslice, 
+			      wslice, &DualEnergyFormalismEta1, 
+			      &DualEnergyFormalismEta2, &GridDimension[2], 
+			      &GridDimension[0], &is_m3, &ie_p3, &js, &je, 
+			      &Gamma, &MinimumPressure);
+  else
+    FORTRAN_NAME(pgas2d)(dslice, eslice, pslice, uslice, vslice, 
+			 wslice, &GridDimension[2], &GridDimension[0], 
+			 &is_m3, &ie_p3, &js, &je, &Gamma, &MinimumPressure);
+
+  /* If requested, compute diffusion and slope flattening coefficients */
+
+  if (PPMDiffusionParameter != 0 || PPMFlatteningParameter != 0)
+    FORTRAN_NAME(calcdiss)(dslice, eslice, uslice, BaryonField[Vel1Num],
+			   BaryonField[Vel2Num], pslice, CellWidthTemp[1],
+			   CellWidthTemp[0], CellWidthTemp[1], 
+			   &GridDimension[2], &GridDimension[0], 
+			   &GridDimension[1], &is, &ie, &js, &je, &k_p1,
+			   &nyz, &dim, &GridDimension[0],
+			   &GridDimension[1], &GridDimension[2],
+			   &dtFixed, &Gamma, &PPMDiffusionParameter,
+			   &PPMFlatteningParameter, diffcoef, flatten);
+
+  /* Compute Eulerian left and right states at zone edges via interpolation */
+
+  FORTRAN_NAME(inteuler)(dslice, pslice, &GravityOn, grslice, geslice, uslice,
+			 vslice, wslice, CellWidthTemp[2], flatten,
+			 &GridDimension[2], &GridDimension[0],
+			 &is, &ie, &js, &je, &DualEnergyFormalism, 
+			 &DualEnergyFormalismEta1, &DualEnergyFormalismEta2,
+			 &PPMSteepeningParameter, &PPMFlatteningParameter,
+			 &dtFixed, &Gamma, &PressureFree, 
+			 dls, drs, pls, prs, gels, gers, uls, urs, vls, vrs,
+			 wls, wrs, &NumberOfColours, colslice, colls, colrs);
+
+  /* Compute (Lagrangian part of the) Riemann problem at each zone boundary */
+
+  FORTRAN_NAME(twoshock)(dls, drs, pls, prs, uls, urs,
+			 &GridDimension[2], &GridDimension[0],
+			 &is, &ie_p1, &js, &je,
+			 &dtFixed, &Gamma, &MinimumPressure, &PressureFree,
+			 pbar, ubar, &GravityOn, grslice,
+			 &DualEnergyFormalism, &DualEnergyFormalismEta1);
+
+  /* Compute Eulerian fluxes and update zone-centered quantities */
+
+  FORTRAN_NAME(euler)(dslice, eslice, grslice, geslice, uslice, vslice, wslice,
+		      CellWidthTemp[2], diffcoef, 
+		      &GridDimension[2], &GridDimension[0], 
+		      &is, &ie, &js, &je, &dtFixed, &Gamma, 
+		      &PPMDiffusionParameter, &GravityOn, &DualEnergyFormalism, 
+		      &DualEnergyFormalismEta1, &DualEnergyFormalismEta2,
+		      dls, drs, pls, prs, gels, gers, uls, urs, vls, vrs,
+		      wls, wrs, pbar, ubar, df, ef, uf, vf, wf, gef,
+		      &NumberOfColours, colslice, colls, colrs, colf);
+
+  /* If necessary, recompute the pressure to correctly set ge and e */
+
+  if (DualEnergyFormalism)
+    FORTRAN_NAME(pgas2d_dual)(dslice, eslice, geslice, pslice, uslice, vslice, 
+			      wslice, &DualEnergyFormalismEta1, 
+			      &DualEnergyFormalismEta2, &GridDimension[2], 
+			      &GridDimension[0], &is_m3, &ie_p3, &js, &je, 
+			      &Gamma, &MinimumPressure);
+
+  /* Check this slice against the list of subgrids (all subgrid
+     quantities are zero-based) */
+
+  int jstart, jend, offset, nfi, lface, rface, lindex, rindex, 
+    fistart, fiend, fjstart, fjend, clindex, crindex;
+  
+  for (n = 0; n < NumberOfSubgrids; n++) {
+
+    fistart = SubgridFluxes[n]->RightFluxStartGlobalIndex[dim][idim] - 
+      GridGlobalStart[idim];
+    fiend = SubgridFluxes[n]->RightFluxEndGlobalIndex[dim][idim] -
+      GridGlobalStart[idim];
+    fjstart = SubgridFluxes[n]->RightFluxStartGlobalIndex[dim][jdim] - 
+      GridGlobalStart[jdim];
+    fjend = SubgridFluxes[n]->RightFluxEndGlobalIndex[dim][jdim] -
+      GridGlobalStart[jdim];
+
+    if (j >= fjstart && j <= fjend) {
+
+      nfi = fiend - fistart + 1;
+      for (i = fistart; i <= fiend; i++) {
+
+	offset = (i-fistart) + (j-fjstart)*nfi;
+
+	lface = SubgridFluxes[n]->LeftFluxStartGlobalIndex[dim][dim] -
+	  GridGlobalStart[dim];
+	lindex = i * GridDimension[dim] + lface;
+
+	rface = SubgridFluxes[n]->RightFluxStartGlobalIndex[dim][dim] -
+	  GridGlobalStart[dim] + 1;
+	rindex = i * GridDimension[dim] + rface;	
+
+	SubgridFluxes[n]->LeftFluxes [DensNum][dim][offset] = df[lindex];
+	SubgridFluxes[n]->RightFluxes[DensNum][dim][offset] = df[rindex];
+	SubgridFluxes[n]->LeftFluxes [TENum][dim][offset]   = ef[lindex];
+	SubgridFluxes[n]->RightFluxes[TENum][dim][offset]   = ef[rindex];
+
+	if (nxz > 1) {
+	  SubgridFluxes[n]->LeftFluxes [Vel1Num][dim][offset] = vf[lindex];
+	  SubgridFluxes[n]->RightFluxes[Vel1Num][dim][offset] = vf[rindex];
+	} // ENDIF x-data
+
+	if (nyz > 1) {
+	  SubgridFluxes[n]->LeftFluxes [Vel2Num][dim][offset] = wf[lindex];
+	  SubgridFluxes[n]->RightFluxes[Vel2Num][dim][offset] = wf[rindex];
+	} // ENDIF y-data
+
+	SubgridFluxes[n]->LeftFluxes [Vel3Num][dim][offset] = uf[lindex];
+	SubgridFluxes[n]->RightFluxes[Vel3Num][dim][offset] = uf[rindex];
+
+	if (DualEnergyFormalism) {
+	  SubgridFluxes[n]->LeftFluxes [GENum][dim][offset] = gef[lindex];
+	  SubgridFluxes[n]->RightFluxes[GENum][dim][offset] = gef[rindex];
+	} // ENDIF DualEnergyFormalism
+
+	for (ncolour = 0; ncolour < NumberOfColours; ncolour++) {
+	  clindex = (i + ncolour * GridDimension[0]) * GridDimension[dim] +
+	    lface;
+	  crindex = (i + ncolour * GridDimension[0]) * GridDimension[dim] +
+	    rface;
+
+	  SubgridFluxes[n]->LeftFluxes [colnum[ncolour]][dim][offset] = 
+	    colf[clindex];
+	  SubgridFluxes[n]->RightFluxes[colnum[ncolour]][dim][offset] = 
+	    colf[crindex];
+	} // ENDFOR ncolour
+
+      } // ENDFOR J
+
+    } // ENDIF k inside
+
+  } // ENDFOR n
+
+  /* Copy from slice to field */
+
+  for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++) {
+    index2 = i * GridDimension[2];
+    for (k = 0; k < GridDimension[2]; k++) {
+      index3 = (k*GridDimension[1] + j)*GridDimension[0] + i;
+      BaryonField[DensNum][index3] = dslice[index2+k];
+      BaryonField[TENum][index3] = eslice[index2+k];
+      BaryonField[Vel1Num][index3] = vslice[index2+k];
+    } // ENDFOR i
+
+    if (GridRank > 1)
+      for (k = 0; k < GridDimension[2]; k++) {
+	index3 = (k*GridDimension[1] + j)*GridDimension[0] + i;
+	BaryonField[Vel2Num][index3] = wslice[index2+k];
+      }
+
+    if (GridRank > 2)
+      for (k = 0; k < GridDimension[2]; k++) {
+	index3 = (k*GridDimension[1] + j)*GridDimension[0] + i;
+	BaryonField[Vel3Num][index3] = uslice[index2+k];
+      }
+
+    if (DualEnergyFormalism)
+      for (k = 0; k < GridDimension[2]; k++) {
+	index3 = (k*GridDimension[1] + j)*GridDimension[0] + i;
+	BaryonField[GENum][index3] = geslice[index2+k];
+      }
+
+    for (n = 0; n < NumberOfColours; n++) {
+      index2 = (n*GridDimension[0] + i) * GridDimension[2];
+      for (k = 0; k < GridDimension[2]; k++) {
+	index3 = (k*GridDimension[1] + j) * GridDimension[0] + i;
+	BaryonField[colnum[n]][index3] = colslice[index2+k];
+      }
+    } // ENDFOR colours
+
+  } // ENDFOR j
+
+  /* Delete all temporary slices */
+
+  delete [] dslice;
+  delete [] eslice;
+  delete [] uslice;
+  delete [] vslice;
+  delete [] wslice;
+  delete [] pslice;
+  if (GravityOn)
+    delete [] grslice;
+  if (DualEnergyFormalism)
+    delete [] geslice;
+  if (NumberOfColours > 0)
+    delete [] colslice;
+
+  delete [] dls;
+  delete [] drs;
+  delete [] flatten;
+  delete [] pbar;
+  delete [] pls;
+  delete [] prs;
+  delete [] ubar;
+  delete [] uls;
+  delete [] urs;
+  delete [] vls;
+  delete [] vrs;
+  delete [] gels;
+  delete [] gers;
+  delete [] wls;
+  delete [] wrs;
+  delete [] diffcoef;
+  delete [] df;
+  delete [] ef;
+  delete [] uf;
+  delete [] vf;
+  delete [] wf;
+  delete [] gef;
+  delete [] colf;
+  delete [] colls;
+  delete [] colrs;
+
+  return SUCCESS;
+
+}

src/enzo/Make.config.objects

         FastSiblingLocatorFinalize.o \
         FastSiblingLocatorInitialize.o \
         FastSiblingLocatorInitializeStaticChainingMesh.o \
-        feuler_sweep.o \
         fft66.o \
         fft90.o \
         ffte4X.o \
         Grid_PoissonSolver.o                    \
         Grid_PoissonSolverCGA.o                 \
         Grid_PoissonSolverTestInitializeGrid.o  \
-        Grid_PPMDirectEuler.o \
         Grid_PrepareBoundaryFluxes.o \
         Grid_PrepareFFT.o \
         Grid_PrepareGreensFunction.o \
         Grid_SolveForPotential.o \
         Grid_SolveHighDensityPrimordialChemistry.o \
         Grid_SolveHydroEquations.o \
+	Grid_SolvePPM_DE.o \
         Grid_SolveRadiativeCooling.o \
         Grid_SolveRateAndCoolEquations.o \
         Grid_SolveRateEquations.o \
         Grid_ZeldovichPancakeInitializeGrid.o \
         Grid_ZeroSolutionUnderSubgrid.o \
         Grid_ZeusSolver.o \
+	Grid_xEulerSweep.o \
+	Grid_yEulerSweep.o \
+	Grid_zEulerSweep.o \
         Group_ReadAllData.o \
         Group_ReadDataHierarchy.o \
         Group_WriteAllData.o \

src/enzo/euler_sweep.h

+extern "C" void FORTRAN_NAME(pgas2d_dual)(
+	        float *dslice, float *eslice, float *geslice, float *pslice,
+		float *uslice, float *vslice, float *wslice, float *eta1, 
+		float *eta2, int *idim, int *jdim, int *i1, int *i2, 
+		int *j1, int *j2, float *gamma, float *pmin);
+
+extern "C" void FORTRAN_NAME(pgas2d)(
+	        float *dslice, float *eslice, float *pslice,
+		float *uslice, float *vslice, float *wslice, int *idim, 
+		int *jdim, int *i1, int *i2, int *j1, int *j2, float *gamma, 
+		float *pmin);
+
+extern "C" void FORTRAN_NAME(calcdiss)(
+	        float *dslice, float *eslice, float *uslice, float *v, 
+		float *w, float *pslice, float *dx, float *dy, float *dz, 
+		int *idim, int *jdim, int *kdim, int *i1, int *i2, int *j1, 
+		int *j2, int *k, int *nzz, int *idir, int *dimx, int *dimy, 
+		int *dimz, float *dt, float *gamma, int *idiff, int *iflatten, 
+		float *diffcoef, float *flatten);
+
+extern "C" void FORTRAN_NAME(inteuler)(
+		float *dslice, float *pslice, int *gravity, float *grslice, 
+		float *geslice, float *uslice, float *vslice, float *wslice, 
+		float *dxi, float *flatten, int *idim, int *jdim, int *i1, 
+		int *i2, int *j1, int *j2, int *idual, float *eta1, 
+		float *eta2, int *isteep, int *iflatten, float *dt, 
+		float *gamma, int *ipresfree, float *dls, float *drs, 
+		float *pls, float *prs, float *gels, float *gers, float *uls, 
+		float *urs, float *vls, float *vrs, float *wls, float *wrs, 
+		int *ncolor, float *colslice, float *colls, float *colrs);
+
+extern "C" void FORTRAN_NAME(twoshock)(
+		float *dls, float *drs, float *pls, float *prs, 
+		float *uls, float *urs, int *idim, int *jdim, int *i1,
+		int *i2, int *j1, int *j2, float *dt, float *gamma, 
+		float *pmin, int *ipresfree, float *pbar, float *ubar,
+		int *gravity, float *grslice, int *idual, float *eta1);
+
+extern "C" void FORTRAN_NAME(euler)(
+		float *dslice, float *pslice, float *grslice, 
+		float *geslice, float *uslice, float *vslice, float *wslice, 
+		float *dx, float *diffcoef, int *idim, int *jdim, int *i1, 
+		int *i2, int *j1, int *j2, float *dt, float *gamma, int *idiff, 
+		int *gravity, int *idual, float *eta1, float *eta2, 
+		float *dls, float *drs, float *pls, float *prs, float *gels, 
+		float *gers, float *uls, float *urs, float *vls, float *vrs, 
+		float *wls, float *wrs, float *pbar, float *ubar, float *df,
+		float *ef, float *uf, float *vf, float *wf, float *gef,
+		int *ncolor, float *colslice, float *colls, float *colrs,
+		float *colf);