cell.cpp

Go to the documentation of this file.

 
/* ===============================================================================================================================
 
   This file is part of CoastalME, the Coastal Modelling Environment.
 
   CoastalME is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 3 of the License, or (at your option) any later
version.
 
   This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave,
Cambridge, MA 02139, USA.
 
===============================================================================================================================*/
#include <vector>
using std::vector;
 
#include <assert.h>
 
#include "cell.h"
#include "cme.h"
#include "simulation.h"

CGeomCell::CGeomCell()
    : m_bInContiguousSea(false), m_bInContiguousFlood(false),
      m_bIsInActiveZone(false), m_bCoastline(false), m_bCliff(false),
      m_bFloodLine(false), m_bWaveFlood(false),
      // m_bCheckCell(false),
      m_bCheckFloodCell(false), m_bShadowBoundary(false),
      m_bPossibleCoastStartCell(false), m_bPossibleFloodStartCell(false),
      m_nBoundingBoxEdge(NO_DIRECTION), m_nPolygonID(INT_NODATA),
      m_nPolygonCoastID(INT_NODATA), m_nCoastlineNormal(INT_NODATA),
      m_nShadowZoneNumber(0), m_nDownDriftZoneNumber(0), m_dLocalConsSlope(0),
      m_dBasementElevation(0), m_dSlope(0), m_dSeaDepth(0), m_dTotSeaDepth(0),
      m_dWaveHeight(0), m_dTotWaveHeight(0), m_dWaveAngle(DBL_NODATA),
      m_dTotWaveAngle(DBL_NODATA), m_dDeepWaterWaveHeight(DBL_NODATA),
      m_dDeepWaterWaveAngle(DBL_NODATA), m_dBeachProtectionFactor(DBL_NODATA),
      m_dSuspendedSediment(0), m_dTotSuspendedSediment(0),
      m_dPotentialPlatformErosionThisIter(0), m_dTotPotentialPlatformErosion(0),
      m_dActualPlatformErosionThisIter(0), m_dTotActualPlatformErosion(0),
      m_dCliffCollapseFineThisIter(0), m_dCliffCollapseSandThisIter(0),
      m_dCliffCollapseCoarseThisIter(0), m_dTotFineCliffCollapse(0),
      m_dTotSandCliffCollapse(0), m_dTotCoarseCliffCollapse(0),
      m_dTalusSandDepositionThisIter(0), m_dTotTalusSandDeposition(0),
      m_dTalusCoarseDepositionThisIter(0), m_dTotTalusCoarseDeposition(0),
      m_dPotentialBeachErosionThisIter(0), m_dTotPotentialBeachErosion(0),
      m_dActualBeachErosionThisIter(0), m_dTotActualBeachErosion(0),
      m_dBeachDepositionThisIter(0), m_dTotBeachDeposition(0), m_dUnconsD50(0),
      m_dInterventionHeight(0) {
  m_Landform.SetLFCategory(LF_NONE);
}

CGeomCell::~CGeomCell(void) {}

void CGeomCell::SetBoundingBoxEdge(int const nDirection) {
  m_nBoundingBoxEdge = nDirection;
}

int CGeomCell::nGetBoundingBoxEdge(void) const { return m_nBoundingBoxEdge; }

bool CGeomCell::bIsBoundingBoxEdge(void) const {
  return (m_nBoundingBoxEdge != NO_DIRECTION);
}

void CGeomCell::SetInContiguousSea(void) { m_bInContiguousSea = true; }

bool CGeomCell::bIsInContiguousSea(void) const { return m_bInContiguousSea; }

void CGeomCell::SetInContiguousFlood(void) { m_bInContiguousFlood = true; }

void CGeomCell::UnSetInContiguousFlood(void) { m_bInContiguousFlood = false; }

void CGeomCell::SetFloodBySetupSurge(void) { m_bFloodBySetupSurge = true; }

bool CGeomCell::bIsFloodBySetupSurge(void) const {
  return m_bFloodBySetupSurge;
}

void CGeomCell::SetFloodBySetupSurgeRunup(void) {
  m_bFloodBySetupSurgeRunup = true;
}

bool CGeomCell::bIsFloodBySetupSurgeRunup(void) const {
  return m_bFloodBySetupSurgeRunup;
}

bool CGeomCell::bIsInContiguousSeaArea(void) const {
  return m_bInContiguousFlood;
}

void CGeomCell::SetInActiveZone(bool const bFlag) { m_bIsInActiveZone = bFlag; }

bool CGeomCell::bIsInActiveZone(void) const { return m_bIsInActiveZone; }

void CGeomCell::SetShadowZoneBoundary(void) { m_bShadowBoundary = true; }

bool CGeomCell::bIsShadowZoneBoundary(void) const { return m_bShadowBoundary; }

void CGeomCell::SetPossibleCoastStartCell(void) {
  m_bPossibleCoastStartCell = true;
}

bool CGeomCell::bIsPossibleCoastStartCell(void) const {
  return m_bPossibleCoastStartCell;
}

void CGeomCell::SetPossibleFloodStartCell(void) {
  m_bPossibleFloodStartCell = true;
}

bool CGeomCell::bIsPossibleFloodStartCell(void) const {
  return m_bPossibleFloodStartCell;
}

bool CGeomCell::bPotentialPlatformErosion(void) const {
  return (m_dPotentialPlatformErosionThisIter > 0);
}
 
// bool CGeomCell::bActualPlatformErosion(void) const
// {
// return (m_dActualPlatformErosionThisIter > 0);
// }

void CGeomCell::SetAsCoastline(bool const bNewFlag) { m_bCoastline = bNewFlag; }

bool CGeomCell::bIsCoastline(void) const { return m_bCoastline; }

void CGeomCell::SetAsCliff(bool const bNewFlag) { m_bCliff = bNewFlag; }

bool CGeomCell::bIsCliff(void) const { return m_bCliff; }

void CGeomCell::SetAsFloodLine(bool const bNewFlag) { m_bFloodLine = bNewFlag; }

bool CGeomCell::bIsFloodLine(void) const { return m_bFloodLine; }

void CGeomCell::SetProfileID(int const nNormal) {
  m_nCoastlineNormal = nNormal;
}

int CGeomCell::nGetProfileID(void) const { return m_nCoastlineNormal; }

bool CGeomCell::bIsProfile(void) const {
  if (m_nCoastlineNormal == INT_NODATA)
    return false;
 
  return true;
}

void CGeomCell::SetPolygonID(int const nPolyID) { m_nPolygonID = nPolyID; }

int CGeomCell::nGetPolygonID(void) const { return m_nPolygonID; }

void CGeomCell::SetPolygonCoastID(int const nPolyCoastID) {
  m_nPolygonCoastID = nPolyCoastID;
}

int CGeomCell::nGetPolygonCoastID(void) const { return m_nPolygonCoastID; }

void CGeomCell::SetCoastAndPolygonID(int const nPolyCoastID,
                                     int const nPolyID) {
  m_nPolygonCoastID = nPolyCoastID;
  m_nPolygonID = nPolyID;
}

void CGeomCell::SetShadowZoneNumber(int const nCode) {
  m_nShadowZoneNumber = nCode;
}

int CGeomCell::nGetShadowZoneNumber(void) const { return m_nShadowZoneNumber; }

bool CGeomCell::bIsinThisShadowZone(int const nZone) const {
  if (m_nShadowZoneNumber == nZone)
    return true;
 
  return false;
}

bool CGeomCell::bIsinAnyShadowZone(void) const {
  if (m_nShadowZoneNumber != 0)
    return true;
 
  return false;
}

void CGeomCell::SetWaveFlood(void) { m_bWaveFlood = true; }
 
// void CGeomCell::SetWaveSetup(int const dWaveSetup)
// {
// m_dWaveSetup = dWaveSetup;
// }
 
// void CGeomCell::SetStormSurge(int const dStormSurge)
// {
// m_dStormSurge = dStormSurge;
// }
 
// void CGeomCell::SetRunUp(int const dRunUp)
// {
// m_dRunUp = dRunUp;
// }
 
// void CGeomCell::SetTotLevel(void) const
// {
// m_dTotLevel = m_dWaveSetup + m_dStormSurge + m_dRunUp +
// m_pGrid->pGetSim()->CSimulation::dGetThisIterSWL();
// }
 
// int CGeomCell::nGetTotLevel(void) const
// {
// return m_dTotLevel;
// }

bool CGeomCell::bIsElevLessThanWaterLevel(void) const {
  return ((m_VdAllHorizonTopElev.back() + m_dInterventionHeight) <
          (m_pGrid->pGetSim()->dGetThisIterTotWaterLevel() +
           m_pGrid->pGetSim()->dGetThisIterSWL()));
}
 
// //! Set this cell as checked TODO What is this used for?
// void CGeomCell::SetCheckCell(void)
// {
// m_bCheckCell = true;
// }
 
// //! Returns true if this cell is checked, false otherwise TODO 007 What is
// this used for? bool CGeomCell::bIsCellCheck(void) const
// {
// return m_bCheckCell;
// }

void CGeomCell::SetCheckFloodCell(void) { m_bCheckFloodCell = true; }

void CGeomCell::UnSetCheckFloodCell(void) { m_bCheckFloodCell = false; }

bool CGeomCell::bIsCellFloodCheck(void) const { return m_bCheckFloodCell; }

void CGeomCell::SetDownDriftZoneNumber(int const nCode) {
  m_nDownDriftZoneNumber = nCode;
}

int CGeomCell::nGetDownDriftZoneNumber(void) const {
  return m_nDownDriftZoneNumber;
}

CRWCellLandform *CGeomCell::pGetLandform(void) { return &m_Landform; }

void CGeomCell::SetLocalConsSlope(double const dNewSlope) {
  m_dLocalConsSlope = dNewSlope;
}

double CGeomCell::dGetLocalConsSlope(void) const { return m_dLocalConsSlope; }

void CGeomCell::SetBasementElev(double const dNewElev) {
  m_dBasementElevation = dNewElev;
}

double CGeomCell::dGetBasementElev(void) const {
  return (m_dBasementElevation);
}

void CGeomCell::SetSlope(double const dNewSlope) { m_dSlope = dNewSlope; }

double CGeomCell::dGetSlope(void) const { return (m_dSlope); }

bool CGeomCell::bBasementElevIsMissingValue(void) const {
  if (bFPIsEqual(m_dBasementElevation,
                 m_pGrid->pGetSim()->CSimulation::dGetMissingValue(),
                 TOLERANCE))
    return true;
 
  return false;
}

double CGeomCell::dGetSeaDepth(void) const { return (m_dSeaDepth); }

double CGeomCell::dGetTotSeaDepth(void) const { return (m_dTotSeaDepth); }

void CGeomCell::SetSuspendedSediment(double const dNewSedDepth) {
  // Note no checks here to see if new equiv depth is sensible (e.g.
  // non-negative)
  m_dSuspendedSediment = dNewSedDepth;
  m_dTotSuspendedSediment += dNewSedDepth;
}

void CGeomCell::AddSuspendedSediment(double const dIncSedDepth) {
  // Note no checks here to see if increment equiv depth is sensible (e.g.
  // non-negative)
  m_dSuspendedSediment += dIncSedDepth;
  m_dTotSuspendedSediment += dIncSedDepth;
}

double CGeomCell::dGetSuspendedSediment(void) const {
  return (m_dSuspendedSediment);
}

double CGeomCell::dGetTotSuspendedSediment(void) const {
  return (m_dTotSuspendedSediment);
}

int CGeomCell::nGetTopNonZeroLayerAboveBasement(void) const {
  if (m_VLayerAboveBasement.empty())
    return INT_NODATA;
 
  int nTop = static_cast<int>(m_VLayerAboveBasement.size()) - 1;
 
  while (m_VLayerAboveBasement[nTop].dGetTotalThickness() <= 0) {
    if (--nTop < 0)
      return NO_NONZERO_THICKNESS_LAYERS;
  }
 
  return nTop;
}

int CGeomCell::nGetTopLayerAboveBasement(void) const {
  if (m_VLayerAboveBasement.empty())
    return INT_NODATA;
 
  return static_cast<int>(m_VLayerAboveBasement.size()) - 1;
}

double CGeomCell::dGetConsSedTopForLayerAboveBasement(int const nLayer) const {
  // Note no check to see if nLayer < m_VLayerAboveBasement.size()
  double dTopElev = m_dBasementElevation;
 
  for (int n = 0; n < nLayer; n++) {
    dTopElev += m_VLayerAboveBasement[n].dGetUnconsolidatedThickness();
    dTopElev += m_VLayerAboveBasement[n].dGetConsolidatedThickness();
  }
 
  dTopElev += m_VLayerAboveBasement[nLayer].dGetConsolidatedThickness();
 
  return dTopElev;
}

CRWCellLayer *CGeomCell::pGetLayerAboveBasement(int const nLayer) {
  // TODO 055 No check that nLayer < size()
  return &m_VLayerAboveBasement[nLayer];
}
 
// //! Returns the volume-equivalent elevation of the sediment's top surface for
// this cell (i.e. if there is a cliff notch, then lower the elevation by the
// notch's volume) double CGeomCell::dGetVolEquivSedTopElev(void) const
// {
// double dTopElev = m_dBasementElevation;
// for (unsigned int n = 0; n < m_VLayerAboveBasement.size(); n++)
// {
// dTopElev += (m_VLayerAboveBasement[n].dGetUnconsolidatedThickness() -
// m_VLayerAboveBasement[n].dGetNotchUnconsolidatedLost()); dTopElev +=
// (m_VLayerAboveBasement[n].dGetConsolidatedThickness() -
// m_VLayerAboveBasement[n].dGetNotchConsolidatedLost());
// }
//
// return dTopElev;
// }

double CGeomCell::dGetSedimentTopElev(void) const {
  return m_VdAllHorizonTopElev.back();
}

double CGeomCell::dGetSedimentPlusInterventionTopElev(void) const {
  return m_VdAllHorizonTopElev.back() + m_dInterventionHeight;
}

double CGeomCell::dGetOverallTopElev(void) const {
  return m_VdAllHorizonTopElev.back() + m_dInterventionHeight + m_dSeaDepth;
}

bool CGeomCell::bIsInundated(void) const {
  return ((m_VdAllHorizonTopElev.back() + m_dInterventionHeight) <
          m_pGrid->pGetSim()->CSimulation::dGetThisIterSWL());
}

double CGeomCell::dGetThisIterSWL(void) const {
  return m_pGrid->pGetSim()->CSimulation::dGetThisIterSWL();
}

double CGeomCell::dGetThisIterTotWaterLevel(void) const {
  return m_pGrid->pGetSim()->CSimulation::dGetThisIterTotWaterLevel();
}
 
// //! Returns true if the elevation of the sediment top surface for this cell
// is greater than or equal to the grid's this-timestep still water elevation.
// Also returns true if the cell has unconsolidated sediment on it and the
// elevation of the sediment top surface, minus a tolerance value, is less than
// the grid's this-timestep still water elevation bool
// CGeomCell::bIsSeaIncBeach(void) const
// {
// if (m_bInContiguousSea)
//       // Sea
// return true;
//
// double
// dWaterLevel = m_pGrid->pGetSim()->CSimulation::dGetThisIterSWL(),
// dSedTop = m_VdAllHorizonTopElev.back();
//
//    // Beach
// if ((m_VLayerAboveBasement.back().dGetUnconsolidatedThickness() > 0) &&
// ((dSedTop - m_pGrid->pGetSim()->CSimulation::dGetMaxBeachElevAboveSWL()) <
// dWaterLevel)) return true;
//
// return false;
// }

double CGeomCell::dGetTotConsFineThickConsiderNotch(void) const {
  double dTotThick = 0;
 
  for (unsigned int n = 0; n < m_VLayerAboveBasement.size(); n++) {
    CRWCellLayer m_Layer = m_VLayerAboveBasement[n];
    double dLayerThick = m_Layer.dGetFineConsolidatedThickness();
    double dNotchEquiv =
        m_Layer.pGetConsolidatedSediment()->dGetNotchFineLost();
 
    dTotThick += (dLayerThick - dNotchEquiv);
  }
 
  return dTotThick;
}

double CGeomCell::dGetTotUnconsFine(void) const {
  double dTotThick = 0;
 
  for (unsigned int n = 0; n < m_VLayerAboveBasement.size(); n++)
    dTotThick += m_VLayerAboveBasement[n].dGetFineUnconsolidatedThickness();
 
  return dTotThick;
}

double CGeomCell::dGetTotConsSandThickConsiderNotch(void) const {
  double dTotThick = 0;
 
  for (unsigned int n = 0; n < m_VLayerAboveBasement.size(); n++) {
    CRWCellLayer m_Layer = m_VLayerAboveBasement[n];
    double dLayerThick = m_Layer.dGetSandConsolidatedThickness();
    double dNotchEquiv =
        m_Layer.pGetConsolidatedSediment()->dGetNotchSandLost();
 
    dTotThick += (dLayerThick - dNotchEquiv);
  }
 
  return dTotThick;
}

double CGeomCell::dGetTotUnconsSand(void) const {
  double dTotThick = 0;
 
  for (unsigned int n = 0; n < m_VLayerAboveBasement.size(); n++)
    dTotThick += m_VLayerAboveBasement[n].dGetSandUnconsolidatedThickness();
 
  return dTotThick;
}

double CGeomCell::dGetTotConsCoarseThickConsiderNotch(void) const {
  double dTotThick = 0;
 
  for (unsigned int n = 0; n < m_VLayerAboveBasement.size(); n++) {
    CRWCellLayer m_Layer = m_VLayerAboveBasement[n];
    double dLayerThick = m_Layer.dGetCoarseConsolidatedThickness();
    double dNotchEquiv =
        m_Layer.pGetConsolidatedSediment()->dGetNotchCoarseLost();
 
    dTotThick += (dLayerThick - dNotchEquiv);
  }
 
  return dTotThick;
}

double CGeomCell::dGetTotUnconsCoarse(void) const {
  double dTotThick = 0;
 
  for (unsigned int n = 0; n < m_VLayerAboveBasement.size(); n++)
    dTotThick += m_VLayerAboveBasement[n].dGetCoarseUnconsolidatedThickness();
 
  return dTotThick;
}

double CGeomCell::dGetTotConsThickness(void) const {
  double dTotThick = 0;
 
  for (unsigned int n = 0; n < m_VLayerAboveBasement.size(); n++)
    dTotThick += m_VLayerAboveBasement[n].dGetConsolidatedThickness();
 
  return dTotThick;
}

double CGeomCell::dGetTotUnconsThickness(void) const {
  double dTotThick = 0;
 
  for (unsigned int n = 0; n < m_VLayerAboveBasement.size(); n++)
    dTotThick += m_VLayerAboveBasement[n].dGetUnconsolidatedThickness();
 
  return dTotThick;
}

double CGeomCell::dGetTotAllSedThickness(void) const {
  return (this->dGetTotUnconsThickness() + this->dGetTotConsThickness());
}

void CGeomCell::AppendLayers(int const nLayer) {
  for (int i = 0; i < nLayer; i++)
    m_VLayerAboveBasement.push_back(CRWCellLayer());
}

void CGeomCell::CalcAllLayerElevsAndD50(void) {
  m_VdAllHorizonTopElev.clear();
  m_VdAllHorizonTopElev.push_back(
      m_dBasementElevation); // Elevation of top of the basement
 
  // Calculate the elevation of the top of all other layers
  int m = 0;
 
  for (unsigned int n = 0; n < m_VLayerAboveBasement.size(); n++)
    m_VdAllHorizonTopElev.push_back(
        m_VLayerAboveBasement[n].dGetTotalThickness() +
        m_VdAllHorizonTopElev[m++]); // Elevation of top of layer n
 
  // Now calculate the d50 of the topmost unconsolidated sediment layer with
  // non-zero thickness. If there is no unconsolidated sediment, m_dUnconsD50 is
  // set to DBL_NODATA
  m_dUnconsD50 = DBL_NODATA;
 
  for (int n = static_cast<int>(m_VLayerAboveBasement.size()) - 1; n >= 0;
       n--) {
    double dUnconsThick =
        m_VLayerAboveBasement[n].dGetUnconsolidatedThickness();
 
    if (dUnconsThick > 0) {
      // This is a layer with non-zero thickness of unconsolidated sediment
      CRWCellSediment const *pUnconsSedLayer =
          m_VLayerAboveBasement[n].pGetUnconsolidatedSediment();
      double dFineProp = pUnconsSedLayer->dGetFineDepth() / dUnconsThick;
      double dSandProp = pUnconsSedLayer->dGetSandDepth() / dUnconsThick;
      double dCoarseProp = pUnconsSedLayer->dGetCoarseDepth() / dUnconsThick;
 
      // Calculate d50 for the unconsolidated sediment
      m_dUnconsD50 = (dFineProp * m_pGrid->pGetSim()->dGetD50Fine()) +
                     (dSandProp * m_pGrid->pGetSim()->dGetD50Sand()) +
                     (dCoarseProp * m_pGrid->pGetSim()->dGetD50Coarse());
 
      break;
    }
  }
}

int CGeomCell::nGetLayerAtElev(double const dElev) const {
  if (dElev < m_VdAllHorizonTopElev[0])
    return ELEV_IN_BASEMENT;
 
  for (unsigned int nLayer = 1; nLayer < m_VdAllHorizonTopElev.size();
       nLayer++) {
    if ((m_VLayerAboveBasement[nLayer - 1].dGetTotalThickness() > 0) &&
        (dElev >= m_VdAllHorizonTopElev[nLayer - 1]) &&
        (dElev <= m_VdAllHorizonTopElev[nLayer]))
      return (nLayer - 1);
  }
 
  return ELEV_ABOVE_SEDIMENT_TOP;
}

double CGeomCell::dCalcLayerElev(const int nLayer) {
  // Note no check to see if nLayer < m_VLayerAboveBasement.size()
  double dTopElev = m_dBasementElevation;
 
  for (int n = 0; n <= nLayer; n++)
    dTopElev += m_VLayerAboveBasement[n].dGetTotalThickness();
 
  return dTopElev;
}

void CGeomCell::SetPotentialPlatformErosion(double const dPotentialIn) {
  m_dPotentialPlatformErosionThisIter = dPotentialIn;
  m_dTotPotentialPlatformErosion += dPotentialIn;
}

double CGeomCell::dGetPotentialPlatformErosion(void) const {
  return m_dPotentialPlatformErosionThisIter;
}

double CGeomCell::dGetTotPotentialPlatformErosion(void) const {
  return m_dTotPotentialPlatformErosion;
}

void CGeomCell::SetActualPlatformErosion(double const dThisActualErosion) {
  m_dActualPlatformErosionThisIter = dThisActualErosion;
  m_dTotActualPlatformErosion += dThisActualErosion;
}

double CGeomCell::dGetActualPlatformErosion(void) const {
  return m_dActualPlatformErosionThisIter;
}

double CGeomCell::dGetTotActualPlatformErosion(void) const {
  return m_dTotActualPlatformErosion;
}

void CGeomCell::SetSeaDepth(void) {
  m_dSeaDepth = tMax(m_pGrid->pGetSim()->CSimulation::dGetThisIterSWL() -
                         (m_VdAllHorizonTopElev.back() + m_dInterventionHeight),
                     0.0);
}

void CGeomCell::InitCell(void) {
  m_bInContiguousSea = m_bInContiguousFlood = m_bCoastline = m_bFloodLine =
      m_bIsInActiveZone =
          // m_bEstimated =
      m_bShadowBoundary = m_bPossibleCoastStartCell =
          m_bPossibleFloodStartCell = m_bWaveFlood =
              // m_bCheckCell =
      m_bCheckFloodCell = false;
 
  m_nPolygonID = m_nPolygonCoastID = m_nCoastlineNormal = INT_NODATA;
 
  m_nShadowZoneNumber = m_nDownDriftZoneNumber = 0;
 
  m_dLocalConsSlope = m_dPotentialPlatformErosionThisIter =
      m_dActualPlatformErosionThisIter = m_dCliffCollapseFineThisIter =
          m_dCliffCollapseSandThisIter = m_dCliffCollapseCoarseThisIter =
              m_dTalusSandDepositionThisIter =
                  m_dTalusCoarseDepositionThisIter =
                      m_dPotentialBeachErosionThisIter =
                          m_dActualBeachErosionThisIter =
                              m_dBeachDepositionThisIter = m_dSeaDepth =
                                  m_dWaveHeight = m_dWaveAngle = 0;
 
  m_dBeachProtectionFactor = DBL_NODATA;
 
  // Initialize this-iteration sediment input event values
  int nThisLayer = this->nGetTopNonZeroLayerAboveBasement();
 
  // Safety check
  if ((nThisLayer == NO_NONZERO_THICKNESS_LAYERS) ||
      (nThisLayer == INT_NODATA)) {
    // Uh-oh, problem. So just return TODO 021
    return;
  }
 
  m_VLayerAboveBasement[nThisLayer]
      .pGetUnconsolidatedSediment()
      ->InitThisIterSedimentInputAll();
}

void CGeomCell::SetWaveHeight(double const dWaveHeight) {
  m_dWaveHeight = dWaveHeight;
  m_dTotWaveHeight += dWaveHeight;
 
  // if (m_dWaveHeight != DBL_NODATA)
  // assert(m_dWaveHeight >= 0);
}

double CGeomCell::dGetWaveHeight(void) const { return m_dWaveHeight; }

double CGeomCell::dGetTotWaveHeight(void) const { return m_dTotWaveHeight; }

void CGeomCell::SetWaveAngle(double const dWaveAngle) {
  m_dWaveAngle = dWaveAngle;
  m_dTotWaveAngle += dWaveAngle;
}

double CGeomCell::dGetWaveAngle(void) const { return m_dWaveAngle; }

double CGeomCell::dGetTotWaveAngle(void) const { return m_dTotWaveAngle; }

void CGeomCell::SetCellDeepWaterWaveHeight(double const dWaveHeight) {
  m_dDeepWaterWaveHeight = dWaveHeight;
}

double CGeomCell::dGetCellDeepWaterWaveHeight(void) const {
  return m_dDeepWaterWaveHeight;
}

void CGeomCell::SetCellDeepWaterWaveAngle(double const dWaveAngle) {
  m_dDeepWaterWaveAngle = dWaveAngle;
}

double CGeomCell::dGetCellDeepWaterWaveAngle(void) const {
  return m_dDeepWaterWaveAngle;
}

void CGeomCell::SetCellDeepWaterWavePeriod(double const dWavePeriod) {
  m_dDeepWaterWavePeriod = dWavePeriod;
}

double CGeomCell::dGetCellDeepWaterWavePeriod(void) const {
  return m_dDeepWaterWavePeriod;
}

void CGeomCell::SetWaveValuesToDeepWaterWaveValues(void) {
  m_dWaveHeight = m_dDeepWaterWaveHeight;
  m_dWaveAngle = m_dDeepWaterWaveAngle;
  m_dWavePeriod = m_dDeepWaterWavePeriod;
}
 
// Sets this cell's beach protection factor
void CGeomCell::SetBeachProtectionFactor(double const dFactor) {
  m_dBeachProtectionFactor = dFactor;
}

double CGeomCell::dGetBeachProtectionFactor(void) const {
  return m_dBeachProtectionFactor;
}

void CGeomCell::IncrCliffCollapseErosion(double const dFineDepth,
                                         double const dSandDepth,
                                         double const dCoarseDepth) {
  m_dCliffCollapseFineThisIter += dFineDepth;
  m_dCliffCollapseSandThisIter += dSandDepth;
  m_dCliffCollapseCoarseThisIter += dCoarseDepth;
 
  m_dTotFineCliffCollapse += dFineDepth;
  m_dTotSandCliffCollapse += dSandDepth;
  m_dTotCoarseCliffCollapse += dCoarseDepth;
}

double CGeomCell::dGetThisIterCliffCollapseErosionFine(void) const {
  return m_dCliffCollapseFineThisIter;
}

double CGeomCell::dGetThisIterCliffCollapseErosionSand(void) const {
  return m_dCliffCollapseSandThisIter;
}

double CGeomCell::dGetThisIterCliffCollapseErosionCoarse(void) const {
  return m_dCliffCollapseCoarseThisIter;
}

double CGeomCell::dGetTotCliffCollapseFine(void) const {
  return m_dTotFineCliffCollapse;
}

double CGeomCell::dGetTotCliffCollapseSand(void) const {
  return m_dTotSandCliffCollapse;
}

double CGeomCell::dGetTotCliffCollapseCoarse(void) const {
  return m_dTotCoarseCliffCollapse;
}

void CGeomCell::AddSandTalusDeposition(double const dDepth) {
  m_dTalusSandDepositionThisIter += dDepth;
  m_dTotTalusSandDeposition += dDepth;
}

void CGeomCell::AddCoarseTalusDeposition(double const dDepth) {
  m_dTalusCoarseDepositionThisIter += dDepth;
  m_dTotTalusCoarseDeposition += dDepth;
}

double CGeomCell::dGetThisIterCliffCollapseSandTalusDeposition(void) const {
  return m_dTalusSandDepositionThisIter;
}

double CGeomCell::dGetThisIterCliffCollapseCoarseTalusDeposition(void) const {
  return m_dTalusCoarseDepositionThisIter;
}

double CGeomCell::dGetTotSandTalusDeposition(void) const {
  return m_dTotTalusSandDeposition;
}

double CGeomCell::dGetTotCoarseTalusDeposition(void) const {
  return m_dTotTalusCoarseDeposition;
}

void CGeomCell::SetPotentialBeachErosion(double const dPotentialIn) {
  m_dPotentialBeachErosionThisIter = dPotentialIn;
  m_dTotPotentialBeachErosion += dPotentialIn;
}

double CGeomCell::dGetPotentialBeachErosion(void) const {
  return m_dPotentialBeachErosionThisIter;
}

double CGeomCell::dGetTotPotentialBeachErosion(void) const {
  return m_dTotPotentialBeachErosion;
}

void CGeomCell::SetActualBeachErosion(double const dThisActualErosion) {
  m_dActualBeachErosionThisIter = dThisActualErosion;
  m_dTotActualBeachErosion += dThisActualErosion;
}

double CGeomCell::dGetActualBeachErosion(void) const {
  return m_dActualBeachErosionThisIter;
}

double CGeomCell::dGetTotActualBeachErosion(void) const {
  return m_dTotActualBeachErosion;
}
 
// //! Returns true if there has been actual beach erosion this timestep
// bool CGeomCell::bActualBeachErosionThisIter(void) const
// {
// return (m_dActualBeachErosionThisIter > 0 ? true : false);
// }

void CGeomCell::IncrBeachDeposition(double const dThisDeposition) {
  m_dBeachDepositionThisIter += dThisDeposition;
  m_dTotBeachDeposition += dThisDeposition;
}

double CGeomCell::dGetBeachDeposition(void) const {
  return m_dBeachDepositionThisIter;
}

double CGeomCell::dGetTotBeachDeposition(void) const {
  return m_dTotBeachDeposition;
}
 
// //! Returns true if there has been beach deposition this timestep
// bool CGeomCell::bBeachDepositionThisIter(void) const
// {
// return (m_dBeachDepositionThisIter > 0 ? true : false);
// }

bool CGeomCell::bBeachErosionOrDepositionThisIter(void) const {
  if ((m_dActualBeachErosionThisIter > 0) || (m_dBeachDepositionThisIter > 0))
    return true;
 
  return false;
}

double CGeomCell::dGetUnconsD50(void) const { return m_dUnconsD50; }

void CGeomCell::SetInterventionClass(int const nSubCatCode) {
  if (nSubCatCode != LF_NONE) {
    this->m_Landform.SetLFCategory(LF_CAT_INTERVENTION);
 
    if (nSubCatCode == IO_INTERVENTION_STRUCT)
      this->m_Landform.SetLFSubCategory(LF_SUBCAT_INTERVENTION_STRUCT);
 
    else if (nSubCatCode == IO_INTERVENTION_NON_STRUCT)
      this->m_Landform.SetLFSubCategory(LF_SUBCAT_INTERVENTION_NON_STRUCT);
  }
}

int CGeomCell::nGetInterventionClass(void) const {
  int nTmp = INT_NODATA;
 
  if (this->m_Landform.nGetLFCategory() == LF_CAT_INTERVENTION) {
    if (this->m_Landform.nGetLFSubCategory() == LF_SUBCAT_INTERVENTION_STRUCT)
      nTmp = IO_INTERVENTION_STRUCT;
 
    else if (this->m_Landform.nGetLFSubCategory() ==
             LF_SUBCAT_INTERVENTION_NON_STRUCT)
      nTmp = IO_INTERVENTION_NON_STRUCT;
  }
 
  return nTmp;
}

void CGeomCell::SetInterventionHeight(double const dHeight) {
  m_dInterventionHeight = dHeight;
}

double CGeomCell::dGetInterventionHeight(void) const {
  return m_dInterventionHeight;
}

double CGeomCell::dGetInterventionTopElev(void) const {
  return m_VdAllHorizonTopElev.back() + m_dInterventionHeight;
}