#include "MHDManager.h"
#include "MHD.h"

/***********************************************************************/

/// Open an MHD HDF file for writting
/**
 * Creates a file named "[prefix]_mhd_[stepIdentification].[suffix]".  
 *
 *  If mjd > 0, assume UTIO... and write the MJD to the HDF file.
 *
 * @example Using UTIO naming, "MyRun_mhd_2008-01-01T01-01-00.hdf"
 * @example Using standard I/O naming: "MyRun_mhd_1234567.hdf"
 * 
 * @param prefix - prefix of filename
 * @param stepIdentification - unique string identifying a particular timestep (i.e. an 0-padded integer or a date/time string)
 * @param suffix - file extension (i.e. "hdf", "dmp", etc.)
 * @param mjd - Modified Julian Date.  If (mjd < 0), it is not stored to the file.
 *
 * @seealso DateTime class for more information on Modified Julian Dates
 */
bool MHDManager::open_write_hdf(char *prefix,  char *stepIdentification, char *suffix, double *mjd)
{
#ifdef DEBUG_MODE_ON
  printf("DEBUG:  process (%d) in MHDManager.C::open_write_hdf(...)\n", RTIME.local_proc);
#endif  
  int istatus = MHD_OPEN_WRITE_HDF(prefix, stepIdentification, suffix, mjd, &currentStep, &simTime);
  return true;
}

/***********************************************************************/

/// Write MHD variables to file.
/** 
 * @param mhd - MHD object that contains the correct pointers to the FORTRAN data structures.
 */
bool MHDManager::write_data(MHD *mhd, const Average<FloatARRAY> &avg )
{
  return write_data(mhd, NULL, avg);
}

/***********************************************************************/

/// Write MHD variables to file, including accumulate arrays.
/** 
 * @param mhd - MHD object that contains the correct pointers to the FORTRAN data structures.
 * @param innerMag - Inner Magnetosphere coupling object.  Set to NULL if there is no data to write.
 */
bool MHDManager::write_data(MHD *mhd, IM_Interface const * const innerMag, 
			    const Average<FloatARRAY> &avg)
{
#ifdef DEBUG_MODE_ON
  printf("DEBUG:  process (%d) in MHDManager.C::write_data(...)\n", RTIME.local_proc);
#endif  
  int no2, ni, nj, nk, nii_0, njj_0, nkk_0;
  int *nii = &nii_0, *njj = &njj_0, *nkk = &nkk_0;

  int istatus = 0;  // everything O.K.

  //  int C_HDFDump_first( FLoaT *, FLoaT *, FLoaT *);
  //  int C_HDFDump_var(char *, FLoaT *, int *, int *, int *, int *, FLoaT *, int *);

  no2 = mhd->NO2; ni = mhd->NI; nj = mhd->NJ; nk = mhd->NK;;
  //Partitioning_Type pppDummy = mhd->MHD_P;
  Range Iall(-no2+1,ni+no2); // Fortran indexing (1 - NI) active
  Range Jall(-no2+1,nj+no2); 
  Range Kall(-no2+1,nk+no2);
  Range I(1,ni);
  Range J(1,nj);
  Range K(1,nk);
  Range Ip1(1,ni+1);
  Range Jp1(1,nj+1);
  Range Kp1(1,nk+1);
  Range Imid(ni/2,ni/2);
  Range Jmid(nj/2,nj/2);
  Range Izero(1,1);
  Range Jzero(1,1);
  
  //   Guard partition matches MHD_P -- needed for data transfer
  Range Master_range(0,0);
  Partitioning_Type Guard(Master_range);
  Guard.partitionAlongAxis(0,TRUE,no2);
  Guard.partitionAlongAxis(1,TRUE,no2);
  Guard.partitionAlongAxis(2,FALSE,no2);
  //  noGuard strips guard cells to make indexing easy for dump
  
  FloatARRAY D_VAR(ni+1,nj+1,nk+1,Guard);
  D_VAR.setBase(1);
  
  //  Grid variables
  int grid = -1; int not_grid = 0;
  
  D_VAR(Ip1,Jp1,Kp1) = mhd->X(Ip1,Jp1,Kp1);  
  *nii=ni+1;  *njj=nj+1;   *nkk=nk+1; int *npad = &no2;
  
  istatus = HDFDUMP_VAR(fcharred("X_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&grid,fcharred("cm\0")); 

  D_VAR(Ip1,Jp1,Kp1) = mhd->Y(Ip1,Jp1,Kp1);
  istatus = HDFDUMP_VAR(fcharred("Y_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&grid,fcharred("cm\0")); 

  D_VAR(Ip1,Jp1,Kp1) = mhd->Z(Ip1,Jp1,Kp1);
  istatus = HDFDUMP_VAR(fcharred("Z_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&grid,fcharred("cm\0")); 

// Density
// Rho
    D_VAR(Ip1,Jp1,Kp1) = mhd->Rho(Ip1,Jp1,Kp1);  
   istatus = HDFDUMP_VAR(fcharred("rho_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("gm/cm^3\0"));

// Vx
	D_VAR(Ip1,Jp1,Kp1) = mhd->Vx(Ip1,Jp1,Kp1);
   istatus = HDFDUMP_VAR(fcharred("vx_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("cm/sec\0"));

// Vy
	D_VAR(Ip1,Jp1,Kp1) = mhd->Vy(Ip1,Jp1,Kp1); 
   istatus = HDFDUMP_VAR(fcharred("vy_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("cm/sec\0"));

// Vz
	D_VAR(Ip1,Jp1,Kp1) = mhd->Vz(Ip1,Jp1,Kp1); 
   istatus = HDFDUMP_VAR(fcharred("vz_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("cm/sec\0"));

// Sound Speed
	D_VAR(Ip1,Jp1,Kp1) = mhd->Csound(Ip1,Jp1,Kp1);
   istatus = HDFDUMP_VAR(fcharred("c_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("cm/sec\0"));

// magnetic field variables need to add back in the subtracted stuff
   
	FIXBOUT();  // adds in subtracted field to bx,by,bz
	AVE2FLUX(); // turns it into a flux over the face

// Bx
	D_VAR(Ip1,Jp1,Kp1) = mhd->Bx(Ip1,Jp1,Kp1);
   istatus = HDFDUMP_VAR(fcharred("bx_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("Gauss\0"));

// By
	D_VAR(Ip1,Jp1,Kp1) = mhd->By(Ip1,Jp1,Kp1);
   istatus = HDFDUMP_VAR(fcharred("by_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("Gauss\0"));

// Bz
	D_VAR(Ip1,Jp1,Kp1) = mhd->Bz(Ip1,Jp1,Kp1);
   istatus = HDFDUMP_VAR(fcharred("bz_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("Gauss\0"));

// Bi
	D_VAR(Ip1,Jp1,Kp1) = mhd->Bi(Ip1,Jp1,Kp1);
   istatus = HDFDUMP_VAR(fcharred("bi_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("Gauss-cm^2\0"));

// Bj
	D_VAR(Ip1,Jp1,Kp1) = mhd->Bj(Ip1,Jp1,Kp1);
   istatus = HDFDUMP_VAR(fcharred("bj_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("Gauss-cm^2\0"));

// Bk
	D_VAR(Ip1,Jp1,Kp1) = mhd->Bk(Ip1,Jp1,Kp1);
   istatus = HDFDUMP_VAR(fcharred("bk_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("Gauss-cm^2\0"));
//   now subtract them back off

	FLUX2AVE();
	RESETB();

// Ei
	D_VAR(Ip1,Jp1,Kp1) = mhd->Ei(Ip1,Jp1,Kp1);
   istatus = HDFDUMP_VAR(fcharred("ei_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("[cm/sec*Gauss]*cm\0"));

// Ej
	D_VAR(Ip1,Jp1,Kp1) = mhd->Ej(Ip1,Jp1,Kp1);
   istatus = HDFDUMP_VAR(fcharred("ej_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("[cm/sec*Gauss]*cm\0"));

// Ek
	D_VAR(Ip1,Jp1,Kp1) = mhd->Ek(Ip1,Jp1,Kp1);
   istatus = HDFDUMP_VAR(fcharred("ek_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("[cm/sec*Gauss]*cm\0"));
   

   // write variables stored in Average object
   map<int, string> avgVars = avg.getVariableIds();
   typedef map<int, string> mapType;

   for (mapType::const_iterator it = avgVars.begin(); it != avgVars.end(); ++it){
     Communication_Manager::Sync();
     // it-> first is the "key"
     D_VAR = avg.getVariable(it->first);
     // it-> second is the "value", mapped to variable identifier
     string varName = it->second + "\0";
     istatus = HDFDUMP_VAR(fcharred(const_cast<char *>(varName.c_str())), D_VAR.getDataPointer(), nii,njj,nkk,npad,&currentStep,&simTime,&not_grid,fcharred("Units of source variable\0"));
   }

   // Write variables stored in Inner Magnetosphere coupling object.
   bool return_value = true;
   if (innerMag != NULL){
     return_value = write_accumulate_and_bleed_data(mhd, innerMag);
   }

   return return_value;
}

/// Write MHD accumulate and bleed arrays to file.  Used for LFM-RCM restart.
/** 
 * @param innerMag - Inner Magnetosphere coupling interface
 */
bool MHDManager::write_accumulate_and_bleed_data(MHD *mhd, IM_Interface const * const innerMag)
{
  int iStatus;
  int not_grid = 0;

  floatArray D_VAR;

  ////////////////////////////////////
  // Accumulate Arrays              //
  //    dimensions: nip1,njp1,nkp1  //
  //    ghost cells: standard LFM   //
  ////////////////////////////////////
  int nip1 = mhd->NI + 1;
  int njp1 = mhd->NJ + 1;
  int nkp1 = mhd->NK + 1;
  int npad = mhd->NO2;
  D_VAR.redim(nip1, njp1, nkp1);
  Partitioning_Type Guard(Range(0,0));
  Guard.partitionAlongAxis(0, TRUE, npad);
  Guard.partitionAlongAxis(1, TRUE, npad);
  Guard.partitionAlongAxis(2, FALSE, npad);
  D_VAR.partition(Guard);

  // Accumulate Density
  D_VAR.setBase(0);
  D_VAR = innerMag->getAccumulateDensity();
  D_VAR.setBase(1);
  iStatus = HDFDUMP_VAR(fcharred("accum_rho\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &currentStep, &simTime, &not_grid, fcharred("g/cm^3\0"));

  // Accumulate Pressure
  D_VAR.setBase(0);
  D_VAR = innerMag->getAccumulatePressure();
  D_VAR.setBase(1);
  iStatus = HDFDUMP_VAR(fcharred("accum_p\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &currentStep, &simTime, &not_grid, fcharred("g/(cm s^2)\0"));

  // Accumulate Bi
  D_VAR.setBase(0);
  D_VAR = innerMag->getAccumulateBi();
  D_VAR.setBase(1);
  iStatus = HDFDUMP_VAR(fcharred("accum_Bi\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &currentStep, &simTime, &not_grid, fcharred("Gauss\0"));

  // Accumulate Bj
  D_VAR.setBase(0);
  D_VAR = innerMag->getAccumulateBj();
  D_VAR.setBase(1);
  iStatus = HDFDUMP_VAR(fcharred("accum_Bj\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &currentStep, &simTime, &not_grid, fcharred("Gauss\0"));

  // Accumulate Bk
  D_VAR.setBase(0);
  D_VAR = innerMag->getAccumulateBk();
  D_VAR.setBase(1);
  iStatus = HDFDUMP_VAR(fcharred("accum_Bk\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &currentStep, &simTime, &not_grid, fcharred("Gauss\0"));

  // Accumulate V
  D_VAR.setBase(0);
  D_VAR = innerMag->getAccumulateV();
  D_VAR.setBase(1);
  iStatus = HDFDUMP_VAR(fcharred("accum_V\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &currentStep, &simTime, &not_grid, fcharred("cm/s\0"));

#ifdef SAVE_ACCUMULATE_TO_HDF
  //FIXME: This bit of code in the preprocessor flag
  // "SAVE_ACCUMULATE_TO_HDF" should be configurable by a variable to
  // the INPUT1.xml file.  This code is LFM-RCM specific.
  //
  // Why does this code exist?  See #81: https://oberon.bu.edu/bugs/issues/show/81

  // Update the accum_b{x,y,z} values.
  innerMag->subtractDipole();

  // Accumulate Bx (dipole removed @ cell center)
  // note:  this is only valid at ni,nj,nk, but we write data at all indices for consistency.
  D_VAR.setBase(0);
  D_VAR = innerMag->getAccumulateBx();
  D_VAR.setBase(1);
  iStatus = HDFDUMP_VAR(fcharred("accum_Bx\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &currentStep, &simTime, &not_grid, fcharred("Gauss\0"));

  // Accumulate By (dipole removed @ cell center)
  // note:  this is only valid at ni,nj,nk, but we write data at all indices for consistency.
  D_VAR.setBase(0);
  D_VAR = innerMag->getAccumulateBy();
  D_VAR.setBase(1);
  iStatus = HDFDUMP_VAR(fcharred("accum_By\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &currentStep, &simTime, &not_grid, fcharred("Gauss\0"));

  // Accumulate Bz (dipole removed @ cell center)
  // note:  this is only valid at ni,nj,nk, but we write data at all indices for consistency.
  D_VAR.setBase(0);
  D_VAR = innerMag->getAccumulateBz();
  D_VAR.setBase(1);
  iStatus = HDFDUMP_VAR(fcharred("accum_Bz\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &currentStep, &simTime, &not_grid, fcharred("Gauss\0"));
#endif


  ///////////////////////////////////
  // Bleed Arrays                  //
  //    dimensions: ni,nj,nk)      //
  //    ghost cells: standard LFM  //
  ///////////////////////////////////
  int ni = mhd->NI;
  int nj = mhd->NJ;
  int nk = mhd->NK;
  D_VAR.redim(ni,nj,nk);
  D_VAR.partition(Guard);

  // Bleed Density
  D_VAR.setBase(0);
  D_VAR = innerMag->getBleedDensity();
  D_VAR.setBase(1);
  iStatus = HDFDUMP_VAR(fcharred("bleed_rho\0"), D_VAR.getDataPointer(), &ni,&nj,&nk, &npad, &currentStep, &simTime, &not_grid, fcharred("g/cm^3\0"));

  // Bleed Pressure
  D_VAR.setBase(0);
  D_VAR = innerMag->getBleedPressure();
  D_VAR.setBase(1);
  iStatus = HDFDUMP_VAR(fcharred("bleed_p\0"), D_VAR.getDataPointer(), &ni,&nj,&nk, &npad, &currentStep, &simTime, &not_grid, fcharred("g/(cm s^2)\0"));

  // Bleed Mask
  intArray int_D_VAR(ni,nj,nk, Guard);
  int_D_VAR = innerMag->getBleedMask();
  D_VAR.setBase(0);
  D_VAR = int_D_VAR.convertTo_floatArray();
  D_VAR.setBase(1);
  iStatus = HDFDUMP_VAR(fcharred("bleed_msk\0"), D_VAR.getDataPointer(), &ni,&nj,&nk, &npad, &currentStep, &simTime, &not_grid, fcharred("boolean\0"));

  // Success if we've made it this far!
  return true;
}

/***********************************************************************/

/// Open an MHD HDF file for reading
/**
 * Reads data from a file named "[prefix]_mhd_[stepIdentification].[suffix]".  
 *
 * @example Using UTIO naming, "MyRun_mhd_2008-01-01T01-01-00.hdf"
 * @example Using standard I/O naming: "MyRun_mhd_1234567.hdf"
 * 
 * @param prefix - prefix of filename
 * @param stepIdentification - unique string (i.e. an 0-padded integer or a date/time string) identifying a particular timestep to read.
 * @param suffix - file extension (i.e. "hdf", "dmp", etc.)
 * @param mjd - Modified Julian Date obtained from the file.  
 *
 * @seealso DateTime class for more information on Modified Julian Dates
 */
bool MHDManager::open_read_hdf(char *prefix,  char *stepIdentification, char *suffix, double *mjd)
{
#ifdef DEBUG_MODE_ON
  printf("DEBUG:  process (%d) in MHDManager.C::open_read_hdf(...)\n", RTIME.local_proc);
#endif  
  int istatus = MHD_OPEN_READ_HDF(prefix, stepIdentification, suffix, mjd, &currentStep, &simTime);

  return true;
}

/***********************************************************************/

/// Read MHD variables from file.
/** 
 * @param mhd - MHD object that contains the correct pointers to the FORTRAN data structures.
 *
 * @fixme - Make sure that currentStep and simTime are read in properly!
 */
bool MHDManager::read_data(MHD *mhd)
{
#ifdef DEBUG_MODE_ON
  printf("DEBUG:  process (%d) in MHDManager.C::read_data(...)\n", RTIME.local_proc);
#endif  
  int no2, ni, nj, nk, nii_0, njj_0, nkk_0;
  int *nii = &nii_0, *njj = &njj_0, *nkk = &nkk_0;
  int *iStep = &currentStep;  /* We will set currentStep while reading data */
  double *Time = &simTime;     /* We will set simTime while reading data*/
  //  int HDFTAKE_FIRST(fcharred(char *, FLoaT *, FLoaT *, FLoaT *, int *, FLoaT *, int *, int *);
  //  int C_HDFTake_var(char *, FLoaT *, int *, int *, int *, int *, FLoaT *, int *);

  int istatus = 0;  // everything O.K.

  no2 = mhd->NO2; ni = mhd->NI; nj = mhd->NJ; nk = mhd->NK;
  Range Iall(-no2+1,ni+no2); // Fortran indexing (1 - NI) active
  Range Jall(-no2+1,nj+no2); 
  Range Kall(-no2+1,nk+no2);
  Range I(1,ni);
  Range J(1,nj);
  Range K(1,nk);
  Range Ip1(1,ni+1);
  Range Jp1(1,nj+1);
  Range Kp1(1,nk+1);  

  //   Guard partition matches MHD_P -- needed for data transfer
  Range Master_range(0,0);
  Partitioning_Type Guard(Master_range);
  Guard.partitionAlongAxis(0,TRUE,no2);
  Guard.partitionAlongAxis(1,TRUE,no2);
  Guard.partitionAlongAxis(2,FALSE,no2);
  FloatARRAY D_VAR(ni+1,nj+1,nk+1,Guard);
  D_VAR.setBase(1);

  //  Grid variables
  int grid = -1; int not_grid = 0;
  *nii=ni+1;  *njj=nj+1;   *nkk=nk+1; int *npad = &no2;

  istatus = HDFTAKE_VAR(fcharred("X_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&grid); 
    //    FLoaT *DDV = D_VAR.getDataPointer();
           mhd->X(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

  istatus = HDFTAKE_VAR(fcharred("Y_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&grid); 
      mhd->Y(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

  istatus = HDFTAKE_VAR(fcharred("Z_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&grid); 
      mhd->Z(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

// Density
   istatus = HDFTAKE_VAR(fcharred("rho_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&not_grid);
              mhd->Rho(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

// Vx
   istatus = HDFTAKE_VAR(fcharred("vx_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&not_grid);
         	 mhd->Vx(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

// Vy
   istatus = HDFTAKE_VAR(fcharred("vy_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&not_grid);
         	 mhd->Vy(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

// Vz
   istatus = HDFTAKE_VAR(fcharred("vz_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&not_grid);
         	 mhd->Vz(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

// Sound Speed
   istatus = HDFTAKE_VAR(fcharred("c_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&not_grid);
         	 mhd->Csound(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

// Bx
   istatus = HDFTAKE_VAR(fcharred("bx_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&not_grid);
         	 mhd->Bx(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

// By
   istatus = HDFTAKE_VAR(fcharred("by_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&not_grid);
         	 mhd->By(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

// Bz
   istatus = HDFTAKE_VAR(fcharred("bz_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&not_grid);
         	 mhd->Bz(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);


// Bi
   istatus = HDFTAKE_VAR(fcharred("bi_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&not_grid);
         	 mhd->Bi(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

// Bj
   istatus = HDFTAKE_VAR(fcharred("bj_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&not_grid);
         	 mhd->Bj(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

// Bk
   istatus = HDFTAKE_VAR(fcharred("bk_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&not_grid);
         	 mhd->Bk(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

// Ei
   istatus = HDFTAKE_VAR(fcharred("ei_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&not_grid);
         	 mhd->Ei(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

// Ej
   istatus = HDFTAKE_VAR(fcharred("ej_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad,iStep,Time,&not_grid);
         	 mhd->Ej(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

// Ek
   istatus = HDFTAKE_VAR(fcharred("ek_\0"),D_VAR.getDataPointer(),nii,njj,nkk,npad, iStep,Time,&not_grid);
         	 mhd->Ek(Ip1,Jp1,Kp1) = D_VAR(Ip1,Jp1,Kp1);

   return true;
}

/***********************************************************************/

/// Read Inner Magnetosphere coupling variables into IM_Interface object
/** 
 * @param mhd - MHD object that contains the correct pointers to the FORTRAN data structures.
 * @param innerMag - Inner Magnetosphere coupling object that contains
 * the correct pointers to the accumulate & bleed arrays that we
 * should read data into.
 */
bool MHDManager::read_accumulate_and_bleed_data(MHD *mhd, IM_Interface * const innerMag)
{
#ifdef DEBUG_MODE_ON
  cout << "DEBUG: process (" << RTIME.local_proc << ") at " << __FILE__ 
       <<  " (Line " << __LINE__ << "):  " << __FUNCTION__ << endl;
#endif

  int iStatus;
  int iStep;
  double simTime;
  int not_grid = 0;

  floatArray D_VAR;

  ////////////////////////////////////
  // Accumulate Arrays              //
  //    dimensions: nip1,njp1,nkp1  //
  //    ghost cells: standard LFM   //
  ////////////////////////////////////
  int nip1 = mhd->NI + 1;
  int njp1 = mhd->NJ + 1;
  int nkp1 = mhd->NK + 1;
  int npad = mhd->NO2;
  D_VAR.redim(nip1, njp1, nkp1);
  Partitioning_Type Guard(Range(0,0));
  Guard.partitionAlongAxis(0, TRUE, npad);
  Guard.partitionAlongAxis(1, TRUE, npad);
  Guard.partitionAlongAxis(2, FALSE, npad);
  D_VAR.partition(Guard);

  // FIXME:  Seems like HDFTAKE_VAR doesn't set iStatus, so how do we tell if there's an error?
  // Kludge fix:  print warning message telling this user what _might_ be going on.
  if ( RTIME.local_proc == 0){
    cout << "If the code dies here, there was an error reading 'accum_*' arrays from file!  "
	 << "Is this a LFM-RCM restart?  Does the hdf restart file contain " 
	 << "bleed & accumulate data?" << endl;
    cout << "See file \"" <<  __FILE__ << "\", Function \"" << __FUNCTION__ << ", Line  " << __LINE__ << endl;
  }

  // Accumulate Density
  D_VAR.setBase(1);
  iStatus = HDFTAKE_VAR(fcharred("accum_rho\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &iStep, &simTime, &not_grid);
  D_VAR.setBase(0);
  innerMag->setAccumulateDensity(D_VAR);

  // Accumulate Pressure
  D_VAR.setBase(1);
  iStatus = HDFTAKE_VAR(fcharred("accum_p\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &iStep, &simTime, &not_grid);
  D_VAR.setBase(0);
  innerMag->setAccumulatePressure(D_VAR);

  // Accumulate Bi
  D_VAR.setBase(1);
  iStatus = HDFTAKE_VAR(fcharred("accum_Bi\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &iStep, &simTime, &not_grid);
  D_VAR.setBase(0);
  innerMag->setAccumulateBi(D_VAR);

  // Accumulate Bj
  D_VAR.setBase(1);
  iStatus = HDFTAKE_VAR(fcharred("accum_Bj\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &iStep, &simTime, &not_grid);
  D_VAR.setBase(0);
  innerMag->setAccumulateBj(D_VAR);

  // Accumulate Bk
  D_VAR.setBase(1);
  iStatus = HDFTAKE_VAR(fcharred("accum_Bk\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &iStep, &simTime, &not_grid);
  D_VAR.setBase(0);
  innerMag->setAccumulateBk(D_VAR);

  // Accumulate V
  D_VAR.setBase(1);
  iStatus = HDFTAKE_VAR(fcharred("accum_V\0"), D_VAR.getDataPointer(), &nip1,&njp1,&nkp1, &npad, &iStep, &simTime, &not_grid);
  D_VAR.setBase(0);
  innerMag->setAccumulateV(D_VAR);

  ///////////////////////////////////
  // Bleed Arrays                  //
  //    dimensions: ni,nj,nk)      //
  //    ghost cells: standard LFM  //
  ///////////////////////////////////
  int ni = mhd->NI;
  int nj = mhd->NJ;
  int nk = mhd->NK;
  D_VAR.redim(ni,nj,nk);
  D_VAR.partition(Guard);

  // FIXME:  HDFTAKE_VAR doesn't set iStatus, so how do we tell if there's an error?
  // Kludge fix:  print warning message telling this user what _might_ be going on.
  if ( RTIME.local_proc == 0){
    cout << "If the code dies here, there is likely was an error reading 'bleed_*' arrays from file!  "
	 << "Is this a LFM-RCM restart?  Does the hdf restart file contain " 
	 << "bleed & accumulate data?" << endl;
    cout << "See file \"" <<  __FILE__ << "\", Function \"" << __FUNCTION__ << ", Line  " << __LINE__ << endl;
  }

  // Bleed Density
  D_VAR.setBase(1);
  iStatus = HDFTAKE_VAR(fcharred("bleed_rho\0"), D_VAR.getDataPointer(), &ni,&nj,&nk, &npad, &iStep, &simTime, &not_grid);
  D_VAR.setBase(0);
  innerMag->setBleedDensity(D_VAR);

  // Bleed Pressure
  D_VAR.setBase(1);
  iStatus = HDFTAKE_VAR(fcharred("bleed_p\0"), D_VAR.getDataPointer(), &ni,&nj,&nk, &npad, &iStep, &simTime, &not_grid);
  D_VAR.setBase(0);
  innerMag->setBleedPressure(D_VAR);

  // Bleed Mask
  D_VAR.setBase(1);
  iStatus = HDFTAKE_VAR(fcharred("bleed_msk\0"), D_VAR.getDataPointer(), &ni,&nj,&nk, &npad, &iStep, &simTime, &not_grid);
  D_VAR.setBase(0);
  intArray int_D_VAR(ni,nj,nk, Guard);
  int_D_VAR = D_VAR.convertTo_intArray();
  innerMag->setBleedMask(int_D_VAR);

  // Success if we've made it this far!
  return true;
}
/***********************************************************************/

/// Close access to the HDF file.
bool MHDManager::close_hdf(void)
{
#ifdef DEBUG_MODE_ON
  printf("DEBUG:  process (%d) in MHDManager.C::close_hdf(...)\n", RTIME.local_proc);
#endif  
  int istatus;
  istatus = HDFDUMP_CLOSE();
  return true;
}

/***********************************************************************/