/**
 * \file IO.C
 * \brief IO class definition.
 */

#include "IO.h"
#include <iostream>
#include <fstream>

// Get revision number store in HDF file:
#include "revision.h"

using namespace std;

IO::IO(Params *params)
{
  PARAMS = params;
}

/********************************************************************//**
 *  \author Pete Schmitt (schmitt at ucar.edu)
 *  \since October 2008
 *
 * Open HDF file, add HDF variables, close (save) hdf file
 * 
 * \param[out] tilt - tilt angle
 * \param[in] UT - Date & Time of file to open
 * \param[in] mhdTimeStep - timestep of file to open
 * \param[out] mhdTimeStep - timestep of corresponding UT file, if UT.getYear() == -9999
 * \param[out] varFunction - MIX variables to read
 *
 * \see IO::Open(const intArray &UT, const long &mhdTimeStep)
 * \see IO::ReadAttributes(real &tilt, const intArray &UT, long &mhdTimeStep);
 * \see IO::GetData(realCompositeGridFunction &varFunction);
 * \see IO::Close()
 ************************************************************************/
void IO::Read(real &tilt, const DateTime &UT, long &mhdTimeStep, realCompositeGridFunction &varFunction)
{
  intArray utime(6);
  utime(0) = UT.getYear();
  utime(1) = UT.getMonth();
  utime(2) = UT.getDay();
  utime(3) = UT.getHours();
  utime(4) = UT.getMinutes();
  utime(5) = (int) UT.getSeconds();

  Open(utime, mhdTimeStep, false);
  ReadAttributes(tilt, utime, mhdTimeStep);
  GetData(varFunction);
  Close();
}

/********************************************************************//**
 *  \author Pete Schmitt (schmitt at ucar.edu)
 *  \since October 2008
 *
 * Open HDF file, add HDF variables, close (save) hdf file
 * 
 * \param[out] tilt - tilt angle
 * \param[in] UT - Date & Time of file to open
 * \param[in] mhdTimeStep - timestep of file to open
 * \param[out] mhdTimeStep - timestep of corresponding UT file, if UT.getYear() == -9999
 * \param[out] varFunction - MIX variables to read
 *
 * \see IO::Open(const intArray &UT, const long &mhdTimeStep)
 * \see IO::ReadAttributes(real &tilt, const intArray &UT, long &mhdTimeStep);
 * \see IO::GetData(realCompositeGridFunction &varFunction);
 * \see IO::Close()
 ************************************************************************/
void IO::ReadVariable(const DateTime &UT, realCompositeGridFunction &varFunction, int VARIABLE_ID)
{
  intArray utime(6);
  utime(0) = UT.getYear();
  utime(1) = UT.getMonth();
  utime(2) = UT.getDay();
  utime(3) = UT.getHours();
  utime(4) = UT.getMinutes();
  utime(5) = (int) UT.getSeconds();
  
  long mhdTimeStep;

  Open(utime, mhdTimeStep, false);
  assert( ReadOneVariable(varFunction, VARIABLE_ID) );
  Close();
}

/********************************************************************//**
 *  \author Pete Schmitt (schmitt at ucar.edu)
 *  \since October 2008
 *
 * Open HDF file, add HDF variables, close (save) hdf file
 * 
 * \see IO::Open(const intArray &UT, const long &mhdTimeStep)
 * \see IO::SaveAttributes(const real &tilt, const intArray &UT, const long &mhdTimeStep
 * \see IO::Add(const realCompositeGridFunction &, const char*);
 * \see IO::Close()
 ************************************************************************/
void IO::Write(const real &tilt, const DateTime &UT, const long &mhdTimeStep, const realCompositeGridFunction &varFunction)
{  
  intArray utime(6);
  utime(0) = UT.getYear();
  utime(1) = UT.getMonth();
  utime(2) = UT.getDay();
  utime(3) = UT.getHours();
  utime(4) = UT.getMinutes();
  utime(5) = (int) UT.getSeconds();

  Write(tilt, utime, mhdTimeStep, varFunction);
}
void IO::Write(const real &tilt, const intArray &UT, const long &mhdTimeStep, const realCompositeGridFunction &varFunction)
{
  Open(UT, mhdTimeStep, true);
  SaveAttributes(tilt, UT, mhdTimeStep);
  Add(varFunction);
  Close();
}

/************************************************************************
 * Private member functions
 ***********************************************************************/

/********************************************************************//**
 * \author Peter Schmitt (schmitt (at) ucar.edu)
 * \since 10-2008
 *
 * \param[out] tilt - tilt angle read in from file

 * \param[in] UT - anticipated UT (YYYY-MM-DD HH-MM-SS) in the HDF file.
 *                 If year=-9999, then it's a "TimeStep" run and not UTIO.
 * \param[in] mhdTimeStep - anticipated timestep for file
 * \param[out] mhdTimeStep - timestep corresponding to UT (if UT.getYear() == -9999)
 * 
 * Read tilt angle, UT and time step from HDF file.
 *
 * \note HDF file must already be opened for read using IO::Open(...)
 *
 * \see IO::SaveAttributes
 ************************************************************************/
void IO::ReadAttributes(real &tilt, const intArray &UT, long &mhdTimeStep)
{
  // Read some attributes
  // 1. stuff received from the MHD code
  float fltTilt;
  //int status = SDsetattr(hdfFileID,"Tilt",DFNT_FLOAT32,1,(VOIDP)&fltTilt);
  int status = SDfindattr(hdfFileID,"Tilt");
  status = SDreadattr(hdfFileID,status,(VOIDP)&fltTilt);
  long timestep;
  //status = SDsetattr(hdfFileID,"MHD time step",DFNT_INT32,1,(VOIDP)&timestep); 
  status = SDfindattr(hdfFileID,"MHD time step");
  status = SDreadattr(hdfFileID,status,(VOIDP)&timestep); 
  intArray utime(6);
  //status = SDsetattr(hdfFileID,"UT (Year, Month, Day, Hour, Min, Sec)",DFNT_INT32,6,
  //		     (VOIDP)(utime.getDataPointer()));
  status = SDfindattr(hdfFileID,"UT (Year, Month, Day, Hour, Min, Sec)");
  status = SDreadattr(hdfFileID,status,(VOIDP)(utime.getDataPointer()));

  // See "SaveAttributes" for more details on the HDF attributes in the file...  

  // Make sure either utime or timestep match:
  if ( ( 
	(utime(0) != UT(0)) && (utime(1) != UT(1)) && (utime(2) != UT(2)) &&
	(utime(3) != UT(3)) && (utime(4) != UT(4)) && (utime(5) != UT(5)) 
	 )
       &&
       (timestep != mhdTimeStep) 
       ){
    cerr << "Error reading MIX HDF file:  "
	 << " Neither the UT (YYYY-MM-DD HH-MM-SS) nor the integer time step match expected values!" << endl;
    // kill the program.
    assert("Neither UT nor timestep match expected values" == false);
  }   
}


/********************************************************************//**
 * \author Peter Schmitt (schmitt (at) ucar.edu)
 * \since 10-2008
 *
 * \param[in] varFunction - Overture realCompositeGridFunction that
 *            already has its grid setup.  Note: it's important that
 *            the grid is properly defined before calling this function!
 * 
 * Read all non-grid variables into a realCompositeGridFUnction from an HDF file.
 *
 * \note Since the realCompositeGridFunction is already setup, we dont read any grid data from file!
 *
 * \note HDF file must already be opened for read using IO::Open(...)
 ************************************************************************/
void IO::GetData(realCompositeGridFunction &varFunction)
{
  // dont bother reading the grid information...
  // varFunction should already have that before entering this function!

  MappedGrid *grid = varFunction[0].getMappedGrid();
  Index I1,I2,I3;
  // include the periodic points here (gridIndexRange) to make it compliant with the grid storage
  getIndex((*grid).gridIndexRange(),I1,I2,I3);

  // Read the data:
  realArray varArray;
#if 0
  // "GRID X"
  varArray = (*grid).vertex()(I1,I2,I3,0);
  ReadOneVariable("Grid X",varArray);
  // "GRID Y"
  varArray = (*grid).vertex()(I1,I2,I3,1);
  ReadOneVariable("Grid Y",varArray);

#endif

  assert( ReadOneVariable(varFunction, SOLVER_POTENTIAL_NORTH) );
  assert( ReadOneVariable(varFunction, SOLVER_POTENTIAL_SOUTH) );

  assert( ReadOneVariable(varFunction, SOLVER_FAC_NORTH) );
  assert( ReadOneVariable(varFunction, SOLVER_FAC_SOUTH) );

  assert( ReadOneVariable(varFunction, SOLVER_SIGMAP_NORTH) );
  assert( ReadOneVariable(varFunction, SOLVER_SIGMAP_SOUTH) );

  assert( ReadOneVariable(varFunction, SOLVER_SIGMAH_NORTH) );
  assert( ReadOneVariable(varFunction, SOLVER_SIGMAH_SOUTH) );

  assert( ReadOneVariable(varFunction, SOLVER_AVG_ENG_NORTH) );
  assert( ReadOneVariable(varFunction, SOLVER_AVG_ENG_SOUTH) );

  assert( ReadOneVariable(varFunction, SOLVER_NUM_FLUX_NORTH) );
  assert( ReadOneVariable(varFunction, SOLVER_NUM_FLUX_SOUTH) );

  assert( ReadOneVariable(varFunction, SOLVER_WIND_NORTH) );
  assert( ReadOneVariable(varFunction, SOLVER_WIND_SOUTH) );

  // This is how you'd print at an individual variable:
  //varFunction[0](I1,I2,I3,varID).display("SOLVER_FAC_SOUTH");

  // ... Or how you'd print at the whole grid:
  //varFunction[0].display("real cgf data"); 
}

/********************************************************************//**
 * \author Peter Schmitt (schmitt (at) ucar.edu)
 * \since 10-2008
 *
 * \param[in] cgf - Overture realCompositeGridFunction that already has its grid setup
 * \param[out] cgf - reads data into cgf[0](*,*,*, VARIABLE_ID)
 * \param[in] VARIABLE_ID - integer variable ID that exists within the CGF.
 * \returns bool - returns true if successful
 *
 * Read a variable into a realCompositeGridFUnction from the HDF file.
 *
 * \note HDF file must already be opened for read using IO::Open(...)
 ************************************************************************/
bool IO::ReadOneVariable(realCompositeGridFunction &cgf, int VARIABLE_ID){
  bool returnValue = true;

  MappedGrid *grid = cgf[0].getMappedGrid();
  Index I1,I2,I3;
  // include the periodic points here (gridIndexRange) to make it compliant with the grid storage
  getIndex((*grid).gridIndexRange(),I1,I2,I3);

  doubleArray data_realArray = cgf[0](I1,I2,I3, VARIABLE_ID);
  returnValue = ReadOneVariable(cgf.getName(VARIABLE_ID).c_str(), data_realArray);
  cgf[0](I1,I2,I3,VARIABLE_ID) = data_realArray;

  int nLat = (*grid).vertex()(I1,I2,I3,0).getBound(0) - (*grid).vertex()(I1,I2,I3,0).getBase(0) + 1;
  int nLon = (*grid).vertex()(I1,I2,I3,0).getBound(1) - (*grid).vertex()(I1,I2,I3,0).getBase(1) + 1;

#if 0
  // debugging output:
  cout << "Verify the var.........." << endl;
  for (int i = 0; i < nLon; i++){
    cout << "i=" << i << "\t";
    for (int j=0; j < nLat; j++){
      cout << cgf[0](I1,I2,I3,VARIABLE_ID)(j,i) << " ";
    }
    cout << endl;
  }
#endif
  return returnValue;
}


/********************************************************************//**
 * \author Peter Schmitt (schmitt (at) ucar.edu)
 * \since 10-2008
 *
 * \param[in] variable - String of the variable to open (e.g. "Potential North [V]")
 * \param[out] data - array to save data to
 * \returns bool - succesful if returns true.
 *
 * Read #{variable} from the HDF file.  
 *
 * \note  Data I/O is done with FLOATs.  This function will convert from a floatArray to a doubleArray.
 *
 * \note HDF file must already be opened for read using IO::Open(...)
 ************************************************************************/
bool IO::ReadOneVariable(const char *variable, realArray &data_realArray)
{
#ifndef OV_USE_DOUBLE
#error The MIX file I/O requires double-precision Overture, however Overture does not appear to be compiled with double-precision.
#endif
  //float *data = (data_realArray.convertTo_floatArray()).getDataPointer();

  int32 errorId;                       // HDF error ID
  int32 SDselectId;                    // HDF selection ID
  int32 varId;                         // HDF variable ID
  int32 dataType;                      // Internal HDF datatype of variable
  int32 nAttributes;                   // number of attributes for the data set
  int32 rank;                          // rank of the variable.
  int32 dimensions_i32[MAX_VAR_DIMS];  // used to compute indexEnd
  int32 indexStart[MAX_VAR_DIMS];      // [x,y,z] starting location from where data is read
  int32 indexEnd[MAX_VAR_DIMS];        // [x,y,z] number of values to be read along each dimension
  char dataName[MAX_NC_NAME];          // 


  int i;
  
  for(i=0;i<MAX_VAR_DIMS;i++) {
    dimensions_i32[i] = 0;
    indexStart[i] = 0;
    indexEnd[i] = 0;
  }
  /*Get the info about thie dataset */
  varId = SDnametoindex(hdfFileID,variable);
  if (varId < 0) {
    cerr << "*** LFMpara::getData Unable to open " << variable << "\n";
    return false;
  }
  else {
    SDselectId = SDselect(hdfFileID,varId);
    errorId = SDgetinfo(SDselectId,dataName,&rank,dimensions_i32,&dataType,&nAttributes);
    if (errorId < 0){
      cerr << "*** LFMpara::getData unable to read info for variable " << variable << "\n";
      return false;
    }
    /* Create reterival information arrays and allocate space */
    for(i=0;i<rank;i++) {
      indexStart[i]=0;
      indexEnd[i]=dimensions_i32[rank-i-1];
      indexEnd[i]=dimensions_i32[i];
    }
  }

  // Allocate memory that we'll read in from the file.
  int nLon = indexEnd[0]-indexStart[0];
  int nLat = indexEnd[1]-indexStart[1];
  float *data = new float[(nLon+1)*(nLat+1)];

  if (varId > 0){

    /* Get dataset */
    errorId = SDreaddata(SDselectId,indexStart,NULL,indexEnd,data);
    if (errorId < 0){
      cerr << "*** LFMpara::getData unable to read data for variable " << variable << "\n";
      return false;
    }

    errorId = SDendaccess(SDselectId);
    if (errorId < 0){
      cerr << "*** LFMpara::getData unable to end access for variable " << variable << "\n";
      return false;
    }
  }

  // Copy data from temporary memory into realArray object

  for (int i=0; i < nLon; i++){
    for (int j=0; j < nLat; j++){
      data_realArray(j,i) = data[j+i*nLat];
    }
  }

#if 0
  // Debuging output:
  for (int i=0; i < nLon; i++){
    cout << "i=" << i << "\t";
    for (int j=0; j < nLat; j++){
      data_realArray(j,i) = data[j+i*nLat];
      cout << data_realArray(j,i) << " ";
      //cout << data[j+i*nLat] << " ";
    }
    cout << endl;
  }
#endif  
  // Free the data pointer
  delete data;

  // if we made it this far, success!
  return true;
}


/********************************************************************//**
 *  \author Pete Schmitt (schmitt at ucar.edu)
 *  \since March 2008
 *
 * Open MIX HDF file handle for I/O.
 *
 * if (UT(0) == -9999), assume we're using timestep I/O and name files:
 *       [filename]_mix_00000000.hdf
 * else, use UTIO and name files:
 *       [filename]_mix_YYYY-MM-DDTHH-MM-SSZ.hdf
 *
 * \param[in] UT - UT integer array (YYYY, MM, DD, HH, MM, SS)
 * \param[in] mhdTimeStep - integer timestep of data dump
 * \param[in] openForWrite - Open HDF file for write (true) or read (false)
 ************************************************************************/
void IO::Open(const intArray &UT, const long &mhdTimeStep, const bool &openForWrite)
{
  // Construct output file name
  // Are we using Standard timestep IO?
  if ( UT(0) <= 0 ){
    int nsym=sprintf(hdfName,"%s%s%07d%s",PARAMS->outputFilePrefix,"_mix_",mhdTimeStep,".hdf");
  }
  else{
    // We are using UTIO, name the MIX file according to UT timestamp.
    // i.e. of the form YYYY-MM-DDTHH-MM-SSZ
    int nsym=sprintf(hdfName,"%s%s%04d-%02d-%02dT%02d-%02d-%02dZ%s",PARAMS->outputFilePrefix,"_mix_",
		     UT(0), UT(1), UT(2),
		     UT(3), UT(4), UT(5), ".hdf");
  }

#ifdef DEBUG_MODE_ON
  cout << "MIX/src/IO.C: Opening MIX data file \"" << hdfName << "\"";
#endif

  if (openForWrite){
#ifdef DEBUG_MODE_ON
    cout << " for writting." << endl;
#endif
    hdfFileID = SDstart(hdfName,DFACC_CREATE);
    if (hdfFileID == -1){
      cerr << "ERROR opening \"" << hdfName << "\" for writing." << endl
	   << "      Is the disk full?  Do you have write permissions?." << endl;
      assert("ERROR opening MIX HDF file.  Is the disk full? Do you have write permissions?" == false);
    }
  }
  else{
#ifdef DEBUG_MODE_ON
    cout << " for reading." << endl;
#endif
    hdfFileID = SDstart(hdfName,DFACC_RDONLY);
    if (hdfFileID == -1){
      cerr << "ERROR opening \"" << hdfName << "\" for reading." << endl
	   << "      If you're doing a restart run, make sure that the previous dump exists." << endl;
      assert("ERROR opening MIX HDF restart file" == false);
    }
  }
}

/********************************************************************//**
 * \author Peter Schmitt (schmitt (at) ucar.edu)
 * \since 10-2008
 *
 * \param[in] tilt - tilt angle 
 * \param[in] UT - Simlation UT (YYYY-MM-DD HH-MM-SS) 
 *                 If year=-9999, then it's a "TimeStep" run and not UTIO.
 * \param[in] mhdTimeStep - timestep of simulation
 * 
 * Write HDF meta data and save some attributes (tilt angle, UT and time step) to an HDF file.
 *
 * \note HDF file must already be opened for read using IO::Open(...)
 *
 * \see IO::ReadAttributes
 ************************************************************************/
void IO::SaveAttributes(const real &tilt, const intArray &UT, const long &mhdTimeStep)
{
  // Save some attributes
  // 1. stuff received from the MHD code
  float fltTilt = (float)tilt;
  int status = SDsetattr(hdfFileID,"Tilt",DFNT_FLOAT32,1,(VOIDP)&fltTilt);
  status = SDsetattr(hdfFileID,"MHD time step",DFNT_INT32,1,(VOIDP)&mhdTimeStep); 
  status = SDsetattr(hdfFileID,"UT (Year, Month, Day, Hour, Min, Sec)",DFNT_INT32,6,
		     (VOIDP)(UT.getDataPointer()));
  DateTime utime(UT(0), UT(1), UT(2),
		 UT(3), UT(4), UT(5));
  double mjd = utime.getMJD();
  status = SDsetattr(hdfFileID, "mjd", DFNT_FLOAT64, 1, &mjd);

  // 2. Some useful parameters
  //  int status = SDsetattr(hdfFileID,"Local name",DFNT_CHAR,50,(VOIDP)(params->localName));
  //  status = SDsetattr(hdfFileID,"XJD name",DFNT_CHAR8,50,(VOIDP)(params->xjdName));
  //  status = SDsetattr(hdfFileID,"MIX file name",DFNT_CHAR8,STRLEN,(VOIDP)(params->ionGridFileName));
  //  status = SDsetattr(hdfFileID,"MHD grid file name",DFNT_CHAR8,STRLEN,(VOIDP)(params->mhdGridFileName));
  //  status = SDsetattr(hdfFileID,"MHD grid file type",DFNT_INT32,STRLEN,(VOIDP)(params->mhdGridFileType));
  //  status = SDsetattr(hdfFileID,"ITM grid file name",DFNT_CHAR8,STRLEN,(VOIDP)(params->itmGridFileName));
  //  status = SDsetattr(hdfFileID,"ITM grid file type",DFNT_INT32,1,(VOIDP)(params->itmGridFileType));
  //  status = SDsetattr(hdfFileID,"ITM model type",DFNT_CHAR8,STRLEN,(VOIDP)(params->itmModelType));
  
  status = SDsetattr(hdfFileID,"EUV model type",DFNT_INT32,1,(VOIDP)&(PARAMS->euvModelType));

  float fltAlpha = (float)(PARAMS->alpha);
  status = SDsetattr(hdfFileID,"Conductivity alpha",DFNT_FLOAT32,1,(VOIDP)&fltAlpha);

  float fltBeta = (float)(PARAMS->beta);
  status = SDsetattr(hdfFileID,"Conductivity beta",DFNT_FLOAT32,1,(VOIDP)&fltBeta);

  float fltR = (float)(PARAMS->r);
  status = SDsetattr(hdfFileID,"Conductivity R",DFNT_FLOAT32,1,(VOIDP)&fltR);

  float fltF107 = (float)(PARAMS->f107);
  status = SDsetattr(hdfFileID,"Conductivity F107",DFNT_FLOAT32,1,(VOIDP)&fltF107);

  float fltPedMin = (float)(PARAMS->pedmin);
  status = SDsetattr(hdfFileID,"Conductivity pedmin",DFNT_FLOAT32,1,(VOIDP)&fltPedMin);

  float fltHallMin = (float)(PARAMS->hallmin);
  status = SDsetattr(hdfFileID,"Conductivity hallmin",DFNT_FLOAT32,1,(VOIDP)&fltHallMin);

  float fltPed0 = (float)(PARAMS->ped0);
  status = SDsetattr(hdfFileID,"Conductivity ped0",DFNT_FLOAT32,1,(VOIDP)&fltPed0);

  float fltSigmaRatio = (float)(PARAMS->sigma_ratio);
  status = SDsetattr(hdfFileID,"Conductivity sigma_ratio",DFNT_FLOAT32,1,(VOIDP)&fltSigmaRatio);

  status = SDsetattr(hdfFileID,"Conductivity const_sigma",DFNT_INT32,1,(VOIDP)&(PARAMS->const_sigma));

  /*
   *   Use properties (svn:keywords) to add revision meta data to HDF file.
   *   http://svnbook.red-bean.com/en/1.1/ch07s02.html
   *
   *   Subversion defines five variables:
   *      - Date
   *      - Revision
   *      - Author
   *      - HeadURL
   *      - Id (a combination of the above)
   *
   *   These will auto-updated when the property is set for a source file:
   *
   *      svn propset svn:keywords "Id" [filename]
   *
   *   For now, we use the Id...
   *
   *   Where's $GlobalRev$?
   *      http://svnbook.red-bean.com/en/1.4/svn.advanced.props.special.keywords.html
   *
   *   Set Global Revision via method described here:
   *       http://subversion.tigris.org/faq.html#version-value-in-source
   *
   *   See also 
   *     - LFM-para/src/hdfdump-para.F
   *     - LFM-para/src/iondump-para.F
   */
  std::string revisionInformation = "$Id: IO.C 858 2009-04-21 19:13:35Z schmitt $";
  status = SDsetattr(hdfFileID,"I/O Revision",DFNT_CHAR8,revisionInformation.length(),
                     (VOIDP)(revisionInformation.c_str()));

  
  revisionInformation = REPOSITORY_REVISION_NUMBER.c_str();
  status = SDsetattr(hdfFileID,"Repository Revision",DFNT_CHAR8,revisionInformation.length(),
                     (VOIDP)(revisionInformation.c_str()));

}

/********************************************************************//**
 *  \author Slava Merkin (vgm at bu.edu)
 *  \since 06-2007 
 *
 *  \param[in] varFunction The grid function to be saved in the HDF file
 *
 *  \param[in] strLabel Label for the function
 *
 * This function was written before I have packed all the solver grid
 * functions into one as different components. Now that I have them
 * all together, we can use the other version of this function.
 ************************************************************************/
void IO::Add(const realCompositeGridFunction & varFunction, const char * strLabel)
{
  Index I1,I2,I3;
  getIndex(varFunction[0],I1,I2,I3);   // This will be based on indexRange, i.e. exclude ghost and periodic lines
  realArray varArray;
  varArray = varFunction[0](I1,I2);

  dataHDF.put(varArray,strLabel);
}

/********************************************************************//**
 *  \author Slava Merkin (vgm at bu.edu)
 *  \since 01-2008
 *
 *  \param[in] varFunction The grid function to be saved in the HDF file
 *
 * This function replaces the older one which had to be called
 * separately for different variables. Now they are all packed into
 * one grid function as components, so this function can be called
 * only once from ionosphere.C
 ************************************************************************/
void IO::Add(const realCompositeGridFunction & varFunction)
{
  MappedGrid *grid = varFunction[0].getMappedGrid();
  Index I1,I2,I3;
  getIndex((*grid).gridIndexRange(),I1,I2,I3); // Unlike the
						  // function
					       // above we are going
					       // to include the
					       // periodic points here
					       // (gridIndexRange) to
					       // make it compliant
					       // with the grid
					       // storage

  int32 status, sdID, dim_sizes[2], start[2];
  floatArray data;
  realArray varArray;

  varArray = varFunction[0](I1,I2,I3,SOLVER_POTENTIAL_NORTH);
  start[0] = varArray.getBase(1);
  start[1] = varArray.getBase(0);
  dim_sizes[0] = varArray.getBound(1) - start[1] + 1; // Number of elements along dimension 0
  dim_sizes[1] = varArray.getBound(0) - start[0] + 1; // Number of elements along dimension 1

  // Get grid vertices 
  const realArray & x = (*grid).vertex()(I1,I2,I3,0);
  const realArray & y = (*grid).vertex()(I1,I2,I3,1);

  //  Write grid X array
  sdID = SDcreate(hdfFileID,"Grid X",DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = x.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(x.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write grid Y array
  sdID = SDcreate(hdfFileID,"Grid Y",DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = y.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(y.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write POTENTIAL_NORTH array
  sdID = SDcreate(hdfFileID,varFunction.getName(SOLVER_POTENTIAL_NORTH),DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = varArray.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(varArray.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write POTENTIAL_SOUTH array
  varArray = varFunction[0](I1,I2,I3,SOLVER_POTENTIAL_SOUTH);
  sdID = SDcreate(hdfFileID,varFunction.getName(SOLVER_POTENTIAL_SOUTH),DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = varArray.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(varArray.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write FAC_NORTH array
  varArray = varFunction[0](I1,I2,I3,SOLVER_FAC_NORTH);
  sdID = SDcreate(hdfFileID,varFunction.getName(SOLVER_FAC_NORTH),DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = varArray.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(varArray.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write FAC_SOUTH array
  varArray = varFunction[0](I1,I2,I3,SOLVER_FAC_SOUTH);
  sdID = SDcreate(hdfFileID,varFunction.getName(SOLVER_FAC_SOUTH),DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = varArray.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(varArray.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write SIGMAP_NORTH array
  varArray = varFunction[0](I1,I2,I3,SOLVER_SIGMAP_NORTH);
  sdID = SDcreate(hdfFileID,varFunction.getName(SOLVER_SIGMAP_NORTH),DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = varArray.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(varArray.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write SIGMAP_SOUTH array
  varArray = varFunction[0](I1,I2,I3,SOLVER_SIGMAP_SOUTH);
  sdID = SDcreate(hdfFileID,varFunction.getName(SOLVER_SIGMAP_SOUTH),DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = varArray.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(varArray.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write SIGMAH_NORTH array
  varArray = varFunction[0](I1,I2,I3,SOLVER_SIGMAH_NORTH);
  sdID = SDcreate(hdfFileID,varFunction.getName(SOLVER_SIGMAH_NORTH),DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = varArray.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(varArray.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write SIGMAH_SOUTH array
  varArray = varFunction[0](I1,I2,I3,SOLVER_SIGMAH_SOUTH);
  sdID = SDcreate(hdfFileID,varFunction.getName(SOLVER_SIGMAH_SOUTH),DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = varArray.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(varArray.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write AVG_ENG_NORTH array
  varArray = varFunction[0](I1,I2,I3,SOLVER_AVG_ENG_NORTH);
  sdID = SDcreate(hdfFileID,varFunction.getName(SOLVER_AVG_ENG_NORTH),DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = varArray.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(varArray.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write AVG_ENG_SOUTH array
  varArray = varFunction[0](I1,I2,I3,SOLVER_AVG_ENG_SOUTH);
  sdID = SDcreate(hdfFileID,varFunction.getName(SOLVER_AVG_ENG_SOUTH),DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = varArray.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(varArray.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write NUM_FLUX_NORTH array
  varArray = varFunction[0](I1,I2,I3,SOLVER_NUM_FLUX_NORTH);
  sdID = SDcreate(hdfFileID,varFunction.getName(SOLVER_NUM_FLUX_NORTH),DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = varArray.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(varArray.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write NUM_FLUX_SOUTH array
  varArray = varFunction[0](I1,I2,I3,SOLVER_NUM_FLUX_SOUTH);
  sdID = SDcreate(hdfFileID,varFunction.getName(SOLVER_NUM_FLUX_SOUTH),DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = varArray.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(varArray.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write WIND_NORTH array
  varArray = varFunction[0](I1,I2,I3,SOLVER_WIND_NORTH);
  sdID = SDcreate(hdfFileID,varFunction.getName(SOLVER_WIND_NORTH),DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = varArray.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(varArray.getDataPointer()));
#endif
  status = SDendaccess(sdID);

  //  Write WIND_SOUTH array
  varArray = varFunction[0](I1,I2,I3,SOLVER_WIND_SOUTH);
  sdID = SDcreate(hdfFileID,varFunction.getName(SOLVER_WIND_SOUTH),DFNT_FLOAT32,2,dim_sizes);
#ifdef OV_USE_DOUBLE
  data = varArray.convertTo_floatArray();
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(data.getDataPointer()));
#else
  status = SDwritedata(sdID, start,NULL,dim_sizes,(VOIDP)(varArray.getDataPointer()));
#endif
  status = SDendaccess(sdID);


//   realArray varArray;

//   varArray = varFunction[0](I1,I2,I3,SOLVER_POTENTIAL_NORTH);
//   dataHDF.put(varArray,varFunction.getName(SOLVER_POTENTIAL_NORTH));

//   varArray = varFunction[0](I1,I2,I3,SOLVER_POTENTIAL_SOUTH);
//   dataHDF.put(varArray,varFunction.getName(SOLVER_POTENTIAL_SOUTH));

//   varArray = varFunction[0](I1,I2,I3,SOLVER_FAC_NORTH);
//   dataHDF.put(varArray,varFunction.getName(SOLVER_FAC_NORTH));

//   varArray = varFunction[0](I1,I2,I3,SOLVER_FAC_SOUTH);
//   dataHDF.put(varArray,varFunction.getName(SOLVER_FAC_SOUTH));

//   varArray = varFunction[0](I1,I2,I3,SOLVER_SIGMAP_NORTH);
//   dataHDF.put(varArray,varFunction.getName(SOLVER_SIGMAP_NORTH));

//   varArray = varFunction[0](I1,I2,I3,SOLVER_SIGMAP_SOUTH);
//   dataHDF.put(varArray,varFunction.getName(SOLVER_SIGMAP_SOUTH));

//   varArray = varFunction[0](I1,I2,I3,SOLVER_SIGMAH_NORTH);
//   dataHDF.put(varArray,varFunction.getName(SOLVER_SIGMAH_NORTH));

//   varArray = varFunction[0](I1,I2,I3,SOLVER_SIGMAH_SOUTH);
//   dataHDF.put(varArray,varFunction.getName(SOLVER_SIGMAH_SOUTH));

//   varArray = varFunction[0](I1,I2,I3,SOLVER_AVG_ENG_NORTH);
//   dataHDF.put(varArray,varFunction.getName(SOLVER_AVG_ENG_NORTH));

//   varArray = varFunction[0](I1,I2,I3,SOLVER_AVG_ENG_SOUTH);
//   dataHDF.put(varArray,varFunction.getName(SOLVER_AVG_ENG_SOUTH));

//   varArray = varFunction[0](I1,I2,I3,SOLVER_NUM_FLUX_NORTH);
//   dataHDF.put(varArray,varFunction.getName(SOLVER_NUM_FLUX_NORTH));

//   varArray = varFunction[0](I1,I2,I3,SOLVER_NUM_FLUX_SOUTH);
//   dataHDF.put(varArray,varFunction.getName(SOLVER_NUM_FLUX_SOUTH));

//   varArray = varFunction[0](I1,I2,I3,SOLVER_WIND_NORTH);
//   dataHDF.put(varArray,varFunction.getName(SOLVER_WIND_NORTH));

//   varArray = varFunction[0](I1,I2,I3,SOLVER_WIND_SOUTH);
//   dataHDF.put(varArray,varFunction.getName(SOLVER_WIND_SOUTH));
}

void IO::Close()
{
  int status = SDend(hdfFileID);
  //  dataHDF.unmount();
}