CH_HDF5.H

Go to the documentation of this file.

#ifdef CH_LANG_CC
/*
 *      _______              __
 *     / ___/ /  ___  __ _  / /  ___
 *    / /__/ _ \/ _ \/  V \/ _ \/ _ \
 *    \___/_//_/\___/_/_/_/_.__/\___/
 *    Please refer to Copyright.txt, in Chombo's root directory.
 */
#endif

#ifndef _CH_HDF5_H_
#define _CH_HDF5_H_

#define CHOFFSET(object, member) (int)((char*)&(object.member) - (char*)&object)

#ifdef CH_USE_HDF5  // if you don't have CH_USE_HDF5, then this file is useless

#include <iostream>
using std::cout;
using std::endl;

#ifdef CH_MPI
#include "mpi.h"
#endif

#include "LevelData.H"
#include "HDF5Portable.H"
// hdf5 #defines inline.... duh!
#undef inline
#include <string>
#include <map>
#include "RealVect.H"
#include "CH_Timer.H"
#include "LoadBalance.H"
#include "LayoutIterator.H"
#include "Vector.H"
#include "memtrack.H"
#include "FluxBox.H"
#include "EdgeDataBox.H"
#include "NodeFArrayBox.H"

#ifdef CH_MULTIDIM
#include "BoxTools_ExternC_Mangler.H" // Generated by lib/utils/multidim/mangle_externs.sh
#endif
#include "NamespaceHeader.H"
#ifdef H5_USE_16_API
#define H516
#endif

template <class T>
void read(T& item, Vector<Vector<char> >& a_allocatedBuffers, const Box& box,
          const Interval& comps)
{
  Vector<void*> b(a_allocatedBuffers.size());
  for (int i=0; i<b.size(); ++i) b[i]= &(a_allocatedBuffers[i][0]);
  read(item, b, box, comps);
}

namespace CH_HDF5
{
  enum IOPolicy
  {
    IOPolicyDefault           = 0,
    IOPolicyMultiDimHyperslab = (1<<0),  // HDF5 will internally linearize the
                                         // data
    IOPolicyCollectiveWrite   = (1<<1)   // HDF5 will collectively write data
  };
}

using std::map;

class HDF5Handle;

// CH_HDF5.H
// ============

/// user-friendly function to write out data on a AMR level
/**
   all data is IN data.
   a_handle:  handle open and ready.
   a_data:    data, what else do you want to know ?
   a_dx:      the grid spacing at this level
   a_dt:      the timestep size that was last completed
   a_time:    the time of this level (might not be the same as other levels)
   a_domain:  the problem domain, represented at this level of refinement
   a_refRatio:the refinement of a_level+1 wrt a_level. for vis systems it
   would probably help if you use 1 for a_level==max_level.
*/
template <class T>
int writeLevel(HDF5Handle& a_handle,
               const int& a_level,
               const T& a_data,
               const Real& a_dx,
               const Real& a_dt,
               const Real& a_time,
               const Box& a_domain,
               const int& a_refRatio,
               const IntVect& outputGhost = IntVect::Zero,
               const Interval& comps = Interval());

template <class T>
int readLevel(HDF5Handle& a_handle,
              const int& a_level,
              LevelData<T>& a_data,
              Real& a_dx,
              Real& a_dt,
              Real& a_time,
              Box& a_domain,
              int& a_refRatio,
              const Interval& a_comps = Interval(),
              bool setGhost = false);

template <class T>
int readLevel(HDF5Handle& a_handle,
              const int& a_level,
              LevelData<T>& a_data,
              RealVect& a_dx,
              Real& a_dt,
              Real& a_time,
              Box& a_domain,
              IntVect& a_refRatio,
              const Interval& a_comps = Interval(),
              bool setGhost = false);

// More basic HDF5 functions.  These can be used at a persons own
// discretion.  Refer to the User's Guide for a description of how these
// functions interact with the convenience functions writeLevel, etc.

/// writes BoxLayout to HDF5 file.
/**
    Writes BoxLayout to HDF5 file.  Only one BoxLayout per group is permitted, this operation overwrites
    previous entries.
    This operation assumes boxes are cell-centered.\\
    returns: success:    0\\
    HDF5 error: negative error code.\\
*/
int write(HDF5Handle& a_handle,
          const BoxLayout& a_layout,
          const std::string& name = "boxes");

/// writes a BoxLayoutData<T> to an HDF5 file.
/**
   writes a BoxLayoutData<T> to an HDF5 file.\\
   returns: success:    0\\
   HDF5 error: negative error code.\\
*/
template <class T>
int write(HDF5Handle& a_handle,
          const BoxLayoutData<T>& a_data,
          const std::string& a_name,
          IntVect outputGhost = IntVect::Zero,
          const Interval& comps = Interval(),
          bool  newForm = false);

/// writes a LevelData<T> to an HDF5 file.
/**
   Writes a LevelData<T> to an HDF5 file.
   the DisjointBoxLayout is not written out with the data, the user is required to
   handle that object seperately. (see the "read" interface below).\\
   returns: success:    0\\
   HDF5 error: negative error code.\\
*/
template <class T>
int write(HDF5Handle& a_handle,
          const LevelData<T>& a_data,
          const std::string& a_name,
          const IntVect& outputGhost = IntVect::Zero,
          const Interval& comps = Interval());

/// reads Vector<Box> from location specified by a_handle.
/**
     Reads BoxLayout from the group specified by a_handle.
     Only one BoxLayout per group is permitted, this operation overwrites.
     This operation assumes boxes are cell-centered.\\
     returns: success:    0\\
     HDF5 error: negative error code.\\

     Arg name refers to the name of an HDF5 dataset in which the boxes you want
     are stored in a particular format, which is the one used to output DataLayouts
     to HDF5.  Here is an example of that format (as shown by h5dump):

      DATASET "datalayout" {
         DATATYPE  H5T_COMPOUND {
            H5T_STD_I32LE "lo_i";
            H5T_STD_I32LE "lo_j";
            H5T_STD_I32LE "hi_i";
            H5T_STD_I32LE "hi_j";
         }
         DATASPACE  SIMPLE { ( 2 ) / ( 2 ) }
         DATA {
         (0): {
               0,
               0,
               1,
               0
            },
         (1): {
               2,
               0,
               2,
               1
            }
         }
      }

*/
int read(HDF5Handle& a_handle,
         Vector<Box>& boxes,
         const std::string& name = "boxes");

/// reads the set of Boxes out from the level_* groups of a Chombo HDF5 AMR file
/**
   goes to all groups named level_n for level n = 0 to numLevel and fills the
   vector of boxes.  CH_HDF5.Handle must be set to the root group (the default location when the
   file is just opened.

   A handy *get everything* version of int read(HDF5Handle& a_handle, Vector<Box>& boxes)
*/
int readBoxes(HDF5Handle& a_handle,
              Vector<Vector<Box> >& boxes);

/// FArrayBox-at-a-time read function.  FArrayBox gets redefined in the function.  Reads data field named by a_dataName.
/**  FArrayBox gets redefined in the function.
     it will be sized to the box at the [level,boxNumber] position and have components
     the go from [0,a_components.size()-1]

     some meta-data has to be read again and again internally to do this function, but
     it should make some out-of-core tools possible.
*/

int readFArrayBox(HDF5Handle& a_handle,
                  FArrayBox&  a_fab,
                  int a_level,
                  int a_boxNumber,
                  const Interval& a_components,
                  const std::string& a_dataName = "data" );

/// read BoxLayoutData named a_name from location specified by a_handle.
/**
    Read BoxLayoutData named a_name from location specified by a_handle.  User must supply the correct BoxLayout for this function if redefineData == true.  \\
    returns: success:      0\\
    bad location: 1\\
    HDF5 error:   negative error code.\\
*/
template <class T>
int read(HDF5Handle& a_handle,
         BoxLayoutData<T>& a_data,
         const std::string& a_name,
         const BoxLayout& a_layout,
         const Interval&  a_comps = Interval(),
         bool redefineData = true);

/// read LevelData named a_name from location specified by a_handle.
/**

    Read LevelData named a_name from location specified by a_handle.
    User must supply the correct BoxLayout for this function.

    Arg a_name is significant: the HDF5 group to which a_handle is set is
    is assumed to contain a dataset called <a_name>:datatype=<some integer>,
    a dataset called <a_name>:offsets=<some integer>, and a subgroup named
    <a_name>_attributes.  You will have all these items if you dumped your
    LevelData out using the corresponding write() function defined here.

    If arg redefineData==false, then the user must pass in a valid  LevelData.
    Otherwise, this function figures out how many components and ghost cells there
    are, and allocates the correct amount of space.  The actual FArray data held
    by the LevelData gets filled in here, regardless of redefineData; "redefine"
    alludes to the family of define() functions.

    returns: success:      0\\
    bad location: 1\\
    HDF5 error:   negative error code.\\
*/
template <class T>
int read(HDF5Handle& a_handle,
         LevelData<T>& a_data,
         const std::string& a_name,
         const DisjointBoxLayout& a_layout,
         const Interval& a_comps = Interval(),
         bool redefineData = true);

/// Handle to a particular group in an HDF file.
/**
    HDF5Handle is a handle to a particular group in an HDF file.  Upon
    construction, group is defined to be the root.  All data is
    written and read assuming the native representations for data on
    the architecture it is running on.  When a file is opened, these
    settings are checked and an error is flagged when things don't
    match up.  It is the USER'S responsibility to close() this object
    when it is no longer needed.

*/
class HDF5Handle
{
public:
  ///
  /**
     Enumeration of opening modes for HDF5 files.  \\

     CREATE: file is created if it didn't exist, or an existing file of
     the same name is clobbered.\\

     CREATE_SERIAL: in serial execution this is equivalent to CREATE. In parallel this opens
     a file just on this particular calling processor.  Working with LevelData and this kind of
     file is difficult to get right and most users will not use this mode.  \\

     OPEN_RDONLY: existing file is opened in read-only mode.  If the
     file doesn't already exist then open fails and isOpen() returns
     false.\\

     OPEN_RDWR: existing file is opened in read-write mode.  If the file
     doesn't already exist then open fails and isOpen() returns false.\\

  */
  enum mode
  {
    CREATE,
    CREATE_SERIAL,
    OPEN_RDONLY,
    OPEN_RDWR
  };

  ///    {\bf constructor}

  ///
  /**
     Default constructor.  User must call open() prior to using
     constructed object.
  */
  HDF5Handle();

  ///
  /** Opens file and sets the current group to the root "/" group. \\

  if mode == CREATE, then file is created if it didn't exist, or an
  existing file of the same name is clobbered.\\

  if mode == OPEN_*, then existing file is opened, if the file doesn't
  already exist then open fails and isOpen() returns false.\\

  Writes basic information such as SpaceDim and testReal to a global
  group. The global group is named Chombo_global in the call signature
  where no globalGroupName is passed.\\
  */
  HDF5Handle(
        const std::string& a_filename,
        mode a_mode,
        const char *a_globalGroupName="Chombo_global");

  // HDF5Handle(const std::string& a_filename, mode a_mode);

  ~HDF5Handle();

  /// {\bf File functions}

  ///
  /**
      Opens file and sets the current group of this HDF5Handle to the
      root "/" group.  File that this HDF5Handle previously pointed at is
      NOT closed, that is the users responsibility.\\

      if mode == CREATE, then file is created if it didn't exist, or an
      existing file of the same name is clobbered.\\

      if mode == OPEN_*, then existing file is opened, if the file doesn't
      already exist then open fails and isOpen() returns false.\\

      Writes basic information such as SpaceDim and testReal to a global
      group. The global group is named Chombo_global in the call signature
      where no globalGroupName is passed.\\

      returns:\\
      0  on success\\
      negative number if file open failed (return code from HDF5)\\
      1  file does not appear to contain datacheck info, probably not a Chombo file\\
      2  on data bit size differences between code and file.\\

      aborts on SpaceDim not matching between code and file\\

  */
  int open(
        const std::string& a_filename,
        mode a_mode,
        const char *a_globalGroupName="Chombo_global");

  // int open(const std::string& a_filename, mode a_mode);

  ///
  /**
     A NULL or failed constructed HDF5Handle will return false.
  */
  bool isOpen() const;

  ///
  /**
     Closes the file.  Must be called to close file.  Files are not
     automatically closed.

  */
  void close();

  /// {\bf Group functions}

  ///
  /**
     Sets the current group to be "/level_x" where x=a_level.
  */
  void setGroupToLevel(int a_level);

  ///
  /**
      Set group to users choice, referenced from file root.
      groupAbsPath will look like a Unix file path:
      "/mySpecialData/group1/" "/" is the root group of the
      file. returns a negative value on failure

  */
  int setGroup(const std::string& groupAbsPath);

  ///
  /**
      Add the indicated string to the group path.  For example, if
      getGroup() returns "/foo" then after pushGroup("bar"), getGroup()
      will return "/foo/bar".
      Return value is whatever the internal setGroup returns.
  */
  int pushGroup( const std::string& grp );

  ///
  /**
      Pop off the last element of the group path.  For example, "/foo/bar"
      becomes "/foo".  It's an error to call this function if the group,
      going in, is "/".
      Return value is whatever the internal setGroup returns, when we call
      it to reset the group's path.
  */
  int popGroup();

  ///
  /**
     Returns name of current group.  groupAbsPath will look like a Unix
     file path: "/mySpecialData/group1/" "/" is the root group of the
     file.

  */
  const std::string& getGroup() const;

  HDF5Handle::mode openMode() const {return m_mode;}
  const hid_t& fileID() const;
  const hid_t& groupID() const;
  static hid_t box_id;
  static hid_t intvect_id;
  static hid_t realvect_id;
  static map<std::string, std::string> groups;

private:

  HDF5Handle(const HDF5Handle&);
  HDF5Handle& operator=(const HDF5Handle&);

  hid_t         m_fileID;
  hid_t         m_currentGroupID;
  mode          m_mode;
  bool          m_isOpen;
  std::string   m_filename; // keep around for debugging
  std::string   m_group;
  int           m_level;

  //  static hid_t  file_access;
  static bool   initialized;
  static void   initialize();

};

/// data to be added to HDF5 files.
/**
   HDF5HeaderData is a wrapper for some data maps to be added to HDF5
   files.  instead of an overdose of access functions, the maps are
   made public and they can be manipulated by the user at will.  They
   maintain type safety.

   to add a Real data entry, a user can simply program as follows:

   <PRE>
   Real dx;
   .
   .
   HDF5HeaderData metaData;
   metaData.m_real["dx"] = dx;
   </PRE>

   If "dx" already existed, then it is overwritten, otherwise an entry is
   created and added with the new value;

   To search for entries, the user does the following:
   <PRE>

   HDF5HeaderData metaData;
   HDF5Handle currentStep(filename);
   currentStep.setGroupToLevel(0);
   metaData.readFromFile(currentStep);
   if (metaData.m_intvect.find("ghost") != metaData.m_intvect.end())
   ghost = metaData.m_intvect["ghost"];
   else
   ghost = defaultGhostIntVect;

   </PRE>

   A user can skip the check for existence if they have reason to "know" the
   data will be there.  It is just good coding practice.

   To erase an entry, you can use:
   <PRE>
   metaData.m_real.erase("dx");
   </PRE>

*/
class HDF5HeaderData
{
public:

  ///
  /**
      Writes this HDF5HeaderData's current attribute list to the
      current group in 'file.'  Returns 0 on success, returns the
      error code from HDF5 on failure.

  */
  int writeToFile(HDF5Handle& file) const;

  ///
  /**
      Reads into this HDF5HeaderData's attribute list from file.  Read
      process is add/change, does not remove key-value pairs. Reads
      from current group.  Returns 0 on success, positive number if a
      particular member of group caused an error, negative on general
      error.

  */
  int readFromFile(HDF5Handle& file);

  ///
  void clear();

  ///
  map<std::string, Real>        m_real;

  ///
  map<std::string, int>         m_int;

  ///
  map<std::string, std::string> m_string;

  ///
  map<std::string, IntVect>     m_intvect;

  ///
  map<std::string, Box>         m_box;

  ///
  map<std::string, RealVect>    m_realvect;

  //users should not need these functions in general

  int writeToLocation(hid_t loc_id) const;
  int readFromLocation(hid_t loc_id);

  /// useful for debugging.  dumps contents to std::cout
  void dump() const;

private:
  static herr_t attributeScan(hid_t loc_id, const char *name, void *opdata);
};

extern "C"
{
#ifdef H516
  herr_t HDF5HeaderDataattributeScan(hid_t loc_id, const char *name, void *opdata);
#else
  herr_t HDF5HeaderDataattributeScan(hid_t loc_id, const char *name, const H5A_info_t* info, void *opdata);
#endif
}

std::ostream& operator<<(std::ostream& os, const HDF5HeaderData& data);

/// Methods for writing multiple LevelData to an HDF5 file.
/**
    While the write functions can be used to write a single LevelData, they
    do not support writing intervals from multiple LevelData.  This is often
    the case when you want to add diagnostic information to the output file.
    The purpose of this class is to handle writing multiple LevelData to
    a single file.  The dataset for the LevelData is created during
    construction, when the boxes are written.  Any number of intervals from
    LevelData can then be written.

    Special behavior can be selected with the policy flags:
    <ul>
      <li> CH_HDF5::IOPolicyMultiDimHyperslab - The memory dataspace will be
           set up so that HDF5 linearizes the data.  Using this requires
           T.dataPtr() and contiguous memory (so things like T=FluxBox will
           not work).
      <li> CH_HDF5::IOPolicyCollectiveWrite - The write for each dataset will
           be collective.  The write is always parallel but otherwise will
           be done independently.  May or may not provide a speedup but it
           should work.  Rumor has it the BG requires this.
    </ul>

    With MPI-everywhere on x86 architectures, it is highly likely that
    the default policy CH_HDF5::IOPolicyDefault gives the best performance
    (i.e., use internal linearization and independent I/O).  In some limited
    testing (2 procs, 20 cores each):
      IOPolicyDefault - 6.77
      IOPolicyCollectiveWrite - 33.33
      IOPolicyMultiDimHyperslab - 168.14
      IOPolicyCollectiveWrite | IOPolicyMultiDimHyperslab - 30.94
    which is difficult to make sense of.

    \note
    <ul>
      <li> Users are resposible for matching the total number of components
           used to allocate the dataset match the number of levelData
           added to the dataset
      <li> Only a single datatype is supported per file
      <li> Since the dataset must be created in advance, T can only be
           statically allocatable (T::preAllocatable() == 0)
      <li> API 1.6 is not supported
    </ul>
*/
template <class T>
class WriteMultiData
{

public:

  /// Constructor writes boxes and allocates dataset for LevelData
  WriteMultiData(HDF5Handle&      a_handle,
                 const BoxLayout& a_layout,
                 const int        a_numIntv,
                 const string&    a_name,
                 const int        a_policyFlags = CH_HDF5::IOPolicyDefault,
                 const IntVect&   a_outputGhost = IntVect::Zero,
                 const bool       a_newForm = false);

  /// Destructor
  ~WriteMultiData()
    {
      H5Dclose(m_dataSet);
    }

private:

  // Copy and assignment not allowed
  WriteMultiData(const WriteMultiData&);
  WriteMultiData& operator=(const WriteMultiData);

public:

  /// Write an interval of LevelData to the dataset
  int writeData(const BoxLayoutData<T>& a_data,
                const Interval&         a_intvMem,
                const Interval&         a_intvFile);

protected:

  const int m_policyFlags;            ///< Policies
  const IntVect m_outputGhost;        ///< Number of ghost cells written
  const bool m_newForm;               ///< ?

  char m_dataname[128];               ///< Name for level data dataset
  hid_t m_dataSet;                    ///< Dataset for level data
  Interval m_allIntvFile;             ///< Interval for components in file
  long m_maxBoxPerProc;               ///< Maximum boxes written by any proc

  // Both of these vectors are size 1 in the outermost dimension since
  // only 1 type is allowed.
  Vector<Vector<long long> > m_offsets;
                                      ///< Offset in file for each process to
                                      ///< write to
  Vector<hid_t> m_types;              ///< Type of data written
};

//=============================================================================
//
// end of declarations.
//
//=============================================================================

#if ( H5_VERS_MAJOR == 1 && H5_VERS_MINOR > 6 )
typedef hsize_t ch_offset_t;
#else
#if ( H5_VERS_MAJOR == 1 && H5_VERS_MINOR == 6 && H5_VERS_RELEASE >= 4 )
typedef hsize_t ch_offset_t;
#else
typedef hssize_t ch_offset_t;
#endif
#endif

template<class T>
hid_t H5Type(const T* dummy);

template< >
hid_t H5Type(const int* dummy);

template< >
hid_t H5Type(const long long* dummy);

template< >
hid_t H5Type(const float* dummy);

template< >
hid_t H5Type(const double* dummy);

template< >
hid_t H5Type(const Box* dummy);

template< >
hid_t H5Type(const RealVect* dummy);

template< >
hid_t H5Type(const IntVect* dummy);

template<class T>
hid_t H5Type(const T* dummy)
{
  // no such definition;
  MayDay::Error(" H5Type(const T* dummy)");
  return -4;
}

void createData(hid_t& a_dataset,
                hid_t& a_dataspace,
                HDF5Handle& handle,
                const std::string& name,
                hid_t type,
                hsize_t size);

template <class T>
void createDataset(hid_t& a_dataset,
                   hid_t& a_dataspace,
                   HDF5Handle& handle,
                   const std::string& name,
                   const T* dummy,
                   hsize_t size)
{
  createData(a_dataset, a_dataspace, handle, name, H5Type(dummy), size);

}

void writeDataset(hid_t a_dataset,
                  hid_t a_dataspace,
                  const void* start,
                  ch_offset_t off,
                  hsize_t  count);

void readDataset(hid_t a_dataset,
                 hid_t a_dataspace,
                 void* start,
                 ch_offset_t off,
                 hsize_t  count);

// non-user code used in implementation of communication

struct OffsetBuffer
{
  Vector<int> index;
  Vector<Vector<int> > offsets;
  void operator=(const OffsetBuffer& rhs);
};

ostream& operator<<(ostream& os, const OffsetBuffer& ob);

#include "NamespaceFooter.H"

#include "BaseNamespaceHeader.H"

#include "NamespaceVar.H"

//OffsetBuffer specialization of linearSize
template < >
int linearSize(const CH_XDIR::OffsetBuffer& a_input);

//OffsetBuffer specialization of linearIn
template < >
void linearIn(CH_XDIR::OffsetBuffer& a_outputT, const void* const a_inBuf);

//OffsetBuffer specialization of linearOut
template < >
void linearOut(void* const a_outBuf, const CH_XDIR::OffsetBuffer& a_inputT);

template < > int linearSize(const Vector<CH_XDIR::OffsetBuffer>& a_input);
template < > void linearIn(Vector<CH_XDIR::OffsetBuffer>& a_outputT, const void* const inBuf);
template < > void linearOut(void* const a_outBuf, const Vector<CH_XDIR::OffsetBuffer>& a_inputT);

#include "BaseNamespaceFooter.H"
#include "NamespaceHeader.H"

// First, template specializations for read/write for FArrayBox.

template <>
inline void dataTypes(Vector<hid_t>& a_types, const BaseFab<int>& dummy)
{
  a_types.resize(1);
  a_types[0] = H5T_NATIVE_INT;
}

template <>
inline void dataTypes(Vector<hid_t>& a_types, const BaseFab<char>& dummy)
{
  a_types.resize(1);
  a_types[0] = H5T_NATIVE_CHAR;
}

/* since many compilers choke on the proper template specialization syntax
   I am forced to pass in a dummy specialization argument
*/
template <>
inline void dataTypes(Vector<hid_t>& a_types, const FArrayBox& dummy)
{
  a_types.resize(1);
  a_types[0] = H5T_NATIVE_REAL;
}

template <>
inline void dataTypes(Vector<hid_t>& a_types, const FluxBox& dummy)
{
  a_types.resize(1);
  a_types[0] = H5T_NATIVE_REAL;
}

template <>
inline void dataTypes(Vector<hid_t>& a_types, const EdgeDataBox& dummy)
{
  a_types.resize(1);
  a_types[0] = H5T_NATIVE_REAL;
}

template <>
inline void dataTypes(Vector<hid_t>& a_types, const NodeFArrayBox& dummy)
{
  a_types.resize(1);
  a_types[0] = H5T_NATIVE_REAL;
}


template <>
inline void dataSize(const BaseFab<int>& item, Vector<int>& a_sizes,
                     const Box& box, const Interval& comps)
{
  a_sizes[0] = box.numPts() * comps.size();
}

template <>
inline void dataSize(const BaseFab<char>& item, Vector<int>& a_sizes,
                     const Box& box, const Interval& comps)
{
  a_sizes[0] = box.numPts() * comps.size();
}

template <>
inline void dataSize(const FArrayBox& item, Vector<int>& a_sizes,
                     const Box& box, const Interval& comps)
{
  a_sizes[0] = box.numPts() * comps.size();
}

template <>
inline void dataSize(const FluxBox& item, Vector<int>& a_sizes,
                     const Box& box, const Interval& comps)
{
  int& t = a_sizes[0];
  t = 0;
  for (int dir =0; dir<CH_SPACEDIM; dir++)
    {
      Box edgeBox(surroundingNodes(box,dir));
      t += edgeBox.numPts()*comps.size();
    }
}


template <>
inline void dataSize(const EdgeDataBox& item, Vector<int>& a_sizes,
                     const Box& box, const Interval& comps)
{
  int& t = a_sizes[0];
  t = 0;
  for (int dir =0; dir<CH_SPACEDIM; dir++)
    {
      Box edgeBox(surroundingNodes(box));
      edgeBox.enclosedCells(dir);
      t += edgeBox.numPts()*comps.size();
    }
}


template <>
inline void dataSize(const NodeFArrayBox& item, Vector<int>& a_sizes,
                     const Box& box, const Interval& comps)
{
  Box boxNodes = surroundingNodes(box);
  a_sizes[0] = boxNodes.numPts() * comps.size();
}



template <>
inline const char* name(const FArrayBox& a_dummySpecializationArg)
{
  // Attempt to get rid of warnings on IBM...
  //static const char* name = "FArrayBox";
  const char* name = "FArrayBox";
  return name;
}

template <>
inline const char* name(const BaseFab<int>& a_dummySpecializationArg)
{
  // Attempt to get rid of warnings on IBM...
  //static const char* name = "BaseFab<int>";
  const char* name = "BaseFab<int>";
  return name;
}

template <>
inline const char* name(const BaseFab<char>& a_dummySpecializationArg)
{
  // Attempt to get rid of warnings on IBM...
  //static const char* name = "BaseFab<int>";
  const char* name = "BaseFab<char>";
  return name;
}
//
// now, generic, non-binary portable version of template functions
// for people who just want ot rely on linearIn/linearOut

template <class T>
inline void dataTypes(Vector<hid_t>& a_types, const T& dummy)
{
  a_types.resize(1);
  a_types[0] = H5T_NATIVE_CHAR;
}

template <class T>
inline void dataSize(const T& item, Vector<int>& a_sizes,
                     const Box& box, const Interval& comps)
{
  a_sizes[0] = item.size(box, comps);
}

template <class T>
inline void write(const T& item, Vector<void*>& a_allocatedBuffers,
                  const Box& box, const Interval& comps)
{
  item.linearOut(a_allocatedBuffers[0], box, comps);
}

template <class T>
inline void read(T& item, Vector<void*>& a_allocatedBuffers,
                 const Box& box, const Interval& comps)
{
  item.linearIn(a_allocatedBuffers[0], box, comps);
}

template <class T>
inline const char* name(const T& a_dummySpecializationArg)
{
  static const char* name = "unknown";
  return name;
}

template <class T>
void getOffsets(Vector<Vector<long long> >& offsets, const BoxLayoutData<T>& a_data,
                int types, const Interval& comps, const IntVect& outputGhost)
{
  CH_TIME("getOffsets");
  const BoxLayout& layout = a_data.boxLayout();
  {
    CH_TIMELEAF("offsets.resize");
    offsets.resize(types, Vector<long long>(layout.size()+1));
  //  offsets.resize(layout.size() + 1, Vector<long long>(types));
  }
  for (int t=0; t<types; t++) offsets[t][0] = 0;
  Vector<int> thisSize(types);
  if (T::preAllocatable()==0)
    { // static preAllocatable
      T dummy;
      unsigned int index = 1;
      for (LayoutIterator it(layout.layoutIterator()); it.ok(); ++it)
        {
          Box region = layout[it()];
          region.grow(outputGhost);
          {
            CH_TIMELEAF("dataSize");
            dataSize(dummy, thisSize, region, comps);
          }
          for (unsigned int i=0; i<thisSize.size(); ++i)
            {
              CH_TIMELEAF("offsets calc");
              //offsets[index][i] = offsets[index-1][i] + thisSize[i];
              offsets[i][index] = offsets[i][index-1] + thisSize[i];
            }
          ++index;
        }
    }
  else
    { // symmetric and dynamic preallocatable need two pass I/O
      OffsetBuffer buff;
      //int index = 0;
      for (DataIterator dit(a_data.dataIterator()); dit.ok(); ++dit)
        {
          int index = a_data.boxLayout().index(dit());
          //int index = dit().intCode();
          buff.index.push_back(index);
          Box region = layout[dit()];
          region.grow(outputGhost);
          { 
            CH_TIMELEAF("dataSize");
            dataSize(a_data[dit()], thisSize, region, comps);
          }
          buff.offsets.push_back(thisSize);
        }
      Vector<OffsetBuffer> gathering(numProc());
      {
        CH_TIMELEAF("gather");
        gather(gathering, buff, uniqueProc(SerialTask::compute));
      }
      {
        CH_TIMELEAF("broadcast");
        broadcast(gathering,  uniqueProc(SerialTask::compute));
      }
      //  pout() << gathering<<endl;
      for (int i=0; i<numProc(); ++i)
        {
          OffsetBuffer& offbuf = gathering[i];
          for (int num=0; num<offbuf.index.size(); num++)
            {
              int index = offbuf.index[num];
              for (unsigned int j=0; j<types; ++j)
                {
                  CH_TIMELEAF("offsets calc");
                  //offsets[index+1][j] = offbuf.offsets[num][j];
                  offsets[j][index+1] = offbuf.offsets[num][j];
                }
            }
        }
      for (int i=0; i<layout.size(); i++)
        {
          for (unsigned int j=0; j<types; ++j)
            {
              CH_TIMELEAF("offsets calc");
              //offsets[i+1][j] += offsets[i][j];
              offsets[j][i+1] += offsets[j][i];
            }
        }
    }

  // pout() << offsets<<endl;
}

//==================================================================
//
//  Statically pre-allocatable version only requires BoxLayout, not
//  the data
//
template <class T>
void getOffsets(Vector<Vector<long long> >& a_offsets,
                const BoxLayout             a_layout,
                int                         a_numTypes, 
                const Interval&             a_comps,
                const IntVect&              a_outputGhost)
{
  CH_TIME("getOffsets (prealloc)");
  CH_assert(T::preAllocatable() == 0);

  a_offsets.resize(a_numTypes, Vector<long long>(a_layout.size()+1));
  for (int t = 0; t != a_numTypes; ++t)
    {
      a_offsets[t][0] = 0;
    }
  Vector<int> thisSize(a_numTypes);
  T dummy;
  unsigned int index = 1;
  for (LayoutIterator it(a_layout.layoutIterator()); it.ok(); ++it)
    {
      Box region = a_layout[it()];
      region.grow(a_outputGhost);
      dataSize(dummy, thisSize, region, a_comps);
      for (int i = 0; i != thisSize.size(); ++i)  // Loop over (1) types
        {
          //offsets[index][i] = offsets[index-1][i] + thisSize[i];
          a_offsets[i][index] = a_offsets[i][index-1] + thisSize[i];
        }
      ++index;
    }
}

//==================================================================
//
// Now, linear IO routines for a BoxLayoutData of T
//
template <class T>
int write(HDF5Handle& a_handle, const BoxLayoutData<T>& a_data,
          const std::string& a_name, IntVect outputGhost,
          const Interval& in_comps, bool newForm)
{
  CH_TIME("write_Level");
  int ret = 0;

  Interval comps(in_comps);
  if ( comps.size() == 0) comps = a_data.interval();
  T dummy; // used for preallocatable methods for dumb compilers.
  Vector<hid_t> types;
  dataTypes(types, dummy);

  Vector<Vector<long long> > offsets;
  Vector<long long> bufferCapacity(types.size(), 1);  // noel (was 0)
  Vector<void*> buffers(types.size(), NULL);

  getOffsets(offsets, a_data, types.size(), comps, outputGhost);

  // create datasets collectively.
  hsize_t flatdims[1];
  char dataname[100];
  Vector<hid_t> dataspace(types.size());
  Vector<hid_t> dataset(types.size());

  herr_t err;
  hsize_t count[1];
  ch_offset_t offset[1];
  CH_assert(!(newForm && types.size() != 1));

  for (unsigned int i=0; i<types.size(); ++i)
    {
      flatdims[0] = offsets[i][offsets[i].size()-1];
      if (newForm)
      {
        if (a_name == "M")
        {
          strcpy(dataname, "Mask");
        }
        else
        {
          sprintf(dataname, "%sRegular", a_name.c_str());
        }
      }
      else
      {
        sprintf(dataname, "%s:datatype=%i",a_name.c_str(), i);
      }
      {
        CH_TIME("H5Screate");
        dataspace[i]      = H5Screate_simple(1, flatdims, NULL);
      }
      CH_assert(dataspace[i] >=0);
      {
        CH_TIME("H5Dcreate");
#ifdef H516
        dataset[i]        = H5Dcreate(a_handle.groupID(), dataname,
                                      types[i],
                                      dataspace[i], H5P_DEFAULT);
#else
        dataset[i]        = H5Dcreate2(a_handle.groupID(), dataname,
                                       types[i],
                                       dataspace[i], H5P_DEFAULT,
                                       H5P_DEFAULT, H5P_DEFAULT);
#endif
      }
      CH_assert(dataset[i] >= 0);
    }

  hid_t offsetspace, offsetData;
  for (unsigned int i=0; i<types.size(); ++i)
    {
      flatdims[0] =  offsets[i].size();
      if (newForm)
      {
        if (a_name == "M")
        {
          strcpy(dataname, "MaskOffsets");
        }
        else
        {
          sprintf(dataname, "%sOffsets",a_name.c_str());
        }
      }
      else
      {
        sprintf(dataname, "%s:offsets=%i",a_name.c_str(), i);
      }
      {
        CH_TIME("H5S_H5D_offsets_create");
        offsetspace =  H5Screate_simple(1, flatdims, NULL);
        CH_assert(offsetspace >= 0);
#ifdef H516
        offsetData  =   H5Dcreate(a_handle.groupID(), dataname,
                                  H5T_NATIVE_LLONG, offsetspace, 
                                  H5P_DEFAULT);
#else
        offsetData  =   H5Dcreate2(a_handle.groupID(), dataname,
                                   H5T_NATIVE_LLONG, offsetspace, 
                                   H5P_DEFAULT,
                                   H5P_DEFAULT,H5P_DEFAULT);
#endif
        CH_assert(offsetData >= 0);
      }
      if (procID() == 0)
        {
          CH_TIME("WriteOffsets");
          hid_t memdataspace = H5Screate_simple(1, flatdims, NULL);
          CH_assert(memdataspace >= 0);
          err = H5Dwrite(offsetData, H5T_NATIVE_LLONG, memdataspace, offsetspace,
                         H5P_DEFAULT, &(offsets[i][0]));
          CH_assert(err >= 0);
          H5Sclose(memdataspace);
        }
      {
        CH_TIME("H5S_H5D_offsets_close");
        H5Sclose(offsetspace);
        H5Dclose(offsetData);
      }
    }

  // write BoxLayoutData attributes into Dataset[0]
  if (!newForm)
  {
    CH_TIME("WriteAttributes");
    HDF5HeaderData info;
    info.m_int["comps"] = comps.size();
    info.m_string["objectType"] = name(dummy);
    info.m_intvect["outputGhost"] = outputGhost;
    info.m_intvect["ghost"] = outputGhost;
    std::string group = a_handle.getGroup();
    a_handle.setGroup(group+"/"+a_name+"_attributes");
    info.writeToFile(a_handle);
    a_handle.setGroup(group);
  }

  // collective operations finished, now perform parallel writes
  // to specified hyperslabs.

  Vector<size_t> type_size(types.size());
  for (unsigned int i=0; i<types.size(); ++i)
    {
      type_size[i] = H5Tget_size(types[i]);
    }

  Vector<int> thisSize(types.size());

  // Hooks for aggregated collective buffering (ACB).
  // Devendran, et al. Collective I/O Optimizations for Adaptive Mesh 
  // Refinement Data Writes on Lustre File System, CUG 2016 for more on
  // ACB.
  // Comment out next line if you don't want ACB.
  // Also, comment out the line "#define TRY_MPI_COLLECTIVES_ 0" below if
  // you don't want MPI-IO collective buffering. 
  // In general, you should turn on collective buffering if you use
  // ACB.
#ifdef CH_MPI
//#define AGGREGATE_BOXES_ 1
#endif
#ifdef AGGREGATE_BOXES_
  // pout() << "Turning on aggregated collective buffering (ACB)." << endl;
  // For each type types[i], aggCount[i] is the total number of data points for all 
  // the boxes on this processor.
  Vector<hsize_t> aggCount(types.size(), 0);
  // size of aggregated buffers; initialize to 0 so we can add up
  // the box sizes to compute the size of the buffers
  Vector<long long> aggBufferSize(types.size(), 0); 
  Vector<void*> aggBuffers(types.size(), NULL);
#endif

 // step 1, create buffer big enough to hold the biggest linearized T
 // that I will have to output.
 //  pout()<<"offsets ";
 // for (int i=0; i<offsets[0].size(); i++)
 // {
 //    pout()<<" "<<offsets[0][i];
 // }
  // pout()<<"\n";
  {
    CH_TIME("ComputeBufferCapacity");
    for (DataIterator it = a_data.dataIterator(); it.ok(); ++it)
      {
        unsigned int index = a_data.boxLayout().index(it());
        for (unsigned int i=0; i<types.size(); ++i)
          {
            long long size = (offsets[i][index+1] - offsets[i][index])
              * type_size[i];
            // pout()<<"index:size "<<index<<" "<<size<<"\n";
            CH_assert(size >= 0);
            if (size > bufferCapacity[i]) // grow buffer if necessary.....
              {
                bufferCapacity[i] = size;
              }
#ifdef AGGREGATE_BOXES_
            aggBufferSize[i] += size;
#endif
          }
      }
  }
  // CH_assert(bufferCapacity[0] > 1);
  for (unsigned int i=0; i<types.size(); ++i)
    {
      CH_TIME("mallocMT_buffers");
      buffers[i] = mallocMT(bufferCapacity[i]);
      if (buffers[i] == NULL)
      {
        pout() << " i=" << i
               << " types.size() = " << types.size()
               << " bufferCapacity[i] = " << (int)bufferCapacity[i]
               << endl;
        MayDay::Error("memory error in buffer allocation write");
      }
#ifdef AGGREGATE_BOXES_
      aggBuffers[i] = mallocMT(aggBufferSize[i]);
      if (aggBuffers[i] == NULL)
      {
        MayDay::Error("memory error in aggregate buffer allocation in write");
      }
#endif
    }

#ifdef CH_MPI
//#define TRY_MPI_COLLECTIVES_ 1
#endif
#ifdef TRY_MPI_COLLECTIVES_
  // pout() << "Turned on MPI-IO collective buffering." << endl;
  // MPI collective buffering requires all processes call H5Dwrite collectively.
  // In particular, the processes must issue the same number of H5Dwrite calls,
  // or the application will hang when collective is turned on.
  // In the case where we do a separate write for each box, gather the 
  // maximum number of boxes assigned to any process, so we know how many 
  // H5Dwrites to do.
  DataIterator dit = a_data.dataIterator();
  int maxNumBoxes = dit.size();
  int sendBuf = maxNumBoxes;
  int result = MPI_Allreduce(&sendBuf, &maxNumBoxes, 1, MPI_INT,MPI_MAX, Chombo_MPI::comm);
  if (result != MPI_SUCCESS)
    {
      MayDay::Error("Couldn't collect maximum number of boxes!");
    }
  // Set dataset transfer property to collective mode. This is how we turn
  // on MPI-IO collective in HDF5.
  hid_t DXPL = H5Pcreate(H5P_DATASET_XFER);
  if(!(a_handle.openMode()==HDF5Handle::CREATE_SERIAL)) // can't set MPI collective if file was created for serial IO
    H5Pset_dxpl_mpio(DXPL, H5FD_MPIO_COLLECTIVE);
#endif /*end TRY_MPI_COLLECTIVES_ */

  // Step 2.  actually a) write each of my T objects into the
  // buffer, then b) write that buffered data out to the
  // write position in the data file using hdf5 hyperslab functions.
  {
    CH_TIME("linearize_H5Dwrite");
#ifdef AGGREGATE_BOXES_
    // the first non-empty hyperslab needs to be handled separately from
    // the others, so keep track of the first non-empty hyperslab
    bool isFirstHyperslab = true;

    // bufferLoc keeps track of last location written to in the aggregrated 
    // buffer
    Vector<void*> bufferLoc(types.size(), NULL);
    // Set the bufferLoc to the beginning of the aggregated buffer
    for(unsigned int i=0; i<types.size(); ++i)
      {
        bufferLoc[i] = aggBuffers[i];
      }
#endif
    for (DataIterator it = a_data.dataIterator(); it.ok(); ++it)
      {
        const T& data = a_data[it()];
        unsigned int index = a_data.boxLayout().index(it());
        Box box = a_data.box(it());
        box.grow(outputGhost);
        // First, linearize the box, and put data into the buffer
        {       
          CH_TIMELEAF("linearize");
#ifdef AGGREGATE_BOXES_
          // Under aggregation, we need to be careful to write to
          // the buffer starting at bufferLoc.
          for(unsigned int i=0; i<types.size(); i++)
            {
              data.linearOut(bufferLoc[i], box, comps);
              char* tempLoc = ((char*)bufferLoc[i]) + data.size(box,comps);
              bufferLoc[i] = (void*)tempLoc;
              //              bufferLoc[i] += data.size(box, comps);
            }
#else
          write(data, buffers, box, comps); //write T to buffer
#endif
        }
        // Next select HDF5 hyperslabs to specify where to write in HDF5 file
        // BUG: The hyperslab union generation is wrong if types.size() > 1
        // because isFirstHyperslab isn't independent for different sizes.
        // We don't make a fix because types.size() is almost always 1 (and
        // will probably be set to 1 in future redesigns).
        for (unsigned int i=0; i<types.size(); ++i)
          {
            offset[0] = offsets[i][index];
            count[0] = offsets[i][index+1] - offset[0];
#ifdef AGGREGATE_BOXES_
            // Under aggregation, we may be sending multiple boxes to HDF5.
            // Because boxes on a process are not necessarily written to 
            // contiguous locations in the HDF5 file, we need to specify
            // multiple file locations to HDF5. This is done using a 
            // union of hyperslabs.

            aggCount[i] += count[0];
            if (isFirstHyperslab)
              {
                // If the box is non-empty, select the first hyperslab
                // and set isFirstHyperslab to false.  Otherwise,
                // we haven't encountered the first non-empty 
                // hyperslab, so keep isFirstHyperslab equal to true.
                if(count[0] > 0)
                  {
                    err =  H5Sselect_hyperslab(dataspace[i], H5S_SELECT_SET,
                                             offset, NULL,
                                             count, NULL);
                    CH_assert(err >= 0);
                    isFirstHyperslab = false;
                  }
                else // must explicitly tell HDF5 we are doing an empty write
                  {
                    H5Sselect_none(dataspace[i]);
                  }
              }
            else
              {
                if(count[0] > 0)
                  {
                    // H5S_SELECT_OR creates a union of the selected
                    // hyperslab with the already selected hyperslabs
                    // in dataspace
                    err =  H5Sselect_hyperslab(dataspace[i], H5S_SELECT_OR,
                                               offset, NULL,
                                               count, NULL);
                    CH_assert(err >= 0);
                  }
                // else don't do H5Sselect_none in case it overrides
                // the already existing union of hyperslabs.
              }
#else
            // Without aggregation, we create a simple hyperslab for the box.
            hid_t memdataspace=0;
            if (count[0] > 0)
              {
                err =  H5Sselect_hyperslab(dataspace[i], H5S_SELECT_SET,
                                         offset, NULL,
                                         count, NULL);
                CH_assert(err >= 0);
                memdataspace = H5Screate_simple(1, count, NULL);
                CH_assert(memdataspace >= 0);
              }
            else // must explicitly tell HDF5 we are doing an empty write
              {
                H5Sselect_none(dataspace[i]);
                H5Sselect_none(memdataspace);
              }

            // Write out box if we are NOT performing aggregation.
            // (Under aggregation, one single write call is issued at the
            // end of the function.)
            {
              CH_TIMELEAF("H5Dwrite");
#ifdef TRY_MPI_COLLECTIVES_
              err = H5Dwrite(dataset[i], types[i], memdataspace, dataspace[i],
                             DXPL, buffers[i]);
#else
              err = H5Dwrite(dataset[i], types[i], memdataspace, dataspace[i],
                             H5P_DEFAULT, buffers[i]);
#endif
            }
            CH_assert(err >= 0);
            H5Sclose(memdataspace);
            if (err < 0)
              {
                ret = err;
                pout() << "Before goto cleanup" << endl;
                goto cleanup;
              }
#endif // end of ifdef AGGREGATE_BOXES_ 
          } // end of loop over types
      } // end of loop over data iterator

    // Under aggregation, now we issue one write call to send all boxes to HDF5
    // If aggregation is turned off, but MPI collective is turned on, we may
    // have to do empty writes to make sure all processes call H5Dwrite 
    // the same number of times.
#ifdef AGGREGATE_BOXES_
    for(unsigned int i=0; i<types.size(); ++i)
      {
        if(aggCount[i] > 0)
          {
            hid_t memdataspace = 0;
            memdataspace = H5Screate_simple(1, &(aggCount[i]), NULL);
            CH_assert(memdataspace >= 0);
            {
              CH_TIMELEAF("H5Dwrite");
#ifdef TRY_MPI_COLLECTIVES_
              err = H5Dwrite(dataset[i], types[i], memdataspace, dataspace[i],
                           DXPL, aggBuffers[i]);
#else
              err = H5Dwrite(dataset[i], types[i], memdataspace, dataspace[i],
                                   H5P_DEFAULT, aggBuffers[i]);
#endif
            }
            H5Sclose(memdataspace);
            if (err < 0)
              {
                ret = err;
                pout() << "Error! goto cleanup" << endl;
                goto cleanup;
              }
          }
        //else aggCount[i] is 0, and this processor has no data to write.
        // For MPI collectives, still have to do an empty write call
#ifdef TRY_MPI_COLLECTIVES_
        else
          {
            hid_t memdataspace = 0;
            memdataspace = H5Screate_simple(1, &(aggCount[i]), NULL);
            H5Sselect_none(memdataspace);
            H5Sselect_none(dataspace[i]);
            err = H5Dwrite(dataset[i], types[i], memdataspace, dataspace[i],
                           DXPL, aggBuffers[i]);
            if (err < 0)
              {
                ret = err;
                pout() << "Before goto cleanup" << endl;
                goto cleanup;
              }
          }
#endif
      } // end for loop over types
#else // not using aggregated collective buffering
#ifdef TRY_MPI_COLLECTIVES_
    // MPI collectives expects all processes to make the same number of H5Dwrite calls,
    // or it will hang.  So, call H5Dwrite with empty data
    // First create memdataspace according to example in
    // https://www.hdfgroup.org/ftp/HDF5/examples/parallel/coll_test.c
    hid_t memdataspace = 0;
    // Setting first argument to 1 b/c that's the value used for non-empty writes.
    // (See H5Sselect_hyperslab code above.)
    memdataspace = H5Screate_simple(1, count, NULL);
    H5Sselect_none(memdataspace);
    int nBoxes = a_data.dataIterator().size();
    for(int iwrite = nBoxes; iwrite < maxNumBoxes; iwrite++)
      {
        for (unsigned int i=0; i<types.size(); ++i)
          {
            H5Sselect_none(dataspace[i]);
            // for debugging: buffers has junk data in it (but none of that data should
            // be written out here)
            err = H5Dwrite(dataset[i], types[i], memdataspace, dataspace[i],
                           DXPL, buffers[i]);
            if (err < 0)
              {
                ret = err;
                pout() << "Before goto cleanup" << endl;
                goto cleanup;
              }
          }
      }
    H5Sclose(memdataspace);

#endif // end of #ifdef TRY_MPI_COLLECTIVES_
#endif // end of #ifdef AGGREGATE_BOXES_

  } // end of region for CH_TIME("linearize_H5Dwrite")

#ifdef TRY_MPI_COLLECTIVES_
  H5Pclose(DXPL);
#endif
  
  // OK, clean up data structures
  
 cleanup:
  for (unsigned int i=0; i<types.size(); ++i)
    {
      {
        CH_TIME("freeMT");
        freeMT(buffers[i]);
        #ifdef AGGREGATE_BOXES_
        freeMT(aggBuffers[i]);
        #endif
      }
      {
        CH_TIME("H5Sclose");
        H5Sclose(dataspace[i]);
      }
      {
        CH_TIME("H5Dclose");
        H5Dclose(dataset[i]);
      }
    }
  return ret;
  
}

template <class T>
int write(HDF5Handle& a_handle, const LevelData<T>& a_data,
          const std::string& a_name, const IntVect& outputGhost, const Interval& in_comps)
{
  CH_TIMERS("Write Level");
  CH_TIMER("calc minimum in outputGhost",t1);
  CH_TIMER("writeToFile",t2);
  CH_TIMER("setGroup",t3);
  HDF5HeaderData info;
  info.m_intvect["ghost"] = a_data.ghostVect();
  IntVect og(outputGhost);
  CH_START(t1); 
  og.min(a_data.ghostVect());
  CH_STOP(t1); 
  info.m_intvect["outputGhost"] = og;
  std::string group = a_handle.getGroup();
  a_handle.setGroup(group+"/"+a_name+"_attributes");
  CH_START(t2); 
  info.writeToFile(a_handle);
  CH_STOP(t2); 
  CH_START(t3); 
  a_handle.setGroup(group);
  CH_STOP(t3); 
  return write(a_handle, (const BoxLayoutData<T>&)a_data, a_name, og, in_comps);
}

template <class T>
int read(HDF5Handle& a_handle, LevelData<T>& a_data, const std::string& a_name,
         const DisjointBoxLayout& a_layout, const Interval& a_comps, bool a_redefineData)
{
  if (a_redefineData)
    {
      HDF5HeaderData info;
      std::string group = a_handle.getGroup();
      if (a_handle.setGroup(group+"/"+a_name+"_attributes"))
        {
          std::string message = "error opening "
                                +a_handle.getGroup()+"/"+a_name+"_attributes" ;
          MayDay::Warning(message.c_str());
          return 1;
        }
      info.readFromFile(a_handle);
      a_handle.setGroup(group);
      int ncomp =  info.m_int["comps"];
      IntVect ghost = info.m_intvect["ghost"];
      if (a_comps.end() > 0 && ncomp < a_comps.end())
        {
          MayDay::Error("attempt to read component interval that is not available");
        }
      if (a_comps.size() == 0)
        a_data.define(a_layout, ncomp, ghost);
      else
        a_data.define(a_layout, a_comps.size(), ghost);
    }
  return read(a_handle, (BoxLayoutData<T>&)a_data, a_name, a_layout, a_comps, false);

}

template <class T>
int read(HDF5Handle& a_handle, BoxLayoutData<T>& a_data, const std::string& a_name,
         const BoxLayout& a_layout, const Interval& a_comps, bool a_redefineData)
{
  CH_TIME("read");
  int ret = 0; // return value;

  herr_t err;

  char dataname[100];
  hsize_t count[1];
  ch_offset_t offset[1];
  Vector<Vector<long long> > offsets;

  T dummy;
  Vector<hid_t> types;
  dataTypes(types, dummy);
  Vector<hid_t> dataspace(types.size());
  Vector<hid_t> dataset(types.size());
  offsets.resize(types.size(), Vector<long long>(a_layout.size() +1));

  //Vector<int> bufferCapacity(types.size(), 500);
  //Vector<void*> buffers(types.size(), NULL);
  Vector<Vector<char> > buffers(types.size(), 500);

  for (unsigned int i=0; i<types.size(); ++i)
    {
      sprintf(dataname, "%s:datatype=%i",a_name.c_str(), i);
#ifdef H516
      dataset[i]        = H5Dopen(a_handle.groupID(), dataname);
#else
      dataset[i]        = H5Dopen2(a_handle.groupID(), dataname, H5P_DEFAULT);
#endif
      if (dataset[i] < 0)
      {
        MayDay::Warning("dataset open failure"); return dataset[i];
      }
      dataspace[i]      = H5Dget_space(dataset[i]);
      if (dataspace[i] < 0)
      {
        MayDay::Warning("dataspace open failure"); return dataspace[i];
      }
    }

  hid_t offsetspace, offsetData;
  hsize_t flatdims[1];
  for (unsigned int i=0; i<types.size(); ++i)
    {
      flatdims[0] =  offsets[i].size();
      sprintf(dataname, "%s:offsets=%i",a_name.c_str(), i);
      offsetspace =  H5Screate_simple(1, flatdims, NULL);
      CH_assert(offsetspace >= 0);
#ifdef H516
      offsetData  =   H5Dopen(a_handle.groupID(), dataname);
#else
      offsetData  =   H5Dopen2(a_handle.groupID(), dataname,H5P_DEFAULT);
#endif
      CH_assert(offsetData >= 0);
      hid_t memdataspace = H5Screate_simple(1, flatdims, NULL);
      CH_assert(memdataspace >= 0);
      err = H5Dread(offsetData, H5T_NATIVE_LLONG, memdataspace, offsetspace,
                    H5P_DEFAULT, &(offsets[i][0]));
      CH_assert(err >=0);
      H5Sclose(memdataspace);
      H5Sclose(offsetspace);
      H5Dclose(offsetData);
    }

  HDF5HeaderData info;
  std::string group = a_handle.getGroup();
  if (a_handle.setGroup(a_handle.getGroup()+"/"+a_name+"_attributes"))
    {
      std::string message = "error opening "+a_handle.getGroup()+"/"+a_name ;
      MayDay::Warning(message.c_str());
      return 1;
    }

  info.readFromFile(a_handle);
  a_handle.setGroup(group);
  int ncomps = info.m_int["comps"];
  IntVect outputGhost(IntVect::Zero); // backwards file compatible mode.
  if (info.m_intvect.find("outputGhost") != info.m_intvect.end())
    {
      outputGhost = info.m_intvect["outputGhost"];
    }
  if (ncomps <= 0)
  {
    MayDay::Warning("ncomps <= 0 in read");
    return ncomps;
  }

  if (a_redefineData)
  {
    if (a_comps.size() != 0)
    {
      a_data.define(a_layout, a_comps.size());
    }
    else
    {
      a_data.define(a_layout, ncomps);
    }
  }

  Interval comps(0, ncomps-1);
  //huh?
  //  if (a_comps.size() != 0)  comps = Interval(0, a_comps.size());

  //  getOffsets(offsets, a_data, types.size(), comps);

  Vector<size_t> type_size(types.size());
  for (unsigned int i=0; i<types.size(); ++i)
    {
      type_size[i] = H5Tget_size(types[i]);

    }

  Vector<int> thisSize(types.size());
  DataIterator it = a_data.dataIterator();
  for ( ; it.ok(); ++it)
    {
      CH_TIMELEAF("H5Dread");
      T& data = a_data[it()];
      unsigned int index = a_data.boxLayout().index(it());
      Box box = a_data.box(it());

      for (unsigned int i=0; i<types.size(); ++i)
        {

          offset[0] = offsets[i][index];
          count[0] = offsets[i][index+1] - offset[0];
          if (count[0] > 0)
          {
            size_t size = count[0] * type_size[i];
            while (size > buffers[i].size())
              {
                buffers[i].resize(2*buffers[i].size());
            }
            
            err =  H5Sselect_hyperslab(dataspace[i], H5S_SELECT_SET,
                                       offset, NULL,
                                       count, NULL);
            CH_assert(err >= 0);
            hid_t memdataspace = H5Screate_simple(1, count, NULL);
            CH_assert(memdataspace >= 0);
            err = H5Dread(dataset[i], types[i], memdataspace, dataspace[i],
                          H5P_DEFAULT, &(buffers[i][0]));
            CH_assert(err >= 0);
            H5Sclose(memdataspace);
            if (err < 0)
              {
              ret = err;
              goto cleanup;
              }
          }
        }
      box.grow(outputGhost);
      read(data, buffers,  box, comps);
    }
//   if (it.size()==0)
//     {
//       //  try doing an empty H5Dread operation to make the collective system happier.
//       for (unsigned int i=0; i<types.size(); ++i)
//         {
//        hid_t filespace, memspace;
//        H5Sselect_none(filespace);
//        H5Sselect_none(memspace);
//        err = H5Dread(dataset[i], types[i], memspace, filespace,
//                      H5P_DEFAULT, 0);
//        H5Sclose(filespace);
//        H5Sclose(memspace);
//      }
//     }
                        
 cleanup:
  for (unsigned int i=0; i<types.size(); ++i)
    {
      // freeMT(buffers[i]);
      H5Sclose(dataspace[i]);
      H5Dclose(dataset[i]);
    }
  return ret;
}

template <class T>
int writeLevel(HDF5Handle& a_handle,
               const int&  a_level,
               const T& a_data,
               const Real& a_dx,
               const Real& a_dt,
               const Real& a_time,
               const Box&  a_domain,
               const int&  a_refRatio,
               const IntVect& outputGhost,
               const Interval& comps)
{
  int error;
  char levelName[10];
  std::string currentGroup = a_handle.getGroup();
  sprintf(levelName, "/level_%i",a_level);
  error = a_handle.setGroup(currentGroup + levelName);
  if (error != 0) return 1;

  HDF5HeaderData meta;
  meta.m_real["dx"] = a_dx;
  meta.m_real["dt"] = a_dt;
  meta.m_real["time"] = a_time;
  meta.m_box["prob_domain"] = a_domain;
  meta.m_int["ref_ratio"] = a_refRatio;

  error = meta.writeToFile(a_handle);
  if (error != 0) return 2;

  error = write(a_handle, a_data.boxLayout());
  if (error != 0) return 3;

  error = write(a_handle, a_data, "data", outputGhost, comps);
  if (error != 0) return 4;

  a_handle.setGroup(currentGroup);

  return 0;
}

template <class T>
int writeLevel(HDF5Handle& a_handle,
               const int&  a_level,
               const T& a_data,
               const RealVect& a_dx, // dx for each direction
               const Real& a_dt,
               const Real& a_time,
               const Box&  a_domain,
               const IntVect&  a_refRatios, // ref ratio for each direction
               const IntVect& outputGhost,
               const Interval& comps)
{
  int error;
  char levelName[10];
  std::string currentGroup = a_handle.getGroup();
  sprintf(levelName, "/level_%i",a_level);
  error = a_handle.setGroup(currentGroup + levelName);
  if (error != 0) return 1;

  HDF5HeaderData meta;
  meta.m_realvect["vec_dx"] = a_dx;
  meta.m_real["dt"] = a_dt;
  meta.m_real["time"] = a_time;
  meta.m_box["prob_domain"] = a_domain;
  meta.m_intvect["vec_ref_ratio"] = a_refRatios;

  error = meta.writeToFile(a_handle);
  if (error != 0) return 2;

  error = write(a_handle, a_data.boxLayout());
  if (error != 0) return 3;

  error = write(a_handle, a_data, "data", outputGhost, comps);
  if (error != 0) return 4;

  a_handle.setGroup(currentGroup);

  return 0;
}

template <class T>
int readLevel(HDF5Handle&   a_handle,
              const int&    a_level,
              LevelData<T>& a_data,
              Real& a_dx,
              Real& a_dt,
              Real& a_time,
              Box&  a_domain,
              int&  a_refRatio,
              const Interval&   a_comps,
              bool  setGhost)
{ CH_TIME("readLevel");
  HDF5HeaderData header;
  header.readFromFile(a_handle);
  //unused
  // int nComp = header.m_int["num_components"];

  int error;
  char levelName[10];
  std::string currentGroup = a_handle.getGroup();
  sprintf(levelName, "/level_%i",a_level);
  error = a_handle.setGroup(currentGroup + levelName);
  if (error != 0) return 1;

  HDF5HeaderData meta;
  error = meta.readFromFile(a_handle);
  if (error != 0) return 2;
  a_dx       = meta.m_real["dx"];
  a_dt       = meta.m_real["dt"];
  a_time     = meta.m_real["time"];
  a_domain   = meta.m_box["prob_domain"];
  a_refRatio = meta.m_int["ref_ratio"];
  Vector<Box> boxes;
  { CH_TIME("readLevel read boxes");
  error = read(a_handle, boxes);
  }
  Vector<int> procIDs;
  { CH_TIME("readLevel LoadBalance");
    LoadBalance(procIDs, boxes);
  }

  DisjointBoxLayout layout;
  { CH_TIME("readLevel define layout");
    layout.define(boxes, procIDs, a_domain);
  }

  { CH_TIME("readLevel close layout");
    layout.close();
  }
  if (error != 0) return 3;

  { CH_TIME("readLevel read general");
  error = read(a_handle, a_data, "data", layout, a_comps, true);
  }

  if (error != 0) return 4;

  a_handle.setGroup(currentGroup);

  return 0;
}

template <class T>
int readLevel(HDF5Handle&   a_handle,
              const int&    a_level,
              LevelData<T>& a_data,
              RealVect& a_dx,
              Real& a_dt,
              Real& a_time,
              Box&  a_domain,
              IntVect&  a_refRatio,
              const Interval&   a_comps,
              bool  setGhost)
{
  HDF5HeaderData header;
  header.readFromFile(a_handle);
  //unused
  // int nComp = header.m_int["num_components"];

  int error;
  char levelName[10];
  std::string currentGroup = a_handle.getGroup();
  sprintf(levelName, "/level_%i",a_level);
  error = a_handle.setGroup(currentGroup + levelName);
  if (error != 0) return 1;

  HDF5HeaderData meta;
  error = meta.readFromFile(a_handle);
  if (error != 0) return 2;
  // a_dx       = meta.m_realvect["vec_dx"];
  // Allow for vec_dx to be absent, as it will be in isotropic files.
  if (meta.m_realvect.find("vec_dx") != meta.m_realvect.end())
    {
      a_dx = meta.m_realvect["vec_dx"];
    }
  else
    { // vec_dx is not present, so get dx.
      Real dxScalar = meta.m_real["dx"];
      a_dx = dxScalar * RealVect::Unit;
    }
  a_dt       = meta.m_real["dt"];
  a_time     = meta.m_real["time"];
  a_domain   = meta.m_box["prob_domain"];
  // a_refRatio = meta.m_intvect["vec_ref_ratio"];
  // Allow for vec_ref_ratio to be absent, as it will be in isotropic files.
  if (meta.m_intvect.find("vec_ref_ratio") != meta.m_intvect.end())
    {
      a_refRatio = meta.m_intvect["vec_ref_ratio"];
    }
  else
    { // vec_ref_ratio is not present, so get ref_ratio.
      int refRatioScalar = meta.m_int["ref_ratio"];
      a_refRatio = refRatioScalar * IntVect::Unit;
    }
  Vector<Box> boxes;
  error = read(a_handle, boxes);
  Vector<int> procIDs;
  LoadBalance(procIDs, boxes);

  DisjointBoxLayout layout(boxes, procIDs, a_domain);

  layout.close();
  if (error != 0) return 3;

  error = read(a_handle, a_data, "data", layout, a_comps, true);

  if (error != 0) return 4;

  a_handle.setGroup(currentGroup);

  return 0;
}


/*******************************************************************************
 *
 * Class WriteMultiData: member definitions
 *
 ******************************************************************************/

/*--------------------------------------------------------------------*/
//  Constructor writes boxes and allocates dataset for LevelData
/** Write boxes, processors, attributes and creates a dataset for all
 *  the level data to be written.  The dataset is closed on
 *  destruction.
 *  \param[in]  a_handle
 *                      Chombo HDF5 handle holding open file for
 *                      writing
 *  \param[in]  a_layout
 *                      Box layout for all data to be written
 *  \param[in]  a_numIntv
 *                      Total number of components that will be
 *                      written from all BoxLayoutData or LevelData
 *  \param[in]  a_name  Name of the dataset for the level data.
 *  \param[in]  a_policyFlags
 *                      Either
 *                        CH_HDF5::IOPolicyDefault - Use internal
 *                          T.linearOut(...) to linearize data and
 *                          processes independently write to file,
 *                          but still in parallel.
 *                      or a union (|) of the following flags: 
 *                        CH_HDF5::IOPolicyMultiDimHyperslab - The
 *                          memory dataspace will be set up so that
 *                          HDF5 linearizes the data.  Using this
 *                          requires T.dataPtr() and contiguous memory
 *                          (so things like T=FluxBox will not work).
 *                        CH_HDF5::IOPolicyCollectiveWrite - The write
 *                          for each dataset will be collective.
 *  \param[in]  a_outputGhost
 *                      Number of ghost cells that will be written to
 *                      the data file.  Any data written must have
 *                      this many ghosts
 *  \param[in]  a_newForm
 *                      ?
 *//*-----------------------------------------------------------------*/

template <class T>
WriteMultiData<T>::WriteMultiData(HDF5Handle&      a_handle,
                                  const BoxLayout& a_layout,
                                  const int        a_numIntv,
                                  const string&    a_name,
                                  const int        a_policyFlags,
                                  const IntVect&   a_outputGhost,
                                  const bool       a_newForm)
  :
  m_policyFlags(a_policyFlags),
  m_outputGhost(a_outputGhost),
  m_newForm(a_newForm)
{
  CH_TIME("WriteMultiData::constructor");
  CH_assert(T::preAllocatable() == 0);
#ifdef H516
  CH_assert(false);  // API 1.6 not supported
#endif

//--Write the boxes

  write(a_handle, a_layout);

//--Create the dataset for the level data

  T dummy;  // Used for preallocatable methods for dumb compilers.
  dataTypes<T>(m_types, dummy);
  // Only allow one type
  CH_assert(m_types.size() == 1);

  m_allIntvFile.define(0, a_numIntv-1);
  getOffsets<T>(m_offsets,
                a_layout,
                m_types.size(),
                m_allIntvFile,
                m_outputGhost);

  hsize_t flatdims[1];
  // Create dataset for data
  {
    hid_t dataspace;
    flatdims[0] = m_offsets[0].back();
    if (m_newForm)
      {
        if (a_name == "M")
          {
            strcpy(m_dataname, "Mask");
          }
        else
          {
            sprintf(m_dataname, "%sRegular", a_name.c_str());
          }
      }
    else
      {
        sprintf(m_dataname, "%s:datatype=%i", a_name.c_str(), 0);
      }
    {
      CH_TIME("H5Screate");
      dataspace = H5Screate_simple(1, flatdims, NULL);
    }
    CH_assert(dataspace >=0);
    {
      CH_TIME("H5Dcreate");
      m_dataSet = H5Dcreate2(a_handle.groupID(), m_dataname,
                             m_types[0],
                             dataspace, H5P_DEFAULT,
                             H5P_DEFAULT, H5P_DEFAULT);
    }
    CH_assert(m_dataSet >= 0);
    {
      CH_TIME("H5S_H5D_data_close");
      H5Sclose(dataspace);
    }
  }

  // Create dataset for offsets
  {
    hid_t offsetspace, offsetdataset;
    char offsetname[128];
    flatdims[0] = m_offsets[0].size();  // Number of boxes
    if (m_newForm)
      {
        if (a_name == "M")
          {
            strcpy(offsetname, "MaskOffsets");
          }
        else
          {
            sprintf(offsetname, "%sOffsets", a_name.c_str());
          }
      }
    else
      {
        sprintf(offsetname, "%s:offsets=%i", a_name.c_str(), 0);
      }
    {
      CH_TIME("H5S_H5D_offsets_create");
      offsetspace = H5Screate_simple(1, flatdims, NULL);
      CH_assert(offsetspace >= 0);
      offsetdataset = H5Dcreate2(a_handle.groupID(), offsetname,
                                 H5T_NATIVE_LLONG, offsetspace, 
                                 H5P_DEFAULT,
                                 H5P_DEFAULT,H5P_DEFAULT);
      CH_assert(offsetdataset >= 0);
    }
    if (procID() == 0)
      // Write the offsets
      {
        herr_t err;
        CH_TIME("WriteOffsets");
        hid_t memdataspace = H5Screate_simple(1, flatdims, NULL);
        CH_assert(memdataspace >= 0);
        err = H5Dwrite(offsetdataset, H5T_NATIVE_LLONG, memdataspace,
                       offsetspace, H5P_DEFAULT, &(m_offsets[0][0]));
        CH_assert(err >= 0);
        H5Sclose(memdataspace);
      }
    {
      CH_TIME("H5S_H5D_offsets_close");
      H5Sclose(offsetspace);
      H5Dclose(offsetdataset);
    }
  }

  // Write attributes
  {
    CH_TIME("WriteAttributes");
    HDF5HeaderData info;
    info.m_intvect["ghost"] = m_outputGhost;
    info.m_intvect["outputGhost"] = m_outputGhost;
    if (!m_newForm)
      {
        info.m_int["comps"] = m_allIntvFile.size();
        info.m_string["objectType"] = name(dummy);
      }
    std::string group = a_handle.getGroup();
    a_handle.setGroup(group+"/"+a_name+"_attributes");
    info.writeToFile(a_handle);
    a_handle.setGroup(group);
  }

  // Get the maximum boxes per process
  {
    m_maxBoxPerProc = a_layout.dataIterator().size();
#ifdef CH_MPI
    long myCountBoxes = m_maxBoxPerProc;
    MPI_Allreduce(&myCountBoxes,
                  &m_maxBoxPerProc,
                  1, MPI_LONG, MPI_MAX, Chombo_MPI::comm);
#endif    
  }
}

/*--------------------------------------------------------------------*/
//  Write an interval of LevelData to the dataset
/** Call this as many times as you want to write level data.  User is
 *  resposible for making sure 'a_intvFile' has relevance to
 *  'a_numIntv' specified during construction.
 *  \param[in]  a_data  Level data to write
 *  \param[in]  a_intvMem
 *                      Interval of a_data to read
 *  \param[in]  a_intvFile
 *                      Interval to write to file
 *//*-----------------------------------------------------------------*/

template <class T>
int
WriteMultiData<T>::writeData(const BoxLayoutData<T>& a_data,
                             const Interval&         a_intvMem,
                             const Interval&         a_intvFile)
{
  CH_TIME("WriteMultiData::writeData");
  CH_assert(a_intvMem.size() > 0);
  CH_assert(a_intvMem.size() == a_intvFile.size());
  CH_assert(a_intvMem.begin()  >= a_data.interval().begin() &&
            a_intvMem.end()    <= a_data.interval().end());
  CH_assert(a_intvFile.begin() >= m_allIntvFile.begin() &&
            a_intvFile.end()   <= m_allIntvFile.end());

  int ret = 0;

  hsize_t fileCount[1], memCount[SpaceDim+1];
  ch_offset_t fileOffset[1], memOffset[SpaceDim+1];  // Same type as hsize_t
  herr_t err;

//--Open the file dataspace

  hid_t fileDataSpace = H5Dget_space(m_dataSet);

//--Data transfer policies (independent or collective)

  hid_t DXPL = H5Pcreate(H5P_DATASET_XFER);
#ifdef CH_MPI
  if (m_policyFlags & CH_HDF5::IOPolicyCollectiveWrite)
    {
      H5Pset_dxpl_mpio(DXPL, H5FD_MPIO_COLLECTIVE);
    }
#endif

  if (m_policyFlags & CH_HDF5::IOPolicyMultiDimHyperslab)

//--Perform parallel writes from specified hyperslabs (when defining hyperslabs
//--in memory, remember that HDF5 assumes row-ordering).

    {
      CH_TIME("hyperslab_H5Dwrite");
      DataIterator dit = a_data.dataIterator();
      const void* buffer;
      for (int iBox = 0; iBox != m_maxBoxPerProc; ++iBox)
        {
          hid_t memDataSpace;
          if (dit.ok())
            {
              const T& data = a_data[dit];
              unsigned globalIdx = a_data.boxLayout().index(dit());
              Box dbox = data.box();         // Data box
              Box rbox = a_data.box(dit());  // Read box
              rbox.grow(m_outputGhost);
              // Create the dataspace in memory
              memCount[0] = a_data.nComp();
              D_TERM6(memCount[SpaceDim-0] = dbox.size(0);,
                      memCount[SpaceDim-1] = dbox.size(1);,
                      memCount[SpaceDim-2] = dbox.size(2);,
                      memCount[SpaceDim-3] = dbox.size(3);,
                      memCount[SpaceDim-4] = dbox.size(4);,
                      memCount[SpaceDim-5] = dbox.size(5);)
              memDataSpace = H5Screate_simple(SpaceDim+1, memCount, NULL);
              // Select the hyperslab from the memory dataspace
              memCount[0] = a_intvMem.size();
              D_TERM6(memCount[SpaceDim-0] = rbox.size(0);,
                      memCount[SpaceDim-1] = rbox.size(1);,
                      memCount[SpaceDim-2] = rbox.size(2);,
                      memCount[SpaceDim-3] = rbox.size(3);,
                      memCount[SpaceDim-4] = rbox.size(4);,
                      memCount[SpaceDim-5] = rbox.size(5);)
              memOffset[0] = a_intvMem.begin();
              D_TERM6(
                memOffset[SpaceDim-0] = rbox.smallEnd(0) - dbox.smallEnd(0);,
                memOffset[SpaceDim-1] = rbox.smallEnd(1) - dbox.smallEnd(1);,
                memOffset[SpaceDim-2] = rbox.smallEnd(2) - dbox.smallEnd(2);,
                memOffset[SpaceDim-3] = rbox.smallEnd(3) - dbox.smallEnd(3);,
                memOffset[SpaceDim-4] = rbox.smallEnd(4) - dbox.smallEnd(4);,
                memOffset[SpaceDim-5] = rbox.smallEnd(5) - dbox.smallEnd(5);)
              err = H5Sselect_hyperslab(memDataSpace, H5S_SELECT_SET,
                                        memOffset, NULL, memCount, NULL);
              CH_assert(err >= 0);
              // Create the hyperslab in the file dataspace
              fileOffset[0] = m_offsets[0][globalIdx];
              fileCount[0] = m_offsets[0][globalIdx + 1] - fileOffset[0];
              if (fileCount[0] > 0)  // Else catches more processes than boxes
                {
                  // Revise offsets based on selection of interval
                  const hsize_t dpnts = rbox.numPts();  // Points per comp
                  fileOffset[0] += dpnts*a_intvFile.begin();
                  fileCount[0] = dpnts*a_intvFile.size();
                  CH_assert(fileOffset[0] + fileCount[0] <=
                            m_offsets[0][globalIdx+1]);
                  err = H5Sselect_hyperslab(fileDataSpace, H5S_SELECT_SET,
                                            fileOffset, NULL, fileCount, NULL);
                  CH_assert(err >= 0);
                }
              else  // More processes than boxes
                {
                  H5Sselect_none(memDataSpace);
                  H5Sselect_none(fileDataSpace);
                }
              buffer = data.dataPtr();
              ++dit;
            }
          else  // Does not have a box to participate in collective :(
            {
              std::memset(memCount, 0, (SpaceDim+1)*sizeof(hsize_t));
              memDataSpace = H5Screate_simple(SpaceDim+1, memCount, NULL);
              H5Sselect_none(memDataSpace);
              H5Sselect_none(fileDataSpace);
              buffer = 0;
            }
          // Write
          err = H5Dwrite(m_dataSet, m_types[0], memDataSpace, fileDataSpace,
                         DXPL, buffer);
          CH_assert(err >= 0);
          H5Sclose(memDataSpace);
          if (err < 0)
            {
              ret = err;
              goto cleanup;
            }
        }
    }
  else

//--Perform parallel writes with 1-D hyperslabs and using internal mechanisms
//--for linearization.

    {
      CH_TIME("linearize_H5Dwrite");
      void* linearBuffer;

      // Step 1: create a buffer to linearize T
      {
        CH_TIMELEAF("allocateLinearBuffer");
        long long bufferSize = 0;
        for (DataIterator dit = a_data.dataIterator(); dit.ok(); ++dit)
          {
            unsigned globalIdx = a_data.boxLayout().index(dit());
            {
              const long long allIntvSize =
                (m_offsets[0][globalIdx + 1] - m_offsets[0][globalIdx]);
              CH_assert(allIntvSize >= 0);
              if (allIntvSize > bufferSize)
                {
                  bufferSize = allIntvSize;
                }
            }
          }
        // Get buffer size in bytes, and adjusted for write interval
        bufferSize = ((bufferSize/m_allIntvFile.size())*a_intvMem.size())*
          H5Tget_size(m_types[0]);
        linearBuffer = mallocMT(bufferSize);
        if (linearBuffer == NULL)
          {
            pout() << " bufferCapacity = " << (int)bufferSize << endl;
            MayDay::Error("Memory error in buffer allocation write");
          }
      }

      // Step 2: linearize the data and then write to file
      DataIterator dit = a_data.dataIterator();
      for (int iBox = 0; iBox != m_maxBoxPerProc; ++iBox)
        {
          hid_t memDataSpace;
          if (dit.ok())
            {
              const T& data = a_data[dit];
              unsigned globalIdx = a_data.boxLayout().index(dit());
              Box rbox = a_data.box(dit());  // Read box
              rbox.grow(m_outputGhost);
              {       
                CH_TIMELEAF("linearize");
                data.linearOut(linearBuffer, rbox, a_intvMem);
              }
              const hsize_t dpnts = rbox.numPts();  // Points per comp
              // Create the dataspace in memory
              memCount[0] = dpnts*a_intvMem.size();
              memDataSpace = H5Screate_simple(1, memCount, NULL);
              // Create the hyperslab in the file dataspace
              fileOffset[0] = m_offsets[0][globalIdx];
              fileCount[0] = m_offsets[0][globalIdx + 1] - fileOffset[0];
              if (fileCount[0] > 0)  // Else catches more processes than boxes
                {
                  // Revise offsets based on selection of interval
                  fileOffset[0] += dpnts*a_intvFile.begin();
                  fileCount[0] = dpnts*a_intvFile.size();
                  CH_assert(fileOffset[0] + fileCount[0] <=
                            m_offsets[0][globalIdx+1]);
                  err = H5Sselect_hyperslab(fileDataSpace, H5S_SELECT_SET,
                                            fileOffset, NULL, fileCount, NULL);
                  CH_assert(err >= 0);
                }
              else  // More processes than boxes
                {
                  H5Sselect_none(memDataSpace);
                  H5Sselect_none(fileDataSpace);
                }
              ++dit;
            }
          else  // Does not have a box to participate in collective :(
            {
              memCount[0] = 0;
              memDataSpace = H5Screate_simple(SpaceDim+1, memCount, NULL);
              H5Sselect_none(memDataSpace);
              H5Sselect_none(fileDataSpace);
            }
          // Collective write
          err = H5Dwrite(m_dataSet, m_types[0], memDataSpace, fileDataSpace,
                         DXPL, linearBuffer);
          CH_assert(err >= 0);
          H5Sclose(memDataSpace);
          if (err < 0)
            {
              ret = err;
              freeMT(linearBuffer);
              goto cleanup;
            }
        }
      freeMT(linearBuffer);
    }

  cleanup: ;
  H5Pclose(DXPL);
  H5Sclose(fileDataSpace);
  // m_dataSet closed in destructor
  return ret;
}

#include "NamespaceFooter.H"

#else // CH_USE_HDF5

// this is the only thing needed when HDF is not used
#define HOFFSET(S,M)    (offsetof(S,M))

#endif // CH_USE_HDF5 not defined
#endif  // CH_HDF5_H