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"
#include "cell.h"
#include "cell_layer.h"
#include "cell_sediment.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 const dLayerThick = m_Layer.dGetFineConsolidatedThickness();
      double const 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 const dLayerThick = m_Layer.dGetSandConsolidatedThickness();
      double const 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 const dLayerThick = m_Layer.dGetCoarseConsolidatedThickness();
      double const 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 const 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 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() + 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 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::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;
}