Britton Smith avatar Britton Smith committed e193993 Merge

Reopening conduction branch. Merged week-of-code changes into conduction branch.

Comments (0)

Files changed (190)

 *~
 
 bin/enzo
+bin/inits
 DEPEND
 DEPEND.bak
 auto_show_*.C

doc/manual/source/EnzoParameters.rst

 
 ``HydroMethod`` (external)
     This integer specifies the hydrodynamics method that will be used.
-    Currently implemented are: 0 - PPM DE (a direct-Eulerian version of
-    PPM), 1 - PPM LR (a Lagrange-Remap version of PPM), 2 - ZEUS (a
-    Cartesian, 3D version of Stone & Norman). The PPM LR version is not
-    recommended. Note that if ZEUS is selected, it automatically turns
-    off ``ConservativeInterpolation`` and the ``DualEnergyFormalism`` flags. 3
-    - Runge Kutta third order based MUSCL solvers. 4 - Same as three
-    but including Dedner MHD (Wang & Abel 2008). For 3 and 4 there are
-    the additional parameters ``RiemannSolver`` and ``ReconstructionMethod``
-    you want to set. Default: 0
+    Currently implemented are
+
+    ============ ===========================
+    Hydro method Description
+    ============ ===========================
+    0            PPM DE (a direct-Eulerian version of PPM)
+    1            PPM LR (a Lagrange-Remap version of PPM). **The PPM LR version is not recommended.**
+    2            ZEUS (a Cartesian, 3D version of Stone & Norman). Note that if ZEUS is selected, it automatically turns off ``ConservativeInterpolation`` and the ``DualEnergyFormalism`` flags. 
+    3            Runge Kutta third order based MUSCL solvers. 
+    4            Same as 3 but including Dedner MHD (Wang & Abel 2008). For 3 and 4 there are the additional parameters ``RiemannSolver`` and ``ReconstructionMethod`` you want to set. 
+    ============ ===========================
+
+    Default: 0
+``RiemannSolver`` (external; only if ``HydroMethod`` is 3 or 4)
+    This integer specifies the Riemann solver used by the MUSCL solver. Choice of
+
+    ============== ===========================
+    Riemann solver Description
+    ============== ===========================
+    0              FluxReconstruction
+    1              HLL (Harten-Lax-van Leer) a two-wave, three-state solver with no resolution of contact waves
+    2              Marquina
+    3              LLF (Local Lax-Friedrichs)
+    4              HLLC (Harten-Lax-van Leer with Contact) a three-wave, four-state solver with better resolution of contacts
+    ============== ===========================
+
+    Default: 1 (HLL)
+``ReconstructionMethod`` (external; only if ``HydroMethod`` is 3 or 4)
+    This integer specifies the reconstruction method for the MUSCL solver. Choice of
+
+    ===================== ====================
+    Reconstruction Method Description
+    ===================== ====================
+    0                     PLM (piecewise linear)
+    1                     PPM (piecwise parabolic)
+    2                     CENO
+    3                     WENO3 (Weighted Essentially Non-Oscillating, 3rd order)
+    4                     WENO5 (Weighted Essentially Non-Oscillating, 5th order)
+    ===================== ====================
+
+    Default: 0 (PLM)
 ``Gamma`` (external)
     The ratio of specific heats for an ideal gas (used by all hydro
     methods). If using multiple species (i.e. ``MultiSpecies`` > 0), then
+#=======================================================================
+#
+# FILE:        Makefile.config
+#
+# SUMMARY:     Configurable Makefile for Enzo
+#
+# AUTHOR:      James Bordner (jobordner@ucsd.edu)
+#
+# DATE:        2007-02-21
+#
+# DESCRIPTION
+#              See 'gmake help' for definitive description of targets
+#
+#              Makefile.config includes the following files:
+# 
+#              Make.config.settings   default configuration settings
+#              Make.config.override   local user settings
+#              Make.config.assemble   maps 'config' settings to 'flag' settings
+#              Make.config.targets    configure targets
+#              Make.mach.*            all machine-dependent settings
+#              Make.config.objects    list of object files
+#              DEPEND                 Make-generated dependencies
+#
+#              Make.mach.* should be the only file that one should
+#              edit when porting Enzo to a new machine.
+#
+#              Make.config.override should be the only file that
+#              one should edit when defining local user settings.
+#              Preferably, this is done implicitly through
+#              the available make targets (e.g. "gmake precision-32").
+#              Use "gmake help-config" for a list of all configuration
+#              settings.  These make targets do error-checking; hand-editing 
+#              Make.config.override is more error-prone.
+#
+#
+# MODIFIED:    Elizabeth Tasker (taskere@mcmaster.ca), 2010-07-06
+#
+#	       This makefile can be run from a directory 'patch' that 
+#	       will compile routines in patch in preference to 
+#              those in src/enzo. 
+#
+#=======================================================================
+
+# Use bash since sh on datastar does not recognize ">&" used in dep: target
+
+SHELL    = /bin/bash
+
+TOP_DIR  = ../
+EXE      = enzo
+OUTPUT   = out.compile
+
+ENZO_DIR = ../src/enzo
+MODULES  = 
+VERBOSE  = 0
+
+SVN      = hg
+
+PATCH = ../patch
+VPATH = $(PATCH):$(ENZO_DIR)
+
+#-----------------------------------------------------------------------
+# Make.config.settings is used for setting default values to all compile-time 
+# configuration settings.
+#-----------------------------------------------------------------------
+
+include $(ENZO_DIR)/Make.config.settings
+
+#-----------------------------------------------------------------------
+# Make.config.machine is used for setting which Make.mach.* file to use
+#-----------------------------------------------------------------------
+
+include $(ENZO_DIR)/Make.config.machine
+
+#-----------------------------------------------------------------------
+# Make.config.override is used for overriding the default settings in
+# Make.config.settings.  This was made separate from the default settings 
+# to enable easily interfacing Enzo with a software testing environment 
+# like lcatest.
+#-----------------------------------------------------------------------
+
+MAKE_CONFIG_OVERRIDE = $(ENZO_DIR)/Make.config.override
+
+include $(MAKE_CONFIG_OVERRIDE)
+
+# THIS WAY OF DOING THE ABOVE DOES NOT WORK:
+#
+# include $(ENZO_DIR)/Make.config.override
+
+#-----------------------------------------------------------------------
+# Make.config.assemble takes the settings in the Make.config.settings
+# and Make.config.override, and generates the appropriate make variables
+# required by this makefile.  E.g. $(CXX), $(CXXFLAGS), etc.
+#-----------------------------------------------------------------------
+
+include $(ENZO_DIR)/Make.config.assemble
+
+#-----------------------------------------------------------------------
+# Make.mach.<machine-name> defines all machine-dependent settings.
+#-----------------------------------------------------------------------
+
+-include $(ENZO_DIR)/Make.mach.$(CONFIG_MACHINE)
+-include $(HOME)/.enzo/Make.mach.$(CONFIG_MACHINE)
+
+#=======================================================================
+# OBJECT FILES
+#=======================================================================
+
+include $(ENZO_DIR)/Make.config.objects
+
+#-----------------------------------------------------------------------
+# MAKE ENZO BY DEFAULT
+#-----------------------------------------------------------------------
+
+all: MACHNOTES $(EXE).exe
+
+#-----------------------------------------------------------------------
+# MAKE AN EXECUTABLE
+#-----------------------------------------------------------------------
+
+MACHNOTES: 
+	@echo -e $(MACHINE_NOTES)
+
+%.exe: 	$(MODULES) autogen dep %.o $(OBJS_LIB) MACHNOTES
+	@rm -f ${ENZO_DIR}/$@
+	@echo "Linking enzo executable. Type  cat $(OUTPUT)  in case it fails."
+	@(if [ $(VERBOSE) -eq 0 ]; then \
+		cd $(ENZO_DIR) ; \
+		$(LD) $(LDFLAGS) -o $*.exe $*.o $(OBJS_LIB) $(LIBS) >& $(OUTPUT) ; \
+		cd $(TOP_DIR)/$(PATCH) ; \
+	else \
+		cd $(ENZO_DIR) \
+		$(LD) $(LDFLAGS) -o $*.exe $*.o $(OBJS_LIB) $(LIBS) >> $(OUTPUT) \
+		2>&1 ; \
+		cd $(TOP_DIR)/$(PATCH) ; \
+	fi)  
+	@(if [ -e $(ENZO_DIR)/$@ ] ; then  \
+		echo "Success!" ; \
+		if [ ! -e $(TOP_DIR)/bin ]; then mkdir $(TOP_DIR)/bin; fi; \
+                cp $(ENZO_DIR)/$(EXE).exe $(TOP_DIR)/bin/$(EXE); \
+	else \
+		echo "$(LD) $(LDFLAGS) -o $< $*.o $(OBJS_LIB) $(LIBS)" > temp1; \
+		cat temp1 $(OUTPUT) > temp2 ;\
+		mv -f temp2 $(OUTPUT); \
+		rm -f temp1 temp2; \
+		echo "Failed! See $(OUTPUT) for error messages"; \
+	fi)
+
+#-----------------------------------------------------------------------
+# WRITE ALL COMPILER OUTPUT TO FILE
+#-----------------------------------------------------------------------
+
+.PHONY: verbose
+verbose: VERBOSE = 1
+#verbose:
+#	@rm -fv $(OUTPUT)
+#	@echo "Writing all compiler output to $(OUTPUT)"
+verbose: $(EXE).exe
+
+#-----------------------------------------------------------------------
+# Implicit rules
+#-----------------------------------------------------------------------
+
+.SUFFIXES: .c .C .F .F90 .o .cu
+
+# Inhibit removing any *.o files after compiling
+
+.PRECIOUS: %.o 
+
+#.src.f:
+#	@$(CPP) $(DEFINES) $(CPPFLAGS) $< > $(ENZO_DIR)/$*.f
+.F.o:
+	@echo "Compiling $<"
+	@rm -f $(ENZO_DIR)/$@
+	@(if [ $(VERBOSE) -eq 0 ]; then \
+		$(FC) -c -o $(ENZO_DIR)/$@ $(FFLAGS) $(DEFINES) $(ENZO_DIR)/$*.F >& $(OUTPUT) ; \
+		if [ ! -e $(ENZO_DIR)/$@ ]; then \
+		echo; \
+		echo "$(FC) -c -o $(ENZO_DIR)/$@ $(FFLAGS) $(DEFINES) $(ENZO_DIR)/$*.F"; \
+		echo; \
+	     	$(FC) -c -o $(ENZO_DIR)/$@ $(FFLAGS) $(DEFINES) $(ENZO_DIR)/$*.F; \
+		echo; \
+		exit 1; \
+	fi ; \
+	else \
+		$(FC) -c -o $(ENZO_DIR)/$@ $(FFLAGS) $(DEFINES) $(ENZO_DIR)/$*.f >> $(OUTPUT) 2>&1 ; \
+		if [ ! -e $(ENZO_DIR)/$@ ]; then \
+	     	echo "See $(OUTPUT) for error messages"; \
+	     	exit 1; \
+	fi ; \
+	fi)
+
+#.src90.f90:#
+#	@echo "Compiling $<"
+#	@$(CPP) $(DEFINES) $(CPPFLAGS) $< > $(ENZO_DIR)/$*.f90
+
+.F90.o:
+	@rm -f $(ENZO_DIR)/$@
+	@(if [ $(VERBOSE) -eq 0 ]; then \
+	  $(F90) -c -o $(ENZO_DIR)/$@ $(F90FLAGS) $(DEFINES) $(ENZO_DIR)/$*.F90 >& $(OUTPUT) ; \
+	  if [ ! -e $(ENZO_DIR)/$@ ]; then \
+             echo; \
+             echo "$(F90) -c -o $(ENZO_DIR)/$@ $(DEFINES) $(F90FLAGS) $(ENZO_DIR)/$*.F90"; \
+             echo; \
+             $(F90) -c -o $(ENZO_DIR)/$@ $(F90FLAGS) $(ENZO_DIR)/$*.F90; \
+             echo; \
+             exit 1; \
+	  fi ; \
+	else \
+	  $(F90) -c -o $(ENZO_DIR)/$@ $(F90FLAGS) $(DEFINES) $(ENZO_DIR)/$*.F90 >> $(OUTPUT) 2>&1 ; \
+	  if [ ! -e $(ENZO_DIR)/$@ ]; then \
+	     echo "See $(OUTPUT) for error messages"; \
+	     exit 1; \
+	  fi ; \
+	fi)
+
+.c.o:
+	@rm -f $(ENZO_DIR)/$@
+	@echo "Compiling $<"
+	@(if [ $(VERBOSE) -eq 0 ]; then \
+	  $(CC) -c -o $(ENZO_DIR)/$@ $(DEFINES) $(CFLAGS) $(INCLUDES) $< \
+	    >& $(OUTPUT) ; \
+	  if [ ! -e $(ENZO_DIR)/$@ ]; then \
+             echo; \
+             echo "$(CC) -c -o $(ENZO_DIR)/$@ $(DEFINES) $(CFLAGS) $(INCLUDES) $<"; \
+             echo; \
+             $(CC) -c -o $(ENZO_DIR)/$@ $(DEFINES) $(CFLAGS) $(INCLUDES) $<;\
+             echo; \
+             exit 1; \
+          fi ; \
+	else \
+	   $(CC) -c -o $(ENZO_DIR)/$@ $(DEFINES) $(CFLAGS) $(INCLUDES) $< \
+	    >> $(OUTPUT) 2>&1 ; \
+	  if [ ! -e $(ENZO_DIR)/$@ ]; then \
+	     echo "See $(OUTPUT) for error messages"; \
+	     exit 1; \
+	  fi ; \
+	fi)
+
+.C.o:
+	@rm -f $(ENZO_DIR)/$@
+	@echo "Compiling $<"
+	@(if [ $(VERBOSE) -eq 0 ]; then \
+	  $(CXX) -c -o $(ENZO_DIR)/$@ $(DEFINES) $(CXXFLAGS) $(INCLUDES) $<  \
+	    >& $(OUTPUT) ; \
+	  if [ ! -e $(ENZO_DIR)/$@ ]; then \
+             echo; \
+             echo "$(CXX) -c -o $(ENZO_DIR)/$@ $(DEFINES) $(CXXFLAGS) $(INCLUDES) $<"; \
+             echo; \
+             $(CXX) -c -o $(ENZO_DIR)/$@ $(DEFINES) $(CXXFLAGS) $(INCLUDES) $<;\
+             echo; \
+             exit 1; \
+          fi ; \
+	else \
+	  $(CXX) -c -o $(ENZO_DIR)/$@ $(DEFINES) $(CXXFLAGS) $(INCLUDES) $< \
+	    >> $(OUTPUT) 2>&1 ; \
+	  if [ ! -e $(ENZO_DIR)/$@ ]; then \
+	     echo "See $(OUTPUT) for error messages"; \
+	     exit 1; \
+	  fi ; \
+	fi)
+
+.cu.o: 
+	@rm -f $@
+	@echo "Compiling $<"
+	@(if [ $(VERBOSE) -eq 0 ]; then \
+	  $(CUDACOMPILER) -c -o $@ $(DEFINES) $(CUDACOMPFLAGS) \
+	    $(INCLUDES) $*.cu >& $(OUTPUT) ; \
+	  if [ ! -e $@ ]; then \
+	     echo; \
+             echo "$(CUDACOMPILER) -c -o $@ $(DEFINES) $(CUDACOMPFLAGS) $(INCLUDES) $*.cu"; \
+             echo; \
+             $(CUDACOMPILER)  -c -o $@ $(DEFINES) $(CUDACOMPFLAGS) $(INCLUDES) $*.cu;\
+             echo; \
+             exit 1; \
+          fi ; \
+	else \
+	  $(CUDACOMPILER) -c -o $@ $(DEFINES) $(CUDACOMPFLAGS) \
+	    $(INCLUDES) $*.cu >> $(OUTPUT) 2>&1 ; \
+	  if [ ! -e $@ ]; then \
+	     echo "See $(OUTPUT) for error messages"; \
+	     exit 1; \
+	  fi ; \
+	fi)
+
+#-----------------------------------------------------------------------
+# Generate all make-generated source files
+#-----------------------------------------------------------------------
+
+.PHONY: autogen
+autogen: svn_version.def auto_show_config.C auto_show_flags.C auto_show_version.C
+
+# Force update of auto_show_config.C
+
+.PHONY: auto_show_config.C
+auto_show_config.C:
+	-@$(MAKE) -s show-config  >& temp.show-config
+	-@awk 'BEGIN {print "#include <stdio.h>\nvoid auto_show_config(FILE *fp) {"}; {print "   fprintf (fp,\""$$0"\\n\");"}; END {print "}"}' < temp.show-config > auto_show_config.C
+
+# Force update of auto_show_flags.C
+
+.PHONY: auto_show_flags.C
+auto_show_flags.C:
+	-@$(MAKE) -s show-flags  >& temp.show-flags
+	-@awk 'BEGIN {print "#include <stdio.h>\nvoid auto_show_flags(FILE *fp) {"}; {print "   fprintf (fp,\""$$0"\\n\");"}; END {print "}"}' < temp.show-flags > auto_show_flags.C
+
+# Force update of auto_show_version.C
+
+.PHONY: auto_show_version.C
+auto_show_version.C:
+	-@$(MAKE) -s show-version  >& temp.show-version
+	-@awk 'BEGIN {print "#include <stdio.h>\nvoid auto_show_version(FILE *fp) {"}; {print "   fprintf (fp,\""$$0"\\n\");"}; END {print "}"}' < temp.show-version > auto_show_version.C
+
+#-----------------------------------------------------------------------
+# Generate SVN version source file 'svn_version.def'
+#-----------------------------------------------------------------------
+
+# Force update of svn_version.def
+.PHONY: svn_version.def
+svn_version.def:
+	-@$(SVN) identify -i | \
+           awk '{print "#define ENZO_SVN_REVISION","\""$$1"\""}' > svn_version.def
+	-@$(SVN) identify -b | awk '{print "#define ENZO_SVN_BRANCH","\""$$1"\""};' >> svn_version.def 
+	-@if [ ! -s svn_version.def ]; then \
+	printf "#define ENZO_SVN_REVISION 0\n" > svn_version.def; \
+	printf "#define ENZO_SVN_BRANCH \"\"\n" >> svn_version.def; \
+	fi
+
+#-----------------------------------------------------------------------
+# Generate dependency file
+#-----------------------------------------------------------------------
+
+.PHONY: dep
+dep:
+	@echo "Updating DEPEND"
+	-@(makedepend $(DEFINES) $(INCLUDES) -fDEPEND -o.o -m -- -- *.C *.c *.src *src90 *.h) >& out.make.DEPEND
+
+include $(ENZO_DIR)/DEPEND
+
+#-----------------------------------------------------------------------
+# Radiative transfer module
+#-----------------------------------------------------------------------
+
+#include $(ENZO_DIR)/photons/Make.config.objects
+#
+#.PHONY: photon
+#photon: OBJS_LIB += photons/*.o
+#photon:
+#	@echo "Making radiative transfer module"
+#	+(cd photons/ ; make photon)
+
+#-----------------------------------------------------------------------
+# HELP TARGET
+#-----------------------------------------------------------------------
+
+help:
+	@echo
+	@echo "========================================================================"
+	@echo "   Enzo Makefile Help"
+	@echo "========================================================================"
+	@echo
+	@echo "   gmake                Compile and generate the executable 'enzo.exe'"
+	@echo "   gmake install        Copy the executable to bin/enzo"
+	@echo "   gmake help           Display this help information"
+	@echo "   gmake clean          Remove object files, executable, etc."
+	@echo "   gmake dep            Create make dependencies in DEPEND file"
+	@echo
+	@echo "   gmake get-version    Generate version-related files"
+	@echo "   gmake show-version   Display Subversion branch and revision"
+	@echo
+	@echo "   gmake help-config    Display detailed help on configuration make targets"
+	@echo "   gmake show-config    Display the configuration settings"
+	@echo "   gmake show-flags     Display specific compilation flags"
+	@echo "   gmake default        Reset the configuration to the default values"
+	@echo
+
+#-----------------------------------------------------------------------
+
+clean:
+	-@rm -f $(ENZO_DIR)/*.o $(ENZO_DIR)/uuid/*.o *.mod $(ENZO_DIR)/*.f $(ENZO_DIR)/*.f90 $(ENZO_DIR)/DEPEND.bak $(ENZO_DIR)/*~ $(OUTPUT) $(ENZO_DIR)/*.exe \
+          auto_show*.C svn_version.def $(ENZO_DIR)/hydro_rk/*.o $(ENZO_DIR)/*.oo $(ENZO_DIR)/hydro_rk/*.oo temp.show* $(ENZO_DIR)/uuid/*.oo DEPEND
+	-@touch DEPEND
+
+#-----------------------------------------------------------------------
+# Include configuration targets
+#-----------------------------------------------------------------------
+
+include $(ENZO_DIR)/Make.config.targets
+# DO NOT DELETE

run/Cooling/CoolingTest/CoolingTest_PrimordialCVODE.enzo

+#
+# PROBLEM DEFINITION FILE: Cooling Test
+#
+#  Iterate rate solver without hydro and output cooling time at the end.
+#
+#
+#  problem setup
+#
+# Cooling Test
+ProblemType =  62
+
+#
+#  grid setup
+#
+TopGridRank = 3
+
+# Dim 0 - H number density
+# Dim 1 - metallicity
+# Dim 2 - temperature
+TopGridDimensions = 16 1 32
+
+CoolingTestMinimumHNumberDensity = 1
+CoolingTestMaximumHNumberDensity = 16
+CoolingTestMinimumMetallicity    = -16
+CoolingTestMaximumMetallicity    = -16
+CoolingTestMinimumTemperature    = 10
+CoolingTestMaximumTemperature    = 1e4
+
+CoolingTestResetEnergies         = 0
+
+#
+#  set I/O and stop/start parameters
+#
+StopTime                  = 0.05
+StopCycle                 = 100000
+dtDataDump                = 0.05
+DataDumpDir               = ct_cv_prim_
+DataDumpName              = ct_cv_prim_
+
+#
+#  set hydro parameters
+#
+HydroMethod               = -1
+UseHydro                  = 0    // no hydro
+DualEnergyFormalism       = 1
+SelfGravity = 0
+FluxCorrection = 0
+
+#
+#  set grid refinement parameters
+#
+StaticHierarchy           = 1   // no AMR
+
+#
+#  set some global parameters
+#
+OutputTemperature         = 1
+
+#
+# Units
+#
+DensityUnits              = 1.67e-24   // 1 g cm^-3
+LengthUnits               = 3.0857e+18 // 1 pc in cm
+TimeUnits                 = 3.1557e+13 // 1 Myr in s
+
+#
+# chemistry/cooling
+#
+RadiativeCooling          = 1
+MultiSpecies              = 3
+MetalCooling              = 0          // no metal cooling
+PrimordialChemistrySolver = 2
+
+# Initial species fractions, fiddle at own risk,
+#TestProblemInitialHIFraction      = 0.998
+#TestProblemInitialHIIFraction     = 1e-10
+#TestProblemInitialHeIFraction     = 1.0
+#TestProblemInitialHeIIFraction    = 1.0e-20
+#TestProblemInitialHeIIIIFraction  = 1.0e-20
+#TestProblemInitialHMFraction      = 1.e-20
+#TestProblemInitialH2IFraction     = 1.e-3
+#TestProblemInitialH2IIFraction    = 1.e-20

src/P-GroupFinder/Make.config.objects

 		io_enzo.o \
 		io_hdf.o \
 		io_hdf5.o \
+		h5utilities.o \
 		cmpfunc.o \
 		allvars.o \
 		allocate.o \

src/P-GroupFinder/Makefile

 #=======================================================================
 
 # Override previous defines with the exception of particle id size
-DEFINES = -DUSE_HDF5 $(ASSEMBLE_IDS_DEFINES)
+DEFINES = -DUSE_HDF5 -DH5_USE_16_API $(ASSEMBLE_IDS_DEFINES)
 
 #-----------------------------------------------------------------------
 # Implicit rules
 	-@($(CC) -c -o $@ $(DEFINES) $(CFLAGS) $(INCLUDES) $*.c) >& $(OUTPUT)
 	@(if [ ! -e $@ ]; then \
              echo; \
-             echo "$(CXX) -c -o $@ $(DEFINES) $(CXXFLAGS) $(INCLUDES) $*.c"; \
+             echo "$(CC) -c -o $@ $(DEFINES) $(CXXFLAGS) $(INCLUDES) $*.c"; \
              echo; \
-             $(CXX) -c -o $@ $(DEFINES) $(CXXFLAGS) $(INCLUDES) $*.c;\
+             $(CC) -c -o $@ $(DEFINES) $(CXXFLAGS) $(INCLUDES) $*.c;\
              echo; \
              exit 1; \
           fi)

src/P-GroupFinder/allvars.c

 
  double  Time;
  double  BoxSize;
+ double  RhoCritical0;
  double  leftEdge[3], rightEdge[3];
 
  double  SearchRadius;

src/P-GroupFinder/allvars.h

 #include "nrsrc/nrutil.h"
 #include "ngbtree.h"
 
+/* HDF5 definitions */
+
+#define HDF5_FILE_I4 H5T_STD_I32BE
+#define HDF5_FILE_I8 H5T_STD_I64BE
+#define HDF5_FILE_R4 H5T_IEEE_F32BE
+#define HDF5_FILE_R8 H5T_IEEE_F64BE
+#define HDF5_FILE_B8 H5T_STD_B8BE
+
+#define HDF5_I4  H5T_NATIVE_INT
+#define HDF5_I8  H5T_NATIVE_LLONG
+#define HDF5_R4  H5T_NATIVE_FLOAT
+#define HDF5_R8  H5T_NATIVE_DOUBLE
+#define HDF5_R16 H5T_NATIVE_LDOUBLE
+
 /* Definitions for controlling the integer type for particle IDs
    (8-byte needed for >2 billion particle simulations) */
 
 #define PISYM "lld"
 #endif
 
+#ifdef CONFIG_PFLOAT_16
+#define PSYM "Lf"
+#define GSYM "g"
+#define GOUTSYM ".21Lg"
+#define HDF5_PREC HDF5_R16
+#define HDF5_FILE_PREC HDF5_R16
+#endif
+
 #define  KERNEL_TABLE 10000
 #define  PI               3.1415927
 #define  GRAVITY     6.672e-8
 extern int     ThisTask, NTask;
 
 extern double  Time;
+extern double  RhoCritical0;
 extern double  BoxSize;
 extern double  leftEdge[3], rightEdge[3];
 

src/P-GroupFinder/h5utilities.c

+/*********************************************************************
+ *
+ *  file    : h5utilities.cpp
+ *
+ *  Project : Visualization of Adaptive Mesh Refinement Data
+ *
+ *  Company : Zuse Institute Berlin
+ *            All rights reserved.
+ *
+ *  Author  : Ralf Kaehler                           
+ *
+ *  Date    : 26.01.2006
+ *
+ *********************************************************************/
+
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include <string.h>
+
+#include <hdf5.h>
+
+
+// need the following ifdef to make standalone version of carpet2amira 
+// independent of hxhdf5 and hence of amira kernel libs ...
+//#define CARPET2AMIRA_MAKESTANDALONE
+//#ifndef CARPET2AMIRA_MAKESTANDALONE
+//#include <hxhdf5/hxhdf5.h>
+//#endif
+
+#include "h5utilities.h"
+
+
+int checkErr(int arg, const char* name) {
+  int callAssert = 0;
+  if (arg<0) {
+    fprintf(stderr,"Error accessing %s. \n",name); 
+    assert(!callAssert);
+    return 1;
+  }
+  return 0;
+}
+
+
+
+hid_t openH5File(const char* filename)
+{
+
+  hid_t fileId = -1;
+
+  int err = 0;
+
+  // need the following ifdef to make standalone version of carpet2amira 
+  // independent of hxhdf5 and hence of amira kernel libs ...
+
+  //#ifdef CARPET2AMIRA_MAKESTANDALONE
+#define DOIT
+#ifdef DOIT
+      
+  // try to open HDF5 files
+  fileId = H5Fopen (filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+
+#else
+
+  hid_t prop = H5Pcreate (H5P_FILE_ACCESS);
+
+  if (prop < 0) { err = 1; }
+
+  if (!err) {
+    if (hx5set_driver (filename, prop) < 0) { err = 2; }
+  }
+
+  if (!err) {
+    fileId = H5Fopen (filename, H5F_ACC_RDONLY, prop);
+    if (fileId<0) {
+      err = 3;
+    }
+  }
+
+  if (prop >= 0) {
+    H5Pclose (prop);
+  }
+
+#endif
+#undef DOIT
+  
+  return fileId;
+
+}
+
+
+
+int readAttribute( const hid_t loc_id, 
+		   const hid_t type_id, 
+		   const char* attrName, 
+		   void* result)
+{
+
+  hid_t attr;
+  int err = 0;
+
+  err |= checkErr ( attr = H5Aopen_name(loc_id,attrName),attrName );
+  err |= checkErr ( H5Aread(attr,type_id,result),attrName );
+  err |= checkErr ( H5Aclose(attr),attrName );
+  
+  return err;
+
+}
+
+
+int writeScalarAttribute( hid_t       loc_id, 
+			  hid_t       type_id, 
+			  const char* attrName,
+			  const void* buffer)
+{
+
+  hid_t attrs, attr;
+  
+  int err = 0;
+
+  err |= checkErr (attrs = H5Screate (H5S_SCALAR), attrName );                                    
+  err |= checkErr (attr  = H5Acreate (loc_id, attrName, type_id, attrs, H5P_DEFAULT), attrName ); 
+  err |= checkErr (H5Awrite (attr,  type_id, buffer), attrName  );                                
+  err |= checkErr (H5Sclose (attrs), attrName );                                                  
+  err |= checkErr (H5Aclose (attr), attrName );                                                   
+
+  return (err);
+
+}
+
+
+int writeArrayAttribute( hid_t         loc_id, 
+			 hid_t         type_id, 
+			 const hsize_t dims, 
+			 const char*   attrName,
+			 const void*   buffer  )
+{
+
+  hid_t dataspace_id, attribute_id;
+  int   err = 0;
+
+ 
+
+  err |= checkErr ( dataspace_id = H5Screate_simple(1, &dims, NULL), attrName );
+  err |= checkErr ( attribute_id = H5Acreate(loc_id, attrName, type_id, dataspace_id, H5P_DEFAULT), attrName );
+  err |= checkErr ( H5Awrite(attribute_id, type_id, buffer), attrName );
+  err |= checkErr ( H5Aclose(attribute_id), attrName );
+  err |= checkErr ( H5Sclose(dataspace_id), attrName );
+
+  return err;
+
+}
+
+int writeArrayDataset( hid_t         loc_id, 
+			 hid_t         type_id, 
+             int ndims,
+			 hsize_t *dims, 
+			 const char*   attrName,
+			 const void*   buffer  )
+{
+
+  hid_t dataspace_id, dataset_id;
+  int   err = 0;
+
+  err |= checkErr ( dataspace_id = H5Screate_simple(ndims, dims, NULL), attrName );
+  err |= checkErr ( dataset_id = H5Dcreate(loc_id, attrName, type_id, dataspace_id, H5P_DEFAULT), attrName );
+  err |= checkErr ( H5Dwrite(dataset_id, type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer), attrName );
+  err |= checkErr ( H5Dclose(dataset_id), attrName );
+  err |= checkErr ( H5Sclose(dataspace_id), attrName );
+
+  return err;
+
+}
+
+
+int writeStringAttribute(hid_t loc_id, const char* name, const char* data)
+{
+
+  int err = 0;
+
+  hid_t	  dataspace, attr, datatype;
+  int str_length = strlen(data);
+  if (str_length == 0) str_length = 1;
+
+  err |= checkErr( dataspace = H5Screate(H5S_SCALAR), data) ;
+  err |= checkErr( datatype  = H5Tcopy(H5T_C_S1), data);
+  err |= checkErr( H5Tset_size(datatype,str_length), data);
+  err |= checkErr( attr = H5Acreate(loc_id, name, datatype, dataspace, H5P_DEFAULT), data);
+
+  err |= checkErr( H5Awrite(attr, datatype, data), data); 
+  err |= checkErr( H5Sclose(dataspace), data); 
+  err |= checkErr( H5Aclose(attr), data); 
+
+  return err;
+
+}
+
+
+
+#ifdef UNUSED
+nativeTypeId h5Type2NativeType (const hid_t h5type) 
+{
+
+  assert(h5type>0);
+
+  const H5T_class_t h5class = H5Tget_class(h5type);
+  const size_t      h5prec  = H5Tget_precision (h5type);
+  const H5T_sign_t  h5sign  = H5Tget_sign(h5type);
+
+  if (h5class==H5T_FLOAT) {
+
+    if (h5prec == 32) { return NTID_float; }
+    else if (h5prec == 64) { return NTID_double; }
+    else { assert(0); }
+
+  }
+  else if (h5class == H5T_INTEGER) {
+
+    if (h5sign == H5T_SGN_2) {
+      if (h5prec == 16) { return NTID_int16; }
+      else if (h5prec == 32) { return NTID_int32; }
+      else { assert(0); }
+    }
+    else if (h5sign == H5T_SGN_NONE) {
+
+      if (h5prec == 8) { return NDTI_uint8; }
+      else if (h5prec == 16) { return NTID_uint16; }
+      else { assert(0); }
+
+    }
+    else {
+      assert(0);
+    }
+
+  }
+  else {
+    assert(0);
+  }
+
+  assert(0);
+
+  return NTID_unknown;
+
+};
+#endif /* UNUSED */
+

src/P-GroupFinder/h5utilities.h

+/*********************************************************************
+ *
+ *  file    : h5utilities.h
+ *
+ *  Project : Visualization of Adaptive Mesh Refinement Data
+ *
+ *  Company : Zuse Institute Berlin
+ *            All rights reserved.
+ *
+ *  Author  : Ralf Kaehler                           
+ *
+ *  Date    : 26.01.2006
+ *
+ *********************************************************************/
+
+#ifndef  _HDF5_UTILITIES_
+#define  _HDF5_UTILITIES_
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <hdf5.h>
+
+//typedef int nativeTypeId;
+//const nativeTypeId NTID_unknown = -1, NTID_float = 0, NTID_double = 1, 
+//  NTID_int16 = 2, NTID_int32 = 3, NDTI_uint8 = 4, NTID_uint16 = 5;
+
+////////////////////////////////////////////////////
+///// HDF5 wrapper functionality
+////////////////////////////////////////////////////
+
+int checkErr( int err, const char* name);
+
+   
+int readAttribute( const hid_t loc_id, 
+		   const hid_t type_id, 
+		   const char* attrName, 
+		   void*       result);
+   
+
+int writeScalarAttribute( hid_t       loc_id, 
+			  hid_t       type_id, 
+			  const char* attrName,
+			  const void* buffer  );
+
+
+int writeArrayAttribute( hid_t         loc_id, 
+			 hid_t         type_id, 
+			 const hsize_t dims, 
+			 const char*   attrName,
+			 const void*   buffer  );
+
+int writeArrayDataset( hid_t         loc_id, 
+			 hid_t         type_id, 
+             int ndims,
+			 hsize_t *dims, 
+			 const char*   attrName,
+			 const void*   buffer  );
+
+
+int writeStringAttribute( hid_t loc_id, 
+                          const char* datasetname, 
+                          const char* data);
+
+hid_t openH5File(const char* filename);
+
+//nativeTypeId h5Type2NativeType (const hid_t h5type); 
+
+#endif

src/P-GroupFinder/io_enzo.c

     Omega * (3.0 * pow(HUBBLE*HubbleConstantNow, 2) / 8.0 / PI / GRAVITY) *
     pow(BoxSize/HubbleConstantNow*CM_PER_MPC,3) /
     (TopGrid[0] * TopGrid[1] * TopGrid[2]);
+  // Critical density in units of Msun / kpc^3
+  RhoCritical0 = 1.4775867e31 * 
+    ((3 * pow(100 * HubbleConstantNow / 3.086e19, 2)) / (8 * M_PI * GRAVITY));
+
 
   EnzoTimeUnit = CM_PER_MPC * (BoxSize/HubbleConstantNow/(1+initialRedshift)) / 
     EnzoVelocityUnit;
   for(i=0; i<files; i++) {
 
     if (i % (files/20) == 0 && ThisTask == 0) {
-      fprintf(stdout, "Read %d out of %d files.\n", i, files);
+      fprintf(stdout, "Read %d out of %d grids.\n", i, files);
       fflush(stdout);
     }
 
       MPI_Finalize();
       exit(1);
     }
-    if ((group_id = H5Gopen(h5_file, buf, H5P_DEFAULT)) < 0) {
+    if ((group_id = H5Gopen(h5_file, buf)) < 0) {
       fprintf(stderr, "enzoCountLoadParticles: error opening %s in %s\n", 
 	      buf, filename[i]);
       MPI_Finalize();

src/P-GroupFinder/io_hdf5.c

   hsize_t freal_size;
   float 	*temp  = NULL;
 
-  if ((data_id = H5Dopen(group_id, label, H5P_DEFAULT)) < 0) {
+  if ((data_id = H5Dopen(group_id, label)) < 0) {
     fprintf(stderr, "ReadParticleField: cannot read %s\n", label);
     fprintf(stderr, "GROUP_ID = %d, DATA_ID = %d\n", group_id, data_id);
     MPI_Finalize();
   hsize_t freal_size;
   float 	*temp  = NULL;
 
-  if ((data_id = H5Dopen(group_id, label, H5P_DEFAULT)) < 0) {
+  if ((data_id = H5Dopen(group_id, label)) < 0) {
     fprintf(stderr, "ReadParticleField: cannot read %s\n", label);
     fprintf(stderr, "GROUP_ID = %d, DATA_ID = %d\n", group_id, data_id);
     MPI_Finalize();
   hid_t data_id, data_type_id;
   PINT 	*temp  = NULL;
 
-  if ((data_id = H5Dopen(group_id, label, H5P_DEFAULT)) < 0) {
+  if ((data_id = H5Dopen(group_id, label)) < 0) {
     fprintf(stderr, "ReadParticleField: cannot read %s\n", label);
     fprintf(stderr, "GROUP_ID = %d, DATA_ID = %d\n", group_id, data_id);
     MPI_Finalize();

src/P-GroupFinder/main.c

 #include <assert.h>
 #include <math.h>
 #include <mpi.h>
+#include <hdf5.h>
 
 #include <sys/stat.h>
 #include <sys/types.h>
 #include "allvars.h"
 #include "nrsrc/nrutil.h"
 #include "proto.h"
+#include "h5utilities.h"
 
 double LinkLength = 0.1;   /* in terms of mean interparticle seperation */
 int    GroupMinLen= 50;    /*  store only groups in the catalogue 
 {
   char input_fname[200];
   char cataloguetxt[200];
+  char scataloguetxt[200];
   char catalogue_fname[200];
   char subhalos_fname[200];
   char particles_fname[200];
+  char particles_fname5[200];
+  char sparticles_fname5[200];
   char parttypes_fname[200];
   char partids_fname[200];
   char subprop_fname[200];
   MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
   MPI_Comm_size(MPI_COMM_WORLD, &NTask);
 
-  if(NTask<=1)
-    {
-      if(ThisTask==0)
-	fprintf(stdout, "Number of processors MUST be a larger than 1.\n");
-      MPI_Finalize(); 
-      exit(1);
-    }
+//  if(NTask<=1)
+//    {
+//      if(ThisTask==0)
+//	fprintf(stdout, "Number of processors MUST be a larger than 1.\n");
+//      MPI_Finalize(); 
+//      exit(1);
+//    }
 
-  if(NTask&1)
+  if(NTask&1 && NTask>1)
     {
       if(ThisTask==0)
 	fprintf(stdout, "Number of processors MUST be a multiple of 2.\n");
 	  fprintf(stderr,"<path>      (path)\n");
 	  fprintf(stderr,"<basename>  (basename of snapshot files)\n");
 	  fprintf(stderr,"<num>       (number of snapshot)\n");
-	  fprintf(stderr,"<format>    (0 = GADGET, 1 = Enzo/HDF4)\n");
+	  fprintf(stderr,"<format>    (0 = GADGET, 1 = Enzo/HDF)\n");
 	  fprintf(stderr,"\n\n");
 	}
       MPI_Finalize(); 
 
   if (formatType != GADGET && formatType != ENZO) {
     fprintf(stderr, "error: please specify a correct format\n"
-	    "0 = GADGET, 1 = Enzo/HDF4\n");
+	    "0 = GADGET, 1 = Enzo/HDF\n");
     MPI_Finalize();
     exit(1);
   }
     sprintf(partids_fname, "%s/groups/groups_%4.4d.ids", path, Snapshot);
     sprintf(catalogue_fname, "%s/groups/groups_%4.4d.fofcat", path, Snapshot);
     sprintf(cataloguetxt, "%s/groups/groups_%4.4d.dat", path, Snapshot);
+    sprintf(scataloguetxt, "%s/groups/subgroups_%4.4d.dat", path, Snapshot);
+    sprintf(particles_fname5, "%s/groups/particles_%4.4d.h5", path, Snapshot);
+    sprintf(sparticles_fname5, "%s/groups/subparticles_%4.4d.h5", path, Snapshot);
     sprintf(subhalos_fname,  "%s/groups/groups_%4.4d.subcat", path, Snapshot);
     sprintf(subprop_fname,   "%s/groups/groups_%4.4d.subprop", path, Snapshot);
     sprintf(fofprop_fname,   "%s/groups/groups_%4.4d.fofprop", path, Snapshot);
   /*  adjust_sfr();*/
 
   marking();
- 
-  exchange_shadow();
+
+  if (NTask>1)
+    exchange_shadow();
 
   init_coarse_grid();
   
   link_local_slab(); 
     
-  do
-    {
-      find_minids();
-    }
-  while(link_accross()>0);
+  if (NTask>1)
+    do
+      {
+	find_minids();
+      }
+    while(link_accross()>0);
   
   find_minids();
   
-  stitch_together(); 
+  if (NTask>1)
+    stitch_together(); 
 
   compile_group_catalogue();
 
-#ifdef FOF_ONLY
-  save_groups(particles_fname, catalogue_fname, parttypes_fname, 
-	      partids_fname);
-#else
+  save_groups(particles_fname, particles_fname5, catalogue_fname, 
+	      parttypes_fname, partids_fname, cataloguetxt);
+#ifndef FOF_ONLY
   subfind(particles_fname, catalogue_fname, subhalos_fname, 
-	  parttypes_fname, partids_fname, subprop_fname, fofprop_fname);
-  write_ascii_catalog(catalogue_fname, fofprop_fname, cataloguetxt);
+	  parttypes_fname, partids_fname, subprop_fname, fofprop_fname,
+	  sparticles_fname5, scataloguetxt);
+  //write_ascii_catalog(catalogue_fname, fofprop_fname, cataloguetxt);
 #endif
 
   MPI_Barrier(MPI_COMM_WORLD);
 
 
 
-void save_groups(char *particles_fname, char *catalogue_fname,
-		 char *parttypes_fname, char *partids_fname)
+void save_groups(char *particles_fname, char *particles_fname5, 
+		 char *catalogue_fname, char *parttypes_fname, 
+		 char *partids_fname, char *cataloguetxt)
 {
-  FILE   *fd, *fdtypes, *fdids;
-  int    i, gr, offset;
+  FILE   *fd, *fdtypes, *fdids, *fdtxt;
+  int    i, gr, offset, index, dim;
   int    ntot;
   int    head, len;
   char   ctype;
-  float cm[3], mtot, mgas, mstars, sfr, mcloud;	  
+  char   halo_name[200];
+  double *temp;
+  PINT *TempPINT;
+  float cm[3], cmv[3], AM[3], vrms, spin, mtot, mgas, mstars, 
+    sfr, mcloud, mvir, rvir;
   struct particle_data *Pbuf;
 
+  hid_t  file_id, dset_id, dspace_id, group_id;
+  hsize_t hdims[2];
+  herr_t status;
+
   if(ThisTask==0)
     {
       printf("writing group catalogue...\n");
       fwrite(&ntot, sizeof(int), 1, fdids);
     }
 
+  if (ThisTask == 0) {
+
+    fprintf(stdout, "Saving halo ASCII list to %s\n", cataloguetxt);
+    fprintf(stdout, "Saving halo HDF5 particle list to %s\n", particles_fname5);
+    if ((fdtxt = fopen(cataloguetxt, "w")) == NULL) {
+      printf("can't open file `%s`\n", cataloguetxt );
+      MPI_Abort(MPI_COMM_WORLD, 1); exit(1);
+    }
+
+    // Write header
+
+    float redshift = 1.0 / Time - 1.0;
+    fprintf(fdtxt, "# Scale Factor = %f\n", Time);
+    fprintf(fdtxt, "# Redshift = %f\n", redshift);
+    fprintf(fdtxt, "# Number of halos = %d\n", NgroupsAll);
+    fprintf(fdtxt, "#\n");
+    fprintf(fdtxt, "# Column 1.  Center of mass (x)\n");
+    fprintf(fdtxt, "# Column 2.  Center of mass (y)\n");
+    fprintf(fdtxt, "# Column 3.  Center of mass (z)\n");
+    fprintf(fdtxt, "# Column 4.  Halo number\n");
+    fprintf(fdtxt, "# Column 5.  Number of particles\n");
+    fprintf(fdtxt, "# Column 6.  Halo mass [solar masses]\n");
+    fprintf(fdtxt, "# Column 7.  Virial mass [solar masses]\n");
+    fprintf(fdtxt, "# Column 8.  Stellar mass [solar masses]\n");
+    fprintf(fdtxt, "# Column 9.  Virial radius (r200) [kpc]\n");
+    fprintf(fdtxt, "# Column 10. Mean x-velocity [km/s]\n");
+    fprintf(fdtxt, "# Column 11. Mean y-velocity [km/s]\n");
+    fprintf(fdtxt, "# Column 12. Mean z-velocity [km/s]\n");
+    fprintf(fdtxt, "# Column 13. Velocity dispersion [km/s]\n");
+    fprintf(fdtxt, "# Column 14. Mean x-angular momentum [Mpc * km/s]\n");
+    fprintf(fdtxt, "# Column 15. Mean y-angular momentum [Mpc * km/s]\n");
+    fprintf(fdtxt, "# Column 16. Mean z-angular momentum [Mpc * km/s]\n");
+    fprintf(fdtxt, "# Column 17. Spin parameter\n");
+    fprintf(fdtxt, "#\n");
+    fprintf(fdtxt, "# datavar lines are for partiview.  Ignore them if you're not partiview.\n");
+    fprintf(fdtxt, "#\n");
+    fprintf(fdtxt, "datavar 0 halo_number\n");
+    fprintf(fdtxt, "datavar 1 number_of_particles\n");
+    fprintf(fdtxt, "datavar 2 halo_mass\n");
+    fprintf(fdtxt, "datavar 3 virial_mass\n");
+    fprintf(fdtxt, "datavar 4 stellar_mass\n");
+    fprintf(fdtxt, "datavar 5 virial_radius\n");
+    fprintf(fdtxt, "datavar 6 x_velocity\n");
+    fprintf(fdtxt, "datavar 7 y_velocity\n");
+    fprintf(fdtxt, "datavar 8 z_velocity\n");
+    fprintf(fdtxt, "datavar 9 velocity_dispersion\n");
+    fprintf(fdtxt, "datavar 10 x_angular_momentum\n");
+    fprintf(fdtxt, "datavar 11 y_angular_momentum\n");
+    fprintf(fdtxt, "datavar 12 z_angular_momentum\n");
+    fprintf(fdtxt, "datavar 13 spin\n");
+    fprintf(fdtxt, "\n");
+
+    file_id = H5Fcreate(particles_fname5, H5F_ACC_TRUNC, H5P_DEFAULT, 
+			H5P_DEFAULT);
+    group_id = H5Gcreate(file_id, "/Parameters", 0);
+    writeScalarAttribute(group_id, H5T_NATIVE_FLOAT, "Redshift", &redshift);
+    writeScalarAttribute(group_id, H5T_NATIVE_DOUBLE, "Scale Factor", &Time);
+    writeScalarAttribute(group_id, H5T_NATIVE_INT, "Number of groups", &NgroupsAll);
+    H5Gclose(group_id);
+
+  } // ENDIF Task 0
+  
 
   for(gr=NgroupsAll-1; gr>=0; gr--)
     {
 	      fwrite(&ctype, sizeof(char), 1, fdtypes);
 	    }
 
-	  get_properties(Pbuf, len, &cm[0], &mtot, &mgas, &mstars, &sfr, &mcloud);
+	  get_properties(Pbuf, len, &cm[0], &mtot, &mgas, &mstars, &sfr, &mcloud,
+			 0, &cmv[0], &mvir, &rvir, AM, &vrms, &spin);
+
+	  /* Write to ASCII catalog */
+
+	  fprintf(fdtxt, "%12.6f %12.6f %12.6f %12d %12d %12.6g %12.6g %12.6g "
+		  "%12.6g %12.6g %12.6g %12.6g %12.6g %12.6g %12.6g %12.6g %12.6g\n",
+		  cm[0], cm[1], cm[2], NgroupsAll-1-gr, len, 
+		  mtot, mvir, mstars, rvir, cmv[0], cmv[1], cmv[2], vrms, 
+		  AM[0], AM[1], AM[2], spin);
+
+	  /* Write to HDF5 particle list */
+
+	  temp = (double*) malloc(3*len*sizeof(double));
+	  TempPINT = (PINT*) malloc(len*sizeof(PINT));
+
+	  index = 0;
+	  for (dim = 0; dim < 3; dim++)
+	    for (i = 0; i < len; i++, index++)
+	      temp[index] = Pbuf[i].Pos[dim] / BoxSize;
+	  for (i = 0; i < len; i++)
+	    TempPINT[i] = Pbuf[i].PartID;
+
+	  sprintf(halo_name, "Halo%8.8d", NgroupsAll-1-gr);
+	  group_id = H5Gcreate(file_id, halo_name, 0);
+	  writeScalarAttribute(group_id, H5T_NATIVE_FLOAT, "Total Mass", &mtot);
+	  writeScalarAttribute(group_id, H5T_NATIVE_FLOAT, "Stellar Mass", &mstars);
+	  writeScalarAttribute(group_id, H5T_NATIVE_FLOAT, "Spin parameter", &spin);
+	  writeScalarAttribute(group_id, H5T_NATIVE_FLOAT, "Velocity dispersion", &vrms);
+	  writeArrayAttribute(group_id, H5T_NATIVE_FLOAT, 3, "Center of mass", cm);
+	  writeArrayAttribute(group_id, H5T_NATIVE_FLOAT, 3, "Mean velocity [km/s]", cmv);
+	  writeArrayAttribute(group_id, H5T_NATIVE_FLOAT, 3, "Angular momentum [Mpc * km/s]", AM);
+
+	  hdims[0] = 3;
+	  hdims[1] = (hsize_t) len;
+	  dspace_id = H5Screate_simple(2, hdims, NULL);
+	  dset_id = H5Dcreate(group_id, "Particle Position", H5T_NATIVE_DOUBLE, dspace_id,
+			      H5P_DEFAULT);
+	  H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, (void*) temp);
+	  H5Sclose(dspace_id);
+	  H5Dclose(dset_id);
+	
+	  hdims[0] = (hsize_t) len;
+	  hdims[1] = 1;
+	  dspace_id = H5Screate_simple(1, hdims, NULL);
+	  dset_id = H5Dcreate(group_id, "Particle ID", HDF5_PINT, dspace_id,
+			      H5P_DEFAULT);
+	  H5Dwrite(dset_id, HDF5_PINT, H5S_ALL, H5S_ALL, H5P_DEFAULT, (void*) TempPINT);
+	  H5Sclose(dspace_id);
+	  H5Dclose(dset_id);
+
+	  H5Gclose(group_id);
+
+	  free(temp);
+	  free(TempPINT);
 
 	  free(Pbuf);
 	}
 
   if(ThisTask==0)
     {
+      fclose(fdtxt);
       fclose(fd);
       fclose(fdids);
       fclose(fdtypes);
+      H5Fclose(file_id);
       printf("done.\n");
     }
 }
   int i, imin, imax, pp, nlocal, nrecv;
   struct particle_data *localbuf;
 
+  if (NTask == 1) {
+    // No communication required.  Just created an array of particles
+    // from the linked list.
+
+    i = 0;
+    pp = Head[minid - Noffset[ThisTask]];
+    do {
+      buf[i++] = P[pp];
+    } while (pp = Next[pp]);
+
+    return len;
+  } // ENDIF
+
   MPI_Bcast(&minid,  1, MPI_INT, dest, MPI_COMM_WORLD);
   MPI_Bcast(&len,    1, MPI_INT, dest, MPI_COMM_WORLD);
 
 	  }
     }
 
-  MPI_Allreduce(&Ngroups,  &NgroupsAll, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-  MPI_Allreduce(&nbound,   &Nbound, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+  if (NTask == 1) {
+    NgroupsAll = Ngroups;
+    Nbound = nbound;
+  }
+  else {
+    MPI_Allreduce(&Ngroups,  &NgroupsAll, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&nbound,   &Nbound, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+  }
 
   if(ThisTask==0)
     {
 
   NgroupsList= mymalloc(sizeof(int)*NTask);
 
-  MPI_Allgather(&Ngroups, 1, MPI_INT, NgroupsList, 1, MPI_INT, MPI_COMM_WORLD);
+  if (NTask == 1)
+    NgroupsList[0] = Ngroups;
+  else
+    MPI_Allgather(&Ngroups, 1, MPI_INT, NgroupsList, 1, MPI_INT, MPI_COMM_WORLD);
 
 
   GroupDatAll= mymalloc(NgroupsAll*sizeof(struct gr_data)); 
 	      CM[3*i+2] / RootBoxSize[2] + leftEdge[2]);
     fclose(out);
 
+    free(Npart_g);
+    free(Mass);
+    free(CM);
+
   } // ENDIF root task
 
 }

src/P-GroupFinder/properties.c

 
 
 
-void get_properties(struct particle_data *p, int len, float *pcm, float *pmtot, float *pmgas, float *pmstars, float *pmsfr, float *pmcold)
+void get_properties(struct particle_data *p, int len, float *pcm, float *pmtot, 
+		    float *pmgas, float *pmstars, float *pmsfr, float *pmcold,
+		    int subgroup, float *pcmv, float *pmvir, float *prvir,
+		    float *pL, float *pvrms, float *pspin)
 {
-  int i,k;
-  double s[3];
-  double mtot, mgas, mstars, sfr, mcold;
-  
+  int i,k,dim, irvir;
+  double s[3], sv[3], L[3], delx[3], delv[3], del, vrms, spin, mvir, rvir;
+  double mtot, mgas, mstars, sfr, mcold, menc, rho, factor, rho200, r3;
+  float *radius;
+  int *pindex;
 
-  
-  s[0]= s[1]= s[2]= 0;
+  for (i = 0; i < 3; i++) {
+    s[i] = 0.0;
+    sv[i] = 0.0;
+    L[i] = 0.0;
+  }
+    
   mtot= mgas= mstars= mcold= 0;
   sfr= 0;
 
   for(i=0; i<len; i++)
     {
-      for(k=0; k<3; k++)
+      for(k=0; k<3; k++) {
 	s[k]+= p[i].Mass * periodic(p[i].Pos[k] - p[0].Pos[k]);
+	sv[k] += p[i].Mass * p[i].Vel[k];
+      }
 
       mtot+= p[i].Mass;
       switch(p[i].Type)
 	  mstars+=  p[i].Mfs;
 	  sfr+=  p[i].Sfr;
 	  break;
-	case 4:
-	  mstars+= p[i].Mass; break;
+	case 2:
+	case 5:
+	case 7:
+	  mstars+= p[i].Mass; 
+	  break;
 	}
     }
 
     {
       s[k]= s[k]/mtot;
       s[k]= periodic_wrap(s[k] + p[0].Pos[k]); 
+      sv[k] = sv[k] / mtot;
     }
 
-  for(k=0; k<3; k++)
+  for(k=0; k<3; k++) {
     pcm[k]= s[k];
+    pcmv[k] = sv[k];
+  }
+
+  /* For groups, find the virial radius (r200) and calculate virial mass */
+
+  // Sort by radius and search for an enclosed density of 200 times
+  // the critical density.  Search outside-in.
+
+  pindex = (int*) malloc(len*sizeof(int));
+
+  if (subgroup == 0) {
+    radius = (float*) malloc(len*sizeof(float));
+    for (i = 0; i < len; i++) {
+      radius[i] = 0;
+      for (dim = 0; dim < 3; dim++) {
+	del = p[i].Pos[dim] - pcm[dim];
+	radius[i] += del*del;
+      }
+      radius[i] = sqrt(radius[i]);
+    }
+    indexx(len, radius-1, pindex-1);
+
+    // Convert one-based (from indexx) to zero-based.
+    for (i = 0; i < len; i++)
+      pindex[i]--;
+
+    // Convert to Msun from 1e10 Msun and pre-compute the (4PI/3)
+    // factor.  Rho will be in units of Msun / kpc^3, as is rho_crit.
+    // Radius is comoving.
+    factor = 1e10 / (4*M_PI/3.0);
+    rho200 = 200 * RhoCritical0;
+
+    menc = mtot;
+    for (i = len-1; i >= 0; i--) {
+      menc -= p[pindex[i]].Mass;
+      r3 = radius[pindex[i]] * radius[pindex[i]] * radius[pindex[i]];
+      rho = factor * menc / max(r3, 1e-20);
+      if (rho > rho200) 
+	break;
+    }
+
+    irvir = max(i, 0);
+    rvir = radius[pindex[irvir]] * Time;  // comoving -> proper
+    mvir = menc + p[pindex[irvir]].Mass;  // Add the last particle removed
+
+    free(radius);
+
+  } // ENDIF !subgroup
+
+  else {
+
+    // Skip finding the overdense sphere for subgroups
+
+    irvir = len-1;
+    rvir = 0;
+    mvir = mtot;
+    for (i = 0; i < len; i++)
+      pindex[i] = i;
+
+  } // ENDELSE !subgroup
+
+  /* Calculate other quantities :: RMS velocity, angular momentum (Mpc * km/s),
+     spin parameter */
+
+  for (i = 0; i <= irvir; i++) {
+    k = pindex[i];
+    for (dim = 0; dim < 3; dim++) {
+      delx[dim] = p[k].Pos[dim] - pcm[dim];
+      delv[dim] = p[k].Vel[dim] - pcmv[dim];
+      vrms += p[k].Mass * delv[dim] * delv[dim];
+    }
+
+    L[0] += p[k].Mass * ( delv[1]*delx[2] - delv[2]*delx[1]);
+    L[1] += p[k].Mass * (-delv[0]*delx[2] + delv[2]*delx[0]);
+    L[2] += p[k].Mass * ( delv[0]*delx[1] - delv[1]*delx[0]);
+
+  } // ENDFOR particles
+
+  /* Divide out weight (m_vir) */
+  // 1e3 to convert from kpc to Mpc, comoving -> proper
+  for (dim = 0; dim < 3; dim++)
+    L[dim] /= mvir * 1e3 / Time;
+  vrms /= mvir;
+  vrms = sqrt(vrms);
+
+  // Convert to solar masses, 0->1 domain
+  mvir *= 1e10;
+  mtot *= 1e10;
+  mstars *= 1e10;
+  for (dim = 0; dim < 3; dim++)
+    pcm[dim] /= BoxSize;
+
+  // Spin parameter
+  float ang_mom = 0;
+  float SpinUnits = (CM_PER_MPC * 1.0e5 * 1.0e5) / (GRAVITY * SOLAR_MASS);
+  for (dim = 0; dim < 3; dim++)
+    ang_mom += L[dim] * L[dim];
+  ang_mom = sqrt(ang_mom);
+  spin = SpinUnits * ang_mom * vrms / mvir;
+
+  free(pindex);
   
   *pmtot=   mtot;
   *pmgas=   mgas;
   *pmstars= mstars;
   *pmsfr=   sfr;
   *pmcold=  mcold;
+  *pvrms = vrms;
+  *pspin = spin;
+  *pmvir = mvir;
+  *prvir = rvir;
+  for (dim = 0; dim < 3; dim++)
+    pL[dim] = L[dim];
+
 }

src/P-GroupFinder/proto.h

 void   find_minids(void);
 void   find_subgroups(void);
 int    get_particles(int dest, int minid, int len, struct particle_data *buf);
-void   get_properties(struct particle_data *p, int len, float *pcm, float *pmtot, float *pmgas, float *pmstars, float *sfr, float *mcloud);
+void get_properties(struct particle_data *p, int len, float *pcm, float *pmtot, 
+		    float *pmgas, float *pmstars, float *pmsfr, float *pmcold,
+		    int subgroup, float *pcmv, float *pmvir, float *prvir,
+		    float *pL, float *pvrms, float *pspin);
 int    get_slab(int index);
 void   iindexx(unsigned int n, int arr[], unsigned int indx[]);
 void   indexx(unsigned long n, float arr[], int indx[]);
 void   order_subgroups_by_potential(void);
 double periodic_wrap(double x);
 double periodic(double x);
-void   save_groups(char *particles_fname, char *catalogue_fname, char *parttypes_fname, char *partids_fname);
+void save_groups(char *particles_fname, char *particles_fname5, 
+		 char *catalogue_fname, char *parttypes_fname, 
+		 char *partids_fname, char *cataloguetxt);
 float  selectb(unsigned long k, unsigned long n, float arr[], int ind[]);
 void   sort2_flt_int(unsigned long n, float arr[], int brr[]);
 void   sort2_int(unsigned long n, int arr[], int brr[]);
 void   set_sph_kernel(void);
 void   stitch_together(void);
-void   subfind(char *particles_fname, char *catalogue_fname, char *subhalo_fname, char *parttypes_fname, char *partids_fname, char *subprop_fname, char *prop_fname);
+void subfind(char *particles_fname, char *catalogue_fname, 
+	     char *subhalo_fname, char *parttypes_fname, char *partids_fname, 
+	     char *subprop_fname, char *prop_fname, char *sparticles_fname5,
+	     char *scataloguetxt);
 int    unbind(int head, int len);
 void   unbind_node(int k);
 void   walk_tree_and_unbind(void);

src/P-GroupFinder/subfind.c

 #include <string.h>
 #include <math.h>
 #include <unistd.h>
+#include "h5utilities.h"
 
 
 #include "allvars.h"
 
 void subfind(char *particles_fname, char *catalogue_fname, 
 	     char *subhalo_fname, char *parttypes_fname, char *partids_fname, 
-	     char *subprop_fname, char *prop_fname)
+	     char *subprop_fname, char *prop_fname, char *sparticles_fname5,
+	     char *scataloguetxt)
 {
   MPI_Status status;
   FILE   *fd, *fdpart, *fdlen, *fdoffset, *fdparent, *fdtypes, *fdids;
   FILE   *fdsubcenter, *fdsubmtot, *fdsubmgas, *fdsubmstars, *fdsubsfr, *fdsubmcloud;
   FILE   *fdcenter, *fdmtot, *fdmgas, *fdmstars, *fdsfr, *fdmcloud;
+  FILE   *fdtxt;
   char   buflen[500], bufoffset[500], bufparent[500], bufcount[500], command[2000];
   char   bufsubcenter[500], bufsubmtot[500], bufsubmgas[500], bufsubmstars[500], bufsubsfr[500], bufsubmcloud[500], commandsubprop[2000];
   char   bufcenter[500], bufmtot[500], bufmgas[500], bufmstars[500], bufsfr[500], bufmcloud[500], commandprop[2000];
-  int    i, k, gr, task, head, len, nsubs, offset, start=0, NSubGroupsAll=0;
-  int    parent, ntot;
+  int    i, k, dim, gr, task, head, len, nsubs, offset, start=0, NSubGroupsAll=0;
+  int    parent, ntot, _idx;
   char   ctype;
-  float  cm[3], mtot, mgas, mstars, sfr, mcloud;
+  float  cm[3], cmv[3], AM[3], mtot, mgas, mstars, sfr, mcloud, redshift, spin, vrms;
+  float  mvir, rvir;
   float  corner[3];
   struct particle_data *Pbuf, *partbuf;
   int    *sublen, *suboffset, *bufsublen, *bufsuboffset;
+  int    *fsuboffset, *fbufsuboffset;
+
+  char   halo_name[200];
+  PINT    *TempPINT;
+  double *temp;
+  float  *msub, *bufmsub;
+
+  hid_t  file_id, dset_id, dspace_id, group_id;
+  hsize_t hdims[2];
 
   sprintf(buflen,    "%s.len",    subhalo_fname);
   sprintf(bufoffset, "%s.offset", subhalo_fname);
 	}
 
       fwrite(&NgroupsAll, sizeof(int), 1, fdcenter);
-    }
+
+      fprintf(stdout, "Saving subhalo list to %s\n", scataloguetxt);
+      fprintf(stdout, "Saving (sub)halo particle list to %s\n", 
+	      sparticles_fname5);
+
+      if ((fdtxt = fopen(scataloguetxt, "w")) == NULL) {
+	printf("can't open file `%s`\n", scataloguetxt);
+	MPI_Abort(MPI_COMM_WORLD, 1); exit(1);
+      }
+
+      // Write header
+
+      redshift = 1.0 / Time - 1.0;
+      fprintf(fdtxt, "# Scale factor = %lg\n", Time);
+      fprintf(fdtxt, "# Redshift = %f\n", redshift);
+      //fprintf(fdtxt, "# Number of subhalos = %"ISYM"\n", AllVars.NgroupsAll);
+      fprintf(fdtxt, "#\n");
+      fprintf(fdtxt, "# Column 1.  Center of mass (x)\n");
+      fprintf(fdtxt, "# Column 2.  Center of mass (y)\n");
+      fprintf(fdtxt, "# Column 3.  Center of mass (z)\n");
+      fprintf(fdtxt, "# Column 4.  Subhalo number\n");
+      fprintf(fdtxt, "# Column 5.  Parent halo number\n");
+      fprintf(fdtxt, "# Column 6.  First particle in halo particle list\n");
+      fprintf(fdtxt, "#            --> All subgroup particles are consecutively listed in\n");
+      fprintf(fdtxt, "#                particle list (if written)\n");
+      fprintf(fdtxt, "# Column 7.  Number of particles\n");
+      fprintf(fdtxt, "# Column 8.  Halo mass [solar masses]\n");
+      fprintf(fdtxt, "# Column 9.  Stellar mass [solar masses]\n");
+      fprintf(fdtxt, "# Column 10. Mean x-velocity [km/s]\n");
+      fprintf(fdtxt, "# Column 11. Mean y-velocity [km/s]\n");
+      fprintf(fdtxt, "# Column 12. Mean z-velocity [km/s]\n");
+      fprintf(fdtxt, "# Column 13. Velocity dispersion [km/s]\n");
+      fprintf(fdtxt, "# Column 14. Mean x-angular momentum [Mpc * km/s]\n");
+      fprintf(fdtxt, "# Column 15. Mean y-angular momentum [Mpc * km/s]\n");
+      fprintf(fdtxt, "# Column 16. Mean z-angular momentum [Mpc * km/s]\n");
+      fprintf(fdtxt, "# Column 17. Spin parameter\n");
+      fprintf(fdtxt, "#\n");
+      fprintf(fdtxt, "# datavar lines are for partiview.  Ignore them if you're not partiview.\n");
+      fprintf(fdtxt, "#\n");
+      fprintf(fdtxt, "datavar 1 subhalo_number\n");
+      fprintf(fdtxt, "datavar 2 parent_number\n");
+      fprintf(fdtxt, "datavar 3 particle_offset\n");
+      fprintf(fdtxt, "datavar 4 number_of_particles\n");
+      fprintf(fdtxt, "datavar 5 halo_mass\n");
+      fprintf(fdtxt, "datavar 6 stellar_mass\n");
+      fprintf(fdtxt, "datavar 7 x_velocity\n");
+      fprintf(fdtxt, "datavar 8 y_velocity\n");
+      fprintf(fdtxt, "datavar 9 z_velocity\n");
+      fprintf(fdtxt, "datavar 10 velocity_dispersion\n");
+      fprintf(fdtxt, "datavar 11 x_angular_momentum\n");
+      fprintf(fdtxt, "datavar 12 y_angular_momentum\n");
+      fprintf(fdtxt, "datavar 13 z_angular_momentum\n");
+      fprintf(fdtxt, "datavar 14 spin\n");
+      fprintf(fdtxt, "\n");
+
+      file_id = H5Fcreate(sparticles_fname5, H5F_ACC_TRUNC, H5P_DEFAULT, 
+			  H5P_DEFAULT);
+      group_id = H5Gcreate(file_id, "/Parameters", 0);
+      writeScalarAttribute(group_id, H5T_NATIVE_FLOAT, "Redshift", &redshift);
+      writeScalarAttribute(group_id, H5T_NATIVE_DOUBLE, "Scale Factor", &Time);
+      writeScalarAttribute(group_id, H5T_NATIVE_INT, "Number of groups", &NgroupsAll);
+      H5Gclose(group_id);
+
+    } // ENDIF TASK 0
 
 
   if(ThisTask==0)
 	      len=  GroupDatAll[gr-task].Len;
 	      sublen=    mymalloc(len/DesLinkNgb*sizeof(int));
 	      suboffset= mymalloc(len/DesLinkNgb*sizeof(int));
+	      fsuboffset= mymalloc(len/DesLinkNgb*sizeof(int));
+	      msub= mymalloc(len/DesLinkNgb*sizeof(float));
 	      Pbuf=      mymalloc(sizeof(struct particle_data)*len);
 	    }
 	  get_particles(task, head, len, Pbuf);
 
       for(task=0; task<NTask && (gr-task)>=0; task++)
 	{
+	  parent = NgroupsAll-(gr-task)-1;
 	  if(ThisTask==0)
 	    {
 	      if(task==0)
 		{
 		  for(i=0; i<nsubs; i++)
 		    {
-		      get_properties(Pbuf+suboffset[i], sublen[i], cm, &mtot, &mgas, &mstars, &sfr, &mcloud);
+		      get_properties(Pbuf+suboffset[i], sublen[i], cm, 
+				     &mtot, &mgas, &mstars, &sfr, &mcloud,
+				     1, &cmv[0], &mvir, &rvir, AM, &vrms,
+				     &spin);
+		      msub[i] = mtot;
 		      fwrite(cm,      sizeof(float), 3, fdsubcenter);
 		      fwrite(&mtot,   sizeof(float), 1, fdsubmtot);
 		      fwrite(&mgas,   sizeof(float), 1, fdsubmgas);
 		      fwrite(&mstars, sizeof(float), 1, fdsubmstars);
 		      fwrite(&sfr,    sizeof(float), 1, fdsubsfr);
 		      fwrite(&mcloud,    sizeof(float), 1, fdsubmcloud);
+
+		      /* Write to ASCII halo catalog */
+
+		      fprintf(fdtxt, "%12.6g %12.6g %12.6g %12d %12d %12d %12d "
+			      "%12.6g %12.6g %12.6g %12.6g %12.6g %12.6g %12.6g "
+			      "%12.6g %12.6g %12.6g\n",
+			      cm[0], cm[1], cm[2], i, parent, suboffset[i], sublen[i], 
+			      mtot, mstars, cmv[0], cmv[1], cmv[2], vrms, AM[0], 
+			      AM[1], AM[2], spin);
+
 		    }
 
-		  get_properties(Pbuf, GroupDatAll[gr-task].Len, cm, &mtot, &mgas, &mstars, &sfr, &mcloud);
+		  for(i=0; i<nsubs; i++) {
+		    fsuboffset[i] = suboffset[i];
+		    suboffset[i]+= start;
+		  }
+		  start+= GroupDatAll[gr-task].Len;
+
+		  len = GroupDatAll[gr-task].Len;
+		  temp = (double*) malloc(3*len*sizeof(double));
+		  TempPINT = (PINT*) malloc(len*sizeof(PINT));
+		  _idx = 0;
+		  for (dim = 0; dim < 3; dim++)
+		    for (i = 0; i < len; i++, _idx++)
+		      temp[_idx] = Pbuf[i].Pos[dim] / BoxSize;
+		  for (i = 0; i < len; i++)
+		    TempPINT[i] = Pbuf[i].PartID;
+
+		  get_properties(Pbuf, GroupDatAll[gr-task].Len, cm, &mtot, 
+				 &mgas, &mstars, &sfr, &mcloud, 1,
+				 &cmv[0], &mvir, &rvir, AM, &vrms, &spin);
 		  fwrite(cm,      sizeof(float), 3, fdcenter);
 		  fwrite(&mtot,   sizeof(float), 1, fdmtot);
 		  fwrite(&mgas,   sizeof(float), 1, fdmgas);
 		  fwrite(&sfr,    sizeof(float), 1, fdsfr);
 		  fwrite(&mcloud,    sizeof(float), 1, fdmcloud);
 
-		  for(i=0; i<nsubs; i++)
-		    suboffset[i]+= start;
-		  start+= GroupDatAll[gr-task].Len;
+		  /* Write to HDF5 particle list */
+
+		  sprintf(halo_name, "Halo%8.8d", parent);
+		  group_id = H5Gcreate(file_id, halo_name, 0);
+		  writeScalarAttribute(group_id, H5T_NATIVE_INT, "NumberOfSubhalos", &nsubs);
+		  writeScalarAttribute(group_id, H5T_NATIVE_FLOAT, "Total Mass", &mtot);
+		  writeScalarAttribute(group_id, H5T_NATIVE_FLOAT, "Stellar Mass", &mstars);
+		  writeScalarAttribute(group_id, H5T_NATIVE_FLOAT, "Spin parameter", &spin);
+		  writeScalarAttribute(group_id, H5T_NATIVE_FLOAT, "Velocity dispersion", &vrms);
+		  writeArrayAttribute(group_id, H5T_NATIVE_FLOAT, 3, "Center of mass", cm);
+		  writeArrayAttribute(group_id, H5T_NATIVE_FLOAT, 3, "Mean velocity [km/s]", cmv);
+		  writeArrayAttribute(group_id, H5T_NATIVE_FLOAT, 3, "Angular momentum [Mpc * km/s]", AM);
+
+		  if (nsubs > 0) {
+		    hdims[0] = (hsize_t) nsubs;
+		    hdims[1] = 1;
+		    dspace_id = H5Screate_simple(1, hdims, NULL);
+		    dset_id = H5Dcreate(group_id, "Subhalo Mass", H5T_NATIVE_FLOAT, dspace_id,
+					H5P_DEFAULT);
+		    H5Dwrite(dset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, 
+			     (void*) msub);
+		    H5Sclose(dspace_id);
+		    H5Dclose(dset_id);
+
+		    dspace_id = H5Screate_simple(1, hdims, NULL);
+		    dset_id = H5Dcreate(group_id, "Subhalo Size", H5T_NATIVE_INT, dspace_id,
+					H5P_DEFAULT);
+		    H5Dwrite(dset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, 
+			     (void*) sublen);
+		    H5Sclose(dspace_id);
+		    H5Dclose(dset_id);
+
+		    dspace_id = H5Screate_simple(1, hdims, NULL);
+		    dset_id = H5Dcreate(group_id, "Subhalo Offset", H5T_NATIVE_INT, dspace_id,
+					H5P_DEFAULT);
+		    H5Dwrite(dset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, 
+			     (void*) fsuboffset);
+		    H5Sclose(dspace_id);
+		    H5Dclose(dset_id);
+
+		  } // ENDIF nsubs>0
+
+		  hdims[0] = 3;
+		  hdims[1] = (hsize_t) len;
+		  dspace_id = H5Screate_simple(2, hdims, NULL);
+		  dset_id = H5Dcreate(group_id, "Particle Position", H5T_NATIVE_DOUBLE,
+				      dspace_id, H5P_DEFAULT);
+		  H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, 
+			   H5P_DEFAULT, (void*) temp);
+		  H5Sclose(dspace_id);
+		  H5Dclose(dset_id);
+	
+		  hdims[0] = (hsize_t) len;
+		  hdims[1] = 1;
+		  dspace_id = H5Screate_simple(1, hdims, NULL);
+		  dset_id = H5Dcreate(group_id, "Particle ID", HDF5_PINT, dspace_id,
+				      H5P_DEFAULT);
+		  H5Dwrite(dset_id, HDF5_PINT, H5S_ALL, H5S_ALL, H5P_DEFAULT, 
+			   (void*) TempPINT);
+		  H5Sclose(dspace_id);
+		  H5Dclose(dset_id);
+
+		  H5Gclose(group_id);
+
+		  free(temp);
+		  free(TempPINT);
+
+		  /* END HDF5 writing */
 
 		  fwrite(sublen, sizeof(int), nsubs, fdlen);
 		  fwrite(suboffset, sizeof(int), nsubs, fdoffset);
 		  MPI_Recv(&nsubs, 1, MPI_INT, task, task, MPI_COMM_WORLD, &status);
 		  if(nsubs)
 		    {
+		      bufmsub=      mymalloc(nsubs*sizeof(float));
 		      bufsublen=    mymalloc(nsubs*sizeof(int));
 		      bufsuboffset= mymalloc(nsubs*sizeof(int));
+		      fbufsuboffset= mymalloc(nsubs*sizeof(int));
 		      MPI_Recv(bufsublen,    nsubs, MPI_INT, task, task, MPI_COMM_WORLD, &status);
 		      MPI_Recv(bufsuboffset, nsubs, MPI_INT, task, task, MPI_COMM_WORLD, &status);
-		      for(i=0; i<nsubs; i++)
+		      for(i=0; i<nsubs; i++) {
+			fbufsuboffset[i] = bufsuboffset[i];
 			bufsuboffset[i]+= start;
+		      }
 		      fwrite(bufsublen,    sizeof(int), nsubs, fdlen);
 		      fwrite(bufsuboffset, sizeof(int), nsubs, fdoffset);
 		    }
 		  
 		  for(i=0; i<nsubs; i++)
 		    {
-		      get_properties(partbuf+bufsuboffset[i]-start, bufsublen[i], cm, &mtot, &mgas, &mstars, &sfr, &mcloud);
+		      get_properties(partbuf+bufsuboffset[i]-start, bufsublen[i], cm, &mtot, 
+				     &mgas, &mstars, &sfr, &mcloud,
+				     1, cmv, &mvir, &rvir, AM, &vrms, &spin);
+		      bufmsub[i] = mtot;
 		      fwrite(cm,      sizeof(float), 3, fdsubcenter);
 		      fwrite(&mtot,   sizeof(float), 1, fdsubmtot);
 		      fwrite(&mgas,   sizeof(float), 1, fdsubmgas);
 		      fwrite(&mstars, sizeof(float), 1, fdsubmstars);
 		      fwrite(&sfr,    sizeof(float), 1, fdsubsfr);
 		      fwrite(&mcloud,    sizeof(float), 1, fdsubmcloud);
+
+		      /* Write to ASCII halo catalog */
+
+		      fprintf(fdtxt, "%12.6g %12.6g %12.6g %12d %12d %12d %12d "
+			      "%12.6g %12.6g %12.6g %12.6g %12.6g %12.6g %12.6g "
+			      "%12.6g %12.6g %12.6g\n",
+			      cm[0], cm[1], cm[2], i, parent, suboffset[i], sublen[i], 
+			      mtot, mstars, cmv[0], cmv[1], cmv[2], vrms, AM[0], 
+			      AM[1], AM[2], spin);
+
 		    }
 
-		  get_properties(partbuf, GroupDatAll[gr-task].Len, cm, &mtot, &mgas, &mstars, &sfr, &mcloud);
+		  get_properties(partbuf, GroupDatAll[gr-task].Len, cm, &mtot, 
+				 &mgas, &mstars, &sfr, &mcloud, 1,
+				 cmv, &mvir, &rvir, AM, &vrms, &spin);
 		  fwrite(cm,      sizeof(float), 3, fdcenter);
 		  fwrite(&mtot,   sizeof(float), 1, fdmtot);
 		  fwrite(&mgas,   sizeof(float), 1, fdmgas);
 		  fwrite(&sfr,    sizeof(float), 1, fdsfr);
 		  fwrite(&mcloud,    sizeof(float), 1, fdmcloud);
 
+		  /* Write to HDF5 particle list */
+
+		  len = GroupDatAll[gr-task].Len;
+		  parent = NgroupsAll-(gr-task)-1;
+		  temp = (double*) malloc(3*len*sizeof(double));
+		  TempPINT = (PINT*) malloc(len*sizeof(PINT));
+		  _idx = 0;
+		  for (dim = 0; dim < 3; dim++)
+		    for (i = 0; i < len; i++, _idx++)
+		      temp[_idx] = Pbuf[i].Pos[dim] / BoxSize;
+		  for (i = 0; i < len; i++)
+		    TempPINT[i] = Pbuf[i].PartID;
+
+		  sprintf(halo_name, "Halo%8.8d", parent);