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 <assert.h>
#include "cme.h"
#include "cell.h"

CGeomCell::CGeomCell()
    : m_bInContiguousSea(false),
      m_bInContiguousFlood(false),        // TODO 007 What is this?
      m_bIsInActiveZone(false),
      m_bCoastline(false),
      m_bFloodLine(false),
      m_bWaveFlood(false),
      // m_bCheckCell(false),             // TODO 007 What is this?
      m_bCheckFloodCell(false),           // TODO 007 What is this?
      m_bShadowBoundary(false),
      m_bPossibleCoastStartCell(false),
      m_bPossibleFloodStartCell(false),   // TODO 007 What is this?
      m_nBoundingBoxEdge(NO_DIRECTION),
      m_nPolygonID(INT_NODATA),
      m_nCoastlineNormal(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_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::bIsInContiguousFlood(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;
}
 
// //! TODO 007 What is this for? Returns a flag which shows whether this cell has been flagged as a possible start- or end-point for a coastline
// 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::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::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;
}
 
// Set this cell as flooded by swl + surge + setup + runup
// 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);
}

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 =           // TODO 007 What is this?
   m_bCoastline =
   m_bFloodLine =
   m_bIsInActiveZone =
   // m_bEstimated =
   m_bShadowBoundary =
   m_bPossibleCoastStartCell =      // TODO 007 What is this?
   m_bPossibleFloodStartCell =
   m_bWaveFlood =
   // m_bCheckCell =                // TODO 007 What is this?
   m_bCheckFloodCell = false;
 
   m_nPolygonID =
   m_nCoastlineNormal = INT_NODATA;
 
   m_nShadowZoneNumber =
   m_nDownDriftZoneNumber = 0;
 
   m_dLocalConsSlope =
   m_dPotentialPlatformErosionThisIter =
   m_dActualPlatformErosionThisIter =
   m_dCliffCollapseFineThisIter =
   m_dCliffCollapseSandThisIter =
   m_dCliffCollapseCoarseThisIter =
   m_dTalusSandDepositionThisIter =
   m_dTotTalusSandDeposition = 
   m_dTalusCoarseDepositionThisIter = 
   m_dTotTalusCoarseDeposition = 
   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;
}