#ifdef DIVTEST
//#define BOUNDS_CHECK
#include "site-stuff.h"

#include "help.h"
/// "help" common block defined in help.inc; help struct defined in help.h
extern HELPER HELP;

/*  run-time definition only needed if MHD.h not used
#include "run-time.h"

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

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

#include <stdio.h>

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

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

#ifdef LC_UNSC
#define    RESETB    resetb_
#define    FIXBOUT   fixbout_
#define    TEST_ARRAYS test_arrays_
#define    FORT_DIVCHECK  fort_divcheck_
#define    ADD_BZERO   add_bzero_
#define    SUBTRACT_BZERO   subtract_bzero_
#define  FDIVCHECK fdivcheck_
#elif LC_NOUNSC
#define    RESETB    resetb
#define    FIXBOUT   fixbout
#define    TEST_ARRAYS test_arrays
#define    FORT_DIVCHECK  fort_divcheck
#define  FDIVCHECK fdivcheck
#define    ADD_BZERO   add_bzero
#define    SUBTRACT_BZERO   subtract_bzero
#endif 

/// Defined in fixB.F
extern "C" void FIXBOUT();
/// Defined in fixB.F
extern "C" void RESETB();
/// Defined in test-arrays.F
extern "C" void TEST_ARRAYS();
/// Defined in test-arrays.F
extern "C" void FORT_DIVCHECK(FLoaT *, FLoaT *, int *);
/// Defined in test-arrays.F
extern "C" void FDIVCHECK(int *, int *, int*);
/// Defined in test-arrays.F
extern "C" void ADD_BZERO();
/// Defined in test-arrays.F
extern "C" void SUBTRACT_BZERO();

/// Compute maximum divergence of B
/**
 *  Compute divergence (DV) of B and return divmax, the maximum value
 *  of the divergence
 *
 *  @param mhd @param FIXB - boolean.  If true, multiply the magnetic
 *  field B[i,j,k] by Face[i,j,k] before checking the div.
 *
 *  @return divmax
 */
float DivTest(MHD *mhd,int FIXB)
{
#ifdef DEBUG_MODE_ON
  printf("DEBUG:  in THEREST.C::(DivTest)\n");
#endif
  
  int NI = mhd->NI;  
  int NK = mhd->NK; 
  int NJ = mhd->NJ;
  Range I(1,NI); 
  Range J(1,NJ);   
  Range K(1,NK);
  Range II(11,15);
  int mype = Communication_Manager::My_Process_Number;  
  int zero = TRUE; int flux = TRUE; int i2 = FALSE;
  
  FloatARRAY DV(NI+1,NJ+1,NK+1); 
  DV.setBase(1); 
  DV.partition(mhd->MHD_P);
  
  FloatARRAY AV(NI+1,NJ+1,NK+1); 
  AV.setBase(1); 
  AV.partition(mhd->MHD_P);
  
  if (FIXB != FALSE) 
  {
    // FIXBOUT();
    mhd->Bi *= mhd->Face_i;
    mhd->Bj *= mhd->Face_j;
    mhd->Bk *= mhd->Face_k;
  }

  TEST_ARRAYS();
  /*   
fflush(NULL); Communication_Manager::Sync();
 mhd->Bi(II,II,II).display("BI ");
fflush(NULL); Communication_Manager::Sync();
 mhd->Bj(II,II,II).display("BJ");
fflush(NULL); Communication_Manager::Sync(); 
 mhd->Bk(II,II,II).display("BK ");
fflush(NULL); Communication_Manager::Sync(); 
 */  
  FDIVCHECK(&zero, &flux, &i2);
  ADD_BZERO();
  
  DV = 0.0;
  DV(I,J,K) =  (mhd->Bi(I+1,J,K) - mhd->Bi(I,J,K));
  DV(I,J,K) += (mhd->Bj(I,J+1,K)- mhd->Bj(I,J,K)); 
  DV(I,J,K) += (mhd->Bk(I,J,K+1)- mhd->Bk(I,J,K));
  //DV(II,II,II).display("DV ");
  //fflush(NULL);

  Communication_Manager::Sync(); 
 
  AV(I,J,K)  =  abs(mhd->Bi(I+1,J,K)) + abs(mhd->Bi(I,J,K));
  AV(I,J,K)  +=	abs(mhd->Bj(I,J+1,K)) + abs(mhd->Bj(I,J,K));
  AV(I,J,K)  += abs(mhd->Bk(I,J,K+1)) + abs(mhd->Bk(I,J,K))+1.e-6*mhd->Face_k(I,J,K);
  // AV(II,II,II).display("AV ");
  //fflush(NULL);
  
  SUBTRACT_BZERO();  

  Communication_Manager::Sync(); 
  
  DV = abs(DV)/AV;
  DV.updateGhostBoundaries();
  TEST_ARRAYS();
  
  //FDIVCHECK(&zero, &flux, &i2);
  
  intArray BadPt(NI+1,NJ+1,NK+1); BadPt.setBase(1); BadPt.partition(mhd->MHD_P);
  
  BadPt = DV >= 1.e-1; 
  int num1 = sum(BadPt); 
  int maxMask1 = sum(BadPt(II,II,II));
  
  BadPt = DV >= 1.e-2; 
  int num2 = sum(BadPt); 
  int maxMask2 = sum(BadPt(II,II,II));
  
  // BadPt(II,II,II).display("Mask ");
  
  BadPt = DV >= 1.e-3; 
  int num3 = sum(BadPt); 
  int maxMask3 = sum(BadPt(II,II,II));
  
  BadPt = DV >= 1.e-4; 
  int num4 = sum(BadPt); 
  int maxMask4 = sum(BadPt(II,II,II));
  
  BadPt = DV >= 1.e-5; 
  int num5 = sum(BadPt); 
  int maxMask5 = sum(BadPt(II,II,II));
  
  // if ( mype == 0 ) printf("number of bad pts: %d %d %d %d %d \n",maxMask1,maxMask2,maxMask3,maxMask4,maxMask5);
  // printf("mype = %d ::number of bad pts: %d %d %d %d %d \n",mype,num1,num2,num3,num4,num5);
  float divmax = max(DV(I,J,K)); 

 /*printf("divmax %g\n",divmax);
 // Optimization_Manager::setOptimizedScalarIndexing(On);
 for (int k=1; k<=NK; k++) {
 for (int j=1; j<=NJ; j++) {
 for (int i=1; i<=NI; i++) {
   if ( DV(i,j,k) >= 0.99*divmax) {
   printf("mype=%d, i=%d,j=%d,k=%d, div=%e\n",RTIME.local_proc,i,j,k,DV(i,j,k));
 }}}} 
 //Optimization_Manager::setOptimizedScalarIndexing(Off);
 */
  
  if (FIXB != FALSE)
  {
    mhd->Bi /= mhd->Face_i;
    mhd->Bj /= mhd->Face_j;
    mhd->Bk /= mhd->Face_k;
    // RESETB();
  }
  // TEST_ARRAYS();
  //FORT_DIVCHECK(DV.getDataPointer(),AV.getDataPointer(),BadPt.getDataPointer());
  
  return divmax;
}
#endif