#include "site-stuff.h"

#include <stdio.h>

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

#include <A++.h>
#include "MHD.h"

#if LC_UNSC
#define  BZZ     bzz_
#define  IBOUND  ibound_
#define  JBOUND  jbound_
#define  KBOUND  kbound_
#define  EDGE_FIX      edge_fix_
#define  IBOUND_FOR  ibound_for_
#define  JBOUND_FOR  jbound_for_
#define  KBOUND_FOR  kbound_for_
#elif LC_NOUNSC
#define  BZZ     bzz
#define  IBOUND  ibound
#define  JBOUND  jbound
#define  KBOUND  kbound
#define  EDGE_FIX      edge_fix
#define  IBOUND_FOR  ibound_for
#define  JBOUND_FOR  jbound_for
#define  KBOUND_FOR  kbound_for
#elif UC_UNSC
#define  BZZ     BZZ_
#define  IBOUND  IBOUND_
#define  JBOUND  JBOUND_
#define  KBOUND  KBOUND_
#define  EDGE_FIX      EDGE_FIX_
#define  IBOUND_FOR  IBOUND_FOR_
#define  JBOUND_FOR  JBOUND_FOR_
#define  KBOUND_FOR  KBOUND_FOR
#else
#define  BZZ     BZZ
#define  IBOUND  IBOUND
#define  JBOUND  JBOUND
#define  KBOUND  KBOUND
#define  EDGE_FIX      EDGE_FIX
#define  IBOUND_FOR  IBOUND_FOR
#define  JBOUND_FOR  JBOUND_FOR
#define  KBOUND_FOR  KBOUND_FOR

#endif

/// Defined in ionosphere.F
extern "C" void BZZ(FLoaT *, FLoaT *, FLoaT *, FLoaT *, FLoaT *, FLoaT *);
/// Defined in bound.F
extern "C" void IBOUND(FLoaT *, FLoaT*, int *, int *, int *);
/// Defined in bound.F
extern "C" void JBOUND(FLoaT *, FLoaT*, int *, int *, int *);
/// Defined in bound.F
extern "C" void KBOUND(FLoaT *, FLoaT*, int *, int *, int *);
/// Defined in test-arrays.F
extern "C" void TEST_ARRAYS();
/// Defined in second.c
extern "C" float dtime_(FLoaT *);
/// Defined in bound.F
extern "C" void EDGE_FIX();
/// Defined in boundaries.F
extern "C" void IBOUND_FOR();
/// Defined in boundaries.F
extern "C" void JBOUND_FOR();
/// Defined in boundaries.F
extern "C" void KBOUND_FOR();

/**********
 * FIXME: *
 **********
 * These values should go in a header or something.... not here! */
//     dummy solar wind values 
#define Rho_SW  1.
#define Vx_SW  0.
#define Vy_SW  0.
#define Vz_SW  0.
#define Csound_SW  1.291
#define Pressure_SW  1.
#define Bx_SW  0.
#define By_SW  0.
#define Bz_SW  0.


/// Boundary
/**
 * Update the boundary conditions.
 */
void MHD::boundary()
{
#ifdef DEBUG_MODE_ON
  printf("DEBUG:  process (%d) in Boundary.C::MHD::boundary()\n", RTIME.local_proc);
#endif  
  FLoaT *bx, *by, *bz, *bi, *bj, *bk;
  
  //
  //     get cell-centered b from fluxes
  //
  
  bx = Bx.getDataPointer();
  by = By.getDataPointer();
  bz = Bz.getDataPointer();
  bi = Bi.getDataPointer();
  bj = Bj.getDataPointer();
  bk = Bk.getDataPointer();
  BZZ(bx,by,bz,bi,bj,bk);
  
#ifdef DEBUG
  int mype = Communication_Manager::My_Process_Number;
  printf("^^^^^^^^^^ in boundary %d\n",mype); fflush(stdout);
#endif

  // 
  //      get internal boundaries set
  //

  Rho.updateGhostBoundaries();
  Vx.updateGhostBoundaries();
  Vy.updateGhostBoundaries();
  Vz.updateGhostBoundaries();
  Csound.updateGhostBoundaries();
  Pressure.updateGhostBoundaries();
  Bx.updateGhostBoundaries();
  By.updateGhostBoundaries();
  Bz.updateGhostBoundaries();
  Bi.updateGhostBoundaries();
  Bj.updateGhostBoundaries();
  Bk.updateGhostBoundaries();

#ifdef PPLUSPLUS   
  int uno=1; 
  int zip = 0;

  // boundaries done using P++
  
  Range I(1,NI+1); // Fortran indexing (1 - NI) active
  Range J(1,NJ+1); 
  Range K(1,NK+1);
  // FloatARRAY BB(I,J,K); 
  FloatARRAY BB(NI+1,NJ+1,NK+1); BB.setBase(1);
  BB.partition(MHD_P);
  /* 
     BB is a dummy boundary array, which we load with the boundary
     guard cell values.  The guard cells are hidden from the C++
     layer and the external bounaries are manipulated by using fortran
     on the data segments of the P++ arrays. BB is laid out
     identically with the variable arrays, so that fortran can deal
     only with local data. BB gets the guard cell data in the first
     NO2 cells of the visible P++ array.  (note, this means the
     variable arrays must be partitioned at least NO2 deep.)
     Transfers from other processors are done by using P++ assignments
     from variable arrays to BB.
  */ 
	   
  // Do the I boundaries first, the rest are symmetry type

  float timer;
  
#ifdef DEBUG
  printf("in boundary before BB %d %f\n", mype, timer,);
#endif

  //  solar wind boundary
  Range I_1(1,NO2);
  Range I_2(NI-NO2+1,NI);
  
  BB(I_2,J,K) = Rho_SW; BB(I_1,J,K) = Rho(I_1,J,K);

#ifdef DEBUG
  TEST_ARRAYS();
#endif

  IBOUND(BB.getDataPointer(),Rho.getDataPointer(),&zip,&zip,&zip);
  BB(I_2,J,K)  = Vx_SW; 
  BB(I_1,J,K) = Vx(I_1,J,K);
  IBOUND(BB.getDataPointer(),Vx.getDataPointer(),&zip,&zip,&zip);

  BB(I_2,J,K)  = Vy_SW; 
  BB(I_1,J,K) = Vy(I_1,J,K);
  IBOUND(BB.getDataPointer(),Vy.getDataPointer(),&zip,&zip,&zip);
  
  BB(I_2,J,K)  = Vz_SW; 
  BB(I_1,J,K) = Vz(I_1,J,K);
  IBOUND(BB.getDataPointer(),Vz.getDataPointer(),&zip,&zip,&zip);
  
  BB(I_2,J,K) = Csound_SW; 
  BB(I_1,J,K) = Csound(I_1,J,K);
  IBOUND(BB.getDataPointer(),Csound.getDataPointer(),&zip,&zip,&zip);
  
  BB(I_2,J,K) = Pressure_SW; 
  BB(I_1,J,K) = Pressure(I_1,J,K);
  IBOUND(BB.getDataPointer(),Pressure.getDataPointer(),&zip,&zip,&zip);
  
  BB(I_2,J,K)  = Bx_SW; 
  BB(I_1,J,K) = Bx(I_1,J,K);
  IBOUND(BB.getDataPointer(),Bx.getDataPointer(),&zip,&zip,&zip);
  
  BB(I_2,J,K)  = By_SW; 
  BB(I_1,J,K) = By(I_1,J,K);
  IBOUND(BB.getDataPointer(),By.getDataPointer(),&zip,&zip,&zip);
  
  BB(I_2,J,K)  = Bz_SW; 
  BB(I_1,J,K) = Bz(I_1,J,K);
  IBOUND(BB.getDataPointer(),Bz.getDataPointer(),&zip,&zip,&zip);
  
  BB(I_2,J,K)  = Bi(I_2+1,J,K); 
  BB(I_1,J,K) = Bi(I_1,J,K);
  IBOUND(BB.getDataPointer(),Bi.getDataPointer(),&uno,&zip,&zip);
  
  BB(I_2,J,K)  = Bj(I_2,J,K); 
  BB(I_1,J,K) = Bj(I_1,J,K);
  IBOUND(BB.getDataPointer(),Bj.getDataPointer(),&zip,&uno,&zip);
  
  BB(I_2,J,K)  = Bk(I_2,J,K); 
  BB(I_1,J,K) = Bk(I_1,J,K);
  IBOUND(BB.getDataPointer(),Bk.getDataPointer(),&zip,&zip,&uno);
  
#ifdef DEBUG
  printf("^^^^^^^^^^ after last ibound %d\n",mype); fflush(stdout);
  printf("in boundary after i %d %f\n", mype, timer);
#endif

  //   Do the J boundaries
  Range Kl(1,NK/2);
  Range Kh(NK/2+1,NK);
  Range J_1(1,NO2);
  Range J_2(NJ-NO2+1,NJ);
  
  BB(I,J_1,Kl) = Rho(I,J_1,Kh); BB(I,J_1,Kh) = Rho(I,J_1,Kl);
  BB(I,J_2,Kl) = Rho(I,J_2,Kh); BB(I,J_2,Kh) = Rho(I,J_2,Kl);
  JBOUND(BB.getDataPointer(),Rho.getDataPointer(),&zip,&zip,&zip);
  
  BB(I,J_1,Kl) = Vx(I,J_1,Kh); BB(I,J_1,Kh) = Vx(I,J_1,Kl);
  BB(I,J_2,Kl) = Vx(I,J_2,Kh); BB(I,J_2,Kh) = Vx(I,J_2,Kl);
  JBOUND(BB.getDataPointer(),Vx.getDataPointer(),&zip,&zip,&zip);
 
  BB(I,J_1,Kl) = Vy(I,J_1,Kh); BB(I,J_1,Kh) = Vy(I,J_1,Kl);
  BB(I,J_2,Kl) = Vy(I,J_2,Kh); BB(I,J_2,Kh) = Vy(I,J_2,Kl);
  JBOUND(BB.getDataPointer(),Vy.getDataPointer(),&zip,&zip,&zip);
  
  BB(I,J_1,Kl) = Vz(I,J_1,Kh); BB(I,J_1,Kh) = Vz(I,J_1,Kl);
  BB(I,J_2,Kl) = Vz(I,J_2,Kh); BB(I,J_2,Kh) = Vz(I,J_2,Kl);
  JBOUND(BB.getDataPointer(),Vz.getDataPointer(),&zip,&zip,&zip);
  
  BB(I,J_1,Kl) = Bx(I,J_1,Kh); BB(I,J_1,Kh) = Bx(I,J_1,Kl);
  BB(I,J_2,Kl) = Bx(I,J_2,Kh); BB(I,J_2,Kh) = Bx(I,J_2,Kl);
  JBOUND(BB.getDataPointer(),Bx.getDataPointer(),&zip,&zip,&zip);
  
  BB(I,J_1,Kl) = By(I,J_1,Kh); BB(I,J_1,Kh) = By(I,J_1,Kl);
  BB(I,J_2,Kl) = By(I,J_2,Kh); BB(I,J_2,Kh) = By(I,J_2,Kl);
  JBOUND(BB.getDataPointer(),By.getDataPointer(),&zip,&zip,&zip);
  
  BB(I,J_1,Kl) = Bz(I,J_1,Kh); BB(I,J_1,Kh) = Bz(I,J_1,Kl);
  BB(I,J_2,Kl) = Bz(I,J_2,Kh); BB(I,J_2,Kh) = Bz(I,J_2,Kl);
  JBOUND(BB.getDataPointer(),Bz.getDataPointer(),&zip,&zip,&zip);
  
  BB(I,J_1,Kl) = Csound(I,J_1,Kh); BB(I,J_1,Kh) = Csound(I,J_1,Kl);
  BB(I,J_2,Kl) = Csound(I,J_2,Kh); BB(I,J_2,Kh) = Csound(I,J_2,Kl);   
  JBOUND(BB.getDataPointer(),Csound.getDataPointer(),&zip,&zip,&zip);
  
  BB(I,J_1,Kl) = Pressure(I,J_1,Kh); BB(I,J_1,Kh) = Pressure(I,J_1,Kl);
  BB(I,J_2,Kl) = Pressure(I,J_2,Kh); BB(I,J_2,Kh) = Pressure(I,J_2,Kl);
  JBOUND(BB.getDataPointer(),Pressure.getDataPointer(),&zip,&zip,&zip);
  
  //face centered fields
  BB(I,J_1,Kl) = Bi(I,J_1,Kh); BB(I,J_1,Kh) = Bi(I,J_1,Kl);
  BB(I,J_2,Kl) = Bi(I,J_2,Kh); BB(I,J_2,Kh) = Bi(I,J_2,Kl);
  JBOUND(BB.getDataPointer(),Bi.getDataPointer(),&uno,&zip,&zip);

  BB(I,J_1,Kl) = Bk(I,J_1,Kh); BB(I,J_1,Kh) = Bk(I,J_1,Kl);
  BB(I,J_2,Kl) = Bk(I,J_2,Kh); BB(I,J_2,Kh) = Bk(I,J_2,Kl);
  BB(I,J,NK+1) = BB(I,J,1);
  JBOUND(BB.getDataPointer(),Bk.getDataPointer(),&zip,&zip,&uno);

  BB(I,J_1,Kl) = Bj(I,J_1,Kh); BB(I,J_1,Kh) = Bj(I,J_1,Kl);
  BB(I,J_2,Kl) = Bj(I,J_2,Kh); BB(I,J_2,Kh) = Bj(I,J_2,Kl);
  JBOUND(BB.getDataPointer(),Bj.getDataPointer(),&zip,&uno,&zip);

#ifdef DEBUG
  printf("^^^^^^^^^^ after last jbound %d\n",mype); fflush(stdout);  
  printf("in boundary after j %d %f\n", mype, timer;
#endif
  
  //  Now the purely periodic K boundaries
  Range K_1(1,NO2);   
  Range K_2(NK-NO2+1,NK);

  Optimization_Manager::setOptimizedScalarIndexing(On);
  BB(I,J,K_1) = Rho(I,J,K_2); BB(I,J,K_2) = Rho(I,J,K_1);
  KBOUND(BB.getDataPointer(),Rho.getDataPointer(),&zip,&zip,&zip);

  BB(I,J,K_1) = Vx(I,J,K_2); BB(I,J,K_2) = Vx(I,J,K_1);
  KBOUND(BB.getDataPointer(),Vx.getDataPointer(),&zip,&zip,&zip);

  BB(I,J,K_1) = Vy(I,J,K_2); BB(I,J,K_2) = Vy(I,J,K_1);
  KBOUND(BB.getDataPointer(),Vy.getDataPointer(),&zip,&zip,&zip);

  BB(I,J,K_1) = Vz(I,J,K_2); BB(I,J,K_2) = Vz(I,J,K_1);
  KBOUND(BB.getDataPointer(),Vz.getDataPointer(),&zip,&zip,&zip);

  BB(I,J,K_1) = Csound(I,J,K_2); BB(I,J,K_2) = Csound(I,J,K_1);
  KBOUND(BB.getDataPointer(),Csound.getDataPointer(),&zip,&zip,&zip);

  BB(I,J,K_1) = Pressure(I,J,K_2); BB(I,J,K_2) = Pressure(I,J,K_1);
  KBOUND(BB.getDataPointer(),Pressure.getDataPointer(),&zip,&zip,&zip);

  BB(I,J,K_1) = Bx(I,J,K_2); BB(I,J,K_2) = Bx(I,J,K_1);
  KBOUND(BB.getDataPointer(),Bx.getDataPointer(),&zip,&zip,&zip);

  BB(I,J,K_1) = By(I,J,K_2); BB(I,J,K_2) = By(I,J,K_1);
  KBOUND(BB.getDataPointer(),By.getDataPointer(),&zip,&zip,&zip);

  BB(I,J,K_1) = Bz(I,J,K_2); BB(I,J,K_2) = Bz(I,J,K_1);
  KBOUND(BB.getDataPointer(),Bz.getDataPointer(),&zip,&zip,&zip);

  // face centered fields
  BB(I,J,K_1) = Bi(I,J,K_2); BB(I,J,K_2) = Bi(I,J,K_1);
  KBOUND(BB.getDataPointer(),Bi.getDataPointer(),&uno,&zip,&zip);

  BB(I,J,K_1) = Bj(I,J,K_2); BB(I,J,K_2) = Bj(I,J,K_1);
  KBOUND(BB.getDataPointer(),Bj.getDataPointer(),&uno,&zip,&zip);

  BB(I,J,K_1) = Bk(I,J,K_2); BB(I,J,K_2) = Bk(I,J,K_1);
  KBOUND(BB.getDataPointer(),Bk.getDataPointer(),&zip,&zip,&uno);

  Optimization_Manager::setOptimizedScalarIndexing(Off);

#ifdef DEBUG
  printf("^^^^^^^^^^ after last kbound %d\n",mype); fflush(stdout);
#endif

  //  Do those pesky corners for the velocities
  EDGE_FIX();
  
  // update the interior ghost cells
#ifdef DEBUG
  printf("in boundary before update ghost %d %f\n", mype, timer);
#endif
  
  TEST_ARRAYS();
  
#else       //   end of P++ stuff for boundaries
  IBOUND_FOR();
  JBOUND_FOR();
  KBOUND_FOR();
#endif     
  
  TEST_ARRAYS();
  
#ifdef DEBUG
  float timer;
  printf("in boundary after update ghost %d %f\n", mype, timer);
  TEST_ARRAYS();
  printf("<><><><><><><> after ghost %d\n",mype); fflush(stdout);
#endif
} /* void MHD::boundary() */