Automated merge with http://hg.spacepope.org/week-of-code/

--HG--
branch : week-of-code

src/enzo/Grid_Group_ReadGrid.C

@@ -335,8 +335,12 @@ int grid::Group_ReadGrid(FILE *fptr, int GridID, HDF5_hid_t file_id,
 
       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 ){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 ){my_exit(EXIT_FAILURE);}

src/enzo/Grid_ReadGrid.C

@@ -43,9 +43,8 @@ int ReadListOfInts(FILE *fptr, int N, int nums[]);
 // extern int ParticleTypeInFile; // declared and set in ReadParameterFile
 
 #ifdef USE_HDF4
-int ReadField(float *temp, int Dims[], int Rank, char *name,
-	      char *field_name);
-static int32 sd_id, sds_index; // HDF4 (SD) handlers                                               
+int ReadField(float *temp, int Dims[], int Rank, char *name, char *field_name);
+static Eint32 sd_id, sds_index; // HDF4 (SD) handlers                                               
 #endif 
 
 int grid::ReadGrid(FILE *fptr, int GridID, 
@@ -82,9 +81,9 @@ int grid::ReadGrid(FILE *fptr, int GridID,
                                     "metallicity_fraction", "alpha_fraction"};
 
 #ifdef USE_HDF4
-    int32 TempIntArray2[MAX_DIMENSION];
-    int32 sds_id, num_type2, attributes, TempInt;  
-    sds_index = 0;  // start at first SDS                                                                            
+  Eint32 TempIntArray2[MAX_DIMENSION];
+  Eint32 sds_id, num_type2, attributes, TempInt;  
+  sds_index = 0;  // start at first SDS                                                                            
 #endif 
 
 #ifdef IO_LOG
@@ -472,10 +471,10 @@ int grid::ReadGrid(FILE *fptr, int GridID,
 
 #ifdef USE_HDF4
       if (!TryHDF5) {
-	float *temp = new float[active_size];
+      TempIntArray2[0] = Eint32(NumberOfParticles);
 
-	int32 HDFDataType = (sizeof(Eflt) == 4) ? DFNT_FLOAT32 : DFNT_FLOAT64;
-	if (sizeof(Eflt) == 16) HDFDataType = DFNT_FLOAT128;
+      Eint32 HDFDataType = (sizeof(FLOAT) == 4) ? DFNT_FLOAT32 : DFNT_FLOAT64;  // DFNT_* are defined by HDF4
+      if (sizeof(FLOAT) == 16) HDFDataType = DFNT_FLOAT128;
 
 	/* Read dims.
 	   If Rank != 1, we may have just read some other field SDS.  If so,
@@ -506,6 +505,18 @@ int grid::ReadGrid(FILE *fptr, int GridID,
 	  return FAIL;  
 	}
 
+      /* Create a temporary buffer (32 bit or twice the size for 64). */
+
+      float32 *temp = NULL;
+      if (num_type2 != HDFDataType) {
+	if (num_type2 == DFNT_FLOAT64)
+	  temp = new float32[NumberOfParticles*2];
+	if (num_type2 == DFNT_FLOAT128)
+	  temp = new float32[NumberOfParticles*4];
+      }
+      if (temp == NULL)
+	temp = new float32[NumberOfParticles];
+
 	/* Read ParticlePosition (use temporary buffer). */ 
 
 	for (dim = 0; dim < GridRank; dim++) {
@@ -535,13 +546,13 @@ int grid::ReadGrid(FILE *fptr, int GridID,
 
 	    if (num_type2 == DFNT_FLOAT32)
 	      for (i = 0; i < NumberOfParticles; i++)
-		ParticlePosition[dim][i] = Eflt(temp[i]);
+		ParticlePosition[dim][i] = FLOAT(temp[i]);
 	    if (num_type2 == DFNT_FLOAT64)
 	      for (i = 0; i < NumberOfParticles; i++)
-		ParticlePosition[dim][i] = Eflt(temp64[i]);
+		ParticlePosition[dim][i] = FLOAT(temp64[i]);
 	    if (num_type2 == DFNT_FLOAT128)
 	      for (i = 0; i < NumberOfParticles; i++)
-		ParticlePosition[dim][i] = Eflt(temp128[i]);
+		ParticlePosition[dim][i] = FLOAT(temp128[i]);
 	  }
 	} // end: loop over dims
 

src/enzo/Grid_StarParticleHandler.C

@@ -50,7 +50,7 @@ extern "C" void FORTRAN_NAME(star_maker1)(int *nx, int *ny, int *nz,
              int *nmax, FLOAT *xstart, FLOAT *ystart, FLOAT *zstart,
      		 int *ibuff, hydro_method *imethod,
              float *odthresh, float *massff, float *smthrest, int *level,
-                 int *np,
+		 int *np,
              FLOAT *xp, FLOAT *yp, FLOAT *zp, float *up, float *vp, float *wp,
              float *mp, float *tdp, float *tcp);
 #endif /* STAR1 */
@@ -65,7 +65,7 @@ extern "C" void FORTRAN_NAME(star_maker2)(int *nx, int *ny, int *nz,
      		 int *ibuff,
              int *imetal, hydro_method *imethod, float *mintdyn,
              float *odthresh, float *massff, float *smthrest, int *level,
-                 int *np,
+		 int *np, 
              FLOAT *xp, FLOAT *yp, FLOAT *zp, float *up, float *vp, float *wp,
              float *mp, float *tdp, float *tcp, float *metalf);
 
@@ -80,7 +80,7 @@ extern "C" void FORTRAN_NAME(star_maker3)(int *nx, int *ny, int *nz,
      		 int *ibuff,
              int *imetal, hydro_method *imethod, float *mintdyn,
              float *odthresh, float *massff, float *smthrest, int *level,
-                 int *np,
+		 int *np, 
              FLOAT *xp, FLOAT *yp, FLOAT *zp, float *up, float *vp, float *wp,
              float *mp, float *tdp, float *tcp, float *metalf);
 
@@ -94,7 +94,7 @@ extern "C" void FORTRAN_NAME(star_maker4)(int *nx, int *ny, int *nz,
      		 int *ibuff, 
              int *imetal, hydro_method *imethod, float *mintdyn,
              float *odthresh, float *massff, float *smthrest, int *level,
-                 int *np,
+	         int *np, 
              FLOAT *xp, FLOAT *yp, FLOAT *zp, float *up, float *vp, float *wp,
              float *mp, float *tdp, float *tcp, float *metalf);
 
@@ -108,7 +108,7 @@ extern "C" void FORTRAN_NAME(star_maker4)(int *nx, int *ny, int *nz,
      		 int *ibuff,
              int *imetal, hydro_method *imethod, float *mintdyn,
              float *odthresh, float *massff, float *smthrest, int *level,
-                 int *np,
+		 int *np,
              FLOAT *xp, FLOAT *yp, FLOAT *zp, float *up, float *vp, float *wp,
              float *mp, float *tdp, float *tcp, float *metalf);
 
@@ -213,7 +213,7 @@ int sink_maker(int *nx, int *ny, int *nz, int *size,
              float *d1, float *x1, float *v1, float *t1,
              int *nmax, FLOAT *xstart, FLOAT *ystart, FLOAT *zstart, 
      		 int *ibuff, 
-	     int *imethod, float *massthresh, int *level, int *np,
+	       int *imethod, float *massthresh, int *level, int *np, 
              FLOAT *xp, FLOAT *yp, FLOAT *zp, 
 	     float *up, float *vp, float *wp, float *mp, 
 		 float *tcp, float *tdp, int *type,
@@ -430,6 +430,7 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level)
 
     /* Call FORTRAN routine to do the actual work. */
 
+    int NumberOfNewParticlesSoFar = 0;
     int NumberOfNewParticles = 0;
 
     if (debug && NumberOfNewParticles > 0) {
@@ -464,6 +465,8 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level)
 
       //---- MODIFIED SF ALGORITHM ("STANDARD VERSION")
 
+      NumberOfNewParticlesSoFar = NumberOfNewParticles;
+
       FORTRAN_NAME(star_maker2)(
        GridDimension, GridDimension+1, GridDimension+2,
        BaryonField[DensNum], dmfield, temperature, BaryonField[Vel1Num],
@@ -475,21 +478,23 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level)
           CellLeftEdge[2], &GhostZones,
        &MetallicityField, &HydroMethod, &StarMakerMinimumDynamicalTime,
        &StarMakerOverDensityThreshold, &StarMakerMassEfficiency,
-          &StarMakerMinimumMass, &level, &NumberOfNewParticles,
+       &StarMakerMinimumMass, &level, &NumberOfNewParticles, 
        tg->ParticlePosition[0], tg->ParticlePosition[1],
           tg->ParticlePosition[2],
        tg->ParticleVelocity[0], tg->ParticleVelocity[1],
           tg->ParticleVelocity[2],
        tg->ParticleMass, tg->ParticleAttribute[1], tg->ParticleAttribute[0],
           tg->ParticleAttribute[2]);
 
-      for (i = 0; i < NumberOfNewParticles; i++)
-          tg->ParticleType[i] = NormalStarType;
+      for (i = NumberOfNewParticlesSoFar; i < NumberOfNewParticles; i++)
+          tg->ParticleType[i] = PARTICLE_TYPE_STAR;
     } 
 
     if (STARMAKE_METHOD(UNIGRID_STAR)) {
 
       //---- UNIGRID ALGORITHM (NO JEANS MASS)
+
+      NumberOfNewParticlesSoFar = NumberOfNewParticles;
 
       FORTRAN_NAME(star_maker3)(
        GridDimension, GridDimension+1, GridDimension+2,
@@ -503,22 +508,24 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level)
           CellLeftEdge[2], &GhostZones,
        &MetallicityField, &HydroMethod, &StarMakerMinimumDynamicalTime,
        &StarMakerOverDensityThreshold, &StarMakerMassEfficiency,
-          &StarMakerMinimumMass, &level, &NumberOfNewParticles,
+       &StarMakerMinimumMass, &level, &NumberOfNewParticles, 
        tg->ParticlePosition[0], tg->ParticlePosition[1],
           tg->ParticlePosition[2],
        tg->ParticleVelocity[0], tg->ParticleVelocity[1],
           tg->ParticleVelocity[2],
        tg->ParticleMass, tg->ParticleAttribute[1], tg->ParticleAttribute[0],
           tg->ParticleAttribute[2]);
 
-      for (i = 0; i < NumberOfNewParticles; i++)
-          tg->ParticleType[i] = NormalStarType;
+      for (i = NumberOfNewParticlesSoFar; i < NumberOfNewParticles; i++)
+          tg->ParticleType[i] = PARTICLE_TYPE_STAR;
     }
 
     if (STARMAKE_METHOD(KRAVTSOV_STAR)) {
 
       //---- KRAVTSOV SF ALGORITHM
 
+      NumberOfNewParticlesSoFar = NumberOfNewParticles;
+
       FORTRAN_NAME(star_maker4)(
        GridDimension, GridDimension+1, GridDimension+2,
        BaryonField[DensNum], dmfield, temperature, BaryonField[Vel1Num],
@@ -530,22 +537,24 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level)
           CellLeftEdge[2], &GhostZones, 
        &MetallicityField, &HydroMethod, &StarMakerMinimumDynamicalTime, 
        &StarMakerOverDensityThreshold, &StarMakerMassEfficiency,
-          &StarMakerMinimumMass, &level, &NumberOfNewParticles, 
+       &StarMakerMinimumMass, &level, &NumberOfNewParticles,
        tg->ParticlePosition[0], tg->ParticlePosition[1], 
           tg->ParticlePosition[2], 
        tg->ParticleVelocity[0], tg->ParticleVelocity[1], 
           tg->ParticleVelocity[2], 
        tg->ParticleMass, tg->ParticleAttribute[1], tg->ParticleAttribute[0],
           tg->ParticleAttribute[2]);
 
-      for (i = 0; i < NumberOfNewParticles; i++)
-          tg->ParticleType[i] = NormalStarType;
+      for (i = NumberOfNewParticlesSoFar; i < NumberOfNewParticles; i++)
+          tg->ParticleType[i] = PARTICLE_TYPE_STAR;
     }
 
     if (STARMAKE_METHOD(POP3_STAR)) {
 
       //---- POPULATION III (SINGLE STAR)
 
+      NumberOfNewParticlesSoFar = NumberOfNewParticles;
+
       FORTRAN_NAME(pop3_maker)
 	(GridDimension, GridDimension+1, GridDimension+2, BaryonField[DensNum], 
 	 dmfield, h2field, temperature, BaryonField[Vel1Num], 
@@ -573,6 +582,8 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level)
       // Convert into a parameter!
       float StarClusterLifeTime = 10e6;  // yr
 
+      NumberOfNewParticlesSoFar = NumberOfNewParticles;
+
       FORTRAN_NAME(cluster_maker)
 	(GridDimension, GridDimension+1, GridDimension+2, 
 	 StarClusterRegionLeftEdge, StarClusterRegionRightEdge,
@@ -597,7 +608,9 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level)
 
     if (STARMAKE_METHOD(INSTANT_STAR)) {
 
-      //---- MODIFIED SF ALGORITHM (NO-JEANS MASS, NO dt DEPENDENCE)
+      //---- MODIFIED SF ALGORITHM (NO-JEANS MASS, NO dt DEPENDENCE, NO stochastic SF)
+
+      NumberOfNewParticlesSoFar = NumberOfNewParticles;
 
       FORTRAN_NAME(star_maker7)(
        GridDimension, GridDimension+1, GridDimension+2,
@@ -610,13 +623,16 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level)
           CellLeftEdge[2], &GhostZones,
        &MetallicityField, &HydroMethod, &StarMakerMinimumDynamicalTime,
        &StarMakerOverDensityThreshold, &StarMakerMassEfficiency,
-          &StarMakerMinimumMass, &level, &NumberOfNewParticles,
+       &StarMakerMinimumMass, &level, &NumberOfNewParticles, 
        tg->ParticlePosition[0], tg->ParticlePosition[1],
           tg->ParticlePosition[2],
        tg->ParticleVelocity[0], tg->ParticleVelocity[1],
           tg->ParticleVelocity[2],
        tg->ParticleMass, tg->ParticleAttribute[1], tg->ParticleAttribute[0],
           tg->ParticleAttribute[2]);
+
+      for (i = NumberOfNewParticlesSoFar; i < NumberOfNewParticles; i++)
+          tg->ParticleType[i] = PARTICLE_TYPE_STAR;
     } 
 
     /* This creates sink particles which suck up mass off the grid. */

src/enzo/ZeusSource.C

@@ -333,7 +333,7 @@ int ZeusSource(float *d, float *e, float *u, float *v, float *w, float *p,
 	  fprintf(stderr, "grav_y=%"GSYM",%"GSYM",%"GSYM"\n", gr_yacc[IDX(i,j-1,k)],gr_yacc[IDX(i,j,k)], gr_yacc[IDX(i,j+1,k)]);
 	  fprintf(stderr, "grav_z=%"GSYM",%"GSYM",%"GSYM"\n", gr_zacc[IDX(i,j,k-1)],gr_zacc[IDX(i,j,k)], gr_zacc[IDX(i,j,k+1)]);
 	  for (n=is; n <= ie; n++)
-	    fprintf(stderr, "d,e,y,v,w=%"GSYM",%"GSYM",%"GSYM",%"GSYM",%"GSYM"\n", d[IDX(n,j,k)],e[IDX(n,j,k)],u[IDX(n,j,k)],v[IDX(n,j,k)],w[IDX(n,j,k)]);
+	    fprintf(stderr, "d,e,u,v,w=%"GSYM",%"GSYM",%"GSYM",%"GSYM",%"GSYM"\n", d[IDX(n,j,k)],e[IDX(n,j,k)],u[IDX(n,j,k)],v[IDX(n,j,k)],w[IDX(n,j,k)]);
 	  ENZO_FAIL("");
 	}
 

src/enzo/Zeus_yTransport.C

@@ -90,7 +90,7 @@ int Zeus_yTransport(float *d, float *e, float *u, float *v, float *w,
       for (j = js-2; j <= je+3; j++) {
 	f1[IDX(i,j,k)] = dstar[j]*ueff[j];
 	if (fabs(v[IDX(i,j,k)]) > 0.5*dy[j]/dt) {
-	  printf("yt problem 1: u=%"GSYM",%"GSYM",%"GSYM"  dstar=%"GSYM",%"GSYM",%"GSYM"  i,j,k=%"ISYM",%"ISYM",%"ISYM"\n",
+	  printf("yt problem 1: v=%"GSYM",%"GSYM",%"GSYM"  dstar=%"GSYM",%"GSYM",%"GSYM"  i,j,k=%"ISYM",%"ISYM",%"ISYM"\n",
                  v[IDX(i,j-1,k)],v[IDX(i,j,k)],v[IDX(i,j+1,k)], 
 		 dstar[j-1], dstar[j],dstar[j+1],i,j,k);
 	  printf("  d=%"GSYM",%"GSYM",%"GSYM"  dy=%"GSYM",%"GSYM"   dt=%"GSYM"\n",

src/enzo/star_maker2.src

@@ -10,7 +10,7 @@ c
      &                      d1, x1, v1, t1,
      &                      nmax, xstart, ystart, zstart, ibuff, 
      &                      imetal, imethod, mintdyn,
-     &                      odthresh, masseff, smthresh, level, np,
+     &                      odthresh, masseff, smthresh, level, np, 
      &                      xp, yp, zp, up, vp, wp,
      &                      mp, tdp, tcp, metalf)
 c
@@ -101,7 +101,7 @@ c
       parameter (pi=3.14159265, G=6.67e-8, sndspdC=1.3095e8,
      &           msolar=1.989e33)
 c
-      ii = 0
+      ii = np
 c
 c  for each zone, : "star" particle is created if answers to all the
 c  following questions are affirmative:

src/enzo/star_maker3.src

@@ -11,7 +11,7 @@ c
      &                      d1, x1, v1, t1,
      &                      nmax, xstart, ystart, zstart, ibuff, 
      &                      imetal, imethod, mintdyn,
-     &                      odthresh, masseff, smthresh, level, np,
+     &                      odthresh, masseff, smthresh, level, np, 
      &                      xp, yp, zp, up, vp, wp,
      &                      mp, tdp, tcp, metalf)
 c
@@ -102,7 +102,7 @@ c
       parameter (pi=3.14159265, G=6.67e-8, sndspdC=1.3095e8,
      &           msolar=1.989e33)
 c
-      ii = 0
+      ii = np
 
 !     print*,'star_maker3: imetal is:',imetal
 

src/enzo/star_maker4.src

@@ -95,7 +95,7 @@ c
       parameter (pi=3.14159265, G=6.67e-8, sndspdC=1.3095e8,
      &           msolar=1.989e33,mproton=1.67e-24)
 c
-      ii = 0
+      ii = np
 c
 c  calculate density threshold.  odthresh is in proper particles per cc
 c  and (d1/mproton) gives the mean density of the universe in 

src/enzo/star_maker7.src

@@ -107,7 +107,7 @@ c
       parameter (pi=3.14159265, G=6.67e-8, sndspdC=1.3095e8,
      &           msolar=1.989e33)
 c
-      ii = 0
+      ii = np
 c
 c  for each zone, : "star" particle is created if answers to all the
 c  following questions are affirmative: