Commits

Matthew Turk committed 9c18eba

Intermediate commit of removing IO routines. Working on compilation next.

Comments (0)

Files changed (23)

 
    ~grid();
 
-/* Read grid data from a file (returns: success/failure) */
-
-   int ReadGrid(FILE *main_file_pointer, int GridID, char DataFilename[], 
-		int ReadText = TRUE, int ReadData = TRUE);
-
 /* Read grid data from a group file (returns: success/failure) */
 
-#ifndef NEW_GRID_IO
-   int Group_ReadGrid(FILE *fptr, int GridID, HDF5_hid_t file_id, 
+   int ReadGrid(FILE *fptr, int GridID, HDF5_hid_t file_id, 
 		      char DataFilename[],
 		      int ReadText, int ReadData, bool ReadParticlesOnly=false);
-#else
-   int Group_ReadGrid(FILE *main_file_pointer, int GridID, HDF5_hid_t file_id, 
-		      char DataFilename[],
-		      int ReadText = TRUE, int ReadData = TRUE, 
-		      bool ReadParticlesOnly=false, int ReadEverything = FALSE);
-#endif
    
 /* Get field or particle data based on name or integer 
    defined in typedefs.h. Details are in Grid_CreateFieldArray.C. */
 
 /* Write grid data to a group file (returns: success/failure) */
 
-#ifndef NEW_GRID_IO
-   int Group_WriteGrid(FILE *fptr, char *base_name, int grid_id, HDF5_hid_t file_id);
-#else
-   int Group_WriteGrid(FILE *main_file_pointer, char *base_name, int grid_id, HDF5_hid_t file_id, int WriteEverything = FALSE);
-#endif
+   int WriteGrid(FILE *main_file_pointer, char *base_name, int grid_id, HDF5_hid_t file_id, int WriteEverything = FALSE);
 
    int WriteAllFluxes(hid_t grid_node);
    int WriteFluxGroup(hid_t top_group, fluxes *fluxgroup);
 /* Interpolate to specified time and write grid data to a group file
    (returns: success/failure). */
 
-   int Group_WriteGridInterpolate(FLOAT WriteTime, FILE *main_file_pointer,
+   int WriteGridInterpolate(FLOAT WriteTime, FILE *main_file_pointer,
                             char *base_name, int grid_id, HDF5_hid_t file_id);
 
    int ComputeVectorAnalysisFields(field_type fx, field_type fy, field_type fz,

src/enzo/Grid_Group_ReadGrid.C

-/***********************************************************************
-/
-/  GRID CLASS (READ GRID)
-/
-/  written by: Greg Bryan
-/  date:       November, 1994
-/  modified1:  Robert Harkness, July 2002
-/  modified2:  Alexei Kritsuk, Jan 2004   a trick for RandomForcing //AK
-/  modified3:  Robert Harkness, Jan 2007 for HDF5 memory buffering
-/  modified4:  Robert Harkness, April 2008
-/  modified5:  Michael Kuhlen, October 2010, HDF5 hierarchy
-/
-/  PURPOSE:
-/
-************************************************************************/
- 
-//  Input a grid from file pointer fpt
- 
-#include <hdf5.h>
-#include <string.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <unistd.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"
-void my_exit(int status);
- 
-#ifdef PROTO /* Remove troublesome HDF PROTO declaration. */
-#undef PROTO
-#endif
- 
-// HDF5 function prototypes
- 
-
- 
-// function prototypes
- 
-int ReadListOfFloats(FILE *fptr, int N, FLOAT floats[]);
-int ReadListOfInts(FILE *fptr, int N, int nums[]);
- 
-static int GridReadDataGridCounter = 0;
- 
- 
-#ifndef NEW_GRID_IO
-int grid::Group_ReadGrid(FILE *fptr, int GridID, HDF5_hid_t file_id, 
-			 char DataFilename[],
-			 int ReadText, int ReadData, bool ReadParticlesOnly)
-{
- 
-  int i, j, k, field, size, active_size;
-  char name[MAX_LINE_LENGTH], dummy[MAX_LINE_LENGTH];
-  char logname[MAX_LINE_LENGTH];
-  char procfilename[MAX_LINE_LENGTH];
- 
-  char id[MAX_GROUP_TAG_SIZE];
-  char pid[MAX_TASK_TAG_SIZE];
-  char gpid[MAX_TASK_TAG_SIZE];
- 
-  int ActiveDim[MAX_DIMENSION];
- 
-  FILE *log_fptr;
- 
-  hid_t       group_id, dset_id;
-  hid_t       float_type_id, FLOAT_type_id;
-  hid_t       file_type_id, FILE_type_id;
-  hid_t       file_dsp_id;
-  hid_t       num_type;
- 
-  hsize_t     OutDims[MAX_DIMENSION];
-  hsize_t     TempIntArray[MAX_DIMENSION];
- 
-  herr_t      h5_status;
-  herr_t      h5_error = -1;
- 
-  int         num_size;
- 
-  char *ParticlePositionLabel[] =
-    {"particle_position_x", "particle_position_y", "particle_position_z"};
-  char *ParticleVelocityLabel[] =
-    {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"};
-#ifdef WINDS
-  char *ParticleAttributeLabel[] =
-    {"creation_time", "dynamical_time", "metallicity_fraction", "particle_jet_x", 
-     "particle_jet_y", "particle_jet_z", "typeia_fraction"};
-#else
-  char *ParticleAttributeLabel[] = 
-    {"creation_time", "dynamical_time", "metallicity_fraction", "typeia_fraction"};
-#endif
- 
-#ifdef IO_LOG
-  int         io_log = 1;
-#else
-  int         io_log = 0;
-#endif
-
-#ifdef IO_64
-#define io_type float64
-#else
-#define io_type float32
-#endif
- 
-  if(ReadText && HierarchyFileInputFormat == 1){
-
-    /* Read general grid class data */
-
-    /* make sure quantities defined at least for 3d */
- 
-    for (int dim = GridRank; dim < 3; dim++) {
-      GridDimension[dim] = 1;
-      GridStartIndex[dim] = 0;
-      GridEndIndex[dim] = 0;
-    }
-    if (fscanf(fptr, "GridRank = %"ISYM"\n", &GridRank) != 1) {
-            ENZO_FAIL("Error reading GridRank.");
-    }
- 
-    if (fscanf(fptr, "GridDimension = ") != 0) {
-            ENZO_FAIL("Error reading GridDimension(0).");
-    }
- 
-    if (ReadListOfInts(fptr, GridRank, GridDimension) == FAIL) {
-            ENZO_FAIL("Error reading GridDimension(1).");
-    }
- 
-    fscanf(fptr, "GridStartIndex = ");
- 
-    if (ReadListOfInts(fptr, GridRank, GridStartIndex) == FAIL) {
-            ENZO_FAIL("Error reading GridStartIndex.");
-    }
- 
-    fscanf(fptr, "GridEndIndex = ");
- 
-    if (ReadListOfInts(fptr, GridRank, GridEndIndex) == FAIL) {
-            ENZO_FAIL("Error reading GridEndIndex.");
-    }
- 
-    fscanf(fptr, "GridLeftEdge = ");
- 
-    if (ReadListOfFloats(fptr, GridRank, GridLeftEdge) == FAIL) {
-            ENZO_FAIL("Error reading GridLeftEdge.");
-    }
- 
-    fscanf(fptr, "GridRightEdge = ");
- 
-    if (ReadListOfFloats(fptr, GridRank, GridRightEdge) == FAIL) {
-            ENZO_FAIL("Error reading GridRightEdge.");
-    }
- 
-    if (fscanf(fptr, "Time = %"PSYM"\n", &Time) != 1) {
-            ENZO_FAIL("Error reading Time.");
-    }
- 
-    if (fscanf(fptr, "SubgridsAreStatic = %"ISYM"\n", &SubgridsAreStatic) != 1) {
-            ENZO_FAIL("Error reading SubgridsAreStatic.");
-    }
- 
-    /* Read baryon field quantities. */
- 
-    if (fscanf(fptr, "NumberOfBaryonFields = %"ISYM"\n",
-	       &NumberOfBaryonFields) != 1) {
-            ENZO_FAIL("Error reading NumberOfBaryonFields.");
-    }
-    if (NumberOfBaryonFields > 0) {
-
-      if (NumberOfBaryonFields >= MAX_NUMBER_OF_BARYON_FIELDS) {
-	ENZO_VFAIL("NumberOfBaryonFields (%"ISYM") exceeds "
-	       "MAX_NUMBER_OF_BARYON_FIELDS (%"ISYM").\n", 
-	       NumberOfBaryonFields, MAX_NUMBER_OF_BARYON_FIELDS)
-      }
-
-      fscanf(fptr, "FieldType = ");
- 
-      if (ReadListOfInts(fptr, NumberOfBaryonFields, FieldType) == FAIL) {
-		ENZO_FAIL("Error reading FieldType.");
-      }
- 
-      fgetpos(fptr, &BaryonFileNamePosition); //AK
- 
-      if (fscanf(fptr, "BaryonFileName = %s\n", procfilename) != 1) {
-		ENZO_FAIL("Error reading BaryonFileName.");
-      }
- 
-      fscanf(fptr, "CourantSafetyNumber    = %"FSYM"\n", &CourantSafetyNumber);
-      fscanf(fptr, "PPMFlatteningParameter = %"ISYM"\n", &PPMFlatteningParameter);
-      fscanf(fptr, "PPMDiffusionParameter  = %"ISYM"\n", &PPMDiffusionParameter);
-      fscanf(fptr, "PPMSteepeningParameter = %"ISYM"\n", &PPMSteepeningParameter);
-    }
-
-    /* 3) Read particle info */
- 
-    if (fscanf(fptr, "NumberOfParticles = %"ISYM"\n", &NumberOfParticles) != 1) {
-            ENZO_FAIL("error reading NumberOfParticles.");
-    }
- 
-    if (NumberOfParticles > 0) {
- 
-      /* Read particle file name. */
-    
-      if (fscanf(fptr, "ParticleFileName = %s\n", procfilename) != 1) {
-		ENZO_FAIL("Error reading ParticleFileName.");
-      }
-    }
- 
-    /* 4) Read gravity info */
- 
-    if (SelfGravity)
-      if (fscanf(fptr, "GravityBoundaryType = %"ISYM"\n",&GravityBoundaryType) != 1) {
-		ENZO_FAIL("Error reading GravityBoundaryType.");
-      }
-
-    // If HierarchyFile has different Ghostzones (which should be a parameter not a macro ...)
-    // (useful in a restart with different hydro/mhd solvers) 
-    int ghosts =DEFAULT_GHOST_ZONES;
-    if (GridStartIndex[0] != ghosts)  {
-	if (GridID < 2)
-     fprintf(stderr,"Grid_Group_ReadGrid: Adjusting Ghostzones which in the hierarchy file did not match the selected HydroMethod.\n");
-      for (int dim=0; dim < GridRank; dim++) {
-	GridDimension[dim]  = GridEndIndex[dim]-GridStartIndex[dim]+1+2*ghosts;
-	GridStartIndex[dim] = ghosts;
-	GridEndIndex[dim]   = GridStartIndex[dim]+GridDimension[dim]-1-2*ghosts;
-	 if (GridID < 2) fprintf(stderr, "dim: GridStart,GridEnd,GridDim:  %i: %i %i %i\n",
-				  dim, GridStartIndex[dim], GridEndIndex[dim], GridDimension[dim]);
-      }
-    }
-
-
-  } // (if (ReadText) )
-
-  // if HDF5 Hierarchy file, then copy DataFilename (read in
-  // Grid::ReadHierarchyInformationHDF5.C) to procfilename
-  if (HierarchyFileInputFormat % 2 == 0) {
-    strcpy(procfilename, DataFilename);
-  }
- 
-
-  int ii = sizeof(io_type);
-
-  switch(ii)
-    {
- 
-    case 4:
-      float_type_id = HDF5_R4;
-      break;
- 
-    case 8:
-      float_type_id = HDF5_R8;
-      break;
- 
-    default:
-      float_type_id = HDF5_R4;
- 
-    }
- 
-  int jj = sizeof(FLOAT);
- 
-  switch(jj)
-    {
- 
-    case 4:
-      FLOAT_type_id = HDF5_R4;
-      FILE_type_id = HDF5_FILE_R4;
-      break;
- 
-    case 8:
-      FLOAT_type_id = HDF5_R8;
-      FILE_type_id = HDF5_FILE_R8;
-      break;
- 
-    case 16:
-      FLOAT_type_id = HDF5_R16;
-      FILE_type_id = H5Tcopy(HDF5_FILE_B8);
-      H5Tset_size(FILE_type_id,16);
-      break;
- 
-    default:
-      printf("INCORRECT FLOAT DEFINITION\n");
- 
-    }
- 
-  sprintf(id, "%"GROUP_TAG_FORMAT""ISYM, GridID);
- 
-  sprintf(pid, "%"TASK_TAG_FORMAT""ISYM, MyProcessorNumber);
- 
-  sprintf(gpid, "%"TASK_TAG_FORMAT""ISYM, ProcessorNumber);
- 
-  strcpy(name, "/Grid");
-  strcat(name, id);
- 
-  if (NumberOfBaryonFields > 0 && ReadData && !ReadParticlesOnly &&
-      (MyProcessorNumber == ProcessorNumber)) {
-
-    strcpy(logname, procfilename);
-    strcat(logname, ".in_log");
-    if (io_log) log_fptr = fopen(logname, "a");
-	
-#ifndef SINGLE_HDF5_OPEN_ON_INPUT
-    if (io_log) fprintf(log_fptr, "H5Fopen with Name %s\n", procfilename);
-    file_id = H5Fopen(procfilename,  H5F_ACC_RDONLY, H5P_DEFAULT);
-    if (io_log) fprintf(log_fptr, "H5Fopen id: %"ISYM"\n", file_id);
-    if( file_id == h5_error ){ENZO_FAIL("Error in IO");}
-
-#endif
- 
-    if (io_log) fprintf(log_fptr, "H5Gopen with Name %s\n", name);
-
-    group_id = H5Gopen(file_id, name);
-    if( group_id == h5_error )
-        ENZO_VFAIL("Error in IO (%s)", name)
- 
-    /* fill in ActiveDim for dims up to 3d */
- 
-    for (int dim = 0; dim < 3; dim++)
-      ActiveDim[dim] = GridEndIndex[dim] - GridStartIndex[dim] +1;
- 
-    /* check dimensions of HDF file against this grid
-       (note: we don't bother to check the coordinate arrays)  */
- 
-    size = 1;
-    active_size = 1;
- 
-    for (int dim = 0; dim < GridRank; dim++) {
-      size *= GridDimension[dim];
-      active_size *= ActiveDim[dim];
-    }
- 
-    //  CAUTION - are the coordinates reversed?
- 
-    for (int dim = 0; dim < GridRank; dim++) {
-      OutDims[GridRank-dim-1] = ActiveDim[dim];
-      if (io_log) fprintf(log_fptr, "Outdims %"ISYM"\n", (int) OutDims[GridRank-dim-1]);
-    }
- 
-    /* allocate temporary space */
- 
-    io_type *temp = new io_type[active_size];
- 
-    /* loop over fields, reading each one */
- 
-    for (field = 0; field < NumberOfBaryonFields; field++) {
- 
-      /* get data into temporary array */
- 
-      file_dsp_id = H5Screate_simple((Eint32) GridRank, OutDims, NULL);
-      if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-      if( file_dsp_id == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      if (io_log) fprintf(log_fptr, "H5Dopen with Name = %s\n", DataLabel[field]);
- 
-      dset_id =  H5Dopen(group_id, DataLabel[field]);
-      if (io_log) fprintf(log_fptr, "H5Dopen id: %"ISYM"\n", dset_id);
-      //      if( dset_id == h5_error ){my_exit(EXIT_FAILURE);}
-       if( dset_id == h5_error ){
-	 fprintf(stderr, "NumberOfBaryonFields = %d", field);
-	 my_exit(EXIT_FAILURE);
-       }
-
-      h5_status = H5Dread(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp);
-      if (io_log) fprintf(log_fptr, "H5Dread: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      h5_status = H5Sclose(file_dsp_id);
-      if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      h5_status = H5Dclose(dset_id);
-      if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      /* copy active region into whole grid */
- 
-      BaryonField[field] = new float[size];
- 
-      for (i = 0; i < size; i++)
-	BaryonField[field][i] = 0;
- 
-      for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++)
-	for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++)
-	  for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++)
-	    BaryonField[field][i + j*GridDimension[0] +
-			       k*GridDimension[0]*GridDimension[1]] =
-	      float(temp[(i-GridStartIndex[0])                         +
-	                 (j-GridStartIndex[1])*ActiveDim[0]            +
-	                 (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ]);
- 
-    } // end: loop over fields
- 
-
-    if (HydroMethod == MHD_RK) { // This is the MHD with Dedner divergence cleaning that needs an extra field
-   
-      int activesize = 1;
-      for (int dim = 0; dim < GridRank; dim++)
-	activesize *= (GridDimension[dim]-2*DEFAULT_GHOST_ZONES);
-      
-      if (divB == NULL) 
-	divB = new float[activesize];
-      
-      /* if we restart from a different solvers output without a PhiField create here and set to zero */
-      int PhiNum; 
-      if ((PhiNum = FindField(PhiField, FieldType, NumberOfBaryonFields)) < 0) {
-	fprintf(stderr, "Starting with Dedner MHD method with no Phi field. \n");
-	fprintf(stderr, "Adding it in Grid_ReadGrid.C \n");
-	char *PhiName = "Phi";
-	PhiNum = NumberOfBaryonFields;
-	int PhiToAdd = PhiField;
-	this->AddFields(&PhiToAdd, 1);
-	DataLabel[PhiNum] = PhiName;
-      } else { 
-	if (0) 
-	  for (int n = 0; n < size; n++)
-	    BaryonField[PhiNum][n] = 0.;
-      }
-
-      /* if we restart from a different solvers output without a Phi_pField 
-	 and yet want to use the divergence cleaning, create here and set to zero */
-      if (UseDivergenceCleaning) {
-	int Phi_pNum; 
-	if ((Phi_pNum = FindField(Phi_pField, FieldType, NumberOfBaryonFields)) < 0) {
-	  fprintf(stderr, "Want to use divergence cleaning with no Phi_p field. \n");
-	  fprintf(stderr, "Adding it in Grid_ReadGrid.C \n");
-	  char *Phi_pName = "Phi_p";
-	  Phi_pNum = NumberOfBaryonFields;
-	  int Phi_pToAdd = Phi_pField;
-	  this->AddFields(&Phi_pToAdd, 1);
-	  DataLabel[Phi_pNum] = Phi_pName;
-	}
-      }
-
-      for (int dim = 0; dim < 3; dim++)
-	if (gradPhi[dim] == NULL)
-	  gradPhi[dim] = new float[activesize];
-      
-      for (int dim = GridRank; dim < 3; dim++)
-	for (int n = 0; n < activesize; n++)
-	  gradPhi[dim][n] = 0.0;
-      
-    } /* if HydroMethod == MHD */
-
-
-    delete [] temp;
- 
-  }  // end:   if (NumberOfBaryonFields > 0 && ReadData &&
-  //      (MyProcessorNumber == ProcessorNumber)) {
-
-  /* Compute Flux quantities */
-
-  this->PrepareGridDerivedQuantities();
- 
- 
-  if (NumberOfParticles > 0 && ReadData &&
-      (MyProcessorNumber == ProcessorNumber)) {
-  
- 
-    /* Open file if not already done (note: particle name must = grid name). */
- 
-    if (NumberOfBaryonFields == 0 || ReadParticlesOnly) {
- 
-      strcpy(logname, procfilename);
-      strcat(logname, ".in_log");
-      if (io_log) log_fptr = fopen(logname, "a");
-
-
-#ifndef SINGLE_HDF5_OPEN_ON_INPUT 
-
-      if (io_log) fprintf(log_fptr, "H5Fopen with Name %s\n", procfilename);
-      file_id = H5Fopen(procfilename, H5F_ACC_RDONLY, H5P_DEFAULT);
-      if (io_log) fprintf(log_fptr, "H5Fopen id: %"ISYM"\n", file_id);
-      if( file_id == h5_error ){ENZO_FAIL("Error in IO");}
-
-#endif
- 
-      if (io_log) fprintf(log_fptr, "H5Gopen with Name %s\n", name);
- 
-      group_id = H5Gopen(file_id, name);
-      if( group_id == h5_error ){ENZO_FAIL("Error in IO");}
- 
-    } // end: if (NumberOfBaryonFields == 0)
- 
-    /* Allocate room for particles. */
- 
-    this->AllocateNewParticles(NumberOfParticles);
- 
-    TempIntArray[0] = NumberOfParticles;
- 
-    /* Create a temporary buffer (32 bit or twice the size for 64). */
- 
-    io_type *temp = NULL;
- 
-    jj = sizeof(FLOAT);
- 
-    switch(jj)
-      {
- 
-      case 4:
-        temp = new io_type[NumberOfParticles];
-        break;
- 
-      case 8:
-        temp = new io_type[NumberOfParticles*2];
-        break;
- 
-      case 16:
-        temp = new io_type[NumberOfParticles*4];
-        break;
- 
-      default:
-        printf("INCORRECT FLOAT DEFINITION\n");
- 
-      }
- 
-    if (temp == NULL)
-      temp = new io_type[NumberOfParticles];
- 
-    /* Read ParticlePosition (use temporary buffer). */
- 
-    for (int dim = 0; dim < GridRank; dim++) {
- 
-      file_dsp_id = H5Screate_simple((Eint32) 1, TempIntArray, NULL);
-      if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-      if( file_dsp_id == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      if (io_log) fprintf(log_fptr,"H5Dopen with Name = %s\n", ParticlePositionLabel[dim]);
- 
-      dset_id =  H5Dopen(group_id, ParticlePositionLabel[dim]);
-      if (io_log) fprintf(log_fptr, "H5Dopen id: %"ISYM"\n", dset_id);
-      if( dset_id == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      num_type = H5Dget_type(dset_id);
-      num_size = H5Tget_size(num_type);
- 
-      if (sizeof(FLOAT) == 16)
-	{
- 
-	  //                                 NOTE: for 128bits this must be FILE_type_id and NOT FLOAT_type_id!
-	  h5_status = H5Dread(dset_id, FILE_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) ParticlePosition[dim]);
-	  if (io_log) fprintf(log_fptr, "H5Dread: %"ISYM"\n", h5_status);
-	  if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-	}
-      else
-	{
- 
-	  h5_status = H5Dread(dset_id, FLOAT_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) ParticlePosition[dim]);
-	  if (io_log) fprintf(log_fptr, "H5Dread: %"ISYM"\n", h5_status);
-	  if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-	}
- 
-      h5_status = H5Sclose(file_dsp_id);
-      if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      h5_status = H5Dclose(dset_id);
-      if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-    }
- 
- 
-    /* Read ParticleVelocity. */
- 
-    for (int dim = 0; dim < GridRank; dim++) {
- 
-      file_dsp_id = H5Screate_simple((Eint32) 1, TempIntArray, NULL);
-      if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-      if( file_dsp_id == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      if (io_log) fprintf(log_fptr, "H5Dopen with Name = %s\n", ParticleVelocityLabel[dim]);
- 
-      dset_id =  H5Dopen(group_id, ParticleVelocityLabel[dim]);
-      if (io_log) fprintf(log_fptr, "H5Dopen id: %"ISYM"\n", dset_id);
-      if( dset_id == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      h5_status = H5Dread(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp);
-      if (io_log) fprintf(log_fptr, "H5Dread: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      h5_status = H5Sclose(file_dsp_id);
-      if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      h5_status = H5Dclose(dset_id);
-      if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      for (i = 0; i < NumberOfParticles; i++)
-	ParticleVelocity[dim][i] = float(temp[i]);
-    }
- 
- 
-    /* Read ParticleMass into temporary buffer and Copy to ParticleMass. */
- 
-    file_dsp_id = H5Screate_simple((Eint32) 1, TempIntArray, NULL);
-    if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-    if( file_dsp_id == h5_error ){ENZO_FAIL("Error in IO");}
- 
-    if (io_log) fprintf(log_fptr,"H5Dopen with Name = particle_mass\n");
- 
-    dset_id =  H5Dopen(group_id, "particle_mass");
-    if (io_log) fprintf(log_fptr, "H5Dopen id: %"ISYM"\n", dset_id);
-    if( dset_id == h5_error ){ENZO_FAIL("Error in IO");}
- 
-    h5_status = H5Dread(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp);
-    if (io_log) fprintf(log_fptr, "H5Dread: %"ISYM"\n", h5_status);
-    if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-    h5_status = H5Sclose(file_dsp_id);
-    if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-    if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-    h5_status = H5Dclose(dset_id);
-    if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-    if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-    for (i = 0; i < NumberOfParticles; i++)
-      ParticleMass[i] = float(temp[i]);
- 
-    /* Read ParticleNumber into temporary buffer and Copy to ParticleNumber. */
- 
-    PINT *tempPINT = new PINT[NumberOfParticles];
- 
-    file_dsp_id = H5Screate_simple((Eint32) 1, TempIntArray, NULL);
-    if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-    if( file_dsp_id == h5_error){ENZO_FAIL("Error in IO");}
- 
-    if (io_log) fprintf(log_fptr,"H5Dopen  with Name = particle_index\n");
- 
-    dset_id =  H5Dopen(group_id, "particle_index");
-    if (io_log) fprintf(log_fptr, "H5Dopen id: %"ISYM"\n", dset_id);
-    if( dset_id == h5_error ){ENZO_FAIL("Error in IO");}
- 
-    h5_status = H5Dread(dset_id, HDF5_PINT, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) tempPINT);
-    if (io_log) fprintf(log_fptr, "H5Dread: %"ISYM"\n", h5_status);
-    if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-    h5_status = H5Sclose(file_dsp_id);
-    if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-    if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-    h5_status = H5Dclose(dset_id);
-    if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-    if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-    for (i = 0; i < NumberOfParticles; i++)
-      ParticleNumber[i] = tempPINT[i];
- 
- 
-    // Read ParticleType if present
-
-    H5E_BEGIN_TRY{
-      dset_id = H5Dopen(group_id, "particle_type");
-    }H5E_END_TRY
- 
-    if (ParticleTypeInFile == TRUE && dset_id != h5_error) {
-
-      H5Dclose(dset_id);
-
-      /* Read ParticleType into temporary buffer and Copy to ParticleType. */
-
-      int *tempint = new int[NumberOfParticles];
- 
-      file_dsp_id = H5Screate_simple((Eint32) 1, TempIntArray, NULL);
-      if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-      if( file_dsp_id == h5_error){ENZO_FAIL("Error in IO");}
- 
-      if (io_log) fprintf(log_fptr,"H5Dopen  with Name = particle_type\n");
- 
-      dset_id =  H5Dopen(group_id, "particle_type");
-      if (io_log) fprintf(log_fptr, "H5Dopen id: %"ISYM"\n", dset_id);
-      if( dset_id == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      h5_status = H5Dread(dset_id, HDF5_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) tempint);
-      if (io_log) fprintf(log_fptr, "H5Dread: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      h5_status = H5Sclose(file_dsp_id);
-      if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      h5_status = H5Dclose(dset_id);
-      if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      for (i = 0; i < NumberOfParticles; i++)
-        ParticleType[i] = tempint[i];
- 
-      int abs_type;
-      for (i = 0; i < NumberOfParticles; i++) {
-	abs_type = ABS(ParticleType[i]);
-        if (abs_type < PARTICLE_TYPE_GAS ||
-            abs_type > NUM_PARTICLE_TYPES-1) {
-          ENZO_VFAIL("file: %s: particle %"ISYM" has unknown type %"ISYM"\n",
-                  name, i, ParticleType[i])
-        }
-      }
-
-      delete [] tempint;
- 
-    } else {
- 
-      /* Otherwise create the type. */
- 
-      for (i = 0; i < NumberOfParticles; i++)
-        ParticleType[i] = ReturnParticleType(i);
- 
-    }
- 
- 
-    /* Read ParticleAttributes. */
-    if (AddParticleAttributes) {
-      for (j = 0; j < NumberOfParticleAttributes; j++) {
-	ParticleAttribute[j] = new float[NumberOfParticles];
-	for (i=0; i < NumberOfParticles; i++)
-	  ParticleAttribute[j][i] = 0;
-      }
-    } else {
-    for (j = 0; j < NumberOfParticleAttributes; j++) {
-
-      H5E_BEGIN_TRY{
-	dset_id = H5Dopen(group_id, ParticleAttributeLabel[j]);
-      }H5E_END_TRY;
-
-      if (dset_id != h5_error) {
-       
-      H5Dclose(dset_id);
- 
-      file_dsp_id = H5Screate_simple((Eint32) 1, TempIntArray, NULL);
-      if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-      if( file_dsp_id == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      if (io_log) fprintf(log_fptr,"H5Dopen with Name = %s\n",ParticleAttributeLabel[j]);
- 
-      dset_id =  H5Dopen(group_id, ParticleAttributeLabel[j]);
-      if (io_log) fprintf(log_fptr, "H5Dopen id: %"ISYM"\n", dset_id);
-      if( dset_id == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      h5_status = H5Dread(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp);
-      if (io_log) fprintf(log_fptr, "H5Dread: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      h5_status = H5Sclose(file_dsp_id);
-      if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      h5_status = H5Dclose(dset_id);
-      if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
- 
-      for (i = 0; i < NumberOfParticles; i++)
-	ParticleAttribute[j][i] = float(temp[i]);
-
-      } // ENDIF dset_id != h5_error
-      else {
-	
-	ParticleAttribute[j] = new float[NumberOfParticles];
-	for (i=0; i < NumberOfParticles; i++)
-	  ParticleAttribute[j][i] = 0;
-
-      } // ENDELSE
- 
-    }
-    } // ENDELSE add particle attributes
- 
-    delete [] temp;
-    delete [] tempPINT;
- 
-
-  } // end: if (NumberOfParticles > 0) && ReadData && (MyProcessorNumber == ProcessorNumber)
- 
-  /* Close file. */
- 
-  if ( (MyProcessorNumber == ProcessorNumber) &&
-       (NumberOfParticles > 0 || 
-	(NumberOfBaryonFields > 0 && !ReadParticlesOnly))
-       && ReadData ){
- 
-    h5_status = H5Gclose(group_id);
-    if (io_log) fprintf(log_fptr, "H5Gclose: %"ISYM"\n", h5_status);
-    if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
-
-#ifndef SINGLE_HDF5_OPEN_ON_INPUT 
-
-    h5_status = H5Fclose(file_id);
-    if (io_log) fprintf(log_fptr, "H5Fclose: %"ISYM"\n", h5_status);
-    if( h5_status == h5_error ){ENZO_FAIL("Error in IO");}
-
-#endif
-  }
- 
-  if (MyProcessorNumber == ProcessorNumber)
-    {
-      if (io_log) fclose(log_fptr);
-
-    }
- 
-  return SUCCESS;
- 
-}
-#endif

src/enzo/Grid_Group_WriteGrid.C

-/***********************************************************************
-/
-/  GRID CLASS (WRITE OUT GRID)
-/
-/  written by: Greg Bryan
-/  date:       November, 1994
-/  modified1:  Robert Harkness, July 2002
-/  modified2:  Robert Harkness, July 2006
-/  modified3:  Robert Harkness, April 2008
-/  modified4:  Michael Kuhlen, October 2010, HDF5 hierarchy
-/
-/  PURPOSE:
-/
-************************************************************************/
- 
-//  Write grid to file pointer fptr
-//     (we assume that the grid is at an appropriate stopping point,
-//      where the Old values aren't required)
- 
-#include <hdf5.h>
-#include <string.h>
-#include <stdio.h>
-#include <stdlib.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"
-void my_exit(int status);
- 
-// HDF5 function prototypes
- 
-
- 
-// function prototypes
- 
-void WriteListOfFloats(FILE *fptr, int N, FLOAT floats[]);
-void WriteListOfInts(FILE *fptr, int N, int nums[]);
-int WriteStringAttr(hid_t dset_id, char *Alabel, char *String, FILE *log_fptr);
-int FindField(int field, int farray[], int numfields);
-
-int GetUnits(float *DensityUnits, float *LengthUnits,
-	     float *TemperatureUnits, float *TimeUnits,
-	     float *VelocityUnits, FLOAT Time);
-
-#ifndef NEW_GRID_IO
-int grid::Group_WriteGrid(FILE *fptr, char *base_name, int grid_id, HDF5_hid_t file_id)
-{
- 
-  int i, j, k, dim, field, size, active_size, ActiveDim[MAX_DIMENSION];
-  int file_status;
-
-  float *temperature, *dust_temperature,
-    *cooling_time;
-
-#ifdef IO_64
-#define io_type float64
-#else
-#define io_type float32
-#endif
- 
-  io_type *temp, *temp_VelAnyl;
- 
-  FILE *log_fptr;
-  FILE *procmap_fptr;
- 
-  hid_t       group_id, dset_id;
-  hid_t       float_type_id, FLOAT_type_id;
-  hid_t       file_type_id, FILE_type_id;
-  hid_t       file_dsp_id;
- 
-  hsize_t     OutDims[MAX_DIMENSION];
-  hsize_t     TempIntArray[1];
- 
-  herr_t      h5_status;
-  herr_t      h5_error = -1;
- 
-  char *ParticlePositionLabel[] =
-     {"particle_position_x", "particle_position_y", "particle_position_z"};
-  char *ParticleVelocityLabel[] =
-     {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"};
-#ifdef WINDS
-  char *ParticleAttributeLabel[] =
-    {"creation_time", "dynamical_time", "metallicity_fraction", "particle_jet_x", 
-     "particle_jet_y", "particle_jet_z", "typeia_fraction"};
-#else
-  char *ParticleAttributeLabel[] = 
-    {"creation_time", "dynamical_time", "metallicity_fraction", "typeia_fraction"};
-#endif
-  char *SmoothedDMLabel[] = {"Dark_Matter_Density", "Velocity_Dispersion",
-			     "Particle_x-velocity", "Particle_y-velocity",
-			     "Particle_z-velocity"};
-  char *GriddedSPLabel[] = {"Star_Particle_Density", "Forming_Stellar_Mass_Density",
-			    "SFR_Density", "Average_creation_time"};
-#ifdef IO_LOG
-  int         io_log = 1;
-#else
-  int         io_log = 0;
-#endif
-
-  /* initialize */
- 
-  char id[MAX_GROUP_TAG_SIZE];
-  sprintf(id, "%"GROUP_TAG_FORMAT""ISYM, grid_id);
- 
-  /* make sure quantities defined at least for 3d */
- 
-  for (dim = GridRank; dim < 3; dim++) {
-    GridDimension[dim] = 1;
-    GridStartIndex[dim] = 0;
-    GridEndIndex[dim] = 0;
-  }
- 
-  for (dim = 0; dim < 3; dim++)
-    ActiveDim[dim] = GridEndIndex[dim] - GridStartIndex[dim] +1;
- 
-  /* ------------------------------------------------------------------- */
-  /* 1) Save general grid class data */
- 
-  if (MyProcessorNumber == ROOT_PROCESSOR && HierarchyFileOutputFormat > 0) {
-
-    fprintf(fptr, "Task              = %"ISYM"\n", ProcessorNumber);
- 
-    fprintf(fptr, "GridRank          = %"ISYM"\n", GridRank);
- 
-    fprintf(fptr, "GridDimension     = ");
-    WriteListOfInts(fptr, GridRank, GridDimension);
- 
-    fprintf(fptr, "GridStartIndex    = ");
-    WriteListOfInts(fptr, GridRank, GridStartIndex);
- 
-    fprintf(fptr, "GridEndIndex      = ");
-    WriteListOfInts(fptr, GridRank, GridEndIndex);
- 
-    fprintf(fptr, "GridLeftEdge      = ");
-    WriteListOfFloats(fptr, GridRank, GridLeftEdge);
- 
-    fprintf(fptr, "GridRightEdge     = ");
-    WriteListOfFloats(fptr, GridRank, GridRightEdge);
- 
-    fprintf(fptr, "Time              = %"GOUTSYM"\n", Time);
- 
-    fprintf(fptr, "SubgridsAreStatic = %"ISYM"\n", SubgridsAreStatic);
- 
-    fprintf(fptr, "NumberOfBaryonFields = %"ISYM"\n", NumberOfBaryonFields);
- 
-  }
- 
-  char pid[MAX_TASK_TAG_SIZE];
-  sprintf(pid, "%"TASK_TAG_FORMAT""ISYM, MyProcessorNumber);
- 
-  char gpid[MAX_TASK_TAG_SIZE];
-  sprintf(gpid, "%"TASK_TAG_FORMAT""ISYM, ProcessorNumber);
- 
-  char *groupfilename = new char[MAX_LINE_LENGTH];
-  strcpy(groupfilename, base_name);
-  strcat(groupfilename, ".cpu");
-  strcat(groupfilename, pid);
- 
-  char *procfilename = new char[MAX_LINE_LENGTH];
-  strcpy(procfilename, base_name);
-  strcat(procfilename, ".cpu");
-  strcat(procfilename, gpid);
- 
-  char *name = new char[MAX_LINE_LENGTH];
-  strcpy(name, "/Grid");
-  strcat(name, id);
- 
-  if (MyProcessorNumber == ProcessorNumber)
-  {
-    char *logname = new char[MAX_LINE_LENGTH];
-    strcpy(logname, groupfilename);
-    strcat(logname, ".log");
-    if (io_log) log_fptr = fopen(logname, "a");
-    delete [] logname;
- 
-    if (io_log) fprintf(log_fptr, "Grid_WriteGrid\n");
-    if (io_log) fprintf(log_fptr, "  ID %"ISYM"  %s\n", grid_id, id);
-    if (io_log) fprintf(log_fptr, "  File %s\n", groupfilename);
-    if (io_log) fprintf(log_fptr, "  Group %s\n", name);
-
-/* 
-    char *procmap = new char[MAX_LINE_LENGTH];
-    strcpy(procmap, base_name);
-    strcat(procmap, ".procmap");
-    procmap_fptr = fopen(procmap, "a");
-    fprintf(procmap_fptr, "%8"ISYM"  %s  Grid%8.8"ISYM"\n", grid_id, procfilename, grid_id);
-    fclose(procmap_fptr);
-    delete [] procmap;
-*/
- 
-  }
- 
-  /* Open HDF file for writing. */
- 
-  if (MyProcessorNumber == ProcessorNumber)
-  {
- 
-//    if (io_log) fprintf(log_fptr, "H5Fopen group file %s\n", groupfilename);
- 
-//    file_id = H5Fopen(groupfilename, H5F_ACC_RDWR, H5P_DEFAULT);
-//      if( file_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-    if (io_log) fprintf(log_fptr,"H5Gcreate with Name = %s\n",name);
- 
-    group_id = H5Gcreate(file_id, name, 0);
-      if (io_log) fprintf(log_fptr, "H5Gcreate id: %"ISYM"\n", group_id);
-      if( group_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-  }
- 
-  /* ------------------------------------------------------------------- */
-  /* 2) save baryon field quantities (including fields). */
- 
-  int ii = sizeof(io_type);
- 
-  switch(ii)
-  {
- 
-    case 4:
-      float_type_id = HDF5_R4;
-      file_type_id = HDF5_FILE_R4;
-      break;
- 
-    case 8:
-      float_type_id = HDF5_R8;
-      file_type_id = HDF5_FILE_R8;
-      break;
- 
-    default:
-      float_type_id = HDF5_R4;
-      file_type_id = HDF5_FILE_R4;
- 
-  }
- 
-  int jj = sizeof(FLOAT);
- 
-  switch(jj)
-  {
- 
-    case 4:
-      FLOAT_type_id = HDF5_R4;
-      FILE_type_id = HDF5_FILE_R4;
-      break;
- 
-    case 8:
-      FLOAT_type_id = HDF5_R8;
-      FILE_type_id = HDF5_FILE_R8;
-      break;
- 
-    case 16:
-      FLOAT_type_id = HDF5_R16;
-      FILE_type_id = H5Tcopy(HDF5_FILE_B8);
-                     H5Tset_size(FILE_type_id,16);
-      break;
- 
-    default:
-      printf("INCORRECT FLOAT DEFINITION\n");
- 
-  }
- 
-  if (NumberOfBaryonFields > 0) {
- 
-    if (MyProcessorNumber == ROOT_PROCESSOR && HierarchyFileOutputFormat > 0) {
- 
-      fprintf(fptr, "FieldType = ");
- 
-      WriteListOfInts(fptr, NumberOfBaryonFields, FieldType);
- 
-      fprintf(fptr, "BaryonFileName = %s\n", procfilename);
- 
-      fprintf(fptr, "CourantSafetyNumber    = %"FSYM"\n", CourantSafetyNumber);
-      fprintf(fptr, "PPMFlatteningParameter = %"ISYM"\n", PPMFlatteningParameter);
-      fprintf(fptr, "PPMDiffusionParameter  = %"ISYM"\n", PPMDiffusionParameter);
-      fprintf(fptr, "PPMSteepeningParameter = %"ISYM"\n", PPMSteepeningParameter);
- 
-    }
- 
-    if (MyProcessorNumber == ProcessorNumber) {
- 
-    /* 2a) Set HDF file dimensions (use FORTRAN ordering). */
- 
-    for (dim = 0; dim < GridRank; dim++)
-      OutDims[GridRank-dim-1] = ActiveDim[dim];
- 
-    /* 2b) Write out co-ordinate values.  Use the centre of each cell. */
- 
-    size = 1;
-    io_type *tempdim[MAX_DIMENSION];
- 
-    for (dim = GridRank-1; dim >= 0; dim--) {
- 
-      /* Compute cell centers and put them in temp. */
- 
-      tempdim[dim] = new io_type[GridDimension[dim]];
-      for (i = 0; i <= GridEndIndex[dim] - GridStartIndex[dim]; i++)
-	tempdim[dim][i] = CellLeftEdge[dim][GridStartIndex[dim] + i] +
-	          0.5 * CellWidth[dim][GridStartIndex[dim] + i];
-      size *= GridDimension[dim];
-    }
- 
-    /* create temporary buffer */
- 
-    temp = new io_type[size];
- 
-    /* 2c) Loop over fields, writing each one. */
- 
-    for (field = 0; field < NumberOfBaryonFields; field++) {
-
-      /* copy active part of field into grid */
- 
-      for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++)
-	for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++)
-	  for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++)
-	    temp[(i-GridStartIndex[0])                           +
-	         (j-GridStartIndex[1])*ActiveDim[0]              +
-	         (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] =
-		       io_type(
-	      BaryonField[field][i + j*GridDimension[0] +
-		                     k*GridDimension[0]*GridDimension[1]]
-                              );
- 
- 
-      file_dsp_id = H5Screate_simple((Eint32) GridRank, OutDims, NULL);
-        if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-        if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      if (io_log) fprintf(log_fptr,"H5Dcreate with Name = %s\n",DataLabel[field]);
- 
-      dset_id =  H5Dcreate(group_id, DataLabel[field], file_type_id, file_dsp_id, H5P_DEFAULT);
-        if (io_log) fprintf(log_fptr, "H5Dcreate id: %"ISYM"\n", dset_id);
-        if( dset_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      /* set datafield name and units, etc. */
- 
-      if ( DataUnits[field] == NULL )
-      {
-        DataUnits[field] = "none";
-      }
- 
-      //printf("OutPut Field2: %d \n", field);
-      WriteStringAttr(dset_id, "Label", DataLabel[field], log_fptr);
-      WriteStringAttr(dset_id, "Units", DataUnits[field], log_fptr);
-      WriteStringAttr(dset_id, "Format", "e10.4", log_fptr);
-      WriteStringAttr(dset_id, "Geometry", "Cartesian", log_fptr);
- 
- 
-      h5_status = H5Dwrite(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp);
-        if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Sclose(file_dsp_id);
-        if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Dclose(dset_id);
-        if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-    }   // end of loop over fields
- 
-    /* If requested, compute and output the temperature field 
-       as well since its such a pain to compute after the fact. */
- 
-    if (OutputTemperature) {
- 
-      /* Allocate field and compute temperature. */
- 
-      temperature = new float[size];
- 
-      if (this->ComputeTemperatureField(temperature) == FAIL) {
-	ENZO_FAIL("Error in grid->ComputeTemperatureField.\n");
-      }
- 
-      /* Copy active part of field into grid */
- 
-      for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++)
-	for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++)
-	  for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++)
-	    temp[(i-GridStartIndex[0])                           +
-	         (j-GridStartIndex[1])*ActiveDim[0]              +
-	         (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] =
-		     io_type(
-		   temperature[(k*GridDimension[1] + j)*GridDimension[0] + i]
-			     );
- 
- 
-      file_dsp_id = H5Screate_simple((Eint32) GridRank, OutDims, NULL);
-        if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-        if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      if (io_log) fprintf(log_fptr,"H5Dcreate with Name = Temperature\n");
- 
-      dset_id = H5Dcreate(group_id, "Temperature", file_type_id, file_dsp_id, H5P_DEFAULT);
-        if (io_log) fprintf(log_fptr, "H5Dcreate id: %"ISYM"\n", dset_id);
-        if( dset_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      if ( DataUnits[field] == NULL )
-      {
-        DataUnits[field] = "none";
-      }
- 
-      WriteStringAttr(dset_id, "Label", "Temperature", log_fptr);
-      WriteStringAttr(dset_id, "Units", "K", log_fptr);
-      WriteStringAttr(dset_id, "Format", "e10.4", log_fptr);
-      WriteStringAttr(dset_id, "Geometry", "Cartesian", log_fptr);
- 
-      h5_status = H5Dwrite(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp);
-        if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Sclose(file_dsp_id);
-        if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Dclose(dset_id);
-        if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-	// If outputing dust temperature, keep temperature field for the calculation.
-	if (!OutputDustTemperature) {
-	  delete [] temperature;
-	}
- 
-    } // end: if (OutputTemperature)
-
-
-    /* If requested, compute and output the dust temperature field 
-       as well since its such a pain to compute after the fact. */
-
-    if (OutputDustTemperature) {
- 
-      /* Get temperature field if we do not already have it. */
-
-      if (!OutputTemperature) {
-	temperature = new float[size];
-
-	if (this->ComputeTemperatureField(temperature) == FAIL) {
-	  ENZO_FAIL("Error in grid->ComputeTemperatureField.\n");
-	}
-      }
-
-      /* Allocate field and compute dust temperature. */
- 
-      dust_temperature = new float[size];
- 
-      if (this->ComputeDustTemperatureField(temperature, 
-					    dust_temperature) == FAIL) {
-	ENZO_FAIL("Error in grid->ComputeDustTemperatureField.\n");
-      }
- 
-      /* Copy active part of field into grid */
- 
-      for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++)
-	for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++)
-	  for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++)
-	    temp[(i-GridStartIndex[0])                           +
-	         (j-GridStartIndex[1])*ActiveDim[0]              +
-	         (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] =
-		     io_type(
-		   dust_temperature[(k*GridDimension[1] + j)*GridDimension[0] + i]
-			     );
- 
- 
-      file_dsp_id = H5Screate_simple((Eint32) GridRank, OutDims, NULL);
-        if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-        if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      if (io_log) fprintf(log_fptr,"H5Dcreate with Name = Dust_Temperature\n");
- 
-      dset_id = H5Dcreate(group_id, "Dust_Temperature", file_type_id, file_dsp_id, H5P_DEFAULT);
-        if (io_log) fprintf(log_fptr, "H5Dcreate id: %"ISYM"\n", dset_id);
-        if( dset_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      if ( DataUnits[field] == NULL )
-      {
-        DataUnits[field] = "none";
-      }
- 
-      WriteStringAttr(dset_id, "Label", "Dust_Temperature", log_fptr);
-      WriteStringAttr(dset_id, "Units", "K", log_fptr);
-      WriteStringAttr(dset_id, "Format", "e10.4", log_fptr);
-      WriteStringAttr(dset_id, "Geometry", "Cartesian", log_fptr);
- 
-      h5_status = H5Dwrite(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp);
-        if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Sclose(file_dsp_id);
-        if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Dclose(dset_id);
-        if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      if (!OutputTemperature) {
-	delete [] temperature;
-      }
-      delete [] dust_temperature;
-    
-    } // end: if (OutputDustTemperature)
-
-
-    if (VelAnyl) {
-      int Vel1Num, Vel2Num, Vel3Num;
-      
-      for (int i=0; i< NumberOfBaryonFields; i++){
-	switch (FieldType[i]){
-	case Velocity1:
-	  Vel1Num=i; break;
-	case Velocity2:
-	  Vel2Num=i; break;
-	case Velocity3:
-	  Vel3Num=i; break;
-	}
-      }
-
-
-      io_type *curlx; io_type *curly;
-
-      
-      if(GridRank==3){
-      curlx = new io_type [size];
-      curly = new io_type [size];}
-
-
-      io_type *curlz = new io_type [size];
-      io_type *div   = new io_type [size];
-
-      FLOAT dx = CellWidth[0][0],
-	dy = CellWidth[1][0], dz;
-      if (GridRank>2)
-	dz = CellWidth[2][0];
-    
-      /* Copy active part of field into grid */
-      int igrid, igridyp1, igridym1, igridzp1, igridzm1;
-      for (int k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) {
-	for (int j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) {
-	  for (int i = GridStartIndex[0]; i <= GridEndIndex[0]; i++) {
-
-	    igrid = (k*GridDimension[1] + j)*GridDimension[0] + i;
-	    igridyp1 = (k*GridDimension[1] + j + 1)*GridDimension[0] + i;
-	    igridym1 = (k*GridDimension[1] + j - 1)*GridDimension[0] + i;
-	    igridzp1 = ((k+1)*GridDimension[1]+j)*GridDimension[0] + i;
-	    igridzm1 = ((k-1)*GridDimension[1]+j)*GridDimension[0] + i;
-
-
-
-	    if (GridRank==3){
-	      
-	    div[(i-GridStartIndex[0])+
-		(j-GridStartIndex[1])*ActiveDim[0]+
-		(k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] =
-              io_type(
-		      (0.5*(BaryonField[Vel1Num][igrid+1]-BaryonField[Vel1Num][igrid-1])/dx +
-		      0.5*(BaryonField[Vel2Num][igridyp1]-BaryonField[Vel2Num][igridym1])/dy +
-			  0.5*(BaryonField[Vel3Num][igridzp1]-BaryonField[Vel3Num][igridzm1])/dz)
-                     );
-	    curlx[(i-GridStartIndex[0])+
-	       (j-GridStartIndex[1])*ActiveDim[0]+ 
-	       (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] =
-	       io_type(
-		      (0.5*(BaryonField[Vel3Num][igridyp1]-BaryonField[Vel3Num][igridym1])/dy -		      
-			  0.5*(BaryonField[Vel2Num][igridzp1]-BaryonField[Vel2Num][igridzm1])/dz)
-		      );
-	    curly[(i-GridStartIndex[0])+
-	       (j-GridStartIndex[1])*ActiveDim[0]+ 
-	       (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] =
-	      io_type(
-		       (0.5*(BaryonField[Vel1Num][igridzp1]-BaryonField[Vel1Num][igridzm1])/dz -		      
-			   0.5*(BaryonField[Vel3Num][igrid+1]-BaryonField[Vel3Num][igrid-1])/dx)
-		       );
-	     curlz[(i-GridStartIndex[0])+
-		(j-GridStartIndex[1])*ActiveDim[0]+ 
-		(k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] =
-	      io_type(
-		      (0.5*(BaryonField[Vel2Num][igrid+1]-BaryonField[Vel2Num][igrid-1])/dx -		      
-			  0.5*(BaryonField[Vel1Num][igridyp1]-BaryonField[Vel1Num][igridym1])/dy)
-		      );
-	    }
-
-	    if (GridRank==2){
-	      
-	    div[(i-GridStartIndex[0])+
-		(j-GridStartIndex[1])*ActiveDim[0]+
-		(k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] =
-              io_type(
-		      (0.5*(BaryonField[Vel1Num][igrid+1]-BaryonField[Vel1Num][igrid-1])/dx +
-		       0.5*(BaryonField[Vel2Num][igridyp1]-BaryonField[Vel2Num][igridym1])/dy 
-		       ));
-	     curlz[(i-GridStartIndex[0])+
-		(j-GridStartIndex[1])*ActiveDim[0]+ 
-		(k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] =
-	      io_type(
-		      (0.5*(BaryonField[Vel2Num][igrid+1]-BaryonField[Vel2Num][igrid-1])/dx -		      
-			  0.5*(BaryonField[Vel1Num][igridyp1]-BaryonField[Vel1Num][igridym1])/dy)
-		      );
-	    }
-
-	    
-
-	  }
-	}
-      }
-
-      /* output */
-      
-      char *DataLabelN[4];
-      if (GridRank==2) {
-	DataLabelN[0]="Velocity_Div";
-	DataLabelN[1]="Velocity_Vorticity";
-      }
-      if (GridRank==3) {
-	DataLabelN[0]="Velocity_Div";
-	DataLabelN[1]="Velocity_Vorticity1";
-	DataLabelN[2]="Velocity_Vorticity2";
-	DataLabelN[3]="Velocity_Vorticity3";
-      }
-
-      
-
-      int tFields=4;
-      if (GridRank==2) tFields=2;
-
-      for (int field=0; field<tFields; field++){
-
-      file_dsp_id = H5Screate_simple((Eint32) GridRank, OutDims, NULL);
-      if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-      if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      if (io_log) fprintf(log_fptr,"H5Dcreate with Name = %s\n",DataLabelN[field]);
- 
-      dset_id =  H5Dcreate(file_id, DataLabelN[field], file_type_id, file_dsp_id, H5P_DEFAULT);
-        if (io_log) fprintf(log_fptr, "H5Dcreate id: %"ISYM"\n", dset_id);
-        if( dset_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      /* set datafield name and units, etc. */
- 
-      if ( DataUnits[field] == NULL )
-      {
-        DataUnits[field] = "none";
-      }
- 
-      
-      WriteStringAttr(dset_id, "Label", DataLabelN[field], log_fptr);
-      WriteStringAttr(dset_id, "Units", NULL, log_fptr);
-      WriteStringAttr(dset_id, "Format", "e10.4", log_fptr);
-      WriteStringAttr(dset_id, "Geometry", "Cartesian", log_fptr);
-
-      switch (field) {
-      case 0:
-	temp_VelAnyl=div;
-	break;
-      case 1:
-	temp_VelAnyl=curlz;
-	break;
-      case 2:
-	temp_VelAnyl=curly;
-	break;
-      case 3:
-	temp_VelAnyl=curlx;
-	break;
-      }
- 
-      h5_status = H5Dwrite(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp_VelAnyl);
-        if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Sclose(file_dsp_id);
-        if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Dclose(dset_id);
-        if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
-
-
-      }
-
-
-
-
-      if(GridRank==3){
-	delete curlx;
-	delete curly;}
-      delete curlz;
-      delete div;
-    }
-
-
-    if (OutputCoolingTime != FALSE) {
- 
-      /* Allocate field and compute cooling time. */
-
-      cooling_time = new float[size];
- 
-      float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1,
-	VelocityUnits = 1, TimeUnits = 1, aUnits = 1;
-
-      GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
-	       &TimeUnits, &VelocityUnits, Time);
-
-      if (this->ComputeCoolingTime(cooling_time) == FAIL) {
-	ENZO_FAIL("Error in grid->ComputeCoolingTime.\n");
-      }
-
-      // Make all cooling time values positive and convert to seconds.
-      for (i = 0;i < size;i++) {
-	cooling_time[i] = fabs(cooling_time[i]) * TimeUnits;
-      }
- 
-      /* Copy active part of field into grid */
- 
-      for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++)
-	for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++)
-	  for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++)
-	    temp[(i-GridStartIndex[0])                           +
-	         (j-GridStartIndex[1])*ActiveDim[0]              +
-	         (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] =
-		     io_type(
-		   cooling_time[(k*GridDimension[1] + j)*GridDimension[0] + i]
-			     );
- 
-      file_dsp_id = H5Screate_simple((Eint32) GridRank, OutDims, NULL);
-        if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-        if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      if (io_log) fprintf(log_fptr,"H5Dcreate with Name = Cooling_Time\n");
- 
-      dset_id = H5Dcreate(group_id, "Cooling_Time", file_type_id, file_dsp_id, H5P_DEFAULT);
-        if (io_log) fprintf(log_fptr, "H5Dcreate id: %"ISYM"\n", dset_id);
-        if( dset_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      if ( DataUnits[field] == NULL )
-      {
-        DataUnits[field] = "none";
-      }
- 
-      WriteStringAttr(dset_id, "Label", "Temperature", log_fptr);
-      WriteStringAttr(dset_id, "Units", "K", log_fptr);
-      WriteStringAttr(dset_id, "Format", "e10.4", log_fptr);
-      WriteStringAttr(dset_id, "Geometry", "Cartesian", log_fptr);
- 
-      h5_status = H5Dwrite(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp);
-        if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Sclose(file_dsp_id);
-        if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Dclose(dset_id);
-        if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      delete [] cooling_time;
- 
-    } // if (OutputCoolingTime)
-
-    /* Make sure that there is a copy of dark matter field to save
-       (and at the right resolution). */
-
-    if (OutputSmoothedDarkMatter == FALSE) {
-    
-
-      if (SelfGravity && NumberOfParticles > 0) {
-	float SaveGravityResolution = GravityResolution;
-	GravityResolution = 1;
-	this->InitializeGravitatingMassFieldParticles(RefineBy);
-	this->ClearGravitatingMassFieldParticles();
-	this->DepositParticlePositions(this, Time,
-				       GRAVITATING_MASS_FIELD_PARTICLES);
-	GravityResolution = SaveGravityResolution;
-      }
-
-      /* If present, write out the GravitatingMassFieldParticles. */
- 
-      if (GravitatingMassFieldParticles != NULL) {
- 
-	/* Set dimensions. */
- 
-	int StartIndex[] = {0,0,0}, EndIndex[] = {0,0,0};
-	for (dim = 0; dim < GridRank; dim++) {
-	  StartIndex[dim] = nint((GridLeftEdge[dim] -
-				  GravitatingMassFieldParticlesLeftEdge[dim])/
-				 GravitatingMassFieldParticlesCellSize);
-	  EndIndex[dim] = nint((GridRightEdge[dim] -
-				GravitatingMassFieldParticlesLeftEdge[dim])/
-			       GravitatingMassFieldParticlesCellSize) - 1;
-//      fprintf(stderr, "%"ISYM" %"ISYM" %10.4"FSYM" %10.4"FSYM" %10.4"FSYM" %10.4"FSYM"\n", StartIndex[dim], EndIndex[dim], GridLeftEdge[dim], GridRightEdge[dim], GravitatingMassFieldParticlesLeftEdge[dim], GravitatingMassFieldParticlesCellSize);
-	}
- 
-	/* Copy active part of field into grid */
- 
-	for (k = StartIndex[2]; k <= EndIndex[2]; k++)
-	  for (j = StartIndex[1]; j <= EndIndex[1]; j++)
-	    for (i = StartIndex[0]; i <= EndIndex[0]; i++)
-	      temp[(i-StartIndex[0])                           +
-		   (j-StartIndex[1])*ActiveDim[0]              +
-		   (k-StartIndex[2])*ActiveDim[0]*ActiveDim[1] ] =
-		io_type(
-			GravitatingMassFieldParticles[ i +
-			j*GravitatingMassFieldParticlesDimension[0] +
-			k*GravitatingMassFieldParticlesDimension[0]*
-			GravitatingMassFieldParticlesDimension[1]]
-			);
- 
- 
-	file_dsp_id = H5Screate_simple((Eint32) GridRank, OutDims, NULL);
-        if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-        if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-	if (io_log) fprintf(log_fptr,"H5Dcreate with Name = Dark_Matter_Density\n");
- 
-	dset_id =  H5Dcreate(group_id, "Dark_Matter_Density", file_type_id, file_dsp_id, H5P_DEFAULT);
-        if (io_log) fprintf(log_fptr, "H5Dcreate id: %"ISYM"\n", dset_id);
-        if( dset_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-	if ( DataUnits[field] == NULL )
-	  {
-	    DataUnits[field] = "none";
-	  }
- 
-	WriteStringAttr(dset_id, "Label", "Dark_Matter_Density", log_fptr);
-	WriteStringAttr(dset_id, "Units", "", log_fptr);
-	WriteStringAttr(dset_id, "Format", "e10.4", log_fptr);
-	WriteStringAttr(dset_id, "Geometry", "Cartesian", log_fptr);
- 
-	h5_status = H5Dwrite(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp);
-        if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-	h5_status = H5Sclose(file_dsp_id);
-        if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-	h5_status = H5Dclose(dset_id);
-        if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-	/* Clean up if we modified the resolution. */
- 
-	if (SelfGravity && GravityResolution != 1)
-	  this->DeleteGravitatingMassFieldParticles();
- 
-      } // end of (if GravitatingMassFieldParticles != NULL)
-
-    } // ENDIF !OutputSmoothedDarkMatter
-
-    delete [] temp;
- 
-    for (dim = 0; dim < GridRank; dim++)
-      delete [] tempdim[dim];
- 
-    /* Write BoundaryFluxes info (why? it's just recreated when the grid
-                                  is read in) */
- 
-   }  // end: if (ProcessorNumber == MyProcessorNumber)
-  } // end: if (NumberOfBaryonFields > 0)
-
-  /* ------------------------------------------------------------------- */
-  /* 3) Save particle quantities smoothed (or gridded) to the grid. */
- 
-  if (OutputSmoothedDarkMatter > 0 && MyProcessorNumber == ProcessorNumber) {
-
-    size = active_size = 1;
-    for (dim = 0; dim < GridRank; dim++) {
-      OutDims[GridRank-dim-1] = ActiveDim[dim];
-      size *= GridDimension[dim];
-      active_size *= ActiveDim[dim];
-    }
- 
-    temp = new io_type[active_size];
-
-    int NumberOfDMFields;
-    switch (OutputSmoothedDarkMatter) {
-    case 1: NumberOfDMFields = 1; break;  // density
-    case 2: NumberOfDMFields = 5; break;  // + rms velocity + 3-velocity
-    } // ENDSWITCH
-      
-    for (field = 0; field < NumberOfDMFields; field++) {
-
-      // Only the active part was calculated, so just copy over.
-      for (i = 0; i < active_size; i++)
-	temp[i] = io_type(InterpolatedField[field][i]);
-
-      file_dsp_id = H5Screate_simple((Eint32) GridRank, OutDims, NULL);
-      if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-      if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      if (io_log) fprintf(log_fptr,"H5Dcreate with Name = %s\n", SmoothedDMLabel[field]);
- 
-      dset_id = H5Dcreate(group_id, SmoothedDMLabel[field], file_type_id, file_dsp_id, H5P_DEFAULT);
-      if (io_log) fprintf(log_fptr, "H5Dcreate id: %"ISYM"\n", dset_id);
-      if( dset_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      WriteStringAttr(dset_id, "Label", SmoothedDMLabel[field], log_fptr);
-      WriteStringAttr(dset_id, "Units", "", log_fptr);
-      WriteStringAttr(dset_id, "Format", "e10.4", log_fptr);
-      WriteStringAttr(dset_id, "Geometry", "Cartesian", log_fptr);
- 
-      h5_status = H5Dwrite(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, 
-			   (VOIDP) temp);
-      if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Sclose(file_dsp_id);
-      if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Dclose(dset_id);
-      if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      delete [] InterpolatedField[field];
-      InterpolatedField[field] = NULL;
-
-    } // ENDFOR field
-
-    delete [] temp;
-      
-  } // ENDIF OutputSmoothedDarkMatter
-
-  if (OutputGriddedStarParticle > 0 && MyProcessorNumber == ProcessorNumber) {
-
-#define NumberOfInterpolatedFieldsForDM 10
-
-    size = active_size = 1;
-    for (dim = 0; dim < GridRank; dim++) {
-      OutDims[GridRank-dim-1] = ActiveDim[dim];
-      size *= GridDimension[dim];
-      active_size *= ActiveDim[dim];
-    }
- 
-    temp = new io_type[active_size];
-
-    // Assign number of fields 
-
-    int NumberOfSPFields;
-    switch (OutputGriddedStarParticle) {
-    case 1: NumberOfSPFields = 1; break;  // particle density
-    case 2: NumberOfSPFields = 4; break;  // + forming star particle density + SFR density, etc. 
-    default: 
-      fprintf(stdout, "Unrecognized value for OutputGriddedStarParticle = %"ISYM"\n",
-	      OutputGriddedStarParticle);
-      fprintf(stdout, "Setting to 1.  Outputting particle density only.\n");
-      OutputGriddedStarParticle = 1;
-      NumberOfSPFields = 1;
-      break;
-    } 
-
-    // Get gridded star particle field
-    if (this->InterpolateStarParticlesToGrid(NumberOfSPFields) == FAIL) {
-      ENZO_FAIL("Error in grid->InterpolateStarParticlesToGrid.\n");
-    }
-
-    for (field = NumberOfInterpolatedFieldsForDM; 
-	 field < NumberOfInterpolatedFieldsForDM+NumberOfSPFields; field++) {
-
-      // Only the active part was calculated, so just copy over.
-      for (i = 0; i < active_size; i++)
-	temp[i] = io_type(InterpolatedField[field][i]);
-
-      file_dsp_id = H5Screate_simple((Eint32) GridRank, OutDims, NULL);
-      if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-      if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      if (io_log) fprintf(log_fptr,"H5Dcreate with Name = %s\n", GriddedSPLabel[field-NumberOfInterpolatedFieldsForDM]);
-
-      dset_id = H5Dcreate(group_id, GriddedSPLabel[field-NumberOfInterpolatedFieldsForDM], file_type_id, file_dsp_id, H5P_DEFAULT);
-      if (io_log) fprintf(log_fptr, "H5Dcreate id: %"ISYM"\n", dset_id);
-      if( dset_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      WriteStringAttr(dset_id, "Label", GriddedSPLabel[field-NumberOfInterpolatedFieldsForDM], log_fptr);
-      WriteStringAttr(dset_id, "Units", "", log_fptr);
-      WriteStringAttr(dset_id, "Format", "e10.4", log_fptr);
-      WriteStringAttr(dset_id, "Geometry", "Cartesian", log_fptr);
- 
-      h5_status = H5Dwrite(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp);
-      if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Sclose(file_dsp_id);
-      if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Dclose(dset_id);
-      if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      delete [] InterpolatedField[field];
-      InterpolatedField[field] = NULL;
-
-    } // ENDFOR field
-
-    delete [] temp;
-      
-  } // ENDIF OutputGriddedStarParticle
- 
-  /* ------------------------------------------------------------------- */
-  /* 4) Save particle quantities. */
- 
-  if (MyProcessorNumber == ROOT_PROCESSOR && HierarchyFileOutputFormat > 0)
-    fprintf(fptr, "NumberOfParticles   = %"ISYM"\n", NumberOfParticles);
-  if (MyProcessorNumber == ROOT_PROCESSOR && 
-      HierarchyFileOutputFormat > 0 && 
-      OutputSmoothedDarkMatter > 0 && NumberOfParticles == 0)
-    fprintf(fptr, "ParticleFileName = %s\n", procfilename);
- 
-  if (NumberOfParticles > 0) {
- 
-    if (MyProcessorNumber == ROOT_PROCESSOR && HierarchyFileOutputFormat > 0)
-      fprintf(fptr, "ParticleFileName = %s\n", procfilename); // must be same as above
- 
-    if (MyProcessorNumber == ProcessorNumber) {
- 
-    /* Sort particles according to their identifier. */
- 
-    if (OutputParticleTypeGrouping)
-      this->SortParticlesByType();
-    else
-      this->SortParticlesByNumber();
- 
-    /* Create a temporary buffer (64 bit). */
- 
-    temp = new io_type[NumberOfParticles];
- 
-    /* Particle positions are not converted to 32 bit first.
-       (128 bit numbers are not supported by HDF so convert to 64). */
- 
-    io_type *temp_pointer = NULL;
-    float128 *long_temp_pointer = NULL;
- 
-    TempIntArray[0] = NumberOfParticles;
-
-    for (dim = 0; dim < GridRank; dim++) {
- 
-      /* Convert to 64 if 128, either just write out. */
- 
-      if (sizeof(FLOAT) == 16) {
-        long_temp_pointer = (float128*) ParticlePosition[dim];
-      }
-      else
-      {
-	temp_pointer = (io_type*) ParticlePosition[dim];
-      }
- 
- 
-      file_dsp_id = H5Screate_simple((Eint32) 1, TempIntArray, NULL);
-        if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-        if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      if (io_log) fprintf(log_fptr,"H5Dcreate with Name = %s\n", ParticlePositionLabel[dim]);
- 
-      dset_id =  H5Dcreate(group_id, ParticlePositionLabel[dim],  FILE_type_id, file_dsp_id, H5P_DEFAULT);
-        if (io_log) fprintf(log_fptr, "H5Dcreate id: %"ISYM"\n", dset_id);
-        if( dset_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      if (sizeof(FLOAT) == 16) {
-//                                  NOTE: for 128bits this must be FILE_type_id and NOT FLOAT_type_id!
-      h5_status = H5Dwrite(dset_id, FILE_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) long_temp_pointer);
-        if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      }
-      else
-      {
- 
-      h5_status = H5Dwrite(dset_id, FLOAT_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp_pointer);
-        if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      }
- 
-      h5_status = H5Sclose(file_dsp_id);
-        if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Dclose(dset_id);
-        if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-    }
- 
-//    if (sizeof(FLOAT) == 16)
-//      delete [] long_temp_pointer;  /* clean up if allocated. */
- 
-    /* Copy particle velocities to temp and write them. */
- 
-    for (dim = 0; dim < GridRank; dim++) {
- 
-      for (i = 0; i < NumberOfParticles; i++)
-	temp[i] = io_type(ParticleVelocity[dim][i]);
- 
-      file_dsp_id = H5Screate_simple((Eint32) 1, TempIntArray, NULL);
-        if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-        if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      if (io_log) fprintf(log_fptr,"H5Dcreate with Name = %s\n",ParticleVelocityLabel[dim]);
- 
-      dset_id =  H5Dcreate(group_id, ParticleVelocityLabel[dim], file_type_id, file_dsp_id, H5P_DEFAULT);
-        if (io_log) fprintf(log_fptr, "H5Dcreate id: %"ISYM"\n", dset_id);
-        if( dset_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Dwrite(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp);
-        if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Sclose(file_dsp_id);
-        if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-      h5_status = H5Dclose(dset_id);
-        if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-        if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-    }
- 
-    /* Copy mass to temp and write it. */
- 
-    for (i = 0; i < NumberOfParticles; i++)
-      temp[i] = io_type(ParticleMass[i]);
- 
- 
-    file_dsp_id = H5Screate_simple((Eint32) 1, TempIntArray, NULL);
-      if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-      if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-    if (io_log) fprintf(log_fptr,"H5Dcreate with Name = particle_mass\n");
- 
-    dset_id =  H5Dcreate(group_id, "particle_mass", file_type_id, file_dsp_id, H5P_DEFAULT);
-      if (io_log) fprintf(log_fptr, "H5Dcreate id: %"ISYM"\n", dset_id);
-      if( dset_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-    h5_status = H5Dwrite(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp);
-      if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-    h5_status = H5Sclose(file_dsp_id);
-      if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-    h5_status = H5Dclose(dset_id);
-      if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
-
-    /* Copy number (ID) to temp and write it. */
- 
-    PINT *tempPINT = new PINT[NumberOfParticles];
- 
-    for (i = 0; i < NumberOfParticles; i++)
-      tempPINT[i] = ParticleNumber[i];
- 
- 
-    file_dsp_id = H5Screate_simple((Eint32) 1, TempIntArray, NULL);
-      if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id);
-      if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-    if (io_log) fprintf(log_fptr,"H5Dcreate with Name = particle_index\n");
- 
-    dset_id =  H5Dcreate(group_id, "particle_index", HDF5_FILE_PINT, file_dsp_id, H5P_DEFAULT);
-      if (io_log) fprintf(log_fptr, "H5Dcreate id: %"ISYM"\n", dset_id);
-      if( dset_id == h5_error ){my_exit(EXIT_FAILURE);}
- 
-    h5_status = H5Dwrite(dset_id, HDF5_PINT, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) tempPINT);
-      if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-    h5_status = H5Sclose(file_dsp_id);
-      if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-    h5_status = H5Dclose(dset_id);
-      if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status);
-      if( h5_status == h5_error ){my_exit(EXIT_FAILURE);}
- 
-    /* Copy type to temp and write it. */
-
-    if (ParticleTypeInFile == TRUE) {
-