AmrIce.H

Go to the documentation of this file.

#ifdef CH_LANG_CC
/*
*      _______              __
*     / ___/ /  ___  __ _  / /  ___
*    / /__/ _ \/ _ \/  V \/ _ \/ _ \
*    \___/_//_/\___/_/_/_/_.__/\___/
*    Please refer to Copyright.txt, in Chombo's root directory.
*/
#endif

#ifndef AMR_ICE_H
#define AMR_ICE_H

#include "LevelData.H"
#include "FArrayBox.H"
#include "LevelSigmaCS.H"
#include "IceVelocitySolver.H"
#include "RealVect.H"
#include "PatchGodunov.H"
#include "AdvectPhysics.H"
#include "SurfaceFlux.H"
#include "BasalFriction.H"
#include "IceThicknessIBC.H"
#include "IceInternalEnergyIBC.H"
#include "CalvingModel.H"
#include "MuCoefficient.H"
#include "CH_HDF5.H"
#include "DomainDiagnosticData.H"
#include "AmrIceBase.H"
#include "ParmParse.H"
#include <sstream>
#ifdef HAVE_PYTHON
#include "PythonInterface.H"
#endif
#include "NamespaceHeader.H"


class AmrIce : public AmrIceBase
{

public:

  AmrIce();


  virtual ~AmrIce();

  void setDefaults();


  void initialize();

  void setDomainSize(const RealVect& a_domainSize)
  { m_domainSize = a_domainSize;}
    


  void setConstitutiveRelation(ConstitutiveRelation* a_constRelPtr) 
  {m_constitutiveRelation = a_constRelPtr->getNewConstitutiveRelation(); }

  //const ConstitutiveRelation* constitutiveRelation() const
  //{return m_constitutiveRelation;}


  void setRateFactor(RateFactor* a_rateFactorPtr)
  {m_rateFactor = a_rateFactorPtr->getNewRateFactor(); }



  void setBasalRateFactor(RateFactor* a_rateFactorPtr)
  {m_basalRateFactor = a_rateFactorPtr->getNewRateFactor(); }
  


  void setBasalFrictionRelation(BasalFrictionRelation* a_basalFrictionRelPtr) 
  {m_basalFrictionRelation = a_basalFrictionRelPtr->getNewBasalFrictionRelation(); }


  void setThicknessBC( IceThicknessIBC* a_thicknessIBC);
  
  void setInternalEnergyBC( IceInternalEnergyIBC* a_internalEnergyIBC);

  void setSurfaceFlux(SurfaceFlux* a_surfaceFluxPtr)
  {
    if (m_surfaceFluxPtr != NULL)
      {
        delete(m_surfaceFluxPtr);
        m_surfaceFluxPtr = NULL;
      }
    m_surfaceFluxPtr = a_surfaceFluxPtr->new_surfaceFlux(); 
  }
  
  void setBasalFlux(SurfaceFlux* a_basalFluxPtr)
  {
    if (m_basalFluxPtr != NULL)
      {
        delete(m_basalFluxPtr);
        m_basalFluxPtr = NULL;
      }
    m_basalFluxPtr = a_basalFluxPtr->new_surfaceFlux(); 
  }

  void setSurfaceHeatBoundaryData(SurfaceFlux* a_surfaceHeatBoundaryDataPtr, bool a_dirichlett, bool a_temperature)
  {
    if (m_surfaceHeatBoundaryDataPtr != NULL)
      {
        delete(m_surfaceHeatBoundaryDataPtr);
        m_surfaceHeatBoundaryDataPtr = NULL;
      }
   
    m_surfaceBoundaryHeatDataDirichlett = a_dirichlett;
    m_surfaceBoundaryHeatDataTemperature = a_temperature;

    if (a_surfaceHeatBoundaryDataPtr != NULL)
      {
        m_surfaceHeatBoundaryDataPtr = a_surfaceHeatBoundaryDataPtr->new_surfaceFlux(); 
      }
  }
  


  void setBasalHeatBoundaryData(SurfaceFlux* a_basalHeatBoundaryDataPtr)
  {
    if (m_basalHeatBoundaryDataPtr != NULL)
      {
        delete(m_basalHeatBoundaryDataPtr);
        m_basalHeatBoundaryDataPtr = NULL;
      }
    if (a_basalHeatBoundaryDataPtr != NULL)
      {
        m_basalHeatBoundaryDataPtr = a_basalHeatBoundaryDataPtr->new_surfaceFlux(); 
      }
  }
  
  


  void setTopographyFlux(SurfaceFlux* a_topographyFluxPtr)
  {
    if (m_topographyFluxPtr != NULL)
      {
        delete(m_topographyFluxPtr);
        m_topographyFluxPtr = NULL;
      }
    m_topographyFluxPtr = a_topographyFluxPtr->new_surfaceFlux(); 
  }

  void setBasalFriction(const BasalFriction* a_basalFrictionPtr)
  {m_basalFrictionPtr = a_basalFrictionPtr->new_basalFriction();}
   
  void setMuCoefficient(const MuCoefficient* a_muCoefficientPtr)
  {m_muCoefficientPtr = a_muCoefficientPtr->new_muCoefficient();}

  void setCalvingModel(CalvingModel* a_calvingModelPtr)
  {
    if (m_calvingModelPtr) delete m_calvingModelPtr;
    m_calvingModelPtr =  a_calvingModelPtr->new_CalvingModel();
  }

  
  void run(Real a_max_time, int a_max_step);

  void timeStep(Real a_dt);

  // compute face-centered half-time ice thickness using ppm
  void computeH_half(Vector<LevelData<FluxBox>* >& a_H_half, Real a_dt);

  // compute thickness fluxes
  void computeThicknessFluxes(Vector<LevelData<FluxBox>* >& a_vectFluxes,
                              const Vector<LevelData<FluxBox>* >& a_H_half,
                              const Vector<LevelData<FluxBox>* >& a_faceVelAdvection);

  // update ice thickness *and* bedrock elevation
  void updateGeometry(Vector<RefCountedPtr<LevelSigmaCS> >& a_vect_coordSys_new, 
                       Vector<RefCountedPtr<LevelSigmaCS> >& a_vectCoords_old, 
                       const Vector<LevelData<FluxBox>* >& a_vectFluxes, 
                       Real a_dt);
  
  

  void regrid();

  Real computeTotalIce() const;


  void getIceThickness(Real* a_data_ptr, int* a_dim_info, 
                       Real* a_dew, Real* a_dns) const;


  bool isDefined() const;

#ifdef CH_USE_HDF5

private:

  HDF5HeaderData m_headerData;

  void writeMetaDataHDF5(HDF5Handle& a_handle) const;

public:
  void setHeader(const std::string& a_key, const std::string& a_val)
  {
    m_headerData.m_string[a_key] = a_val;
  }
  void getHeader(const std::string& a_key, std::string& a_val)
  {
    a_val =  m_headerData.m_string[a_key];
  }
  void setHeader(const std::string& a_key, const int& a_val)
  {
    m_headerData.m_int[a_key] = a_val;
  }
  void getHeader(const std::string& a_key, int& a_val)
  {
    a_val = m_headerData.m_int[a_key];
  }
  void setHeader(const std::string& a_key, const Real& a_val)
  {
    m_headerData.m_real[a_key] = a_val;
  }
  
  void getHeader(const std::string& a_key, Real& a_val)
  {
    a_val =  m_headerData.m_real[a_key];
  }

  void writeAMRHierarchyHDF5(HDF5Handle& a_handle, 
                                     const Vector<DisjointBoxLayout>& a_grids,
                                     const Vector<LevelData<FArrayBox>* > & a_data,
                                     const Vector<string>& a_name,
                                     const Box& a_domain,
                                     const Real& a_dx,
                                     const Real& a_dt,
                                     const Real& a_time,
                                     const Vector<int>& a_ratio,
                             const int& a_numLevels) const;
  
  void writeAMRHierarchyHDF5(const string& filename,
                             const Vector<DisjointBoxLayout>& a_grids,
                             const Vector<LevelData<FArrayBox>* > & a_data,
                             const Vector<string>& a_name,
                             const Box& a_domain,
                             const Real& a_dx,
                             const Real& a_dt,
                             const Real& a_time,
                             const Vector<int>& a_ratio,
                             const int& a_numLevels) const;
 

  void writePlotFile();
  void writeAMRPlotFile();

  void writeCheckpointFile() ;


  void writeCheckpointFile(const string& a_file) ;


  void readCheckpointFile(HDF5Handle& a_handle);

  void restart(const string& a_restart_file);

  



private:

  LevelData<FArrayBox> m_uniform_cf_data;
  int m_cf_level;
  Vector<std::string> m_uniform_cf_data_name;
  Vector<std::string> m_uniform_cf_standard_name;
  Vector<std::string> m_uniform_cf_long_name;
  Vector<std::string> m_uniform_cf_units;
  Vector<std::function<const LevelData<FArrayBox>*(int, LevelData<FArrayBox>&) > >m_cf_field_function;
  // the function prototype include a LevelData<FArrayBox>& that can be used for storage by the implementation
  
  LevelData<FArrayBox> m_cf_tmp_data;
  Vector<Interval> m_cf_field_interval;
  Real m_cf_dt;
  Real m_cf_start_time;
  bool m_plot_style_cf, m_plot_style_amr;
  enum OutputFileNumbering { time_step = 0, time_seconds = 1, time_years = 2, time_yyyymmdd_360 = 3, NUM_PLOT_NUMBERING};
  OutputFileNumbering m_output_file_numbering;

  long outputNumbering(); 

  friend class DomainDiagnosticData;
  DomainDiagnosticData m_cf_domain_diagnostic_data;


  void setOutputOptions(ParmParse& a_pp);

  void flushCFData();

  void accumulateCFData(Real a_dt, bool a_reset = false);

  void initCFData();
public:
  
#endif

  //Real computeTotalGroundedIce() const;

  //Real computeVolumeAboveFlotation() const;



  void computeAreaFraction(LevelData<FArrayBox>& a_area,
                           int a_maskVal, int a_level) const;

  // Real computeArea(int a_maskVal) const;
  
  //Real computeGroundedArea() const;

  //Real computeFloatingArea() const;

  //Real computeFluxOverIce(const Vector<LevelData<FArrayBox>*> a_flux);

  //Real computeDeltaVolumeOverIce() const;

  //Real computeTotalFlux(const Vector<LevelData<FArrayBox>*> a_flux);


  enum velSolverTypes { Picard = 0,
                        JFNK = 1,
                        KnownVelocity = 2,
                        PetscNLSolver = 3,
                        FASMGAMR = 4,
                        Python = 5,
                        InverseVerticallyIntegrated = 6,
                        NUM_SOLVER_TYPES};

  void setTime(Real a_time)
  {
    m_time = a_time - m_offsetTime;
  }

  Real time() const
  {
    return m_time + m_offsetTime;
  };

  Real dt() const
  {
    return m_dt;
  }

  int finestLevel() const
  {
    return m_finest_level;
  }

  const RealVect dx(int a_level) const
  {
    return m_vect_coordSys[a_level]->dx();
  };
  
  const RefCountedPtr<LevelSigmaCS> geometry(int a_level) const 
  {
    return m_vect_coordSys[a_level];
  };
  const DisjointBoxLayout& grids(int a_level) const
  {
    return m_amrGrids[a_level];
  }
  const Vector<DisjointBoxLayout>& grids() const
  {
    return m_amrGrids;
  }
  const Vector<RefCountedPtr<LevelSigmaCS> >& amrGeometry() const
  {
    return m_vect_coordSys;
  }

  void incrementIceThickness(Vector<LevelData<FArrayBox>*> a_thk);
  
  const Vector<int>& refRatios() const { return m_refinement_ratios;}

  const Vector<Real>& amrDx() const {return m_amrDx;}

  const LevelData<FArrayBox>* velocity(int a_level) const
  {
    return m_velocity[a_level];
  };

  const Vector<LevelData<FArrayBox>* >& amrVelocity() 
  {
    return m_velocity;
  }

  SurfaceFlux& surfaceHeatBoundaryData() const
  {
    if (m_surfaceHeatBoundaryDataPtr == NULL)
      {
        CH_assert(m_surfaceHeatBoundaryDataPtr != NULL);
        MayDay::Error("no surface heat boundary data");
      }
    return *m_surfaceHeatBoundaryDataPtr;
  }

  bool surfaceHeatBoundaryDirichlett() const
  {
    return m_surfaceBoundaryHeatDataDirichlett;
  }


  LevelData<FArrayBox>* iceFrac(int a_level) 
  {
    return m_iceFrac[a_level];
  };

  bool surfaceHeatBoundaryTemperature() const
  {
    return m_surfaceBoundaryHeatDataTemperature;
  }

  SurfaceFlux& basalHeatBoundaryData() const
  {
    if (m_basalHeatBoundaryDataPtr == NULL)
      {
        CH_assert(m_basalHeatBoundaryDataPtr != NULL);
        MayDay::Error("no surface heat boundary data");
      }
    return *m_basalHeatBoundaryDataPtr;
  }


  const LevelData<FArrayBox>* iceFrac(int a_level) const
  {
    return m_iceFrac[a_level];
  };
  //** useful for calving models **/
  const LevelData<FArrayBox>& muCoef(int a_level) const
  {
    return *m_cellMuCoef[a_level];
  }
  
  Vector<LevelData<FArrayBox>* >& amrIceFrac() 
  {
    return m_iceFrac;
  }

  const Vector<LevelData<FArrayBox>* >& amrIceFrac() const
  {
    return m_iceFrac;
  }

  const LevelData<FArrayBox>* groundingLineProximity(int a_level) const;

  Real groundingLineProximityScale() const
  {
    return m_groundingLineProximityScale;
  }
  

  const LevelData<FArrayBox>* surfaceThicknessSource(int a_level) const;
  const LevelData<FArrayBox>* basalThicknessSource(int a_level) const;
  const LevelData<FArrayBox>* viscousTensor(int a_level) const;
  const LevelData<FArrayBox>* dragCoefficient(int a_level) const;
  const LevelData<FArrayBox>* viscosityCoefficient(int a_level) const;

  const Vector<LevelData<FluxBox>* >& faceVelocities() {return m_faceVelTotal;}

  const Vector<LevelData<FluxBox>* >& faceVelocitiesLayered() const {return m_layerXYFaceXYVel;} 

  const  Vector<LevelData<FArrayBox>* >& calvedIceThickness() const {return m_calvedIceThickness; }
  const  Vector<LevelData<FArrayBox>* >& removedIceThickness() const {return m_removedIceThickness; }
  const  Vector<LevelData<FArrayBox>* >& addedIceThickness() const {return m_addedIceThickness; }
  
  const  Vector<LevelData<FArrayBox>* >& surfaceInternalEnergy() const {return m_sInternalEnergy; }

  const  Vector<LevelData<FArrayBox>* >& surfaceHeatFlux() const {return m_sHeatFlux; }

protected:

  



  void tagCells(Vector<IntVectSet>& a_tags);

  void tagCellsLevel(IntVectSet& a_tags, int a_level);
 
  void tagCellsInit(Vector<IntVectSet>& a_tags);

  void initGrids(int a_finest_level);
  
  void setupFixedGrids(const std::string& a_gridFile);

  void levelSetup(int a_level, const DisjointBoxLayout& a_grids);
                       
  void initData(Vector<RefCountedPtr<LevelSigmaCS > >& a_vectCoordSys,
                Vector<LevelData<FArrayBox>* >& a_initialVelocity);

  //void initDiagnostics();

  virtual void defineSolver();

  virtual void solveVelocityField(bool a_forceSolve = false, Real a_convergenceMetric = -1.0);

  virtual void defineVelRHS(Vector<LevelData<FArrayBox>* >& a_vectRhs);

  virtual void setBasalFriction(Vector<LevelData<FArrayBox>* >& a_C, Vector<LevelData<FArrayBox>* >& a_C0);

  virtual void setMuCoefficient(Vector<LevelData<FArrayBox>* >& a_cellMuCoef);

  void computeFaceVelocity(Vector<LevelData<FluxBox>* >& a_faceVelAdvection,
                           Vector<LevelData<FluxBox>* >& a_faceVelTotal,
                           Vector<LevelData<FluxBox>* >& a_diffusivity,
                           Vector<LevelData<FluxBox>* >& a_layerXYFaceXYVel,
                           Vector<LevelData<FArrayBox>* >& a_layerSFaceXYVel);


  // /// compute div(vel*H) at a given time
  // void computeDivThicknessFlux(Vector<LevelData<FArrayBox>* >& a_divFlux,
  //                              Vector<LevelData<FluxBox>* >& a_flux,
  //                              Vector<LevelData<FArrayBox>* >& a_thickness,
  //                              Real a_time, Real a_dt);

  // increment phi := phi + dt*dphi
  void incrementWithDivFlux(Vector<LevelData<FArrayBox>* >& a_phi,
                            const Vector<LevelData<FArrayBox>* >& a_dphi,
                            Real a_dt);

  // incremente coordSys with new thickness
  void updateCoordSysWithNewThickness(const Vector<LevelData<FArrayBox>* >& a_thickness);


  void setIceFrac(const LevelData<FArrayBox>& a_thickness, int a_level);



  void updateIceFrac(LevelData<FArrayBox>& a_thickness, int a_level);

  void setIceFracBC(LevelData<FArrayBox>& a_iceFrac, const ProblemDomain& a_domain);
  
  void advectIceFrac(Vector<LevelData<FArrayBox>* >& a_iceFrac,
                     const Vector<LevelData<FluxBox>* >& a_faceVelAdvection,  
                     Real a_dt);


  // /// update real-valued ice fraction through advection from neighboring cells
  // /**
  //    Also increment a_thicknessSource based on the calving rate
  // */
  // // void advectIceFrac(Vector<LevelData<FArrayBox>* >& a_iceFrac,
  // //                      Vector<LevelData<FArrayBox>* >& a_thicknessSource,
  // //                    const Vector<LevelData<FluxBox>* >& a_faceVelAdvection,
  // //                      const Vector<LevelData<FluxBox>* >& a_thicknessFlux,
  // //                    Real a_dt);
  
  // /// update real-valued ice fraction through application of the calving rate
  // /**
  //    Also increment a_thicknessSource based on the calving rate
  // */
  // // void advectIceFracUpstream(Vector<LevelData<FArrayBox>* >& a_iceFrac,
  // //                              Vector<LevelData<FArrayBox>* >& a_thicknessSource,
  // //                              const Vector<LevelData<FluxBox>* >& a_faceVelAdvection,
  // //                              const Vector<LevelData<FArrayBox> *>& a_calvingRate,
  // //                              Real a_subDt, Real a_dt);

  // // /// update real-valued ice fraction through advection from neighboring cell
  // // void advectIceFracDownstream(Vector<LevelData<FArrayBox>* >& a_iceFrac,
  // //                                Vector<LevelData<FArrayBox>* >& a_thicknessSource,
  // //                                const Vector<LevelData<FluxBox>* >& a_faceVelAdvection,
  // //                                const Vector<LevelData<FluxBox>* >& a_thicknessFlux,
  // //                                const Vector<LevelData<FArrayBox> *>& a_calvingRate,
  // //                                Real a_subDt, Real a_dt);

  void updateInvalidIceFrac(Vector<LevelData<FArrayBox> *> a_iceFrac);
  
  Real computeDt();

  Real computeInitialDt();

  //set an entire amr hierarchy to zero
  void setToZero(Vector<LevelData<FArrayBox>*>& a_data);

  void implicitThicknessCorrection(Real a_dt,
                                   const Vector<LevelData<FArrayBox>* >& a_sts,
                                   const Vector<LevelData<FArrayBox>* >& a_bts,
                                   const Vector<LevelData<FArrayBox>* >& a_vts
                                   );
  
  void applyCalvingCriterion(CalvingModel::Stage a_stage);
  
  void eliminateRemoteIce();

  void helmholtzSolve(Vector<LevelData<FArrayBox>* >& a_phi,
                      const Vector<LevelData<FArrayBox>* >& a_rhs, 
                      Real a_alpha, Real a_beta) const;

  void helmholtzSolve(Vector<LevelData<FArrayBox>* >& a_phi,
                      Real a_alpha, Real a_beta) const;
 

#if BISICLES_Z == BISICLES_LAYERED
 

  void computeA(Vector<LevelData<FArrayBox>* >& a_A, 
                        Vector<LevelData<FArrayBox>* >& a_sA,
                        Vector<LevelData<FArrayBox>* >& a_bA,
                        const Vector<LevelData<FArrayBox>* >& a_internalEnergy, 
                        const Vector<LevelData<FArrayBox>* >& a_sInternalEnergy,
                        const Vector<LevelData<FArrayBox>* >& a_bInternalEnergy,
                         const Vector<RefCountedPtr<LevelSigmaCS> >& a_coordSys) const;

  //compute the face- and layer- centered internal energy (a_layerEH_half)
  //and thickness (a_layerH_half) at time a_time + 1/2 * a_dt
  void computeInternalEnergyHalf(Vector<LevelData<FluxBox>* >& a_layerEH_half, 
                                 Vector<LevelData<FluxBox>* >& a_layerH_half,
                                 const Vector<LevelData<FluxBox>* >& a_layerXYFaceXYVel, 
                                 const Real a_dt, const Real a_time);

 
  
  void updateInternalEnergy(Vector<LevelData<FluxBox>* >& a_layerTH_half, 
                            Vector<LevelData<FluxBox>* >& a_layerH_half,
                            const Vector<LevelData<FluxBox>* >& a_layerXYFaceXYVel,
                            const Vector<LevelData<FArrayBox>* >& a_layerSFaceXYVel,
                            const Real a_dt, const Real a_time, 
                            Vector<RefCountedPtr<LevelSigmaCS> >& a_coordSysNew,
                            Vector<RefCountedPtr<LevelSigmaCS> >& a_coordSysOld,
                            const Vector<LevelData<FArrayBox>*>& a_surfaceThicknessSource,
                            const Vector<LevelData<FArrayBox>*>& a_basalThicknessSource,
                            const Vector<LevelData<FArrayBox>*>& a_volumeThicknessSource);
  void updateTemperature();
#endif


  static void postInterpolationReFloat(LevelData<FArrayBox>& a_H,
                                       const LevelData<FArrayBox>& a_coarseH,
                                       const LevelData<FArrayBox>& a_coarseBed,
                                       const DisjointBoxLayout a_newDBL,
                                       const ProblemDomain& a_domain,
                                       int a_refRatio,
                                       Real a_seaLevel, 
                                       Real a_waterDensity, 
                                       Real a_iceDensity);

  // diagnostic routine -- compute discharge
  void computeDischarge(const Vector<LevelData<FluxBox>* >& a_vectFluxes);

  // diagnostic routine 
  //void endTimestepDiagnostics();



  int m_temporalAccuracy;
  

  int m_num_thickness_ghost;

  int m_solverType;

  int m_maxSolverIterations;

  Real  m_velocity_solver_tolerance;

  int m_velocity_solve_interval;

  Real  m_velocitySolveInitialResidualNorm;
  Real  m_velocitySolveFinalResidualNorm;
  IceVelocitySolver* m_velSolver;
  

  ConstitutiveRelation* m_constitutiveRelation;
  RateFactor* m_rateFactor;
  BasalFrictionRelation* m_basalFrictionRelation;
  RateFactor* m_basalRateFactor;


  bool  m_reset_floating_friction_to_zero;

  GodunovPhysics* m_thicknessPhysPtr;

  Vector<PatchGodunov*> m_thicknessPatchGodVect;

  IceThicknessIBC* m_thicknessIBCPtr;
 
  IceInternalEnergyIBC* m_internalEnergyIBCPtr;

  CalvingModel* m_calvingModelPtr;

  SurfaceFlux* m_surfaceFluxPtr;
  SurfaceFlux* m_basalFluxPtr;
  
  SurfaceFlux* m_surfaceHeatBoundaryDataPtr;
  SurfaceFlux* m_basalHeatBoundaryDataPtr;
  
  SurfaceFlux* m_topographyFluxPtr;
  Vector<LevelData<FArrayBox>*> m_deltaTopography;

  BasalFriction* m_basalFrictionPtr;

  MuCoefficient* m_muCoefficientPtr;
  
  int m_max_level;
  
  int m_finest_level;
  
  int m_finest_timestep_level;
  int finestTimestepLevel() const
  {
    return (0 <= m_finest_timestep_level)?
      std::min(m_finest_timestep_level,m_finest_level):m_finest_level;
  }


  int m_tag_cap;

  int m_block_factor;

  IntVectSet m_tag_subset;

  Vector<IntVectSet> m_vectTagSubset;

  Real m_fill_ratio;

  int m_nesting_radius;

  int m_max_box_size;


  int m_max_base_grid_size;

  int m_regrid_interval;
  int m_n_regrids;

  bool m_doInitialVelSolve;


  bool m_doInitialVelGuess;
  // types of initial guess
  enum initialGuessTypes {SlidingLaw=0, // sliding law with the form Cu = rhs 
                          ConstMu=1, // solve the elliptic equations with a supplied mu
                          Function=2, //use a RealFunction<RealVect> 
                          NUM_INITIAL_GUESS_TYPES};
  int m_initialGuessType;

  //paramamters which specify the linear solver in the case m_initialGuessType = ConstMu
  Real m_initialGuessConstMu; // domain wide viscosity
  int m_initialGuessSolverType; // solver type, defaults to the Picard solver
  RealVect m_initialGuessConstVel; // initial velocity, useful when the basal friction relation is nonlinear


  bool m_interpolate_zb;
  int m_regrid_thickness_interpolation_method;

  // if > 0, then the coefficient of friction at the base
  // of the ice and the RHS will be smoothed.
  Real m_basalLengthScale;
 
  Vector<IntVectSet> m_vectTags; 

  bool m_tagOnGradVel;
  
  Real m_tagging_val;

  bool m_tagOnLapVel;

  bool m_tagOnGroundedLapVel;

  Real m_laplacian_tagging_val;
  Real m_laplacian_tagging_max_basal_friction_coef;
  bool m_tagOnEpsSqr;


  Real m_epsSqr_tagVal;

  bool m_tagOnVelRHS;

  Real m_velRHS_tagVal;

  bool m_tagOndivHgradVel;

  Real m_divHGradVel_tagVal;

  bool m_tagGroundingLine;
  Real m_groundingLineTaggingMinVel;
  Real m_groundingLineTaggingMaxBasalFrictionCoef;
  
  bool m_tagVelDx;
  Real m_velDx_tagVal;
  int m_velDx_tagVal_finestLevelGrounded;
  int m_velDx_tagVal_finestLevelFloating;

  // tag thin ice shelf cavities. Refinining thin cavities 
  // sometimes reveals pinning points
  bool m_tag_thin_cavity;
  Real m_tag_thin_cavity_thickness;

  bool m_tagMargin;
  int m_margin_tagVal_finestLevel;

  bool m_tagAllIce;


  bool m_tagAllIceOnLevel0;
  
  bool m_tagEntireDomain;


  //attempt to eliminate floating ice that is unconnected to grounded ice 
  bool m_eliminate_remote_ice;
  int m_eliminate_remote_ice_max_iter;
  Real m_eliminate_remote_ice_tol;

  //attempt to eliminate unconnected floating ice after regridding
  bool m_eliminate_remote_ice_after_regrid;
  
#ifdef HAVE_PYTHON
  bool m_tagPython;
  PyObject* m_tagPythonModule;
  PyObject* m_tagPythonFunction;
#endif
  // grow tags in all directions by m_tags_grow
  int m_tags_grow; 
  //grow tags in all direction dir by max(0,m_tags_grow_dir[dir] -  m_tags_grow)
  IntVect m_tags_grow_dir; 

  Vector<int>  m_refinement_ratios;

  Vector<Real> m_amrDx;

  Vector<ProblemDomain> m_amrDomains;
  
  // problem domain size
  RealVect m_domainSize;

  Vector<DisjointBoxLayout> m_amrGrids;

  Vector<int> m_covered_level;

  Vector<int> m_num_cells;

  Real m_time; Real m_offsetTime;

  Real m_dt;

  Real m_stable_dt;

  Real m_cfl;

  Real m_initial_cfl;

  Real m_max_dt_grow;

  Real m_velocity_exit;

  // current step
  int m_cur_step;

  // isothermal model ?
  bool m_isothermal;
 
  //depth of water in surface crevasses
  Real m_waterDepth;

  bool m_surfaceBoundaryHeatDataDirichlett;
  bool m_surfaceBoundaryHeatDataTemperature;

  //ice density
  Real m_iceDensity;

  //sea water density
  Real m_seaWaterDensity;

  // acceleration due to gravity
  Real m_gravity;

  // size of unit time (for velocity & time step *only*, other quantities are SI) in seconds
  Real m_seconds_per_unit_time;
  
  Vector<LevelData<FArrayBox>*> m_old_thickness;

 

  Vector<LevelData<FArrayBox>*> m_velocity;

  Vector<LevelData<FArrayBox>* > m_velRHS;

  // coeffient of |u|^m-1 u in the basal traction
  Vector<LevelData<FArrayBox>* > m_velBasalC; 

  //coefficient of mu in the stress tensor
  Vector<LevelData<FArrayBox>* > m_cellMuCoef;
  //Vector<LevelData<FluxBox>* > m_faceMuCoef;

  // thickness diffusion coefficient
  Vector<LevelData<FluxBox>* > m_diffusivity;
  Real m_additionalDiffusivity;
  
  //face velocities
  Vector<LevelData<FluxBox>* > m_faceVelAdvection;
  Vector<LevelData<FluxBox>* > m_faceVelTotal; // m_faceVelTotal = m_faceVelAdvection - m_diffusivity * grad(H)
#if BISICLES_Z == BISICLES_LAYERED
  Vector<LevelData<FArrayBox>* > m_layerSFaceXYVel;
  Vector<LevelData<FArrayBox>* > m_layerSFaceSVel;
  Vector<LevelData<FluxBox>* > m_layerXYFaceXYVel;
#endif

 

  // SigmaCS mapping
  Vector<RefCountedPtr<LevelSigmaCS > > m_vect_coordSys;
  

  Vector<LevelData<FArrayBox>* > m_iceFrac;

  // Cell-centered internalEnergy field
  Vector<LevelData<FArrayBox>* > m_internalEnergy, m_temperature;
  //cell-centered till water depth (accompanies the internal energy)
  Vector<LevelData<FArrayBox>* > m_tillWaterDepth;
  // Cell-centered A field
  Vector<LevelData<FArrayBox>* > m_A;
  // used to decide whether A needs updating (after a regrid, or a internalEnergy update)
  mutable bool m_A_valid;
#if BISICLES_Z == BISICLES_LAYERED
  // InternalEnergy and A at the top and bottom surfaces of 
  // the ice, horizontally cell centred (but at vertical faces)
  Vector<LevelData<FArrayBox>* > m_sInternalEnergy, m_sTemperature ;
  Vector<LevelData<FArrayBox>* > m_bInternalEnergy, m_bTemperature;
  Vector<LevelData<FArrayBox>* > m_sHeatFlux;
  Vector<LevelData<FArrayBox>* > m_bHeatFlux;
  Vector<LevelData<FArrayBox>* > m_sA;
  Vector<LevelData<FArrayBox>* > m_bA;
#elif BISICLES_Z == BISICLES_FULL
  // possibly won't even need surface and basal temperature
  // field, but if we do they will need to be 2D
#endif
  //the thickness source will have three components, from the surface and base
  //and distributed throught the volume.
  //Although they are combined for thickness evolution, they must
  //be kept seperate for the internalEnergy calculation. 
  mutable Vector<LevelData<FArrayBox>* > m_surfaceThicknessSource;
  mutable Vector<LevelData<FArrayBox>* > m_basalThicknessSource;
  mutable Vector<LevelData<FArrayBox>* > m_volumeThicknessSource;

  void setStableSources(FArrayBox& a_sts,FArrayBox& a_bts, FArrayBox& a_vts,
                        const FArrayBox& a_divuh, const BaseFab<int>& a_mask,
                        const Box& a_box) const;
  
  Vector<LevelData<FArrayBox>* > m_divThicknessFlux; // set to -dH/dt
  //thickness of calved ice 
  Vector<LevelData<FArrayBox>* > m_calvedIceThickness;
  Vector<LevelData<FArrayBox>* > m_removedIceThickness;
  Vector<LevelData<FArrayBox>* > m_addedIceThickness;
  Vector<LevelData<FArrayBox>* > m_calvedThicknessSource;


  //cell centered storage for the stress-balance equation coefficients
  mutable Vector<LevelData<FArrayBox>* > m_dragCoef;
  mutable Vector<LevelData<FArrayBox>* > m_viscosityCoefCell;

  //cell centered storage for the viscous tensor components
  mutable Vector<LevelData<FArrayBox>* > m_viscousTensorCell;
  //face centered storage for the viscous tensor components
  mutable Vector<LevelData<FluxBox>* > m_viscousTensorFace;
  // cache validity flag
  mutable bool m_viscousTensor_valid;
  void updateViscousTensor() const;

 
  

  //Scalar field that varies in the ice shelf from 0 at the grounding line
  mutable Vector<LevelData<FArrayBox>* > m_groundingLineProximity;
  // cache validity flag;
  mutable bool m_groundingLineProximity_valid;
  // scale length for grounding line proximity
  Real m_groundingLineProximityScale;
  int m_groundingLineProximityCalcType;
  void updateGroundingLineProximity() const;
  
  bool m_report_discharge;

  Real m_report_time_interval;
  Real m_next_report_time;

  int m_beta_type;

  Real m_betaVal;

  Real m_betaEps;


  RealVect m_basalSlope;


  static int s_verbosity;

  bool m_is_defined;
  
  bool m_do_restart;

  int m_restart_step;

  string m_plot_prefix;
  
  string m_check_prefix;

  bool m_check_overwrite; // overwrite check point files
  bool m_check_exit; // exit immediately after writing checkpoint?


  int m_plot_interval;  
  Real m_plot_time_interval;

  bool m_reduced_plot;

  bool m_write_mask;

  bool m_write_dHDt;

  // if true, include (cell-averaged) flux velocity in plotfiles
  bool m_write_fluxVel;

  // if true, include (cell-averaged) viscous tensor components in plotfiles
  bool m_write_viscousTensor;

  // if true, include base velocity in plotfiles
  bool m_write_baseVel;

  bool m_write_solver_rhs;

  //if true, include internalEnergy in plot files
  bool m_write_internal_energy;

  //if true, include thickness source terms in plot files
  bool m_write_thickness_sources;

  //if true, include fields required for ISMIP6 in plot files
  bool m_write_ismip6;

  bool m_write_map_file;
  

  bool m_write_presolve_plotfiles;

  int m_check_interval;

  // if true, allow thickness to change in time
  bool m_evolve_thickness; 

  // if true, allow ice frac to change in time
  bool m_evolve_ice_frac;
  
  //if true, allow velocity field to change in time
  bool m_evolve_velocity;

  // if true, bedrock topogography evolves rather than surface in grounded ice
  bool m_evolve_topography_fix_surface;

  //if true, keep floating or grounded ice stable;
  bool m_grounded_ice_stable ;
  bool m_floating_ice_stable ;
  //if true, use the basal flux to set dhdt
  bool m_floating_ice_basal_flux_is_dhdt;
  bool m_floating_ice_basal_flux_is_min_dhdt;
  bool m_grounded_ice_basal_flux_is_dhdt;

  bool m_frac_sources;


  enum DiffusionTreatment {NONE,IMPLICIT,EXPLICIT};
  DiffusionTreatment m_diffusionTreatment;
 
  // if true, time step are always an integer power of two
  bool m_timeStepTicks;

  //ised to store a fixed time step > 0 , if desired
  Real m_fixed_dt;

  // if true, report sum of grounded ice (volume) -- default is false
  bool m_reportGroundedIce;




#if BISICLES_Z == BISICLES_LAYERED
  //additional members used in poor-mans multidim mode
public:
  const Vector<Real>& getFaceSigma() const
  {
    return m_faceSigma;
  }
protected:
  void setLayers(const Vector<Real>& a_sigma);


  bool m_sigmaSet;
  
  // Number layers within the ice sheet 
  unsigned int m_nLayers;

  //arrangment of layers within the bedrock
  Vector<Real> m_faceSigma;
  
  //if true, include layer velocities in plot files
  bool m_write_layer_velocities;
  bool m_additionalVelocity;

#endif
private: 
  template <class T>
  void levelAllocate(LevelData<T>** a_ptrPtr, const DisjointBoxLayout& a_grids, 
                     int a_nComp, const IntVect& a_ghost)
  {
    CH_assert(a_ptrPtr != NULL);
    if (*a_ptrPtr != NULL)
      delete *a_ptrPtr;
    *a_ptrPtr = new LevelData<T>(a_grids, a_nComp, a_ghost);
  }


public:

  class Observer {

    friend class AmrIce;

    AmrIce* m_observee;

    // setObservee is private so that only friends (AmrIce) can access
    void setObservee(AmrIce* a_amrIce)
    {
      m_observee = a_amrIce;
    }

    // clearObservee is private so that only friends (AmrIce) can access
    void clearObservee()
    {
      m_observee = NULL;
    }

  public:

    enum Notification {PostInitGrids,PreVelocitySolve,PostVelocitySolve,PostGeometryUpdate,PreCalving,PostCalving};
    
    virtual ~Observer()
    {
      if (m_observee != NULL)
        {
          m_observee->removeObserver(this);
        }
    }

    virtual void notify(const Notification, AmrIce&)=0;



#ifdef CH_USE_HDF5
    virtual void addPlotVars(Vector<std::string>& a_vars){;}
    
    virtual void writePlotData(LevelData<FArrayBox>& a_data, int a_level){;}

    virtual void addCheckVars(Vector<std::string>& a_vars){;}
    
    virtual void writeCheckData(HDF5Handle& a_handle, int a_level){;}

    virtual void readCheckData(HDF5Handle& a_handle, HDF5HeaderData& a_header, int a_level, const DisjointBoxLayout& a_grids){;}
#endif
  };

  void addObserver(Observer* a_observer)
  {
    // add a_observer to the list so long as it isn't there already
    if (std::find(m_observers.begin(), m_observers.end(), a_observer) == m_observers.end())
      {
        m_observers.push_back(a_observer);
        a_observer->setObservee(this);
      }
  }

  void removeObserver(Observer* a_observer)
  {
    // remove observer from the list
    std::vector<Observer*>::iterator it = std::find(m_observers.begin(), m_observers.end(), a_observer);
    if ( it != m_observers.end())
      {
        a_observer->clearObservee();
        m_observers.erase(it);
      }
  }


private:

  std::vector<Observer*> m_observers; //needs to be a std::vector for begin(), end()

  void notifyObservers(const Observer::Notification a_n)
  {
    for (int i = 0; i < m_observers.size(); i++)
      m_observers[i]->notify(a_n, *this);
  }

};  
  
#include "NamespaceFooter.H"

#endif