// Generic function to reduce the minloc,maxloc,avg of a Real over all procs onto rank0
int reduce_avg_min_max_loc(Real value, Real& avg, Real& min, Real& max, int& minloc, int& maxloc)
{
  int eek=0;
#ifdef CH_MPI
  struct
  {
    double val;
    int rank;
  } in, out;
  in.val = value;
  in.rank = procID();

  Real sum;
  eek = MPI_Reduce(&value, &sum, 1, MPI_CH_REAL, MPI_SUM, uniqueProc(SerialTask::compute), Chombo_MPI::comm);
  CH_assert(eek == MPI_SUCCESS);
  avg=sum/(Real)numProc();

  eek = MPI_Reduce(&in, &out, 1, MPI_DOUBLE_INT, MPI_MINLOC, uniqueProc(SerialTask::compute), Chombo_MPI::comm);
  CH_assert(eek == MPI_SUCCESS);
  min = out.val;
  minloc = out.rank;

  eek = MPI_Reduce(&in, &out, 1, MPI_DOUBLE_INT, MPI_MAXLOC, uniqueProc(SerialTask::compute), Chombo_MPI::comm);
  CH_assert(eek == MPI_SUCCESS);
  max = out.val;
  maxloc = out.rank;
#else
  avg=value;
  min=value;
  max=value;
  minloc=0;
  maxloc=0;
#endif
  return eek;
}

void
EBCellFAB::setInvalidData(const Real& a_val,
                          const int& a_comp)
{
  CH_assert(a_comp >= 0);
  CH_assert(a_comp < m_nComp);

  if (m_ebisBox.isAllRegular())
    {
      return;
    }
  else if (m_ebisBox.isAllCovered())
    {
      m_regFAB.setVal(a_val, a_comp);
    }
  else
    {
      for (BoxIterator bit(m_region); bit.ok(); ++bit)
        {
          const IntVect& iv = bit();
          if (m_ebisBox.isCovered(iv))
            {
              m_regFAB(iv, a_comp) = a_val;
            }
        }
      //also set the multivalued cells
      for (IVSIterator ivsit(getMultiCells());
          ivsit.ok(); ++ivsit)
        {
          const IntVect& iv = ivsit();
          m_regFAB(iv, a_comp) = a_val;
        }
      //also set the multivalued cells
    }
}

// this preconditioner first initializes phihat to (IA)phihat = rhshat
// (diagonization of L -- A is the matrix version of L)
// then smooths with a couple of passes of levelGSRB
void VCAMRPoissonOp2::preCond(LevelData<FArrayBox>&       a_phi,
                              const LevelData<FArrayBox>& a_rhs)
{
  CH_TIME("VCAMRPoissonOp2::preCond");

  // diagonal term of this operator in:
  //
  //       alpha * a(i)
  //     + beta  * sum_over_dir (b(i-1/2*e_dir) + b(i+1/2*e_dir)) / (dx*dx)
  //
  // The inverse of this is our initial multiplier.

  int ncomp = a_phi.nComp();

  CH_assert(m_lambda.isDefined());
  CH_assert(a_rhs.nComp()    == ncomp);
  CH_assert(m_bCoef->nComp() == ncomp);

  // Recompute the relaxation coefficient if needed.
  resetLambda();

  // don't need to use a Copier -- plain copy will do
  DataIterator dit = a_phi.dataIterator();
  for (dit.begin(); dit.ok(); ++dit)
    {
      // also need to average and sum face-centered bCoefs to cell-centers
      Box gridBox = a_rhs[dit].box();

      // approximate inverse
      a_phi[dit].copy(a_rhs[dit]);
      a_phi[dit].mult(m_lambda[dit], gridBox, 0, 0, ncomp);
    }

  relax(a_phi, a_rhs, 2);
}

void DenseIntVectSet::coarsen(int iref)
{
  if (iref == 1) return;
  CH_assert(iref >= 1);
  // int refinements = iref/2;
  CH_assert((iref/2)*2 == iref); // check iref for power of 2

  Box newDomain(m_domain);
  newDomain.coarsen(iref);
  DenseIntVectSet newSet(newDomain, false);
  BoxIterator bit(m_domain);
  int count=0;
  for (bit.begin(); bit.ok(); ++bit, ++count)
    {
      if (m_bits[count])
        {
          IntVect iv(bit());
          iv.coarsen(iref);
          long index = newDomain.index(iv);
          newSet.m_bits.setTrue(index);
        }
    }

  *this = newSet;
}

// Set boundary slopes:
//   The boundary slopes in a_dW are already set to one sided difference
//   approximations.  If this function doesn't change them they will be
//   used for the slopes at the boundaries.
void ExplosionIBC::setBdrySlopes(FArrayBox&       a_dW,
                                 const FArrayBox& a_W,
                                 const int&       a_dir,
                                 const Real&      a_time)
{
  CH_assert(m_isFortranCommonSet == true);
  CH_assert(m_isDefined == true);

  // In periodic case, this doesn't do anything
  if (!m_domain.isPeriodic(a_dir))
  {
    Box loBox,hiBox,centerBox,domain;
    int hasLo,hasHi;
    Box slopeBox = a_dW.box();
    slopeBox.grow(a_dir,1);

    // Generate the domain boundary boxes, loBox and hiBox, if there are
    // domain boundarys there
    loHiCenter(loBox,hasLo,hiBox,hasHi,centerBox,domain,
               slopeBox,m_domain,a_dir);

    // Set the boundary slopes if necessary
    if ((hasLo != 0) || (hasHi != 0))
    {
      FORT_SLOPEBCSF(CHF_FRA(a_dW),
                     CHF_CONST_FRA(a_W),
                     CHF_CONST_INT(a_dir),
                     CHF_BOX(loBox),
                     CHF_CONST_INT(hasLo),
                     CHF_BOX(hiBox),
                     CHF_CONST_INT(hasHi));
    }
  }
}

void EBPoissonOp::
relax(LevelData<EBCellFAB>&       a_e,
      const LevelData<EBCellFAB>& a_residual,
      int                         a_iterations)
{
  CH_TIME("EBPoissonOp::relax");

  CH_assert(a_e.ghostVect() == m_ghostCellsPhi);
  CH_assert(a_residual.ghostVect() == m_ghostCellsRHS);

  CH_assert(a_e.nComp() == 1);
  CH_assert(a_residual.nComp() == 1);

  if (m_relaxType == 0)
    {
      for (int i = 0; i < a_iterations; i++)
        {
          levelJacobi(a_e,a_residual);
        }
    }
  else  if (m_relaxType == 1)
    {
      for (int i = 0; i < a_iterations; i++)
        {
          levelMulticolorGS(a_e,a_residual);
        }
    }
  else
    {
      MayDay::Error("EBPoissonOp::relax - invalid relaxation type");
    }
}

void PatchGodunov::updateState(FArrayBox&       a_U,
                               FluxBox&         a_F,
                               Real&            a_maxWaveSpeed,
                               const FArrayBox& a_S,
                               const Real&      a_dt,
                               const Box&       a_box)
{
  CH_assert(isDefined());
  CH_assert(a_box == m_currentBox);

  int numPrim = m_gdnvPhysics->numPrimitives();
  int numFlux = m_gdnvPhysics->numFluxes();

  FluxBox whalf(a_box,numPrim);
  whalf.setVal(0.0);

  a_F.resize(a_box,numFlux);
  a_F.setVal(0.0);

  computeWHalf(whalf, a_U, a_S, a_dt, a_box);

  FArrayBox dU(a_U.box(),a_U.nComp());
  computeUpdate(dU, a_F, a_U, whalf, a_dt, a_box);

  a_U += dU;

  // Get and return the maximum wave speed on this patch/grid
  a_maxWaveSpeed = m_gdnvPhysics->getMaxWaveSpeed(a_U, m_currentBox);
}

// Set boundary fluxes
void RampIBC::primBC(FArrayBox&            a_WGdnv,
                     const FArrayBox&      a_Wextrap,
                     const FArrayBox&      a_W,
                     const int&            a_dir,
                     const Side::LoHiSide& a_side,
                     const Real&           a_time)
{
  CH_assert(m_isFortranCommonSet == true);
  CH_assert(m_isDefined == true);

  // Neither the x or y direction can be periodic
  if ((a_dir == 0 || a_dir == 1) && m_domain.isPeriodic(a_dir))
    {
      MayDay::Error("RampIBC::primBC: Neither the x or y boundaries can be periodic");
    }

  Box boundaryBox;
  getBoundaryFaces(boundaryBox, a_WGdnv.box(), a_dir, a_side);

  if (! boundaryBox.isEmpty() )
    {
      // Set the boundary fluxes
      int lohisign = sign(a_side);
      FORT_RAMPBCF(CHF_FRA(a_WGdnv),
                   CHF_CONST_FRA(a_Wextrap),
                   CHF_CONST_FRA(a_W),
                   CHF_CONST_REAL(a_time),
                   CHF_CONST_INT(lohisign),
                   CHF_CONST_REAL(m_dx),
                   CHF_CONST_INT(a_dir),
                   CHF_BOX(boundaryBox));
    }
}

NewPoissonOp* NewPoissonOpFactory::MGnewOp(const ProblemDomain& a_FineindexSpace,
                                           int   a_depth,
                                           bool  a_homoOnly)
{
  CH_assert(a_depth >= 0 );
  CH_assert(m_bc != NULL);
  NewPoissonOp* newOp = new NewPoissonOp();
  RealVect dx = m_dx;
  ProblemDomain domain = a_FineindexSpace;
  for (int i=0; i<a_depth; i++)
    {
      Box d = domain.domainBox();
      d.coarsen(8);
      d.refine(8);
      if (domain.domainBox() == d)
        {
          dx*=2;
          domain.coarsen(2);
        }
      else
        {
          return NULL;
        }
    }
  newOp->define(dx, domain, m_bc);
  return newOp;
}

// -----------------------------------------------------------------------------
// Interpolate ghosts at CF interface using zeros on coarser grids.
// This version acts on a set of IntVects.
// -----------------------------------------------------------------------------
void interpOnIVSHomo (LevelData<FArrayBox>& a_phif,
                      const DataIndex&      a_index,
                      const int             a_dir,
                      const Side::LoHiSide  a_side,
                      const IntVectSet&     a_interpIVS,
                      const Real            a_fineDxDir,
                      const Real            a_crseDxDir)
{
    CH_TIME("interpOnIVSHomo");

    // Sanity checks
    CH_assert((a_dir >= 0) && (a_dir < SpaceDim));
    CH_assert(a_phif.ghostVect()[a_dir] >= 1);

    IVSIterator fine_ivsit(a_interpIVS);
    FArrayBox& a_phi = a_phif[a_index];
    const int isign = sign(a_side);

    if (a_phi.box().size(a_dir) == 3) {
        // Linear interpolation of fine ghosts assuming
        // all zeros on coarser level.

        // we are in a 1-wide box
        IntVect iv;
        Real pa;
        Real factor = 1.0 - 2.0 * a_fineDxDir / (a_fineDxDir + a_crseDxDir);
        for (fine_ivsit.begin(); fine_ivsit.ok(); ++fine_ivsit) {
            iv = fine_ivsit();
            iv[a_dir] -= isign;
            // Use linear interpolation
            for (int ivar = 0; ivar < a_phif.nComp(); ivar++) {
                pa = a_phi(iv, ivar);
                a_phi(fine_ivsit(), ivar) = factor * pa;
            }
        }
    } else {
        // Quadratic interpolation of fine ghosts assuming
        // all zeros on coarser level.

        // Symbolic reduced version of CF quadratic stencil
        Real pa, pb;
        Real c1 = 2.0*(a_crseDxDir-a_fineDxDir)/(a_crseDxDir+    a_fineDxDir); //first inside point
        Real c2 =    -(a_crseDxDir-a_fineDxDir)/(a_crseDxDir+3.0*a_fineDxDir); // next point inward
        IntVect ivf;
        for (fine_ivsit.begin(); fine_ivsit.ok(); ++fine_ivsit) {
            ivf = fine_ivsit();
            // Use quadratic interpolation
            for (int ivar = 0; ivar < a_phif.nComp(); ++ivar) {
                ivf[a_dir]-=2*isign;
                pa = a_phi(ivf, ivar);
                ivf[a_dir]+=isign;
                pb = a_phi(ivf, ivar);

                ivf[a_dir]+=isign;
                a_phi(fine_ivsit(), ivar) = c1*pb + c2*pa;
            } //end loop over components
        } //end loop over fine intvects
    }
}

// ---------------------------------------------------------
void
NodeQCFI::coarseFineInterp(LevelData<NodeFArrayBox>& a_phiFine,
                           const LevelData<NodeFArrayBox>& a_phiCoarse,
                           bool a_inhomogeneous)
{
  CH_assert(isDefined());
  CH_assert(a_phiFine.nComp() == m_ncomp);
  CH_assert(a_phiCoarse.nComp() == m_ncomp);

  if (m_coarsenings == 1)
    {
      m_qcfi2[0]->coarseFineInterp(a_phiFine, a_phiCoarse);
    }
  else // m_coarsenings >= 2
    {
      // if m_coarsenings == 2:
      // m_qcfi2[1] = new NodeQuadCFInterp2(a_grids, true);
      // m_qcfi2[0] = new NodeQuadCFInterp2(a_grids.coarsen(2), false);
      // m_inter[0] = new LevelData(a_grids.coarsen(2));

      // m_qcfi2[0] interpolates m_inter[0] from a_phiCoarse on all of it.
      // m_qcfi2[1] interpolates a_phiFine from m_inter[0] on interface only.

      m_qcfi2[0]->coarseFineInterp(*m_inter[0], a_phiCoarse);
      m_inter[0]->exchange(m_inter[0]->interval());

      // Set domain and dx for coarsest level refined by 2.
      ProblemDomain domain(m_domainPenultimate);
      Real dx = m_dxPenultimate;

      bool homogeneous = !a_inhomogeneous;
      const DisjointBoxLayout& dbl0 = m_inter[0]->disjointBoxLayout();
      for (DataIterator dit = dbl0.dataIterator(); dit.ok(); ++dit)
        {
          m_bc((*m_inter[0])[dit()],  dbl0.get(dit()), domain,  dx, homogeneous);
        }

      for (int interlev = 1; interlev < m_coarsenings - 1; interlev++)
        {
          m_qcfi2[interlev]->coarseFineInterp(*m_inter[interlev],
                                              *m_inter[interlev-1]);

          m_inter[interlev]->exchange(m_inter[interlev]->interval());

          domain.refine(2);
          dx *= 0.5;
          const DisjointBoxLayout& dbl = m_inter[interlev]->disjointBoxLayout();
          for (DataIterator dit = dbl.dataIterator(); dit.ok(); ++dit)
            {
              m_bc((*m_inter[interlev])[dit()],  dbl.get(dit()), domain,  dx, homogeneous);
            }
        }

      m_qcfi2[m_coarsenings-1]->coarseFineInterp(a_phiFine,
                                                 *m_inter[m_coarsenings-2]);
      // Don't need to set boundary conditions for a_phiFine
      // because the calling function will do that.
    }
}

IntVectSet& IntVectSet::grow(int idir, int igrow)
{
  CH_assert(idir >= 0);
  CH_assert(idir < SpaceDim);
  if (m_isdense) m_dense.grow(idir, igrow);
  else          m_ivs.grow(idir, igrow);
  return *this;
}

// ---------------------------------------------------------
FArrayBox&
FluxBox::getFlux(const int dir)
{
  CH_assert(m_nvar >0);
  CH_assert(dir < SpaceDim);
  CH_assert(m_fluxes[dir] != NULL);

  return *m_fluxes[dir];
}

void
NoFlowAdvectBC::
fluxBC(EBFluxFAB&            a_primGdnv,
       const EBCellFAB&      a_primCenter,
       const EBCellFAB&      a_primExtrap,
       const Side::LoHiSide& a_side,
       const Real&           a_time,
       const EBISBox&        a_ebisBox,
       const DataIndex&      a_dit,
       const Box&            a_box,
       const Box&            a_faceBox,
       const int&            a_dir)
{
  CH_assert(m_isDefined);

  Box FBox = a_faceBox;
  Box cellBox = FBox;
  CH_assert(a_primGdnv[a_dir].nComp()==1);

  // Determine which side and thus shifting directions
  int isign = sign(a_side);
  cellBox.shiftHalf(a_dir,isign);

  // Is there a domain boundary next to this grid
  if (!m_domain.contains(cellBox))
    {
      cellBox &= m_domain;
      // Find the strip of cells next to the domain boundary
      Box bndryBox = adjCellBox(cellBox, a_dir, a_side, 1);

      // Shift things to all line up correctly
      bndryBox.shift(a_dir,-isign);

      IntVectSet ivs(bndryBox);
      for (VoFIterator vofit(ivs, a_ebisBox.getEBGraph()); vofit.ok(); ++vofit)
        {
          const VolIndex& vof = vofit();

          Vector<FaceIndex> bndryFaces = a_ebisBox.getFaces(vof, a_dir, a_side);
          for (int iface= 0; iface < bndryFaces.size(); iface++)
            {
              const FaceIndex& face = bndryFaces[iface];
              //set all fluxes to zero then fix momentum flux
              //solid wall
              if (a_dir == m_velComp)
                {
                  a_primGdnv[a_dir](face, 0) = 0.0;
                }
              else
                {
                  a_primGdnv[a_dir](face, 0) = a_primExtrap(vof, 0);
                }
            }
        }
    }
}

void
EBMGInterp::pwlInterp(LevelData<EBCellFAB>&       a_fineData,
                      const LevelData<EBCellFAB>& a_coarData,
                      const Interval&             a_variables)
{
  CH_TIME("EBMGInterp::pwlInterp");
  CH_assert(a_fineData.ghostVect() == m_ghost);
  CH_assert(a_coarData.ghostVect() == m_ghost);
  CH_assert(m_doLinear); //otherwise stencils have not been defined

  if (m_layoutChanged)
    {
      if (m_coarsenable)
        {
          CH_TIME("EBMGInterp::pwlInterp::coarsenable");
          EBCellFactory ebcellfact(m_buffEBISL);
          LevelData<EBCellFAB> coarsenedFineData(m_buffGrids, m_nComp, m_ghost, ebcellfact);
          a_coarData.copyTo(a_variables, coarsenedFineData, a_variables);
          fillGhostCellsPWC(coarsenedFineData, m_buffEBISL, m_coarDomain);
          for (DataIterator dit = m_fineGrids.dataIterator(); dit.ok(); ++dit)
            {
              //does incrementonly = true
              pwlInterpFAB(a_fineData[dit()],
                           m_buffGrids[dit()],
                           coarsenedFineData[dit()],
                           dit(),
                           a_variables);
            }
        }
      else
        {
          CH_TIME("EBMGInterp::pwlInterp::uncoarsenable");
          EBCellFactory ebcellfact(m_buffEBISL);
          fillGhostCellsPWC((LevelData<EBCellFAB>&)a_coarData, m_coarEBISL, m_coarDomain);
          LevelData<EBCellFAB> refinedCoarseData(m_buffGrids, m_nComp, m_ghost, ebcellfact);
          for (DataIterator dit = m_coarGrids.dataIterator(); dit.ok(); ++dit)
            {
              refinedCoarseData[dit()].setVal(0.);
              pwlInterpFAB(refinedCoarseData[dit()],
                           m_coarGrids[dit()],
                           a_coarData[dit()],
                           dit(),
                           a_variables);
            }

          EBAddOp op;
          refinedCoarseData.copyTo(a_variables, a_fineData, a_variables, m_copierRCtoF, op);
        }
    }
  else
    {
      pwcInterpMG(a_fineData, a_coarData, a_variables);
    }
}

// Compute the speed of sound
void PolytropicPhysics::soundSpeed(FArrayBox& a_speed,
                                   const FArrayBox& a_U,
                                   const Box&       a_box)
{
    CH_assert(isDefined());
    CH_assert(a_U.contains(a_box));
    CH_assert(a_speed.contains(a_box));

    FORT_SOUNDSPEEDF(CHF_CONST_FRA1(a_speed, 0),
                     CHF_CONST_FRA(a_U),
                     CHF_BOX(a_box));
}

void EBPoissonOp::
GSColorAllRegular(LevelData<EBCellFAB>&        a_phi,
                  const LevelData<EBCellFAB>&  a_rhs,
                  const IntVect&               a_color,
                  const Real&                  a_weight,
                  const bool&                  a_homogeneousPhysBC)
{
  CH_TIME("EBPoissonOp::GSColorAllRegular");
  CH_assert(a_rhs.ghostVect() == m_ghostCellsRHS);
  CH_assert(a_phi.ghostVect() == m_ghostCellsPhi);

  int nComps = a_phi.nComp();
  for (DataIterator dit = a_phi.dataIterator(); dit.ok(); ++dit)
    {
      Box dblBox(m_eblg.getDBL().get(dit()));
      BaseFab<Real>& phiFAB       = (a_phi[dit()] ).getSingleValuedFAB();
      const BaseFab<Real>& rhsFAB = (a_rhs[dit()] ).getSingleValuedFAB();

      Box loBox[SpaceDim],hiBox[SpaceDim];
      int hasLo[SpaceDim],hasHi[SpaceDim];

      applyDomainFlux(loBox, hiBox, hasLo, hasHi,
                      dblBox, nComps, phiFAB,
                      a_homogeneousPhysBC, dit(),m_beta);

      IntVect loIV = dblBox.smallEnd();
      IntVect hiIV = dblBox.bigEnd();

      for (int idir = 0; idir < SpaceDim; idir++)
        {
          if (loIV[idir] % 2 != a_color[idir])
            {
              loIV[idir]++;
            }
        }

      if (loIV <= hiIV)
        {
          Box coloredBox(loIV, hiIV);

          for (int comp=0; comp<a_phi.nComp(); comp++)
            {
              FORT_DOALLREGULARMULTICOLOR(CHF_FRA1(phiFAB,comp),
                                          CHF_CONST_FRA1(rhsFAB,comp),
                                          CHF_CONST_REAL(a_weight),
                                          CHF_CONST_REAL(m_alpha),
                                          CHF_CONST_REAL(m_beta),
                                          CHF_CONST_REALVECT(m_dx),
                                          CHF_BOX(coloredBox));
            }
        }
    }
}

void
EBPlanarShockIBC::
initialize(LevelData<EBCellFAB>& a_conState,
           const EBISLayout& a_ebisl) const
{
 CH_assert(m_isDefined);
 CH_assert(m_isFortranCommonSet);

  // Iterator of all grids in this level
  for(DataIterator dit = a_conState.dataIterator(); dit.ok(); ++dit)
    {
      const EBISBox& ebisBox = a_ebisl[dit()];
      // Storage for current grid
      EBCellFAB& conFAB = a_conState[dit()];

      BaseFab<Real>& regConFAB = conFAB.getSingleValuedFAB();
      // Box of current grid
      Box uBox = regConFAB.box();
      uBox &= m_domain;
      // Set up initial condition in this grid
      /**/
      FORT_PLANARSHOCKINIT(CHF_CONST_FRA(regConFAB),
                           CHF_CONST_REAL(m_dx),
                           CHF_BOX(uBox));
      /**/

      //now for the multivalued cells.  Since it is pointwise,
      //the regular calc is correct for all single-valued cells.
      IntVectSet ivs = ebisBox.getMultiCells(uBox);
      for(VoFIterator vofit(ivs, ebisBox.getEBGraph()); vofit.ok(); ++vofit)
        {
          const VolIndex& vof = vofit();
          const IntVect& iv = vof.gridIndex();
          RealVect momentum;
          Real energy, density;
          /**/
          FORT_POINTPLANARSHOCKINIT(CHF_REAL(density),
                                    CHF_REALVECT(momentum),
                                    CHF_REAL(energy),
                                    CHF_CONST_INTVECT(iv),
                                    CHF_CONST_REAL(m_dx));
          /**/
          conFAB(vof, CRHO) = density;
          conFAB(vof, CENG) = energy;
          for(int idir = 0; idir < SpaceDim; idir++)
            {
              conFAB(vof, CMOMX+idir) = momentum[idir];
            }
        }//end loop over multivalued cells
    } //end loop over boxes
}

// Define new AMR level
void AMRLevelPluto::define(AMRLevel*            a_coarserLevelPtr,
                           const ProblemDomain& a_problemDomain,
                           int                  a_level,
                           int                  a_refRatio)
{
  if (s_verbosity >= 3)
  {
    pout() << "AMRLevelPluto::define " << a_level << endl;
  }

  // Call inherited define
  AMRLevel::define(a_coarserLevelPtr,
                   a_problemDomain,
                   a_level,
                   a_refRatio);

  // Get setup information from the next coarser level
  if (a_coarserLevelPtr != NULL)
  {
    AMRLevelPluto* amrGodPtr = dynamic_cast<AMRLevelPluto*>(a_coarserLevelPtr);

    if (amrGodPtr != NULL)
    {
      m_cfl = amrGodPtr->m_cfl;
      m_domainLength = amrGodPtr->m_domainLength;
      m_refineThresh = amrGodPtr->m_refineThresh;
      m_tagBufferSize = amrGodPtr->m_tagBufferSize;
    }
    else
    {
      MayDay::Error("AMRLevelPluto::define: a_coarserLevelPtr is not castable to AMRLevelPluto*");
    }
  }

  // Compute the grid spacing
  m_dx = m_domainLength / (a_problemDomain.domainBox().bigEnd(0)-a_problemDomain.domainBox().smallEnd(0)+1.);

  // Nominally, one layer of ghost cells is maintained permanently and
  // individual computations may create local data with more
  m_numGhost = 1;

  CH_assert(m_patchPluto != NULL);
  CH_assert(isDefined());
  m_patchPluto->define(m_problem_domain,m_dx,m_level,m_numGhost);

  // Get additional information from the patch integrator
  m_numStates      = m_patchPluto->numConserved();
  m_ConsStateNames = m_patchPluto->ConsStateNames();
  m_PrimStateNames = m_patchPluto->PrimStateNames();
}

void LinElastPhysics::consToPrim(FArrayBox&       a_W,
    const FArrayBox& a_U,
    const Box&       a_box)
{
    //JK Our primitives are our conserved variables
    //JK pout() << "LinElastPhysics::consToPrim" << endl;
    CH_assert(isDefined());
    CH_assert(a_U.box().contains(a_box));
    CH_assert(a_W.box().contains(a_box));

    FORT_COPYF(CHF_FRA(a_W),
        CHF_CONST_FRA(a_U),
        CHF_BOX(a_box));
}