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 <cstddef>
 
#include <vector>
using std::vector;
 
#include <assert.h>
 
#include "cell.h"
#include "cme.h"
#include "simulation.h"
#include "cell.h"
#include "cell_layer.h"
#include "cell_sediment.h"
#include "cell_talus.h"

CGeomCell::CGeomCell()
   : m_bInContiguousSea(false),
     m_bInContiguousFlood(false),
     m_bIsInActiveZone(false),
     m_bCliffToe(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_nCoastlineID(INT_NODATA),
     m_nProfileID(INT_NODATA),
     m_nProfileCoastID(INT_NODATA),
     m_nPolygonID(INT_NODATA),
     m_nPolygonCoastID(INT_NODATA),
     m_nShadowZoneNumber(0),
     m_nDownDriftZoneNumber(0),
     m_dLocalConsSlope(0),
     m_dBasementElevation(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_dSandTalusToUnconsThisIter(0),
     m_dTotSandTalusToUncons(0),
     m_dCoarseTalusToUnconsThisIter(0),
     m_dTotCoarseTalusToUncons(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_UNKNOWN);
}

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::bIsInContiguousSeaFlood(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(int const nCoast)
{
   m_nCoastlineID = nCoast;
}

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

int CGeomCell::nGetCoastline(void) const
{
   return m_nCoastlineID;
}

void CGeomCell::SetAsCliffToe(bool const bNewFlag)
{
   m_bCliffToe = bNewFlag;
}

bool CGeomCell::bIsCliffToe(void) const
{
   return m_bCliffToe;
}

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

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

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

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

int CGeomCell::nGetProfileCoastID(void) const
{
   return m_nProfileCoastID;
}

void CGeomCell::SetCoastAndProfileID(int const nProfileCoastID, int const nProfile)
{
   m_nProfileCoastID = nProfileCoastID;
   m_nProfileID = nProfile;
}

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

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;
}

double CGeomCell::dGetTalusDepth(void) const
{
   double dTotTalusDepth = 0;
 
   for (int nLayer = 0; nLayer < static_cast<int>(m_VLayerAboveBasement.size()); nLayer++)
   {
      CRWCellTalus const* pTalus = m_VLayerAboveBasement[nLayer].pGetTalus();
      if (pTalus != NULL)
      {
         // There is some talus on this layer
         double const dThisTalusDepth = pTalus->dGetSandDepth() + pTalus->dGetCoarseDepth();
         dTotTalusDepth += dThisTalusDepth;
      }
   }
 
   return dTotTalusDepth;
}

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::bElevLessThanSWL(void) const
{
   // // Note that m_pGrid->pGetSim()->dGetThisIterTotWaterLevel() is zero, since TODO 007 Finish surge and runup stuff
   // return ((m_VdAllHorizonTopElev.back() + this->dGetTalusDepth() + m_dInterventionHeight) < (m_pGrid->pGetSim()->dGetThisIterTotWaterLevel() + m_pGrid->pGetSim()->dGetThisIterSWL()));
 
   // Will need to change this (see above) once surge and runup stuff is working
   return ((m_VdAllHorizonTopElev.back() + this->dGetTalusDepth() + m_dInterventionHeight) < m_pGrid->pGetSim()->dGetThisIterSWL());
}
 
// //! Set this cell as checked TODO 007 Finish surge and runup stuff
// void CGeomCell::SetCheckCell(void)
// {
// m_bCheckCell = true;
// }
 
// //! Returns true if this cell is checked, false otherwise TODO 007 Finish surge and runup stuff
// 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::dGetConsSedSlope(void) const
{
   return m_dLocalConsSlope;
}

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

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

void CGeomCell::SetSlopeForCliffToe(double const dNewSlope)
{
   m_dSlopeForCliffToe = dNewSlope;
}

double CGeomCell::dGetSlopeForCliffToe(void) const
{
   return (m_dSlopeForCliffToe);
}

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::nGetNumLayers(void) const
{
   return static_cast<int>(m_VLayerAboveBasement.size());
}
 

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::nGetNumOfTopLayerAboveBasement(void) const
{
   if (m_VLayerAboveBasement.empty())
      return INT_NODATA;
 
   return static_cast<int>(m_VLayerAboveBasement.size()) - 1;
}

double CGeomCell::dGetConsSedTopElevForLayerAboveBasement(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].dGetAllUnconsDepth();
      dTopElev += m_VLayerAboveBasement[n].dGetAllConsDepth();
   }
 
   dTopElev += m_VLayerAboveBasement[nLayer].dGetAllConsDepth();
 
   return dTopElev;
}

CRWCellLayer* CGeomCell::pGetLayerAboveBasement(int const nLayer)
{
   // TODO 055 No check that nLayer < size()
   return &m_VLayerAboveBasement[nLayer];
}

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

double CGeomCell::dGetAllSedTopElevIncTalus(void)
{
   return this->dGetAllSedTopElevOmitTalus() + this->dGetTalusDepth();
}

double CGeomCell::dGetConsSedTopElevOmitTalus(void) const
{
   return this->dGetBasementElev() + this->dGetAllConsDepthAllLayers();
}

double CGeomCell::dGetConsSedTopElevIncTalus(void)
{
   return this->dGetConsSedTopElevOmitTalus() + this->dGetTalusDepth();
}

double CGeomCell::dGetTopElevIncSea(void)
{
   return m_VdAllHorizonTopElev.back() + this->dGetTalusDepth() + m_dInterventionHeight + m_dSeaDepth;
}

bool CGeomCell::bIsInundated(void)
{
   return ((m_VdAllHorizonTopElev.back() + m_dInterventionHeight + this->dGetTalusDepth()) < 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();
}

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

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

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

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

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

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

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

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

double CGeomCell::dGetAllSedDepthAllLayers(void) const
{
   return (this->dGetAllUnconsDepthAllLayers() + this->dGetAllConsDepthAllLayers());
}

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, ignoring any talus
 
   // 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 const dUnconsThick = m_VLayerAboveBasement[n].dGetAllUnconsDepth();
 
      if (dUnconsThick > 0)
      {
         // This is a layer with non-zero thickness of unconsolidated sediment
         CRWCellSediment const* pUnconsSedLayer = m_VLayerAboveBasement[n].pGetUnconsolidatedSediment();
         double const dFineProp = pUnconsSedLayer->dGetFineDepth() / dUnconsThick;
         double const dSandProp = pUnconsSedLayer->dGetSandDepth() / dUnconsThick;
         double const 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() + this->dGetTalusDepth() + m_dInterventionHeight), 0.0);
}

void CGeomCell::InitCell(void)
{
   m_bInContiguousSea = false;
   m_bInContiguousFlood = false;
   m_bFloodLine = false;
   m_bIsInActiveZone = false;
   // m_bEstimated = false;
   m_bShadowBoundary = false;
   m_bPossibleCoastStartCell = false;
   m_bPossibleFloodStartCell = false;
   m_bWaveFlood = false;
   // m_bCheckCell = false;
   m_bCheckFloodCell = false;
 
   m_nCoastlineID = INT_NODATA;
   m_nPolygonID = INT_NODATA;
   m_nPolygonCoastID = INT_NODATA;
   m_nProfileID = INT_NODATA;
 
   m_nShadowZoneNumber = 0;
   m_nDownDriftZoneNumber = 0;
 
   m_dLocalConsSlope = 0;
   m_dPotentialPlatformErosionThisIter = 0;
   m_dActualPlatformErosionThisIter = 0;
   m_dCliffCollapseFineThisIter = 0;
   m_dCliffCollapseSandThisIter = 0;
   m_dCliffCollapseCoarseThisIter = 0;
   m_dTalusSandDepositionThisIter = 0;
   m_dTalusCoarseDepositionThisIter = 0;
   m_dSandTalusToUnconsThisIter = 0;
   m_dCoarseTalusToUnconsThisIter = 0;
   m_dPotentialBeachErosionThisIter = 0;
   m_dActualBeachErosionThisIter = 0;
   m_dBeachDepositionThisIter = 0;
   m_dSeaDepth = 0;
   m_dWaveHeight = 0;
   m_dWaveAngle = 0;
 
   m_dBeachProtectionFactor = DBL_NODATA;
 
   // Initialize this-iteration sediment input event values
   int const 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::AddSandTalusToUncons(double const dSandTalus)
{
   m_dSandTalusToUnconsThisIter += dSandTalus;
   m_dTotSandTalusToUncons += dSandTalus;
}

double CGeomCell::dGetThisIterSandTalusToUncons(void)
{
   return m_dSandTalusToUnconsThisIter;
}

double CGeomCell::dGetTotSandTalusToUncons(void)
{
   return m_dTotSandTalusToUncons;
}

void CGeomCell::AddCoarseTalusToUncons(double const dCoarseTalus)
{
   m_dCoarseTalusToUnconsThisIter += dCoarseTalus;
   m_dTotCoarseTalusToUncons += dCoarseTalus;
}

double CGeomCell::dGetThisIterCoarseTalusToUncons(void)
{
   return m_dCoarseTalusToUnconsThisIter;
}

double CGeomCell::dGetTotCoarseTalusToUncons(void)
{
   return m_dTotCoarseTalusToUncons;
}

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;
}

int CGeomCell::nGetInterventionClass(void) const
{
   int const nTmp = this->m_Landform.nGetLFCategory();
   if ((nTmp == LF_INTERVENTION_STRUCT) || (nTmp == LF_INTERVENTION_NON_STRUCT))
      return nTmp;
 
   return INT_NODATA;
}

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;
}