#ifndef MHD_H
#define MHD_H

#include "run-time.h"

/// Run-time common block defined in run-time.inc
extern RTime RTIME;

#include "site-stuff.h"

class ION;;
/// C++ class that handles MHD data and calls the FORTRAN MHD solver
/**
 * Magnetohydrodynamics class
 *
 * A C++ interface to the LFM MHD solver.  The actual MHD solver is
 * implmeneted in FORTRAN.  For more information on the numerics, see
 *    - mhd_fortran_arrays.F
 *    - <a href="http://scholar.google.com/scholar?hl=en&lr=&client=firefox-a&q=%22The+Lyon%E2%80%93Fedder%E2%80%93Mobarry+%28LFM%29+global+MHD+magnetospheric+simulation+code%22&btnG=Search">Lyon, Fedder, Mobarry, 2004.  "The Lyon-Fedder-Mobarry (LFM) global MHD magnetospheric simulation code".  Journal of Atmospheric and Solar-Terrestrial Physics 66, 1333-1350.</a>
 * 
 */
class MHD {
 public:
  MHD(int, int, int, int );   // defined in MHD.C
  MHD(int,int,int,int,Range); // defined in MHD.C
  //FIXME:  Should there be a detructor here???

  void initial(int,int,int);  // defined in MHD.C
  void startup( char *);      // defined in MHD.C
  void dummy_start();         // defined in MHD.C
  void boundary();
  void step(FLoaT *);         // defined in MHD.C
  float courant();

// dimensioning information

  // Accessor functions
  FloatARRAY *getX() { return &X; }
  FloatARRAY *getY() { return &Y; }
  FloatARRAY *getZ() { return &Z; }
  FloatARRAY *getBx() { return &Bx; }
  FloatARRAY *getBy() { return &By; }
  FloatARRAY *getBz() { return &Bz; }
  FloatARRAY *getBi() { return &Bi; }
  FloatARRAY *getBj() { return &Bj; }
  FloatARRAY *getBk() { return &Bk; }
  FloatARRAY *getRho() { return &Rho; }
  FloatARRAY *getCsound() { return &Csound; }
  FloatARRAY *getEi() { return &Ei; }
  FloatARRAY *getEj() { return &Ej; }
  FloatARRAY *getEk() { return &Ek; }
  FloatARRAY *getVx() { return &Vx; }
  FloatARRAY *getVy() { return &Vy; }
  FloatARRAY *getVz() { return &Vz; }
  int getNI() { return NI; }
  int getNJ() { return NJ; }
  int getNK() { return NK; }
  int getNO2() { return NO2; }  

private:

  int NI,NJ,NK,NORDER,NO2;
  int MHD_PROCESSOR;
  Partitioning_Type MHD_P;
  Partitioning_Type Tran_P;
  // grid data
  FloatARRAY X;
  FloatARRAY Y;
  FloatARRAY Z;
  FloatARRAY Volume;
  FloatARRAY Volume_inv;
  FloatARRAY Face_i;
  FloatARRAY Face_j;
  FloatARRAY Face_k;
  FloatARRAY Tran_i;
  FloatARRAY Tran_j;
  FloatARRAY Tran_k;
  FloatARRAY Edge_i;
  FloatARRAY Edge_j;
  FloatARRAY Edge_k;
  FloatARRAY Vec_i;
  FloatARRAY Vec_j;
  FloatARRAY Vec_k;
  FloatARRAY Cell_Length;

  // cell centered variables
  /// cell-centered density ( g/cm^3 )
  FloatARRAY Rho;
  /// x-component of cell-centered velocity ( cm/s )
  FloatARRAY Vx;
  /// y-component of cell-centered velocity ( cm/s )
  FloatARRAY Vy;
  /// z-component of cell-centered velocity ( cm/s )
  FloatARRAY Vz;
  /// cell-centered speed of sound ( cm/s )
  FloatARRAY Csound;
  /// cell-centered pressure ( g/(cm s^2) )
  FloatARRAY Pressure;  
  /// x-component of cell-centered magnetic field ( G )
  FloatARRAY Bx;  
  /// y-component of cell-centered magnetic field ( G )
  FloatARRAY By;  
  /// z-component of cell-centered magnetic field ( G )
  FloatARRAY Bz;  

  // face centered magnetic field variables
  /// i-component of magnetic flux through cell face ( G * cm^2 )
  FloatARRAY Bi;
  /// j-component of magnetic flux through cell face ( G * cm^2 )
  FloatARRAY Bj;
  /// k-component of magnetic flux through cell face ( G * cm^2 )
  FloatARRAY Bk;

  // edge centered electric field variables
  /// i-component of edge-centered electric field
  FloatARRAY Ei;
  /// j-component of edge-centered electric field
  FloatARRAY Ej;
  /// k-component of edge-centered electric field
  FloatARRAY Ek;


  // more arrays for second time level
  // cell centered variables
  /// cell-centered density rho for second time level
  FloatARRAY Rho_1;
  // x-component of cell-centered velocity for second time level
  FloatARRAY Vx_1;  
  // y-component of cell-centered velocity for second time level
  FloatARRAY Vy_1;  
  // z-component of cell-centered velocity for second time level
  FloatARRAY Vz_1;
  /// cell-centered speed of sound for second time level
  FloatARRAY Csound_1;
  /// cell-centered pressure for second time level
  FloatARRAY Pressure_1;  
  /// x-component of cell-centered magnetic field for second time level
  FloatARRAY Bx_1;  
  /// y-component of cell-centered magnetic field for second time level
  FloatARRAY By_1;  
  /// z-component of cell-centered magnetic field for second time level
  FloatARRAY Bz_1;  

  // face centered magnetic field variables
  /// i-component of face-centered magnetic field for second time level
  FloatARRAY Bi_1;
  /// j-component of face-centered magnetic field for second time level
  FloatARRAY Bj_1;
  /// k-component of face-centered magnetic field for second time level
  FloatARRAY Bk_1;

#ifdef DO_RINGAV
  //  variables needed for operation of Ringav
  FloatARRAY Face_av_k;
  FloatARRAY Rho_R;
  FloatARRAY Vx_R;  
  FloatARRAY Vy_R;  
  FloatARRAY Vz_R;
  FloatARRAY Csound_R;
  FloatARRAY Pressure_R;  
  FloatARRAY Bx_R;  
  FloatARRAY By_R;  
  FloatARRAY Bz_R;  
  // face centered magnetic field variables
  FloatARRAY Bi_R;
  FloatARRAY Bj_R;
  FloatARRAY Bk_R;
#endif  // end of RINGAV

// friend functions
#if defined (ION_ON) && !defined(USE_MIX)
  friend class ION;;
  friend class TestSuite; // see LFM-para/test/integration/IO/

  friend void mhd2ion(MHD*, ION*);
  friend void ion2mhd(MHD*, ION*);
  friend void grid2ion(MHD*, ION*);
#endif

  friend class MHDManager;

#ifdef DIVTEST
  friend float DivTest(MHD*,int);
#endif

#if defined(USE_RCM) && !defined(USE_MIX)
  friend class MHD2RCM;
#endif
};

#endif // MHD_H