/**
 *  @file MSPHERE.C
 *  @brief The main LFM/MSPHERE program
 *
 *  @mainpage
 *  @authors John Lyon
 *  @authors Joel Fedder
 *  @authors Clark Mobarry
 *  @since 1985
 *  @date 20 May 2008
 *
 *  <hr width="75%" size="1" />
 *
 * The Lyon-Fedder-Mobarry (LFM) Code is an integrated simulation
 * model for the global magnetosphere-ionosphere system. The heart of
 * the model is a time-dependent, ideal MHD calculation of the state
 * of the magnetosphere. This magnetospheric model is tightly coupled
 * to a realistic model for the polar ionospheres and is driven by
 * solar wind plasma and magnetic field data upwind of the calculation
 * domain. While lacking important physical processes in both the
 * magnetosphere and ionosphere, the model is the simplest,
 * self-consistent first principles model possible. With current
 * computational capabilities, it is also the only type of model
 * capable of providing reasonable quantitative accuracy in simulating
 * the magnetosphere-ionosphere system.
 *
 * The architecture of the LFM is based on a number of important
 * design considerations:
 *
 *   -  TVD transport with the highest possible resolving power.
 *   -  div B = 0 maintained to roundoff.  
 *   -  grid adapted to the problem, packing more resolution where needed.
 *   -  integrated ionospheric model.
 *
 * The outer boundary of the MHD calculation is handled as either
 * inflow or outflow. The solar wind in the vicinity of the Earth is
 * generally highly supersonic, so inflow is appropriate for those
 * directions where the solar wind passes into the computational
 * domain. The actual inflow conditions are usually described from an
 * appropriately time shifted time series of solar wind data observed
 * by satellite. Other regions of the outer boundary use simple
 * outflow boundaries. The computational domain is large enough that
 * flows have re-accelerated back to supersonic by the time they reach
 * the edge of the grid.
 *
 * The low altitude boundary of the LFM is handled by incorporating a
 * tightly coupled model for the ionospheric electric field. The MHD
 * calculation is stopped at some altitude between 2 and 3
 * RE. Field-aligned currents are calculated and mapped along dipole
 * field lines to the ionosphere where they are used as the source
 * term for the height-integrated potential equation. The calculated
 * potential is then mapped back out to the MHD lower boundary where
 * it is used to determine a boundary condition for the velocity and
 * electric field.  The algorithms used in the LFM code are described
 * in some detail in Lyon et al. [2004]. Details of the calculation of
 * the ionospheric conductances can be found in Fedder et al. [1995].
 *
 * The LFM grid is a distorted spherical grid with azimuthal symmetry
 * about the polar axis. The polar axis points in the x-direction of
 * the SM coordinate system (roughly toward the Sun). In the (r-?)
 * planes the grid is distorted to place resolution where chosen by
 * the user. The MHD grid usually covers the domain from about 30 RE
 * upwind to 300 RE downwind of the Earth and roughly 100 RE out to
 * the sides. The complementary ionospheric grid is a mapping of the
 * inner surface of the MHD (magnetospheric) domain to the two polar
 * ionospheres. This typically covers the region from the pole to
 * 45-60 degrees latitude.  * <hr width="75%" size="1" />
 *
 *  For more information, see:
 *    - <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>
*      - Paper that serves as the documentation of the numerics behind the LFM.
*    - <a href="http://www.glue.umd.edu/~shay/papers/coupling.pdf">Goodrich, Sussman, Lyon, Shay & Cassak, 2004.  "The CISM code coupling strategy", .  Journal of Atmospheric and Solar-Terrestrial Physics 66, 1469-1479.</a>
*      - Paper that describes how the different codes are going to be coupled to one another.
*    - <a href="https://wiki.ucar.edu/display/LTR">LTR Wiki</a>
*/


//#define BOUNDS_CHECK
#include "site-stuff.h"
#include "time.h"

#include <mpi.h>

#include "DateTime.h"

#include <sys/times.h>
#ifdef TRAP_FPE
#include <sigfpe.h>
#endif

#include <stdio.h>

#if NEWCC
#include <iostream>
#else
#include <iostream.h>
#endif //NEWCC

#include <A++.h>
#include "param.h"
#include "MHD.h"
#if defined(ION_ON) && !defined(USE_MIX)
#include "ION.h"
#endif
#include "help.h"
extern HELPER HELP;

#define PPLUSPLUS
#include <IC_EndPointSet.h>
 
// make global definitions of structures/common blocks
/// "help" common block defined in help.inc
struct RTime RTIME;

#if defined(ION_ON) && !defined(USE_MIX)
#include "IonIO.h"
#include "UTIonIO.h"
#include "TimestepIonIO.h"
#endif

#include "MHDManager.h"
#include "UTMHD.h"
#include "TimestepMHD.h"

#ifdef USE_MIX
#include "MHDInnerBoundaryInterface.h"
#include "CPL_IC_Interface.h"
#include "CPL_FE_Interface.h"
float getTilt();
#ifdef LC_UNSC
#define  GET_DIP_MOMENT get_dip_moment_
#define  GET_TILT get_tilt_
#elif LC_NOUNSC
#define  GET_DIP_MOMENT get_dip_moment
#define GET_TILT get_tilt
#endif

/// Defined in MSPHERE.C
extern "C" void GET_TILT(float *);

void parse_xjd(const char *, char *, char *, int &);

#endif // USE_MIX

#ifndef DISABLE_RCM
extern "C" void GET_DIP_MOMENT(float *);
#include "IM_IC_Interface.h"
#endif //DISABLE_RCM
#include "rcm-ut.h"



/********************************************************************//**
* Data to keep track of model coupling times
************************************************************************/
struct modelCoupling {
  DateTime start;  /// starting date/time to couple
  DateTime next;   /// next date/time to exchange
  int delta_t;     /// frequency (in seconds) to exchange
  int exchangeNum; /// How many times have we exchanged?
  bool isFirstExchange;
  bool isLastExchange;
};

int calcExchangeNum(const DateTime &start, const DateTime &current, const int &delta_t);
bool isLastExchange(const DateTime &stop, const int &delta_t, const DateTime &current);

#ifdef LC_UNSC
#define DTIME dtime_
#define LOAD_B load_b_
#define EXAM_B exam_b_
#define ION_TEST_ARRAYS ion_test_arrays_
#define TEST_ARRAYS test_arrays_
#elif LC_NOUNSC
#define DTIME dtime
#define LOAD_B load_b
#define EXAM_B exam_b
#define ION_TEST_ARRAYS ion_test_arrays
#define TEST_ARRAYS test_arrays
#endif

/// Defined in second.c
extern "C" float DTIME(float *);
/// Defined in test-arrays.F
extern "C" void TEST_ARRAYS();
/// Defined in test-arrays.F
extern "C" void ION_TEST_ARRAYS();
/// Defined in test-arrays.F
extern "C" void LOAD_B();
/// Defined in test-arrays.F
extern "C" void EXAM_B();

#ifdef USE_DEPRECATED_NAMELISTS
#ifdef LC_UNSC
#define NAMELISTS namelists_
#elif LC_NOUNSC
#define NAMELISTS namelists
#endif /* LC_UNSC */

/// Defined in init-fortran.F
extern "C" void NAMELISTS(char *, char *, int *, int *, int *);

# else /* USE_DEPRECATED_NAMLIESTS */
#ifdef LC_UNSC
#define INITIALIZE initialize_
#elif LC_NOUNSC
#define INITIALIZE initialize
#endif
extern "C" void INITIALIZE(void);

/*
 * FIXME Instead of just passing lots of parameters as arguments to the
 *       xmlInputReader function, we should probably use a "Params" class
 *       similar to "MIX/src/Params.h"
 */
bool xmlInputReader(const char * in_xmlFilename, const int in_ni_bin_d, const int in_nj_bin_d,
                    char *out_HDFBASE_IN, char *out_HDFBASE_OUT,
                    int &out_lion,
		    float &out_time_dump_interval,
		    DateTime &out_startDateTime,
		    DateTime &out_stopDateTime,
		    float &out_ion_exchange_frequency,
		    struct RCMut &innerMagnetosphere,
		    bool &out_use_intercomm);
#endif /* USE_DEPRECATED_NAMELISTS */

/// MSPHERE - LFM Driver
/**
 *  Set up and run the MHD magnetosphere code
 *
 *  If the ionosphere is enabled (in the INPUT1.xml file), use the
 *  simple LFM ionospher solver
 *
 */ 
main(int argc, char **argv)
{
#ifdef DEBUG_MODE_ON
  printf("DEBUG:  process (%d) in MSPHERE.C::main(...)\n", RTIME.local_proc);
#endif

  using namespace std;
  
  int lion;
  char *BASE_IN = (char *)malloc(300);
  char *BASE_OUT = (char *)malloc(300);
  
  // this command outputs "A++ Internal_Index bounds checking: ON" to stdout:
  Index::setBoundsCheck(on);
  
#ifdef PARA_PPP
  int num_proc = NUMBER_OF_PROCESSORS;
  int nmhd     = NUMBER_OF_MHD;
  int ionpe    = 0;
  int proc_first = 0;
  int proc_last  = nmhd-1;
  
  /* The "Optimization_Manager" command outputs alot of information
   * about the P++ Virtual Machine to stdout in the form:
   *
   * *****************************************************
   * P++ Virtual Machine Initialized:
   *   ...
   *   ...
   * *****************************************************
   * This A++/P++ function broadcasts command-line arguments to all the processors and
   * sets up the equivalent of setting up MPI_COMM_WORLD, "MPI_Init", etc.
   */
  Optimization_Manager::Initialize_Virtual_Machine("lfm",num_proc,argc,argv);

  int mype = Communication_Manager::My_Process_Number; 
  RTIME.number_of_processors = num_proc;
  RTIME.local_proc = mype;  

#ifdef USE_MIX
  // Read the XJD file & validate some of its input
  if ((argc < 2) || (argv[1] == NULL) ){
    cerr << "Usage:  " << argv[0] << " <xjd file name>" << endl;
    Communication_Manager::Sync();
    Optimization_Manager::Exit_Virtual_Machine();
  }
  const char *xjdName = argv[1];

  char localName[STRLEN]; // STRLEN is defined in MHDInnerBoundaryInterface.h
  int localTasks;
  
  /* Let's figure out the component name and the anticipated number of processors from the xjd-file */
  parse_xjd(xjdName, "mhd", localName, localTasks);
  cout << "\n********** Parameters parsed from the XJD file ************"<<endl;
  cout << "Local component name: " << localName << endl
       << "Local number of processors: " << localTasks 
       << endl;
  cout << "***********************************************************\n" << endl;
  
  if ( strcmp(localName, "lfm") != 0 ){
    cerr << "mhd component must be named \"lfm\" in the xjd file \"" << xjdName << "\".";
    Communication_Manager::Sync();
    Optimization_Manager::Exit_Virtual_Machine();
  }
    
  if ( localTasks != num_proc ) {
    cout << "Number of processors specified in param.h file is inconsistent with the XJD specification." << endl;
    Communication_Manager::Sync();
    Optimization_Manager::Exit_Virtual_Machine();
  }
#endif //USE_MIX
#endif // PARA_PPP

  /********************************************************************
   *                                                                  *
   *   Read input parameters from INPUT1.xml (or deprecated INPUT1    *
   *   file) and initialize file I/O and timestepping for the MHD     *
   *   and LFM ionosphere                                             *
   *                                                                  *
   *******************************************************************/

  float rdumpInt;
  DateTime startDateTime;
  DateTime stopDateTime;
  float ion_exchange_frequency;
  struct RCMut innerMagnetosphereInput;
  bool exitFlag = false;
  bool use_intercomm = true;

  // this is simpler than doing a broadcast
  for ( int lpe=0; lpe<num_proc; lpe++) {
    Communication_Manager::Sync();   
    if ( mype == lpe ) {
      int ni_bin_default=4; int nj_bin_default=N_J/2;    
#ifdef USE_DEPRECATED_NAMELISTS
      NAMELISTS(BASE_IN,BASE_OUT, &lion, &ni_bin_default, &nj_bin_default);
#else /* USE_DEPRECATED_NAMELISTS */
      
      if ( !xmlInputReader("INPUT1.xml", ni_bin_default, nj_bin_default, BASE_IN, BASE_OUT, lion,
			   rdumpInt, startDateTime, stopDateTime, ion_exchange_frequency, innerMagnetosphereInput, use_intercomm) ){
	std::cerr << "ERROR reading \"INPUT1.xml\"! \n";
	exitFlag = true;
      }
      INITIALIZE();
#endif /* USE_DEPRECATED_NAMELISTS */
    }
  }
  
  if (exitFlag){
    Communication_Manager::Sync();
    Optimization_Manager::Exit_Virtual_Machine();
  }
  
  MHDManager *magnetosphereManager;
#if defined(ION_ON) && !defined(USE_MIX)
  IonIO *ionosphere_IO;
#endif
  
  // DateTimes if we are using UTIO
  DateTime current_ut;

  bool isStandardTimestepIO;

  struct modelCoupling ionCoupling;
  struct modelCoupling innerMagCoupling;

  // Are we using UTIO?

  if ( (NSTART == NSTOP) &&
       (NSTOP == NDUMP) ){
    if (mype == 0){
      std::cout << "---Using UTIO---\n";
    }
    magnetosphereManager = new UTMHD(startDateTime, stopDateTime, rdumpInt);    
#if defined(ION_ON) && !defined(USE_MIX)
    ionosphere_IO = new UTIonIO((UTMHD *) magnetosphereManager);    
#endif

    isStandardTimestepIO = false;
    
    current_ut = ((UTMHD *) magnetosphereManager)->getCurrent();
  
#ifndef DISABLE_RCM
    // Setup Inner Magnetosphere Coupling Information
    innerMagCoupling.start = DateTime(innerMagnetosphereInput.startYear,
				      innerMagnetosphereInput.startMonth,
				      innerMagnetosphereInput.startDay,
				      innerMagnetosphereInput.startHour,
				      innerMagnetosphereInput.startMinute,
				      (double) innerMagnetosphereInput.startSecond);
    innerMagCoupling.delta_t = innerMagnetosphereInput.deltaT;
    innerMagCoupling.exchangeNum = calcExchangeNum(innerMagCoupling.start, current_ut, innerMagCoupling.delta_t);
    innerMagCoupling.next = innerMagCoupling.start;
    innerMagCoupling.next.incrementSeconds( (double) innerMagCoupling.delta_t * innerMagCoupling.exchangeNum );
    innerMagCoupling.isFirstExchange = true;
    innerMagCoupling.isLastExchange = false;
#endif

    // Setup Ionosphere Coupling information
    ionCoupling.start = startDateTime; // first exchange
    ionCoupling.delta_t = ion_exchange_frequency; // number of seconds between exchanges
    ionCoupling.exchangeNum = calcExchangeNum(ionCoupling.start, current_ut, ionCoupling.delta_t);
    ionCoupling.next = ionCoupling.start;     // next time to exchange
    ionCoupling.isLastExchange = false;
  }
  else{
    if (mype == 0){
      std::cout << "---Using standard timestep I/O---\n";
    }
    magnetosphereManager = new TimestepMHD(NSTART, NSTOP, NDUMP);
#if defined(ION_ON) && !defined(USE_MIX)
    ionosphere_IO = new TimestepIonIO((TimestepMHD *) magnetosphereManager);
#endif

    isStandardTimestepIO = true;
    current_ut = DateTime(0);
    ionCoupling.next = DateTime(-1);

#ifndef DISABLE_RCM
    // Why compare against year 9000?  We assuem coupling is turned off when year=9999.
    if ( innerMagCoupling.start < DateTime(9000,0,0, 0,0,0) ){
      cerr << "***** ERROR ***** LFM-RCM coupling requires LFM UTIO, but you have "
	   << " configured the run for standard timestep I/O.  Reconfigure "
	   << " LFM to use UTIO." << endl;
      Communication_Manager::Sync();
      Optimization_Manager::Exit_Virtual_Machine();
    }
#endif
  }
  
  /********************************
   *                              *
   *   Initialize magnetosphere   *
   *                              *
   *******************************/
  
  // constructor for the MHD class
#ifdef PARA_PPP
  Range mhd_procs(0,nmhd-1);
  MHD mhd = MHD(N_I,N_J,N_K,N_ORDER,mhd_procs);
#else
  MHD mhd = MHD(N_I,N_J,N_K,N_ORDER);
#endif
  
  //  initialize the arrays -- use old fortran to get info on how to do this
  
  mhd.initial(nmhd,proc_first,proc_last);
  /*
   *   initial does a number of things:
   *   1. sets up pointers and dimensioning information for the fortran 
   *      arrays used locally
   *      almost all of this done using legacy fortran code
  */
  
  Communication_Manager::Sync();
  
  
#if defined(ION_ON) && !defined(USE_MIX)
  // ionospheric initialization stuff goes here
  
  // constructor for the ION class
  int NION_Shell = 2;
  (N_ORDER/2 - 1) > 2  ? NION_Shell = N_ORDER/2-1 : NION_Shell = 2;
  
  Range ion_procs(ionpe,ionpe);
  ION ion = ION(NION_Shell,N_J,N_K,ion_procs, N_ORDER, NION_I, NION_J);
  
  ion.initial(ionpe);
#endif
  
  Communication_Manager::Sync();

  /****************************************
   *                                      *
   *    Read initial timestep from file   *
   *                                      *
   ***************************************/
  
  magnetosphereManager->read_step(BASE_IN, &mhd);
#if defined(ION_ON) && !defined(USE_MIX)
  ionosphere_IO->read_step(BASE_IN,&ion);
  ION_TEST_ARRAYS();
#endif
  
  Communication_Manager::Sync();
  mhd.startup("unknown.hdf");
  mhd.boundary();
  
#ifdef DIVTEST
  float divmax = DivTest(&mhd,TRUE);
  if ( mype == 0 ){
    std::cout << "maximum divergence =" << divmax << endl;
  }
#endif
  
  Communication_Manager::Sync();
  
#if defined(ION_ON) && !defined(USE_MIX)
  grid2ion(&mhd, &ion);
  ion.startup();
#endif

  /*
   *
   * Initialize InterComm, MIX and RCM
   *
   */
  Communication_Manager::Sync();
#ifdef ION_ON
#ifdef USE_MIX
  Range ioProcessRange(ionpe,ionpe);
  Partitioning_Type ioProcess(ioProcessRange);
  Partitioning_Type INTERFACE_P_NOGUARD(mhd_procs);
  
  int NI = mhd.getNI();
  int NJ = mhd.getNJ();
  int NK = mhd.getNK();

  FloatARRAY *X = mhd.getX();
  FloatARRAY *Y = mhd.getY(); 
  FloatARRAY *Z = mhd.getZ(); 

#if (FloatARRAY == floatArray) 
  doubleArray xDum = (*X).convertTo_doubleArray();
  doubleArray yDum = (*Y).convertTo_doubleArray(); 
  doubleArray zDum = (*Z).convertTo_doubleArray(); 
#elif (FloatARRAY == doubleArray)
  doubleArray xDum = (*X);
  doubleArray yDum = (*Y); 
  doubleArray zDum = (*Z);
#else
  cout << "FloatARRAY must be set to doubleArray or floatArray in site-stuff.h" << endl;
  Optimization_Manager::Exit_Virtual_Machine();
#endif
  xDum.setBase(1);
  yDum.setBase(1);
  zDum.setBase(1);

  xDum.partition(INTERFACE_P_NOGUARD);
  yDum.partition(INTERFACE_P_NOGUARD);
  zDum.partition(INTERFACE_P_NOGUARD);

  int ic_err;
  IC_EndPointSet *epset;
  MHDInnerBoundaryInterface *ionosphere = NULL;
  IM_Interface *innerMagnetosphere = NULL;
  if ( use_intercomm ){
    if ( mype == 0 ) cout << "creating IC_EndPointSet " << xjdName << " " << localName << endl;
    epset = new IC_EndPointSet(xjdName,localName,ic_err);

    // MIX coupler/solver interface
    if ( mype == 0 ) cout << "Allocating coupled model interfaces" << endl;
    ionosphere = new CPL_IC_Interface(epset, xDum, yDum, zDum, 2, NJ, NK);
    
#ifndef DISABLE_RCM
    // RCM inner magnetosphere interface

    // Get the Dipole moment (required by IM_Interface)
    float dipole_moment; 
    GET_DIP_MOMENT(&dipole_moment);
    // Convert from Gauss cm^3 to MKS (T s m^3)
    dipole_moment *= 1.0e-10;

    innerMagnetosphere = new IM_IC_Interface(epset, 
					     xDum.convertTo_floatArray(), yDum.convertTo_floatArray(), zDum.convertTo_floatArray(),
					     X->getPartition(),
					     dipole_moment, GAMMA, 
					     innerMagCoupling.start, innerMagCoupling.delta_t);
#else
    innerMagnetosphere = new IM_IC_Interface(NULL, 
					     xDum.convertTo_floatArray(), yDum.convertTo_floatArray(), zDum.convertTo_floatArray(),
					     X->getPartition(),
					     dipole_moment, GAMMA, 
					     DateTime(9999,0,0,0,0,0), 1.0);
#endif
    
    // We commit the InterComm arrays here because this function
    // commits *all* regions created previously.  In our case, both
    // CPL_IC_Interface and IM_IC_Interface define arrays that need to
    // be committed.
    epset->commitArrays(ic_err);
    epset->printErrorMessage("     Arrays commited. Status",ic_err);
  }
  else {
    ionosphere = new CPL_FE_Interface(xDum, yDum, zDum, 2, NJ, NK);

#ifndef DISABLE_RCM
    if ( innerMagCoupling.start > DateTime(-9999,0,0, 0,0,0) ){
      cerr << "***** ERROR ***** You have requested RCM coupling, which requires the InterComm interface, but you have selected file exchange!  Exiting." << endl;
      Communication_Manager::Sync();
      Optimization_Manager::Exit_Virtual_Machine();
    }
#endif    
  }

  Partitioning_Type p_zero(Range(0,0));

#ifndef DISABLE_RCM
    innerMagnetosphere->Initialize();

    intArray imScalars(IM_Interface::NUMBER_OF_SCALARS, p_zero);
    imScalars(IM_Interface::KILL_SIGNAL) = IM_Interface::KILL_SIGNAL_CONTINUE;
    imScalars(IM_Interface::DELTA_T) = innerMagCoupling.delta_t; // number of seconds between LFM-RCM exchange
#endif

  doubleArray cplScalars(CPL_Interface::NUMBER_OF_SCALARS, p_zero); 
  cplScalars *= 0; // initialize to zero.

  Range all, I(1,NI), J(1,NJ), Jp1(1,NJ+1),  K(1,NK), Kp1(1,NK+1), TWO(1,2);

  FloatARRAY *Bx = mhd.getBx();
  FloatARRAY *By = mhd.getBy(); 
  FloatARRAY *Bz = mhd.getBz(); 
  FloatARRAY *Rho = mhd.getRho(); 
  FloatARRAY *Csound = mhd.getCsound(); 

  FloatARRAY *Ej = mhd.getEj();
  FloatARRAY *Ek = mhd.getEk();
  FloatARRAY *Vx = mhd.getVx();
  FloatARRAY *Vy = mhd.getVy();
  FloatARRAY *Vz = mhd.getVz();
  
  Communication_Manager::Sync();
#else
  mhd2ion(&mhd,&ion);
  if ( mype == ionpe ){
    ion.step(); 
  }
  Communication_Manager::Sync();
  //  ion2mhd(&mhd,&ion); 
#endif //USE_MIX
#endif //ION_ON

  /// Average Average Average ///
  // AKA low-pass filter
  vector <int> variables;
  variables.insert(variables.end(), AverageVariables::BX);
  variables.insert(variables.end(), AverageVariables::BY);
  variables.insert(variables.end(), AverageVariables::BZ);
  variables.insert(variables.end(), AverageVariables::EI);
  variables.insert(variables.end(), AverageVariables::EJ);
  variables.insert(variables.end(), AverageVariables::EK);
  FloatARRAY *tmpX = mhd.getX();
  Average<FloatARRAY> mhdAveraging(tmpX->getLength(0), tmpX->getLength(1), tmpX->getLength(2),
				   tmpX->getPartition(),
				   variables);

#ifndef DISABLE_RCM
  /*
   * If we are coupled to the inner magnetosphere and this is a
   * restart run, read the accumulate & bleed arrays into the
   * innerMagnetosphere IM_Interface object.  This code must go after
   * initializing both the Coupler (MIX) & Inner Magnetosphere (RCM)
   * interfaces and before dumping the initial timestep.
   */
  
  DateTime accumulate_start_ut = innerMagCoupling.start;
  accumulate_start_ut.incrementSeconds(-innerMagCoupling.delta_t);
  if (current_ut > accumulate_start_ut){
    magnetosphereManager->read_accumulate_and_bleed_step(BASE_OUT, &mhd, innerMagnetosphere);
  }
#endif
  
#ifdef MSPHERE_DUMP_FIRST_STEP
  /*
   * Dump the initial time step.  Note: This is necessary after
   * loading the initial time step because
   *   - magnetic field information is set in MSPHERE
   *   - we read in the inner magnetosphere accumulate & bleed arrays
   */
#ifndef DISABLE_RCM
  magnetosphereManager->dump_step(BASE_OUT,&mhd, innerMagnetosphere, mhdAveraging);
#else
  magnetosphereManager->dump_step(BASE_OUT,&mhd, NULL, mhdAveraging);
#endif
#if defined(ION_ON) && !defined(USE_MIX)
  ionosphere_IO->dump_step(BASE_OUT,&ion);
#endif
#endif

  
  /*******************************************************************
   *                                                                 *
   *   ----  Main loop  ------------------------------------------   *
   *                                                                 *
   ******************************************************************/

  // set performance timing information
  double sim_time_0 = TIME;
  double sim_step_0 = magnetosphereManager->currentStep;
  double time_0 = MPI_Wtime();
  double time_i_0 = time_0;
  double time_i_1 = time_0;
  double avg_realDT = 0.0;
  double avg_simDT = 0.0;
  int nAvg = -1;
  
  while (magnetosphereManager->advance() ) {
    time_i_0 = time_i_1;
    time_i_1 = MPI_Wtime();
    
    nAvg++;
    avg_realDT += (double) (time_i_1-time_i_0);
    avg_simDT += DT;
    
    // Output timestep information on head node:
    if ( mype == 0 && (magnetosphereManager->currentStep % 10) == 0 ) {
      std::cout << "Loop # " << magnetosphereManager->currentStep << ";  "
		<< fixed 
		<< "real Time: " << setprecision(2) << (time_i_1 - time_0) << ";  "
		<< "real dt: " << (time_i_1-time_i_0) << ";  "
		<< "sim Time: " << TIME << ";  "
		<< "sim DT: " << setprecision(8) << DT
		<< "\n";
    }
    
    ///////////////////////////////////////////////////////
    // Update the magnetosphere Solar Wind boundary condition
    ///////////////////////////////////////////////////////
    
    /*   	   
	   This is where the SPMD P++ model bites us. If we send the
	   ion calculation down a separate track, none of the P++
	   arrays sync properly because everybody is using
	   MPI_COMM_WORLD.
	   
	   Therefore, we have to bite the bullet and have the
	   ionosphere take as long as it takes, potentially holding up
	   the MHD. We gain a little bit of parallelism by doing the
	   MHD boundary first (MPI calls), then the ionosphere stuff,
	   then the MHD step. The latter two do not use MPI, so the
	   stuff is basically done in parallel.
    */

    mhd.boundary();
    
#ifdef DIVTEST
    float divmax = DivTest(&mhd,TRUE);
    if ( mype == 0 ) {
      std::cout << "maximum divergence =" << divmax << endl;
    }
#endif

    ///////////////////////////////////////////////////////
    // Increment Timestep DT
    ///////////////////////////////////////////////////////
    
    // How big of a timestep are we going to take?
    DT =  mhd.courant();
    
    //LSTEP = lstep;
    LSTEP = magnetosphereManager->currentStep;
    
    // Increment simulation time (TIME = TIME + DT)
    magnetosphereManager->incrementTime(DT);
    TIME = magnetosphereManager->simTime;
    
    if (isStandardTimestepIO == false){
      current_ut = ((UTMHD *) magnetosphereManager)->getCurrent();
    }
    
    Communication_Manager::Sync();
    
    LOAD_B();
    
    ///////////////////////////////////////////////////////
    // Advance the Magnetosphere simulation DT seconds
    ///////////////////////////////////////////////////////
    mhd.step(&DT);
    
    Communication_Manager::Sync();

#if defined(USE_MIX) && !defined(DISABLE_RCM)
    // Accumulate LFM data into Inner Mag accumulate arrays if it's
    // time.  Note: We accummulate IMMEDIATELY after time stepping the
    // MHD.  This will populate the Accumulate arrays on the first
    // iteration of the LFM time loop.
    DateTime accumulate_start = innerMagCoupling.start;
    accumulate_start.incrementSeconds(-innerMagCoupling.delta_t);
    if (current_ut > accumulate_start){
      innerMagnetosphere->Accumulate(DT,
				     *(mhd.getRho()), *(mhd.getCsound()),
				     *(mhd.getBi()), *(mhd.getBj()), *(mhd.getBk()),
				     *(mhd.getVx()), *(mhd.getVy()), *(mhd.getVz()));
    }
#endif

#ifdef ION_ON      
    ///////////////////////////////////////////////////////
    // Advance the Ionosphere (if it's time to)
    ///////////////////////////////////////////////////////
    // NLOOP is specified in run-time.inc and is specified in the INPUT1.xml file
    // NLOOP determines how often the ionosphere & MHD communicate.
    if ( 
	( ( ( magnetosphereManager->currentStep % NLOOP ) == 0) && (isStandardTimestepIO) ) ||
	( ( current_ut >= ionCoupling.next ) && (isStandardTimestepIO == false) )
	 ){
      ionCoupling.exchangeNum++;
#ifdef USE_MIX 

#ifndef DISABLE_RCM
      if (current_ut >= innerMagCoupling.next){
	///////////////////////////////////////////////////////
	// Import data from Inner Magnetosphere (RCM) into MHD:
	//   1.  Density on LFM grid
	//   2.  Pressure on LFM grid
	//   3.  Mask: 0==not part of RCM;  1==part of RCM 
	///////////////////////////////////////////////////////
	if (mype == 0)  std::cout << "IM exchange # " << innerMagCoupling.exchangeNum << "\t"
				  << current_ut.getDateTimeString() << std::endl;

	if (innerMagCoupling.isFirstExchange){
	  // Tell the coupler (MIX) that it needs to export only to the RCM
	  cplScalars(CPL_Interface::EXCHANGE_WITH_IM) = CPL_Interface::CPL_SEND_TO_IM;
	}
	else if ( isLastExchange(stopDateTime, innerMagCoupling.delta_t, innerMagCoupling.next) ){
	  innerMagCoupling.isLastExchange = true;
	  // Update kill signal.
	  imScalars(IM_Interface::KILL_SIGNAL) = IM_Interface::KILL_SIGNAL_LAST_EXCHANGE;
	  // Last exchange.  Tell the Coupler (MIX) that it needs import only from RCM
	  cplScalars(CPL_Interface::EXCHANGE_WITH_IM) = CPL_Interface::CPL_RECV_FROM_IM;	
	}
	else{
	  cplScalars(CPL_Interface::EXCHANGE_WITH_IM) = CPL_Interface::DUPLEX_IM_EXCHANGE;
	}

	imScalars(IM_Interface::YEAR)   = current_ut.getYear();
	imScalars(IM_Interface::MONTH)  = current_ut.getMonth();
	imScalars(IM_Interface::DAY)    = current_ut.getDay();
	imScalars(IM_Interface::HOUR)   = current_ut.getHour();
	imScalars(IM_Interface::MINUTE) = current_ut.getMinute();
	imScalars(IM_Interface::SECOND) = int( current_ut.getSecond() );
	imScalars(IM_Interface::LABEL)  = magnetosphereManager->currentStep;

	innerMagnetosphere->SendScalars("mhd-im_scalars", imScalars);
	if (!innerMagCoupling.isFirstExchange){
	  innerMagnetosphere->Import();
	}
      }
#endif //DISABLE_RCM      

      ///////////////////////////////////////////////////////
      // Export MHD (J_para, rho, c_s) data to MIX
      ///////////////////////////////////////////////////////
      if ( mype == 0){
	std::cout << "LFM-MIX exchange # " << ionCoupling.exchangeNum << "; ";
	if (isStandardTimestepIO)
	  std::cout << "MHD Timestep: " << magnetosphereManager->currentStep << "\n";
	else
	  std::cout << "MHD UT: "
		    << current_ut.getDateTimeString() << endl;
      }
      
#if (FloatARRAY == floatArray) 
      doubleArray bxDum = (*Bx).convertTo_doubleArray();
      doubleArray byDum = (*By).convertTo_doubleArray(); 
      doubleArray bzDum = (*Bz).convertTo_doubleArray(); 
      doubleArray rhoDum = (*Rho).convertTo_doubleArray(); 
      doubleArray csoundDum = (*Csound).convertTo_doubleArray(); 
#elif (FloatARRAY == doubleArray) 
      doubleArray bxDum = (*Bx);
      doubleArray byDum = (*By);
      doubleArray bzDum = (*Bz);
      doubleArray rhoDum = (*Rho);
      doubleArray csoundDum = (*Csound);
#else
#error "FloatARRAY must be set to doubleArray or floatArray in site-stuff.h"
#endif
      bxDum.setBase(1);      bxDum.partition(INTERFACE_P_NOGUARD);
      byDum.setBase(1);      byDum.partition(INTERFACE_P_NOGUARD);
      bzDum.setBase(1);      bzDum.partition(INTERFACE_P_NOGUARD);
      rhoDum.setBase(1);     rhoDum.partition(INTERFACE_P_NOGUARD);
      csoundDum.setBase(1);  csoundDum.partition(INTERFACE_P_NOGUARD);     
      
      cplScalars(CPL_Interface::KILL_SIGNAL) = CPL_Interface::KILL_SIGNAL_CONTINUE;
      cplScalars(CPL_Interface::LABEL) = magnetosphereManager->currentStep; 
      cplScalars(CPL_Interface::TILT) = getTilt();
      
      // Send UT info for the next MHD step
      // i.e. the UT after the ionosphere takes a step.
      if (isStandardTimestepIO){      
	// Send the "second-to-last" signal if it's time.
	if (LSTEP + 1 > NSTOP){
	  cplScalars(CPL_Interface::KILL_SIGNAL) = CPL_Interface::KILL_SIGNAL_LAST_EXCHANGE;
	  ionCoupling.isLastExchange = true;
	}
  
	// Give MIX a bogus invalid year since we're using timesteps.
	cplScalars(CPL_Interface::YEAR) = -9999.;  
      }
      else{
	// UTIO.  Is this the last exchange?
	if (isLastExchange(stopDateTime, ionCoupling.delta_t, ionCoupling.next)){
	  cplScalars(CPL_Interface::KILL_SIGNAL) = CPL_Interface::KILL_SIGNAL_LAST_EXCHANGE;
	  ionCoupling.isLastExchange = true;
	}
	cplScalars(CPL_Interface::YEAR) =   current_ut.getYear();   // Year
	cplScalars(CPL_Interface::MONTH) =  current_ut.getMonth();  // Month
	cplScalars(CPL_Interface::DAY) =    current_ut.getDay();    // Day
	cplScalars(CPL_Interface::HOUR) =   current_ut.getHours();  // Hour
	cplScalars(CPL_Interface::MINUTE) = current_ut.getMinute(); // Min
	cplScalars(CPL_Interface::SECOND) = current_ut.getSecond(); // Sec
      }
      
      ionosphere->sendScalars(cplScalars);

      ionosphere->Export( bxDum(I,J,K), // ionosphere needs all magnetic field
			  byDum(I,J,K), // and only first shells of rho and csound
			  bzDum(I,J,K), 
			  rhoDum(TWO,J,K), 
			  csoundDum(TWO,J,K) );	

      ///////////////////////////////////////////////////////
      // Import (Potential) data from MIX
      ///////////////////////////////////////////////////////
      // Note: MIX dumps its data files at an interval specified in
      //       the MIX.param file after the Import(...)
      ///////////////////////////////////////////////////////
      doubleArray ej_pot, ek_pot, velocity;
      ionosphere->Import(ej_pot, ek_pot, velocity);

      // ej and ek are partitioned in the following way within
      // MHDInnerBoundaryInterface.C: 
      // eField_j.redim(2,nj,nkp1); eField_k.redim(2,njp1,nk);	
#if (FloatARRAY == floatArray) 
      floatArray ejDum(*Ej);
      ejDum.partition(INTERFACE_P_NOGUARD);
      ejDum(1,J,Kp1) = ej_pot(1,all,all).convertTo_floatArray();
      
      floatArray ekDum(*Ek);
      ekDum.partition(INTERFACE_P_NOGUARD);
      ekDum(1,Jp1,K) = ek_pot(1,all,all).convertTo_floatArray(); 
      
      floatArray vxDum(*Vx);
      vxDum.partition(INTERFACE_P_NOGUARD);
      vxDum(1,J,K) = velocity(1,all,all,1).convertTo_floatArray(); 
      
      floatArray vyDum(*Vy);
      vyDum.partition(INTERFACE_P_NOGUARD);
      vyDum(1,J,K) = velocity(1,all,all,2).convertTo_floatArray(); 
      
      floatArray vzDum(*Vz);
      vzDum.partition(INTERFACE_P_NOGUARD);
      vzDum(1,J,K) = velocity(1,all,all,3).convertTo_floatArray(); 
      
#elif (FloatARRAY == doubleArray)       
      doubleArray ejDum(*Ej);
      ejDum.partition(INTERFACE_P_NOGUARD); 
      ejDum(1,J,Kp1) = ej_pot(1,all,all);
      
      doubleArray ekDum(*Ek);
      ekDum.partition(INTERFACE_P_NOGUARD);
      ekDum(1,Jp1,K) = ek_pot(1,all,all); 
      
      doubleArray vxDum(*Vx);
      vxDum.partition(INTERFACE_P_NOGUARD);
      vxDum(1,J,K) = velocity(1,all,all,1); 
      
      doubleArray vyDum(*Vy);
      vyDum.partition(INTERFACE_P_NOGUARD);
      vyDum(1,J,K) = velocity(1,all,all,2); 
      
      doubleArray vzDum(*Vz);
      vzDum.partition(INTERFACE_P_NOGUARD);
      vzDum(1,J,K) = velocity(1,all,all,3); 
#else
#error "FloatARRAY must be set to doubleArray or floatArray in site-stuff.h"
#endif
	
      ejDum.partition(*Ej);
      ekDum.partition(*Ek);
      vxDum.partition(*Vx);
      vyDum.partition(*Vy);
      vzDum.partition(*Vz);
      
      (*Ej)(1,J,Kp1) = ejDum(1,J,Kp1);
      (*Ek)(1,Jp1,K) = ekDum(1,Jp1,K);
      (*Vx)(1,J,K) = vxDum(1,J,K);
      (*Vy)(1,J,K) = vyDum(1,J,K);
      (*Vz)(1,J,K) = vzDum(1,J,K);
      
      if ( (current_ut >= ionCoupling.next) && (isStandardTimestepIO == false) ){
	// prevent floating point error from creeping in by specifying the next dump step as
	// start-time + dt * (n+1)
	ionCoupling.next = ionCoupling.start;
	ionCoupling.next.incrementSeconds((double) ionCoupling.delta_t * ionCoupling.exchangeNum );
      }

      // Assume that MIX Coupler and RCM Inner-Magnetosphere do not
      // exchange at the next LFM step.
      cplScalars(CPL_Interface::EXCHANGE_WITH_IM) = CPL_Interface::NO_IM_EXCHANGE;

#ifndef DISABLE_RCM
      ///////////////////////////////////////////////////////
      // Export data from LFM into Inner Magnetosphere (RCM)
      ///////////////////////////////////////////////////////
      
      // Is this an Inner Mag (RCM) coupling exchange step?
      if (current_ut >= innerMagCoupling.next){
	// Do not export to RCM on the last exchange...
	if (! innerMagCoupling.isLastExchange){
	  innerMagnetosphere->Export();     
	}
      }
#endif //DISABLE_RCM
#else
      // (i.e. not using MIX)
      ///////////////////////////////////////////////////////
      // Advance LFM ionosphere simulation:
      ///////////////////////////////////////////////////////
      mhd2ion(&mhd,&ion);
      if ( mype == ionpe ){
	ion.step(); 
      }	
      ion2mhd(&mhd,&ion); 
#endif
      // If we've gotten this far, we've exchanged at least once.
    } // end of lstep % NLOOP if
#endif // ION_ON    

    Communication_Manager::Sync();
    EXAM_B();

#ifndef DISABLE_RCM
    // Average RCM data (pressure, density) into the LFM when it's time.
    DateTime firstBleedTime = innerMagCoupling.start;
    firstBleedTime.incrementSeconds(innerMagCoupling.delta_t);
    if ( current_ut > firstBleedTime) { 
      innerMagnetosphere->Bleed(innerMagnetosphereInput.bleed_rate, innerMagnetosphereInput.RCM_bleed_type_scheme,
				DT, *(mhd.getRho()), *(mhd.getCsound()),
				*(mhd.getVx()), *(mhd.getVy()), *(mhd.getVz()) );
    }
#endif

    /// Average Average Average ///
    // AKA low-pass filter
    mhdAveraging.setDeltaTime(DT);
    mhdAveraging.increment(*(mhd.getBx()), AverageVariables::BX);
    mhdAveraging.increment(*(mhd.getBy()), AverageVariables::BY);
    mhdAveraging.increment(*(mhd.getBz()), AverageVariables::BZ);
    mhdAveraging.increment(*(mhd.getEi()), AverageVariables::EI);
    mhdAveraging.increment(*(mhd.getEj()), AverageVariables::EJ);
    mhdAveraging.increment(*(mhd.getEk()), AverageVariables::EK);
    
    ///////////////////////////////////////////////////////
    // LFM mhd (and ion) I/O, if necessary
    ///////////////////////////////////////////////////////
    if ( magnetosphereManager->atDump() ){
#ifndef DISABLE_RCM
      magnetosphereManager->dump_step(BASE_OUT,&mhd, innerMagnetosphere, mhdAveraging);
#else
      magnetosphereManager->dump_step(BASE_OUT,&mhd, NULL, mhdAveraging);
#endif
      /// Average Average Average ///
      // AKA low-pass filter
      mhdAveraging.reset();
#if defined (ION_ON) && !defined(USE_MIX)
      ionosphere_IO->dump_step(BASE_OUT,&ion);
#endif
    }

#ifndef DISABLE_RCM
    // When it's time, do some Inner Magnetosphere book keeping:
    //
    // 1.  Reset accumulate arrays.  Note:  This must happen _after_ file I/O!!!
    // 2.  Update coupling time.  Note:  This must happen _after_ resetting accum arrays
    if (current_ut >= innerMagCoupling.next) {
      if (! innerMagCoupling.isLastExchange){ 
	innerMagnetosphere->resetAccumulateArrays();
      }
      // Update the time for next exchange.
      innerMagCoupling.next = innerMagCoupling.start;
      innerMagCoupling.next.incrementSeconds(++innerMagCoupling.exchangeNum*innerMagCoupling.delta_t);
      innerMagCoupling.isFirstExchange = false;      
    }
#endif
  }  //  END MAIN LOOP ---------------------------------
  
  if ( mype == 0 ){
    std::cout << "Average wallclock time for a timestep: " << avg_realDT/(double)nAvg << "\n"
	      << "Average simulatime time for a timestep: " << avg_simDT/(double)nAvg << "\n"
	      << "Total wallclock time: " <<  (time_i_1 - time_0)<< "\n"
	      << "Total simulation time: " << TIME - sim_time_0 << "\n"
	      << "for " << (int) magnetosphereManager->currentStep - sim_step_0 << " steps \n";
  }
  
#ifdef USE_MIX
  cplScalars(CPL_Interface::KILL_SIGNAL) = CPL_Interface::KILL_SIGNAL_SHUTDOWN;
  ionosphere->sendScalars(cplScalars);

#ifndef DISABLE_RCM
  imScalars(IM_Interface::KILL_SIGNAL) = IM_Interface::KILL_SIGNAL_SHUTDOWN;
  innerMagnetosphere->SendScalars("mhd-im_scalars", imScalars);
  
  delete innerMagnetosphere;
  innerMagnetosphere = NULL;
#endif

  delete ionosphere;
  ionosphere = NULL;

  if (use_intercomm){
    delete epset;
    epset = NULL;
  }
#endif
  
#ifdef PARA_PPP
  Communication_Manager::Sync();
  Optimization_Manager::Exit_Virtual_Machine();
#endif  
}



/********************************************************************//**
 * Calculate how many exchanges should have occured from start
 * date/time to current date/time, assuming constant exchange interval
 * delta_t
 ************************************************************************/ 								       
int calcExchangeNum(const DateTime &start, const DateTime &current, const int &delta_t)
{
  if (start < current){
    // Note: We _must_ spread this calculation out among some local
    // variables to prevent floating-point arithmetic underflow.
    double numerator = current.getMJD() - start.getMJD();
    double denominator = (double) delta_t / 86400.0;
    double iteration_number = numerator/denominator;

    // Doing the above arithmetic as a one-liner in the return
    // statement using the following Date/Time objects incorrectly
    // returns iteration 0:
    //    start = DateTime(1995,3,21,  4,0,0)
    //    current = DateTime(1995,3,21, 4,1,0)
    //    delta_t = 60.0 seconds
    return floor( iteration_number );
  }
  else{
    return 0;
  }
}


/********************************************************************//**
 * Is this the final exchange that should occur (i.e. is current +
 * delta_t > stop)?
 ************************************************************************/ 								       
bool isLastExchange(const DateTime &stop, const int &delta_t, const DateTime &current)
{
  DateTime nextExchange = current;
  nextExchange.incrementSeconds((double) delta_t);
  if (nextExchange > stop){
    return true;
  }
  return false;
}


#ifdef USE_MIX
/********************************************************************//**
 * MIX requires the tilt angle at every exchange.
 ************************************************************************/ 								       
float getTilt()
{
#ifdef DEBUG_MODE_ON
  printf("DEBUG:  process (%d) in MSPHERE.C::getTilt(...)\n", RTIME.local_proc);
#endif
  float tilt;
  GET_TILT(&tilt);
  return tilt;
}

#endif