CoastalME (Coastal Modelling Environment)
Simulates the long-term behaviour of complex coastlines
Loading...
Searching...
No Matches
do_beach_within_polygon.cpp
Go to the documentation of this file.
1
10
11/* ==============================================================================================================================
12 This file is part of CoastalME, the Coastal Modelling Environment.
13
14 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.
15
16 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.
17
18 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.
19==============================================================================================================================*/
20#include <assert.h>
21
22#include <cmath>
23
24#include <iostream>
25using std::endl;
26
27#include <algorithm>
28using std::shuffle;
29
30#include "cme.h"
31#include "simulation.h"
32#include "coast.h"
33#include "2d_point.h"
34#include "2di_point.h"
35
36//===============================================================================================================================
38//===============================================================================================================================
39int CSimulation::nDoUnconsErosionOnPolygon(int const nCoast, CGeomCoastPolygon* pPolygon, int const nTexture, double const dErosionTargetOnPolygon, double& dEroded)
40{
41 string strTexture;
42
43 if (nTexture == TEXTURE_FINE)
44 strTexture = "fine";
45
46 else if (nTexture == TEXTURE_SAND)
47 strTexture = "sand";
48
49 else if (nTexture == TEXTURE_COARSE)
50 strTexture = "coarse";
51
52 // Intialise
53 double dStillToErodeOnPolygon = dErosionTargetOnPolygon;
54 dEroded = 0;
55
56 // Get the up-coast and down-coast boundary details
57 int const nUpCoastProfile = pPolygon->nGetUpCoastProfile();
58 CGeomProfile* pUpCoastProfile = m_VCoast[nCoast].pGetProfile(nUpCoastProfile);
59
60 int const nDownCoastProfile = pPolygon->nGetDownCoastProfile();
61 CGeomProfile const* pDownCoastProfile = m_VCoast[nCoast].pGetProfile(nDownCoastProfile);
62
63 // We will use only part of the up-coast boundary profile, seaward as far as the depth of closure. First find the seaward end point of this up-coast part-profile. Note that this does not change as the landwards offset changes
64 int const nIndex = pUpCoastProfile->nGetCellGivenDepth(m_pRasterGrid, m_dDepthOfClosure);
65
66 if (nIndex == INT_NODATA)
67 {
68 LogStream << m_ulIter << ": " << ERR << "while eroding unconsolidated " + strTexture + " sediment on polygon " << pPolygon->nGetPolygonCoastID() << ", could not find the seaward end point of the up-coast profile (" << nUpCoastProfile << ") for depth of closure = " << m_dDepthOfClosure << endl;
69
70 return RTN_ERR_NO_SEAWARD_END_OF_PROFILE_BEACH_EROSION;
71 }
72
73 // The part-profile length is one greater than nIndex, since pPtiGetCellGivenDepth() returns the index of the cell at the depth of closure
74 int const nUpCoastPartProfileLen = nIndex + 1;
75
76 // assert(bIsWithinValidGrid(&PtiUpCoastPartProfileSeawardEnd));
77
78 // Store the cell coordinates of the boundary part-profile in reverse (sea to coast) order so we can append to the coastward end as we move inland (i.e. as nInlandOffset increases)
79 vector<CGeom2DIPoint> PtiVUpCoastPartProfileCell;
80
81 // for (int n = 0; n < nUpCoastPartProfileLen; n++)
82 // PtiVUpCoastPartProfileCell.push_back(*pUpCoastProfile->pPtiGetCellInProfile(nUpCoastPartProfileLen - n - 1));
83 for (int n = nUpCoastPartProfileLen - 1; n >= 0; n--)
84 {
85 CGeom2DIPoint const Pti = *pUpCoastProfile->pPtiGetCellInProfile(n);
86 PtiVUpCoastPartProfileCell.push_back(Pti);
87 }
88
89 int const nUpCoastProfileCoastPoint = pUpCoastProfile->nGetCoastPoint();
90 int const nDownCoastProfileCoastPoint = pDownCoastProfile->nGetCoastPoint();
91 int const nXUpCoastProfileExistingCoastPoint = m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nUpCoastProfileCoastPoint)->nGetX();
92 int const nYUpCoastProfileExistingCoastPoint = m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nUpCoastProfileCoastPoint)->nGetY();
93 int nCoastSegLen;
94
95 // Store the coast point numbers for this polygon so that we can shuffle them
96 vector<int> nVCoastPoint;
97
98 if (nDownCoastProfileCoastPoint == m_VCoast[nCoast].nGetCoastlineSize() - 1)
99 {
100 // This is the final down-coast polygon, so also include the down-coast polygon boundary
101 nCoastSegLen = nDownCoastProfileCoastPoint - nUpCoastProfileCoastPoint + 1;
102
103 for (int nCoastPoint = nUpCoastProfileCoastPoint; nCoastPoint <= nDownCoastProfileCoastPoint; nCoastPoint++)
104 nVCoastPoint.push_back(nCoastPoint);
105 }
106 else
107 {
108 // This is not the final down-coast polygon, so do not include the polygon's down-coast boundary
109 nCoastSegLen = nDownCoastProfileCoastPoint - nUpCoastProfileCoastPoint;
110
111 for (int nCoastPoint = nUpCoastProfileCoastPoint; nCoastPoint < nDownCoastProfileCoastPoint; nCoastPoint++)
112 nVCoastPoint.push_back(nCoastPoint);
113 }
114
115 // Estimate the volume of sediment which is to be eroded from each parallel profile
116 double const dAllTargetPerProfile = dErosionTargetOnPolygon / nCoastSegLen;
117
118 // Shuffle the coast points, this is necessary so that leaving the loop does not create sequence-related artefacts
119 shuffle(nVCoastPoint.begin(), nVCoastPoint.begin() + nCoastSegLen, m_Rand[1]);
120
121 // Traverse the polygon's existing coastline in a DOWN-COAST (i.e. increasing coastpoint indices) sequence, at each coast point fitting a Dean profile which is parallel to the up-coast polygon boundary
122 for (int n = 0; n < nCoastSegLen; n++)
123 {
124 // Get the coast point
125 int const nCoastPoint = nVCoastPoint[n];
126
127 CGeom2DIPoint const PtiCoastPoint = *m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nCoastPoint);
128 int const nCoastX = PtiCoastPoint.nGetX();
129 int const nCoastY = PtiCoastPoint.nGetY();
130
131 // LogStream << m_ulIter << ": nCoastX = " << nCoastX << ", nCoastY = " << nCoastY << ", this is {" << dGridCentroidXToExtCRSX(nCoastX) << ", " << dGridCentroidYToExtCRSY(nCoastY) << "}" << endl;
132
133 // Is the coast cell an intervention structure?
134 if (bIsInterventionCell(nCoastX, nCoastY))
135 {
136 // No erosion possible on this parallel profile, so move on
137 // LogStream << m_ulIter << ": intervention structure at coast point [" << nCoastX << "][" << nCoastY << "] = {" << dGridCentroidXToExtCRSX(nCoastX) << ", " << dGridCentroidYToExtCRSY(nCoastY) << "}, cannot erode this parallel profile" << endl;
138
139 continue;
140 }
141
142 // LogStream << m_ulIter << ": PtiVUpCoastPartProfileCell.back().nGetX() = " << PtiVUpCoastPartProfileCell.back().nGetX() << ", PtiVUpCoastPartProfileCell.back().nGetY() = " << PtiVUpCoastPartProfileCell.back().nGetY() << ", this is {" << dGridCentroidXToExtCRSX(PtiVUpCoastPartProfileCell.back().nGetX()) << ", " << dGridCentroidYToExtCRSY(PtiVUpCoastPartProfileCell.back().nGetY()) << "}" << endl;
143
144 // Not an intervention structure, so calculate the x-y offset between this coast point, and the coast point of the up-coast normal
145 int const nXOffset = nCoastX - PtiVUpCoastPartProfileCell.back().nGetX();
146 int const nYOffset = nCoastY - PtiVUpCoastPartProfileCell.back().nGetY();
147
148 // Get the x-y coords of a profile starting from this coast point and parallel to the up-coast polygon boundary profile (these are in reverse sequence, like the boundary part-profile)
149 vector<CGeom2DIPoint> VPtiParProfile;
150
151 for (int m = 0; m < nUpCoastPartProfileLen; m++)
152 {
153 // TODO 017 Check that each point is within valid grid, do same for other similar places in rest of model
154 CGeom2DIPoint const PtiTmp(PtiVUpCoastPartProfileCell[m].nGetX() + nXOffset, PtiVUpCoastPartProfileCell[m].nGetY() + nYOffset);
155 VPtiParProfile.push_back(PtiTmp);
156 }
157
158 // Get the elevations of the start and end points of the parallel profiles (as we extend the profile inland, the elevation of the new coast point of the Dean profile is set to the elevation of the original coast point)
159 int nParProfEndX = VPtiParProfile[0].nGetX();
160 int nParProfEndY = VPtiParProfile[0].nGetY();
161
162 // Safety check
163 if (! bIsWithinValidGrid(nParProfEndX, nParProfEndY))
164 {
165 // if (m_nLogFileDetail >= LOG_FILE_ALL)
166 // LogStream << WARN << "while eroding unconsolidated " + strTexture + " sediment for coast " << nCoast << " polygon " << pPolygon->nGetPolygonCoastID() << ", hit edge of grid at [" << nParProfEndX << "][" << nParProfEndY << "] for parallel profile from coast point " << nCoastPoint << " at [" << nCoastX << "][" << nCoastY << "]. Constraining this parallel profile at its seaward end" << endl;
167
168 KeepWithinValidGrid(nCoastX, nCoastY, nParProfEndX, nParProfEndY);
169 VPtiParProfile[0].SetX(nParProfEndX);
170 VPtiParProfile[0].SetY(nParProfEndY);
171 }
172
173 bool bHitEdge = false;
174 bool bEndProfile = false;
175 bool bZeroGradient = false;
176 bool bEnoughEroded = false;
177
178 int nParProfLen;
179 int nInlandOffset = -1;
180
181 // Note that we must not consider the extra elevation due to to any talus here, since talus is removed quickly compared with the time required to create a Dean profile
182 double const dParProfCoastElev = m_pRasterGrid->m_Cell[nCoastX][nCoastY].dGetAllSedTopElevOmitTalus();
183 double const dParProfEndElev = m_pRasterGrid->m_Cell[nParProfEndX][nParProfEndY].dGetAllSedTopElevOmitTalus();
184
185 vector<double> VdParProfileDeanElev;
186
187 // These are for saving values for each offset
188 vector<int> VnParProfLenEachOffset;
189 vector<double> VdAmountEachOffset;
190 vector<vector<CGeom2DIPoint>> VVPtiParProfileEachOffset;
191 vector<vector<double>> VVdParProfileDeanElevEachOffset;
192
193 // OK, loop either until we can erode sufficient unconsolidated sediment, or until the landwards-moving parallel profile hits the grid edge
194 while (true)
195 {
196 // Move inland by one cell
197 nInlandOffset++;
198
199 if (nInlandOffset > 0)
200 {
201 if (nInlandOffset > (pUpCoastProfile->nGetNumCellsInProfile() - 1))
202 {
203 if (m_nLogFileDetail >= LOG_FILE_HIGH_DETAIL)
204 LogStream << m_ulIter << ": \tcoast " << nCoast << " reached end of up-coast profile " << nUpCoastProfile << " during down-coast erosion of unconsolidated " + strTexture + " sediment for coast " << nCoast << " polygon " << pPolygon->nGetPolygonCoastID() << " (nCoastPoint = " << nCoastPoint << " nInlandOffset = " << nInlandOffset << ")" << endl;
205
206 bEndProfile = true;
207 break;
208 }
209
210 // Extend the parallel profile by one cell in the coastward direction. First get the coords of the cell that is nInlandOffset cells seaward from the existing up-coast part-profile start point
211 CGeom2DIPoint const PtiUpCoastTmp = *pUpCoastProfile->pPtiGetCellInProfile(nInlandOffset);
212
213 // Then get the offset between this PtiUpCoastTmp cell and the existing up-coast part-profile start point, and use the reverse of this offset to get the coordinates of the cell that extends the existing up-coast part-profile landwards
214 int const nXUpCoastStartOffset = PtiUpCoastTmp.nGetX() - nXUpCoastProfileExistingCoastPoint;
215 int const nYUpCoastStartOffset = PtiUpCoastTmp.nGetY() - nYUpCoastProfileExistingCoastPoint;
216 int const nXUpCoastThisStart = nCoastX - nXUpCoastStartOffset;
217 int const nYUpCoastThisStart = nCoastY - nYUpCoastStartOffset;
218
219 // Is the new landwards point within the raster grid?
220 if (! bIsWithinValidGrid(nXUpCoastThisStart, nYUpCoastThisStart))
221 {
222 // It isn't
223 // LogStream << WARN << "reached edge of grid at [" << nXUpCoastThisStart << "][" << nYUpCoastThisStart << "] = {" << dGridCentroidXToExtCRSX(nXUpCoastThisStart) << ", " << dGridCentroidYToExtCRSY(nYUpCoastThisStart) << "} during DOWN-COAST beach erosion for coast " << nCoast << " polygon " << nPoly << ", nCoastPoint = " << nCoastPoint << ", this is [" << nCoastX << "][" << nCoastY << "] = {" << dGridCentroidXToExtCRSX(nCoastX) << ", " << dGridCentroidYToExtCRSY(nCoastY) << "}, nInlandOffset = " << nInlandOffset << ")" << endl << endl;
224
225 // TODO 018 Need to improve this: at present we just abandon erosion on this coast point and move to another coast point
226 bHitEdge = true;
227 break;
228 }
229
230 // CGeom2DIPoint PtiThisUpCoastStart(nXUpCoastThisStart, nYUpCoastThisStart);
231
232 // Calculate the coordinates of a possible new landwards cell for the parallel profile
233 int nXParNew = nXUpCoastThisStart + nXOffset;
234 int nYParNew = nYUpCoastThisStart + nYOffset;
235
236 // Safety check
237 if (! bIsWithinValidGrid(nXParNew, nYParNew))
238 {
239 // LogStream << WARN << "while eroding beach on coast " << nCoast << " polygon " << nPoly << " (nInlandOffset = " << nInlandOffset << "), outside valid grid at [" << nXParNew << "][" << nYParNew << "] for parallel profile from coast point " << nCoastPoint << " at [" << nCoastX << "][" << nCoastY << "]. Constraining this parallel profile at its landward end, is now [";
240
241 KeepWithinValidGrid(nCoastX, nCoastY, nXParNew, nYParNew);
242
243 // LogStream << "[" << nXParNew << "][" << nYParNew << "] = {" << dGridCentroidXToExtCRSX(nXParNew) << ", " << dGridCentroidYToExtCRSY(nYParNew) << "}" << endl;
244
245 // Is this cell already in the parallel profile?
246 if ((VPtiParProfile.back().nGetX() != nXParNew) || (VPtiParProfile.back().nGetY() != nYParNew))
247 {
248 // It isn't, so append it to the parallel profile
249 CGeom2DIPoint const PtiTmp(nXParNew, nYParNew);
250 VPtiParProfile.push_back(PtiTmp);
251 }
252 }
253
254 else
255 {
256 // No problem, so just append this to the parallel profile
257 CGeom2DIPoint const PtiTmp(nXParNew, nYParNew);
258 VPtiParProfile.push_back(PtiTmp);
259 }
260 }
261
262 nParProfLen = static_cast<int>(VPtiParProfile.size());
263
264 if (nParProfLen < MIN_PARALLEL_PROFILE_SIZE)
265 {
266 // Can't have a meaningful parallel profile with very few points
267 // LogStream << m_ulIter << ": only " << nParProfLen << " points in parallel profile, min is " << MIN_PARALLEL_PROFILE_SIZE << ", abandoning" << endl;
268
269 continue;
270 }
271
272 // for (int m = 0; m < static_cast<int>(VPtiParProfile.size()); m++)
273 // LogStream << "[" << VPtiParProfile[m].nGetX() << "][" << VPtiParProfile[m].nGetY() << "] ";
274 // LogStream << endl;
275
276 // Get the distance between the start and end of the parallel profile, in external CRS units. Note that the parallel profile coordinates are in reverse sequence
277 CGeom2DPoint const PtStart = PtGridCentroidToExt(&VPtiParProfile.back());
278 CGeom2DPoint const PtEnd = PtGridCentroidToExt(&VPtiParProfile[0]);
279
280 // Calculate the length of the parallel profile
281 double dParProfileLen = dGetDistanceBetween(&PtStart, &PtEnd);
282
283 // Calculate the elevation difference between the start and end of the parallel profile
284 double const dElevDiff = dParProfCoastElev - dParProfEndElev;
285
286 if (bFPIsEqual(dElevDiff, 0.0, TOLERANCE))
287 {
288 // Can't have a meaningful Dean profile with a near-zero elevation difference
289 // TODO 019 Need to improve this: at present we just abandon erosion on this coast point and move to another coast point
290 // LogStream << m_ulIter << ": zero gradient on parallel profile, abandoning" << endl;
291
292 bZeroGradient = true;
293 break;
294 }
295
296 // Safety check
297 if (dParProfileLen <= 0)
298 dParProfileLen = 1;
299
300 // Solve for dA so that the existing elevations at the end of the parallel profile, and at the end of a Dean equilibrium profile on that part-normal, are the same
301 double const dParProfA = dElevDiff / pow(dParProfileLen, DEAN_POWER);
302 VdParProfileDeanElev.resize(nParProfLen, 0);
303 double const dInc = dParProfileLen / (nParProfLen - 1);
304
305 // For the parallel profile, calculate the Dean equilibrium profile of the unconsolidated sediment h(y) = A * y^(2/3) where h(y) is the distance below the highest point in the profile at a distance y from the landward start of the profile
306 CalcDeanProfile(&VdParProfileDeanElev, dInc, dParProfCoastElev, dParProfA, false, 0, 0);
307
308 vector<double> dVParProfileNow(nParProfLen, 0);
309 vector<bool> bVProfileValid(nParProfLen, true);
310
311 for (int m = 0; m < nParProfLen; m++)
312 {
313 int const nX = VPtiParProfile[nParProfLen - m - 1].nGetX();
314 int const nY = VPtiParProfile[nParProfLen - m - 1].nGetY();
315
316 // Safety check
317 if (! bIsWithinValidGrid(nX, nY))
318 {
319 // if (m_nLogFileDetail >= LOG_FILE_ALL)
320 // LogStream << WARN << "while constructing parallel profile for erosion of unconsolidated " + strTexture + " sediment on coast " << nCoast << " polygon " << nPoly << ", hit edge of grid at [" << nX << "][" << nY << "] for parallel profile from coast point " << nCoastPoint << " at [" << nCoastX << "][" << nCoastY << "]. Constraining this parallel profile at its seaward end" << endl;
321
322 bVProfileValid[m] = false;
323 continue;
324 }
325
326 // Don't erode intervention cells
327 if (bIsInterventionCell(nX, nY))
328 bVProfileValid[m] = false;
329
330 dVParProfileNow[m] = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus();
331 }
332
333 // Get the total difference in elevation (present profile - Dean profile)
334 double const dParProfTotDiff = dSubtractProfiles(&dVParProfileNow, &VdParProfileDeanElev, &bVProfileValid);
335
336 // // DEBUG CODE -----------------------------------------------------
337 // LogStream << m_ulIter<< ": eroding polygon " << nPoly << " at coast point " << nCoastPoint << ", parallel profile with nInlandOffset = " << nInlandOffset << ", from [" << VPtiParProfile.back().nGetX() << "][" << VPtiParProfile.back().nGetY() << "] = {" << dGridCentroidXToExtCRSX(VPtiParProfile.back().nGetX()) << ", " << dGridCentroidYToExtCRSY(VPtiParProfile.back().nGetY()) << "} to [" << VPtiParProfile[0].nGetX() << "][" << VPtiParProfile[0].nGetY() << "] = {" << dGridCentroidXToExtCRSX(VPtiParProfile[0].nGetX()) << ", " << dGridCentroidYToExtCRSY(VPtiParProfile[0].nGetY()) << "}, nParProfLen = " << nParProfLen << " dParProfileLen = " << dParProfileLen << " dParProfCoastElev = " << dParProfCoastElev << " dParProfEndElev = " << dParProfEndElev << " dParProfA = " << dParProfA << endl;
338
339 // LogStream << "Profile now:" << endl;
340 // for (int n = 0; n < nParProfLen; n++)
341 // {
342 // if (bVProfileValid[n])
343 // LogStream << dVParProfileNow[n] << " ";
344 // else
345 // LogStream << "XXX ";
346 // }
347 // LogStream << endl << endl;;
348 //
349 // LogStream << "Parallel Dean profile for erosion:" << endl;
350 // for (int n = 0; n < nParProfLen; n++)
351 // {
352 // if (bVProfileValid[n])
353 // LogStream << VdParProfileDeanElev[n] << " ";
354 // else
355 // LogStream << "XXX ";
356 // }
357 // LogStream << endl << endl;
358 //
359 // LogStream << "Difference (present profile minus Dean profile):" << endl;
360 // for (int n = 0; n < nParProfLen; n++)
361 // {
362 // if (bVProfileValid[n])
363 // LogStream << dVParProfileNow[n] - VdParProfileDeanElev[n] << " ";
364 // else
365 // LogStream << "XXX ";
366 // }
367 // LogStream << endl << endl;
368 //
369 // LogStream << "dParProfTotDiff = " << dParProfTotDiff << " dAllTargetPerProfile = " << dAllTargetPerProfile << endl;
370 // // DEBUG CODE -----------------------------------------------------
371
372 // So will we be able to erode as much as is needed?
373 if (dParProfTotDiff > dAllTargetPerProfile)
374 {
375 // LogStream << m_ulIter << ": eroding polygon " << nPoly << " at coast point " << nCoastPoint << ", on parallel profile with nInlandOffset = " << nInlandOffset << ", can meet erosion target: dParProfTotDiff = " << dParProfTotDiff << " dAllTargetPerProfile = " << dAllTargetPerProfile << endl;
376
377 bEnoughEroded = true;
378 break;
379 }
380
381 // We have not been able to reach the erosion target for this parallel profile. Now, if we have moved at least MIN_INLAND_OFFSET_UNCONS_EROSION cells inland, and dParProfTotDiff is zero, break out of the loop
382 if ((nInlandOffset >= MIN_INLAND_OFFSET_UNCONS_EROSION) && (bFPIsEqual(dParProfTotDiff, 0.0, TOLERANCE)))
383 {
384 if (m_nLogFileDetail >= LOG_FILE_ALL)
385 LogStream << m_ulIter << ": leaving loop because nInlandOffset (" << nInlandOffset << ") >= MIN_INLAND_OFFSET_UNCONS_EROSION) and dParProfTotDiff = " << dParProfTotDiff << endl;
386
387 break;
388 }
389
390 // Save the amount which can be eroded for this offset
391 VdAmountEachOffset.push_back(dParProfTotDiff);
392 VnParProfLenEachOffset.push_back(nParProfLen);
393 VVPtiParProfileEachOffset.push_back(VPtiParProfile);
394 VVdParProfileDeanElevEachOffset.push_back(VdParProfileDeanElev);
395 }
396
397 // If we hit the edge of the grid, or have a zero gradient on the profile, or this is an end profile, then abandon this profile and do the next parallel profile
398 // TODO 019 TODO 018 Improve this, see above
399 if (bHitEdge || bEndProfile || bZeroGradient)
400 continue;
401
402 // OK we will do some erosion on this parallel profile, set the target for this profile to be the full target amount
403 double dStillToErodeOnProfile = dAllTargetPerProfile;
404
405 if (! bEnoughEroded)
406 {
407 // We have not been able to reach the target for erosion on this parallel profile. So find the offset that gives us the largest erosion amount
408 int nOffsetForLargestPossible = -1;
409 double dLargestPossibleErosion = 0;
410
411 for (unsigned int nn = 0; nn < VdAmountEachOffset.size(); nn++)
412 {
413 if (VdAmountEachOffset[nn] > dLargestPossibleErosion)
414 {
415 dLargestPossibleErosion = VdAmountEachOffset[nn];
416 nOffsetForLargestPossible = nn;
417 }
418 }
419
420 // If every offset gave zero erosion then abandon this profile and do the next parallel profile
421 if (nOffsetForLargestPossible < 0)
422 continue;
423
424 // OK, we have an offset which gives us the largest possible erosion (but less than the full target amount), continue with this
425 nInlandOffset = nOffsetForLargestPossible;
426 dStillToErodeOnProfile = dLargestPossibleErosion;
427 nParProfLen = VnParProfLenEachOffset[nInlandOffset];
428 VPtiParProfile = VVPtiParProfileEachOffset[nInlandOffset];
429 VdParProfileDeanElev = VVdParProfileDeanElevEachOffset[nInlandOffset];
430
431 // LogStream << m_ulIter << ": eroding polygon " << nPoly << " at coast point " << nCoastPoint << ", for parallel profile with nInlandOffset = " << nInlandOffset << ", could not meet erosion target dAllTargetPerProfile = " << dAllTargetPerProfile << ", instead using best possible: nInlandOffset = " << nInlandOffset << " which gives dStillToErodeOnProfile = " << dStillToErodeOnProfile << endl;
432 }
433
434 // This value of nInlandOffset gives us some (tho' maybe not enough) erosion. So do the erosion of this sediment size class, by working along the parallel profile from the landward end (which is inland from the existing coast, if nInlandOffset > 0). Note that dStillToErodeOnProfile and dStillToErodeOnPolygon are changed within nDoParallelProfileUnconsErosion()
435 int const nRet = nDoParallelProfileUnconsErosion(pPolygon, nCoastPoint, nCoastX, nCoastY, nTexture, nInlandOffset, nParProfLen, &VPtiParProfile, &VdParProfileDeanElev, dStillToErodeOnProfile, dStillToErodeOnPolygon, dEroded);
436
437 if (nRet != RTN_OK)
438 return nRet;
439
440 // LogStream << m_ulIter << ": eroding polygon " << nPoly << ": finished at coast point " << nCoastPoint << " dStillToErodeOnProfile = " << dStillToErodeOnProfile << " dStillToErodeOnPolygon = " << dStillToErodeOnPolygon << " dAllTargetPerProfile = " << dAllTargetPerProfile << endl;
441
442 // if (dStillToErodeOnProfile > 0)
443 // {
444 // // LogStream << " X Still to erode on profile = " << dStillToErodeOnProfile << " dStillToErodeOnPolygon WAS = " << dStillToErodeOnPolygon;
445 // dStillToErodeOnPolygon += dStillToErodeOnProfile;
446 // dAllTargetPerProfile += (dStillToErodeOnProfile / (nCoastSegLen - n));
447 // // LogStream << " dStillToErodeOnPolygon NOW = " << dStillToErodeOnPolygon << endl;
448 // }
449
450 if (dStillToErodeOnPolygon <= 0)
451 {
452 // LogStream << m_ulIter << ": YYYYYYYYYYYYYYY in polygon " << nPoly << ", leaving loop because dStillToErodeOnPolygon = " << dStillToErodeOnPolygon << endl;
453 break;
454 }
455 }
456
457 // How much have we been able to erode?
458 // LogStream << m_ulIter << ": nPoly = " << nPoly << " dEroded = " << dEroded << endl;
459
460 // // Is the depth eroded within TOLERANCE of the target depth-equivalent?
461 // if (bFPIsEqual(dEroded, dErosionTargetOnPolygon, TOLERANCE))
462 // {
463 // if (m_nLogFileDetail >= LOG_FILE_MIDDLE_DETAIL)
464 // LogStream << m_ulIter << ": polygon " << nPoly << " actual erosion of " + strTexture + " sediment is approximately equal to estimate, actual = " << dEroded * m_dCellArea << " estimate = " << dErosionTargetOnPolygon * m_dCellArea << endl;
465 // }
466 // else
467 // {
468 // if (dEroded < dErosionTargetOnPolygon)
469 // {
470 // if (m_nLogFileDetail >= LOG_FILE_MIDDLE_DETAIL)
471 // LogStream << m_ulIter << ": polygon " << nPoly << " actual erosion of " + strTexture + " sediment is less than estimate, actual = " << dEroded * m_dCellArea << " estimate = " << dErosionTargetOnPolygon * m_dCellArea << " dStillToErodeOnPolygon = " << dStillToErodeOnPolygon * m_dCellArea << ". This reduces the volume of " + strTexture + " sediment exported from this polygon to " << (dErosionTargetOnPolygon - dStillToErodeOnPolygon) * m_dCellArea << endl;
472 // }
473 // else
474 // {
475 // if (m_nLogFileDetail >= LOG_FILE_MIDDLE_DETAIL)
476 // LogStream << ERR << "polygon " << nPoly << " actual erosion of " + strTexture + " sediment is greater than estimate, actual = " << dEroded * m_dCellArea << " estimate = " << dErosionTargetOnPolygon * m_dCellArea << " dStillToErodeOnPolygon = " << dStillToErodeOnPolygon * m_dCellArea << endl;
477 // }
478 // }
479
480 return RTN_OK;
481}
482
483//===============================================================================================================================
485//===============================================================================================================================
486int CSimulation::nDoParallelProfileUnconsErosion(CGeomCoastPolygon* pPolygon, int const nCoastPoint, int const nCoastX, int const nCoastY, int const nTexture, int const nInlandOffset, int const nParProfLen, vector<CGeom2DIPoint> const* pVPtiParProfile, vector<double> const* pVdParProfileDeanElev, double& dStillToErodeOnProfile, double& dStillToErodeOnPolygon, double& dTotEroded)
487{
488 for (int nDistSeawardFromNewCoast = 0; nDistSeawardFromNewCoast < nParProfLen; nDistSeawardFromNewCoast++)
489 {
490 // Leave the loop if we have eroded enough for this polygon
491 if (dStillToErodeOnPolygon <= 0)
492 {
493 if (m_nLogFileDetail >= LOG_FILE_ALL)
494 LogStream << m_ulIter << ": AAA in polygon " << pPolygon->nGetPolygonCoastID() << " at coast point " << nCoastPoint << " nInlandOffset = " << nInlandOffset << ", leaving loop because enough erosion for polygon, dStillToErodeOnPolygon = " << dStillToErodeOnPolygon << " dStillToErodeOnProfile = " << dStillToErodeOnProfile << endl;
495
496 break;
497 }
498
499 // Leave the loop if we have eroded enough for this parallel profile
500 if (dStillToErodeOnProfile <= 0)
501 {
502 if (m_nLogFileDetail >= LOG_FILE_ALL)
503 LogStream << m_ulIter << ": BBB in polygon " << pPolygon->nGetPolygonCoastID() << " at coast point " << nCoastPoint << " nInlandOffset = " << nInlandOffset << ", leaving loop because enough erosion for profile, dStillToErodeOnPolygon = " << dStillToErodeOnPolygon << " dStillToErodeOnProfile = " << dStillToErodeOnProfile << endl;
504
505 break;
506 }
507
508 CGeom2DIPoint PtiTmp = pVPtiParProfile->at(nParProfLen - nDistSeawardFromNewCoast - 1);
509 int
510 nX = PtiTmp.nGetX(),
511 nY = PtiTmp.nGetY();
512
513 // Safety check
514 if (! bIsWithinValidGrid(nX, nY))
515 {
516 // LogStream << WARN << "04 @@@@ while eroding polygon " << nPoly << ", hit edge of grid at [" << nX << "][" << nY << "] for parallel profile from coast point " << nCoastPoint << ". Constraining this parallel profile at its seaward end" << endl;
517
518 KeepWithinValidGrid(nCoastX, nCoastY, nX, nY);
519 PtiTmp.SetX(nX);
520 PtiTmp.SetY(nY);
521 }
522
523 // Don't do anything to intervention cells
524 if (bIsInterventionCell(nX, nY))
525 continue;
526
527 // Don't do cells twice
528 if (! m_pRasterGrid->m_Cell[nX][nY].bBeachErosionOrDepositionThisIter())
529 {
530 // Get this cell's current elevation
531 double const dThisElevNow = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus();
532
533 // LogStream << "\tnPoly = " << nPoly << ", [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "} nCoastPoint = " << nCoastPoint << " nDistSeawardFromNewCoast = " << nDistSeawardFromNewCoast << " dThisElevNow = " << dThisElevNow << " Dean Elev = " << VdParProfileDeanElev[nDistSeawardFromNewCoast] << endl;
534
535 // Subtract the two elevations
536 double const dElevDiff = dThisElevNow - pVdParProfileDeanElev->at(nDistSeawardFromNewCoast);
537
538 if ((dElevDiff > 0) && (m_pRasterGrid->m_Cell[nX][nY].bIsInContiguousSea()))
539 {
540 // The current elevation is higher than the Dean elevation, so we have possible beach erosion (i.e. if not constrained by availability of unconsolidated sediment) here
541 // LogStream << "\tnPoly = " << nPoly << " doing DOWN-COAST, possible beach erosion = " << dElevDiff << " at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "} nCoastPoint = " << nCoastPoint << " nDistSeawardFromNewCoast = " << nDistSeawardFromNewCoast << " dThisElevNow = " << dThisElevNow << " Dean Elev = " << pVdParProfileDeanElev->at(nDistSeawardFromNewCoast) << endl;
542
543 // Now get the number of the highest layer with non-zero thickness
544 int const nThisLayer = m_pRasterGrid->m_Cell[nX][nY].nGetTopNonZeroLayerAboveBasement();
545
546 // Safety check
547 if (nThisLayer == INT_NODATA)
548 return RTN_ERR_NO_TOP_LAYER;
549
550 if (nThisLayer != NO_NONZERO_THICKNESS_LAYERS)
551 {
552 // We still have at least one layer left with non-zero thickness (i.e. we are not down to basement), and the cell's current elevation is higher than the Dean equilibrium profile elevation. So do some beach erosion
553 double const dToErode = tMin(dElevDiff, dStillToErodeOnProfile, dStillToErodeOnPolygon);
554 double dRemoved = 0;
555
556 // assert(dToErode > 0);
557
558 // Erode this sediment size class
559 ErodeCellBeachSedimentSupplyLimited(nX, nY, nThisLayer, nTexture, dToErode, dRemoved);
560
561 if (dRemoved > 0)
562 {
563 // Update totals for the polygon and the profile
564 dTotEroded += dRemoved;
565 dStillToErodeOnProfile -= dRemoved;
566 dStillToErodeOnPolygon -= dRemoved;
567
568 // Update this-timestep totals
569 m_ulThisIterNumActualBeachErosionCells++;
570
571 // LogStream << m_ulIter << ": in polygon " << nPoly << ", actual beach erosion = " << dTmpTot << " at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "} nCoastPoint = " << nCoastPoint << " nDistSeawardFromNewCoast = " << nDistSeawardFromNewCoast << endl;
572
573 // // DEBUG CODE ==================================================================================================
574 // if (nTexture == TEXTURE_SAND)
575 // {
576 // LogStream << m_ulIter << ": $$$$$$$$$ BEACH 2 Dean profile lower than existing profile, SAND depth eroded on [" << nX << "][" << nY << "] in Poly " << nPoly << " is " << dRemoved * m_dCellArea << endl;
577 // }
578 // else
579 // {
580 // LogStream << m_ulIter << ": $$$$$$$$$ BEACH 2 Dean profile lower than existing profile, COARSE depth eroded on [" << nX << "][" << nY << "] in Poly " << nPoly << " is " << dRemoved * m_dCellArea << endl;
581 // }
582 // // DEBUG CODE ==================================================================================================
583
584 // Store the this-polygon depth of sediment eroded during Dean profile deposition of beach unconsolidated sediment
585 if (nTexture == TEXTURE_SAND)
586 pPolygon->AddBeachSandErodedDeanProfile(dRemoved);
587
588 else
589 pPolygon->AddBeachCoarseErodedDeanProfile(dRemoved);
590 }
591 }
592 }
593
594 else if ((dElevDiff < 0) && (m_pRasterGrid->m_Cell[nX][nY].bIsInContiguousSea()))
595 {
596 if ((nTexture == TEXTURE_SAND) || (nTexture == TEXTURE_COARSE))
597 {
598 // The current elevation is below the Dean elevation, so we have can have unconsolidated sediment deposition here provided that we have some previously-eroded unconsolidated sediment (sand and coarse only) to deposit
599 if (dTotEroded > 0)
600 {
601 double dTotToDeposit = tMin(-dElevDiff, dTotEroded);
602
603 int const nTopLayer = m_pRasterGrid->m_Cell[nX][nY].nGetNumOfTopLayerAboveBasement();
604
605 // Safety check
606 if (nTopLayer == INT_NODATA)
607 return RTN_ERR_NO_TOP_LAYER;
608
609 if (dTotToDeposit > 0)
610 {
611 dTotToDeposit = tMin(dTotToDeposit, dTotEroded);
612
613 if (nTexture == TEXTURE_SAND)
614 {
615 double const dSandNow = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetSandDepth();
616 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->SetSandDepth(dSandNow + dTotToDeposit);
617
618 // Set the changed-this-timestep switch
619 m_bUnconsChangedThisIter[nTopLayer] = true;
620
621 dTotEroded -= dTotToDeposit;
622
623 dStillToErodeOnProfile += dTotToDeposit;
624 dStillToErodeOnPolygon += dTotToDeposit;
625
626 // Update per-timestep totals
627 m_ulThisIterNumBeachDepositionCells++;
628 // m_dThisIterBeachDepositionSand += dTotToDeposit;
629 }
630
631 if (nTexture == TEXTURE_COARSE)
632 {
633 double const dCoarseNow = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetCoarseDepth();
634 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->SetCoarseDepth(dCoarseNow + dTotToDeposit);
635
636 // Set the changed-this-timestep switch
637 m_bUnconsChangedThisIter[nTopLayer] = true;
638
639 dTotEroded -= dTotToDeposit;
640
641 dStillToErodeOnProfile += dTotToDeposit;
642 dStillToErodeOnPolygon += dTotToDeposit;
643
644 // Update per-timestep totals
645 m_ulThisIterNumBeachDepositionCells++;
646 m_dThisIterBeachDepositionCoarse += dTotToDeposit;
647 }
648 }
649
650 // Now update the cell's layer elevations
651 m_pRasterGrid->m_Cell[nX][nY].CalcAllLayerElevsAndD50();
652
653 // Update the cell's sea depth
654 m_pRasterGrid->m_Cell[nX][nY].SetSeaDepth();
655
656 // Update the cell's beach deposition, and total beach deposition, values
657 m_pRasterGrid->m_Cell[nX][nY].IncrBeachDeposition(dTotToDeposit);
658
659 // And set the landform category
660 CRWCellLandform* pLandform = m_pRasterGrid->m_Cell[nX][nY].pGetLandform();
661 int const nCat = pLandform->nGetLFCategory();
662
663 if ((nCat != LF_SEDIMENT_INPUT_UNCONSOLIDATED) && (nCat != LF_SEDIMENT_INPUT_CONSOLIDATED))
664 pLandform->SetLFCategory(LF_DRIFT_BEACH);
665
666 // LogStream << m_ulIter << ": nPoly = " << nPoly << ", beach deposition = " << dTotToDeposit << " at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "} nCoastPoint = " << nCoastPoint << " nDistSeawardFromNewCoast = " << nDistSeawardFromNewCoast << endl;
667 }
668 }
669 }
670 }
671 }
672
673 return RTN_OK;
674}
675
676//===============================================================================================================================
678//===============================================================================================================================
679void CSimulation::ErodeCellBeachSedimentSupplyLimited(int const nX, int const nY, int const nThisLayer, int const nTexture, double const dMaxToErode, double& dRemoved)
680{
681 // Find out how much unconsolidated sediment of this size class we have available on this cell
682 double dExistingAvailable = 0;
683 double dErodibility = 0;
684
685 if (nTexture == TEXTURE_FINE)
686 {
687 dExistingAvailable = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->dGetFineDepth();
688 dErodibility = m_dFineErodibilityNormalized;
689 }
690 else if (nTexture == TEXTURE_SAND)
691 {
692 dExistingAvailable = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->dGetSandDepth();
693 dErodibility = m_dSandErodibilityNormalized;
694 }
695 else if (nTexture == TEXTURE_COARSE)
696 {
697 dExistingAvailable = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->dGetCoarseDepth();
698 dErodibility = m_dCoarseErodibilityNormalized;
699 }
700
701 // Is there any unconsolidated sediment on this cell?
702 if (dExistingAvailable <= 0)
703 {
704 return;
705 }
706
707 // Erode some unconsolidated sediment
708 double const dLowering = dErodibility * dMaxToErode;
709
710 // Make sure we don't get -ve amounts left on the cell
711 dRemoved = tMin(dExistingAvailable, dLowering);
712 double const dRemaining = dExistingAvailable - dRemoved;
713
714 if (nTexture == TEXTURE_FINE)
715 {
716 // Set the value for this layer
717 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->SetFineDepth(dRemaining);
718 }
719 else if (nTexture == TEXTURE_SAND)
720 {
721 // Set the value for this layer
722 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->SetSandDepth(dRemaining);
723 }
724 else if (nTexture == TEXTURE_COARSE)
725 {
726 // Set the value for this layer
727 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->SetCoarseDepth(dRemaining);
728 }
729
730 // And set the changed-this-timestep switch
731 m_bUnconsChangedThisIter[nThisLayer] = true;
732
733 // Set the actual erosion value for this cell
734 m_pRasterGrid->m_Cell[nX][nY].SetActualBeachErosion(dRemoved);
735
736 if (dRemoved > 0)
737 {
738 // Recalculate the elevation of every layer
739 m_pRasterGrid->m_Cell[nX][nY].CalcAllLayerElevsAndD50();
740
741 // And update the cell's sea depth
742 m_pRasterGrid->m_Cell[nX][nY].SetSeaDepth();
743 }
744}
745
746//===============================================================================================================================
748//===============================================================================================================================
749int CSimulation::nDoUnconsDepositionOnPolygon(int const nCoast, CGeomCoastPolygon* pPolygon, int const nTexture, double dTargetToDepositOnPoly, double& dDepositedOnPoly)
750{
751 // // DEBUG CODE #####################
752 // if (m_ulIter == 1)
753 // {
754 // int nPoly = pPolygon->nGetPolygonCoastID();
755 // LogStream << m_ulIter << ": entered nDoUnconsDepositionOnPolygon() nCoast = " << nCoast << " nPoly = " << nPoly << " dTargetToDepositOnPoly = " << dTargetToDepositOnPoly * m_dCellArea << " dDepositedOnPoly = " << dDepositedOnPoly * m_dCellArea << endl;
756 // }
757 // // DEBUG CODE #####################
758
759 string strTexture;
760
761 if (nTexture == TEXTURE_SAND)
762 strTexture = "sand";
763
764 else if (nTexture == TEXTURE_COARSE)
765 strTexture = "coarse";
766
767 // Don't bother with tiny amounts of deposition
768 if (dTargetToDepositOnPoly < SED_ELEV_TOLERANCE)
769 return RTN_OK;
770
771 // Get the grid cell coordinates of this polygon's up-coast and down-coast profiles
772 int const nUpCoastProfile = pPolygon->nGetUpCoastProfile();
773 int const nDownCoastProfile = pPolygon->nGetDownCoastProfile();
774
775 CGeomProfile* pUpCoastProfile = m_VCoast[nCoast].pGetProfile(nUpCoastProfile);
776 CGeomProfile* pDownCoastProfile = m_VCoast[nCoast].pGetProfile(nDownCoastProfile);
777
778 // We are using only part of each profile, seaward as far as the depth of closure. First find the seaward end point of the up-coast part-profile
779 // CGeom2DIPoint PtiUpCoastPartProfileSeawardEnd;
780 // int nIndex = pUpCoastProfile->nGetCellGivenDepth(m_pRasterGrid, m_dDepthOfClosure, &PtiUpCoastPartProfileSeawardEnd);
781 int const nIndex = pUpCoastProfile->nGetCellGivenDepth(m_pRasterGrid, m_dDepthOfClosure);
782
783 if (nIndex == INT_NODATA)
784 {
785 LogStream << m_ulIter << ": " << ERR << "while depositing " + strTexture + " unconsolidated sediment for coast " << nCoast << " polygon " << pPolygon->nGetPolygonCoastID() << ", could not find the seaward end point of the up-coast profile (" << nUpCoastProfile << ") for depth of closure = " << m_dDepthOfClosure << endl;
786
787 return RTN_ERR_NO_SEAWARD_END_OF_PROFILE_UPCOAST_BEACH_DEPOSITION;
788 }
789
790 // The part-profile length is one greater than nIndex, since pPtiGetCellGivenDepth() returns the index of the cell at depth of closure. This will be the number of cells in the Dean profile portion of every parallel profile
791 int const nUpCoastDeanLen = nIndex + 1;
792
793 // assert(bIsWithinValidGrid(&PtiUpCoastPartProfileSeawardEnd));
794
795 // Get the distance between the start and end of the part-profile (the Dean length), in external CRS units
796 // CGeom2DPoint
797 // PtUpCoastProfileStart = *pUpCoastProfile->pPtGetPointInProfile(0),
798 // PtUpCoastProfileEnd = PtGridCentroidToExt(&PtiUpCoastPartProfileSeawardEnd);
799
800 // double dUpCoastDeanLen = dGetDistanceBetween(&PtUpCoastProfileStart, &PtUpCoastProfileEnd);
801
802 int const nUpCoastProfileCoastPoint = pUpCoastProfile->nGetCoastPoint();
803 int const nDownCoastProfileCoastPoint = pDownCoastProfile->nGetCoastPoint();
804 int const nXUpCoastProfileExistingCoastPoint = m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nUpCoastProfileCoastPoint)->nGetX();
805 int const nYUpCoastProfileExistingCoastPoint = m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nUpCoastProfileCoastPoint)->nGetY();
806 int nCoastSegLen;
807
808 // Store the coast point numbers for this polygon so that we can shuffle them
809 vector<int> nVCoastPoint;
810
811 if (nDownCoastProfileCoastPoint == m_VCoast[nCoast].nGetCoastlineSize() - 1)
812 {
813 // This is the final down-coast polygon, so also include the down-coast polygon boundary
814 nCoastSegLen = nDownCoastProfileCoastPoint - nUpCoastProfileCoastPoint + 1;
815
816 for (int nCoastPoint = nUpCoastProfileCoastPoint; nCoastPoint <= nDownCoastProfileCoastPoint; nCoastPoint++)
817 nVCoastPoint.push_back(nCoastPoint);
818 }
819 else
820 {
821 // This is not the final down-coast polygon, so do not include the polygon's down-coast boundary
822 nCoastSegLen = nDownCoastProfileCoastPoint - nUpCoastProfileCoastPoint;
823
824 for (int nCoastPoint = nUpCoastProfileCoastPoint; nCoastPoint < nDownCoastProfileCoastPoint; nCoastPoint++)
825 nVCoastPoint.push_back(nCoastPoint);
826 }
827
828 double dStillToDepositOnPoly = dTargetToDepositOnPoly;
829 double dTargetToDepositOnProfile = dTargetToDepositOnPoly / nCoastSegLen;
830 double dStillToDepositOnProfile; // = dTargetToDepositOnProfile;
831
832 // Shuffle the coast points, this is necessary so that leaving the loop does not create sequence-related artefacts
833 shuffle(nVCoastPoint.begin(), nVCoastPoint.begin() + nCoastSegLen, m_Rand[1]);
834
835 // // Get the volume of sediment which is to be deposited on the polygon and on each parallel profile. Note that if dSandToMoveOnPoly is -ve, then don't do any sand deposition. Similarly if dCoarseToMoveOnPoly is -ve then don't do any coarse deposition
836 // double
837 // dSandToMoveOnPoly = pPolygon->dGetToDoBeachDepositionUnconsSand(), // +ve is deposition, -ve is erosion
838 // dCoarseToMoveOnPoly = pPolygon->dGetToDoBeachDepositionUnconsCoarse(); // +ve is deposition, -ve is erosion
839 //
840 // double
841 // dTmpLower = 0,
842 // dSandRatio = 0,
843 // dCoarseRatio = 0,
844 // dSandTargetToDepositOnPoly = 0,
845 // dCoarseTargetToDepositOnPoly = 0;
846 // if (dSandToMoveOnPoly > 0)
847 // dTmpLower += dSandToMoveOnPoly;
848 // if (dCoarseToMoveOnPoly > 0)
849 // dTmpLower += dCoarseToMoveOnPoly;
850 //
851 // if ((dSandToMoveOnPoly > 0) && (dTmpLower > 0))
852 // dSandRatio = dSandToMoveOnPoly / dTmpLower;
853 // if (dCoarseToMoveOnPoly > 0)
854 // dCoarseRatio = 1 - dSandRatio;
855 //
856 // // LogStream << "####### nPoly = " << nPoly << " dSandRatio = " << dSandRatio << " dCoarseRatio = " << dCoarseRatio << endl;
857 //
858 // dSandTargetToDepositOnPoly = dAllTargetToDepositOnPoly * dSandRatio,
859 // dCoarseTargetToDepositOnPoly = dAllTargetToDepositOnPoly * dCoarseRatio;
860 //
861 // double
862 // dAllTargetPerProfile = dTargetToDepositOnPoly / nCoastSegLen,
863 // dTargetStillToDepositOnProfile = dTargetToDepositOnPoly / nCoastSegLen,
864 // dSandDepositedOnPoly = 0, // Grand total for the whole polygon
865 // dCoarseDepositedOnPoly = 0; // Grand total for the whole polygon
866 //
867 // // LogStream << "####### nPoly = " << nPoly << " dSandToMoveOnPoly = " << dSandToMoveOnPoly << " dCoarseToMoveOnPoly = " << dCoarseToMoveOnPoly << endl;
868 //
869 // // LogStream << "####### nPoly = " << nPoly << " dSandTargetToDepositOnPoly = " << dSandTargetToDepositOnPoly << " dCoarseTargetToDepositOnPoly = " << dCoarseTargetToDepositOnPoly << " dAllTargetToDepositOnPoly = " << dAllTargetToDepositOnPoly << endl;
870
871 // Now traverse the polygon's existing coastline in a random (but broadly DOWN-COAST i.e. increasing coast point indices) sequence, fitting a Dean profile at each coast point
872 // DOWN-COAST (increasing coastpoint indices) ================================================================================
873 // LogStream << "####### DOWN-COAST nPoly = " << nPoly << " nCoastSegLen = " << nCoastSegLen << endl;
874 for (int n = 0; n < nCoastSegLen; n++)
875 {
876 // Pick a random coast point
877 int const nCoastPoint = nVCoastPoint[n];
878 CGeom2DIPoint const PtiCoastPoint = *m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nCoastPoint);
879 int const nCoastX = PtiCoastPoint.nGetX();
880 int const nCoastY = PtiCoastPoint.nGetY();
881
882 // Calculate the x-y offset between this coast point, and the coast point of the up-coast normal
883 int const nXOffset = nCoastX - nXUpCoastProfileExistingCoastPoint;
884 int const nYOffset = nCoastY - nYUpCoastProfileExistingCoastPoint;
885 int nSeawardOffset = -1;
886 // nParProfLen;
887 vector<CGeom2DIPoint> PtiVParProfile;
888 vector<double> VdParProfileDeanElev;
889
890 // OK, loop until we can deposit sufficient unconsolidated sediment on the parallel profile starting at this coast point
891 while (true)
892 {
893 // Move seaward by one cell
894 nSeawardOffset++;
895
896 // And lengthen the parallel profile
897 int nParProfLen = nUpCoastDeanLen + nSeawardOffset;
898
899 if (nParProfLen > (pUpCoastProfile->nGetNumCellsInProfile()))
900 {
901 // We've reached the seaward end of the up-coast profile, need to quit
902 // if (m_nLogFileDetail >= LOG_FILE_MIDDLE_DETAIL)
903 // LogStream << m_ulIter << ": " << WARN << "reached seaward end of up-coast profile during DOWN-COAST deposition of unconsolidated sediment for coast " << nCoast << " polygon " << pPolygon->nGetPolygonCoastID() << " (nCoastPoint = " << nCoastPoint << " nSeawardOffset = " << nSeawardOffset << ")" << endl;
904
905 break;
906 }
907
908 // Get the x-y coords of a profile starting from this coast point and parallel to the up-coast polygon boundary profile (these are in natural sequence, like the boundary part-profile)
909 PtiVParProfile.resize(0);
910
911 for (int m = 0; m < nParProfLen; m++)
912 {
913 CGeom2DIPoint const PtiProf = *pUpCoastProfile->pPtiGetCellInProfile(m);
914 CGeom2DIPoint const PtiTmp(PtiProf.nGetX() + nXOffset, PtiProf.nGetY() + nYOffset);
915 PtiVParProfile.push_back(PtiTmp);
916 }
917
918 // Get the existing elevation of the seaward end of the parallel profile
919 int nSeaEndX = PtiVParProfile.back().nGetX();
920 int nSeaEndY = PtiVParProfile.back().nGetY();
921
922 // Safety check
923 if (! bIsWithinValidGrid(nSeaEndX, nSeaEndY))
924 {
925 // if (m_nLogFileDetail >= LOG_FILE_MIDDLE_DETAIL)
926 // LogStream << WARN << "09 @@@@ while doing DOWN-COAST deposition on polygon " << nPoly << ", hit edge of grid at [" << nSeaEndX << "][" << nSeaEndY << "] for parallel profile from coast point " << nCoastPoint << " at [" << nCoastX << "][" << nCoastY << "]. Constraining this parallel profile at its seaward end" << endl;
927
928 KeepWithinValidGrid(nCoastX, nCoastY, nSeaEndX, nSeaEndY);
929 PtiVParProfile.back().SetX(nSeaEndX);
930 PtiVParProfile.back().SetY(nSeaEndY);
931 }
932
933 double const dParProfEndElev = m_pRasterGrid->m_Cell[nSeaEndX][nSeaEndY].dGetAllSedTopElevOmitTalus();
934
935 // Set the start elevation for the Dean profile just a bit above mean SWL for this timestep (i.e. so that it is a Bruun profile)
936 double const dParProfStartElev = m_dThisIterMeanSWL + m_dDeanProfileStartAboveSWL;
937
938 // Calculate the total length of the parallel profile, including any seaward offset
939 double const dParProfLen = dGetDistanceBetween(&PtiVParProfile.front(), &PtiVParProfile.back());
940
941 // Now calculate the length of the Dean profile-only part i.e. without any seaward offset. The approach used here is approximate but probably OK
942 double dParProfDeanLen = dParProfLen - (nSeawardOffset * m_dCellSide);
943
944 // Safety check
945 if (dParProfDeanLen <= 0)
946 dParProfDeanLen = 1;
947
948 // Solve for dA so that the existing elevations at the end of the parallel profile, and at the end of a Dean equilibrium profile on that part-normal, are the same
949 double const dParProfA = (dParProfStartElev - dParProfEndElev) / pow(dParProfDeanLen, DEAN_POWER);
950
951 nParProfLen = static_cast<int>(PtiVParProfile.size());
952 VdParProfileDeanElev.resize(nParProfLen, 0);
953
954 // for (int m = 0; m < static_cast<int>(PtiVParProfile.size()); m++)
955 // LogStream << "[" << PtiVParProfile[m].nGetX() << "][" << PtiVParProfile[m].nGetY() << "] ";
956 // LogStream << endl;
957
958 double const dInc = dParProfDeanLen / (nParProfLen - nSeawardOffset - 2);
959
960 // The elevation of the coast point in the Dean profile is the same as the elevation of the current coast point, ignoring any talus (since talus is assumed to be removed quickly) TODO 020 Is this correct? Should it be dParProfStartElev?
961 double const dCoastElev = m_pRasterGrid->m_Cell[nCoastX][nCoastY].dGetAllSedTopElevOmitTalus();
962
963 // For this depositing parallel profile, calculate the Dean equilibrium profile of the unconsolidated sediment h(y) = A * y^(2/3) where h(y) is the distance below the highest point in the profile at a distance y from the landward start of the profile
964 CalcDeanProfile(&VdParProfileDeanElev, dInc, dParProfStartElev, dParProfA, true, nSeawardOffset, dCoastElev);
965
966 double dParProfTotDiff = 0;
967
968 for (int m = 0; m < nParProfLen; m++)
969 {
970 CGeom2DIPoint PtiTmp = PtiVParProfile[m];
971 int nX = PtiTmp.nGetX();
972 int nY = PtiTmp.nGetY();
973
974 // Safety check
975 if (! bIsWithinValidGrid(nX, nY))
976 {
977 // LogStream << WARN << "10 @@@@ while doing DOWN-COAST deposition on polygon " << nPoly << ", hit edge of grid at [" << nX << "][" << nY << "] for parallel profile from coast point " << nCoastPoint << " at [" << nCoastX << "][" << nCoastY << "]. Constraining this parallel profile at its seaward end" << endl;
978
979 KeepWithinValidGrid(nCoastX, nCoastY, nX, nY);
980 PtiTmp.SetX(nX);
981 PtiTmp.SetY(nY);
982 }
983
984 // Don't do anything to intervention cells
985 if (bIsInterventionCell(nX, nY))
986 continue;
987
988 // Don't do cells twice
989 if (! m_pRasterGrid->m_Cell[nX][nY].bBeachErosionOrDepositionThisIter())
990 {
991 double const dTmpElev = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus();
992 double const dDiff = VdParProfileDeanElev[m] - dTmpElev;
993
994 dParProfTotDiff += dDiff;
995 }
996 }
997
998 // // DEBUG CODE -----------------------------------------------------
999 // LogStream << endl << "\tFor polygon " << nPoly << " doing DOWN-COAST deposition, nSeawardOffset = " << nSeawardOffset << ", parallel profile from [" << PtiVParProfile[0].nGetX() << "][" << PtiVParProfile[0].nGetY() << "] to [" << PtiVParProfile.back().nGetX() << "][" << PtiVParProfile.back().nGetY() << "], nUpCoastDeanLen = " << nUpCoastDeanLen << " dUpCoastDeanLen = " << dUpCoastDeanLen << " nParProfLen = " << nParProfLen << " dParProfDeanLen = " << dParProfDeanLen << " dInc = " << dInc << " dParProfStartElev = " << dParProfStartElev << " dParProfEndElev = " << dParProfEndElev << " dParProfA = " << dParProfA << endl;
1000 //
1001 // LogStream << "\tExisting profile for deposition = ";
1002 // for (int n = 0; n < nParProfLen; n++)
1003 // {
1004 // CGeom2DIPoint PtiTmp = PtiVParProfile[n];
1005 // int
1006 // nX = PtiTmp.nGetX(),
1007 // nY = PtiTmp.nGetY();
1008 //
1009 // // Safety check
1010 // if (! bIsWithinValidGrid(nX, nY))
1011 // KeepWithinValidGrid(nX, nY);
1012 //
1013 // LogStream << m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus() << " ";
1014 // }
1015 // LogStream << endl;
1016 // LogStream << "\tParallel Dean equilibrium profile for deposition = ";
1017 // for (int n = 0; n < nParProfLen; n++)
1018 // {
1019 // LogStream << VdParProfileDeanElev[n] << " ";
1020 // }
1021 // LogStream << endl;
1022 //
1023 // LogStream << "\tnCoastPoint = " << nCoastPoint << " nSeawardOffset = " << nSeawardOffset << " difference = ";
1024 // for (int n = 0; n < nParProfLen; n++)
1025 // {
1026 // CGeom2DIPoint PtiTmp = PtiVParProfile[n];
1027 // int
1028 // nX = PtiTmp.nGetX(),
1029 // nY = PtiTmp.nGetY();
1030 //
1031 // // Safety check
1032 // if (! bIsWithinValidGrid(nX, nY))
1033 // KeepWithinValidGrid(nX, nY);
1034 //
1035 // double
1036 // dTmpElev = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus(),
1037 // dDiff = VdParProfileDeanElev[n] - dTmpElev;
1038 //
1039 // LogStream << dDiff << " ";
1040 // }
1041 // LogStream << endl;
1042 // // END DEBUG CODE -----------------------------------------------------
1043
1044 // So will we be able to deposit as much as is needed?
1045 if (dParProfTotDiff >= dTargetToDepositOnProfile)
1046 {
1047 // LogStream << " DOWN-COAST nPoly = " << pPolygon->nGetPolygonCoastID() << " nCoastPoint = " << nCoastPoint << " nSeawardOffset = " << nSeawardOffset << " dParProfTotDiff = " << dParProfTotDiff << " dTargetToDepositOnProfile = " << dTargetToDepositOnProfile << " WILL BE ABLE TO DEPOSIT ENOUGH" << endl;
1048
1049 break;
1050 }
1051 }
1052
1053 // assert(dParProfTotDiff > 0);
1054
1055 // OK, this value of nSeawardOffset gives us enough deposition. So start depositing on the parallel profile from the coastward end
1056 double dDepositedOnProfile = 0; // Total for this parallel profile
1057
1058 // LogStream << " DOWN-COAST nPoly = " << nPoly << " nCoastPoint = " << nCoastPoint << " doing deposition on parallel profile, dSandTargetStillToDepositOnProfile = " << dSandTargetStillToDepositOnProfile << " dCoarseTargetToDepositOnProfile = " << dCoarseTargetToDepositOnProfile << endl;
1059
1060 dStillToDepositOnProfile = dTargetToDepositOnProfile;
1061
1062 for (unsigned int nSeawardFromCoast = 0; nSeawardFromCoast < PtiVParProfile.size(); nSeawardFromCoast++)
1063 {
1064 // Move along this parallel profile starting from the coast. Leave the loop if we have deposited enough on this polygon
1065 if (dStillToDepositOnPoly < SED_ELEV_TOLERANCE)
1066 {
1067 // LogStream << m_ulIter << ": nPoly = " << nPoly << " texture = " << strTexture << " DOWN-COAST nCoastPoint = " << nCoastPoint << " nSeawardOffset = " << nSeawardOffset << " leaving loop because enough deposited on polygon, dTargetToDepositOnPoly = " << dTargetToDepositOnPoly << " dDepositedOnPoly = " << dDepositedOnPoly << endl;
1068
1069 break;
1070 }
1071
1072 // Leave the loop if we have deposited enough on this parallel profile
1073 if (dStillToDepositOnProfile < SED_ELEV_TOLERANCE)
1074 {
1075 // LogStream << m_ulIter << ": nPoly = " << nPoly << " texture = " << strTexture << " DOWN-COAST nCoastPoint = " << nCoastPoint << " nSeawardOffset = " << nSeawardOffset << " leaving loop because enough deposited on parallel profile, dTargetStillToDepositOnProfile = " << dTargetStillToDepositOnProfile << " dDepositedOnProfile = " << dDepositedOnProfile << endl;
1076
1077 break;
1078 }
1079
1080 // assert(nSeawardFromCoast < PtiVParProfile.size());
1081 CGeom2DIPoint PtiTmp = PtiVParProfile[nSeawardFromCoast];
1082 int nX = PtiTmp.nGetX();
1083 int nY = PtiTmp.nGetY();
1084
1085 // Safety check
1086 if (! bIsWithinValidGrid(nX, nY))
1087 {
1088 // LogStream << WARN << "11 @@@@ while doing DOWN-COAST deposition on polygon " << nPoly << ", hit edge of grid at [" << nX << "][" << nY << "] for parallel profile from coast point " << nCoastPoint << " at [" << nCoastX << "][" << nCoastY << "]. Constraining this parallel profile at its seaward end" << endl;
1089
1090 KeepWithinValidGrid(nCoastX, nCoastY, nX, nY);
1091 PtiTmp.SetX(nX);
1092 PtiTmp.SetY(nY);
1093 }
1094
1095 // Don't do anything to intervention cells
1096 if (bIsInterventionCell(nX, nY))
1097 continue;
1098
1099 // Don't do cells twice
1100 if (! m_pRasterGrid->m_Cell[nX][nY].bBeachErosionOrDepositionThisIter())
1101 {
1102 // Get this cell's current elevation
1103 double const dThisElevNow = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus();
1104
1105 // LogStream << "\tnPoly = " << nPoly << ", [" << nX << "][" << nY << "] nCoastPoint = " << nCoastPoint << " nSeawardFromCoast = " << nSeawardFromCoast << " dThisElevNow = " << dThisElevNow << " Dean Elev = " << VdParProfileDeanElev[nSeawardFromCoast] << endl;
1106
1107 // Subtract the two elevations
1108 double const dElevDiff = VdParProfileDeanElev[nSeawardFromCoast] - dThisElevNow;
1109 CRWCellLandform* pLandform = m_pRasterGrid->m_Cell[nX][nY].pGetLandform();
1110
1111 if (dElevDiff > SED_ELEV_TOLERANCE)
1112 {
1113 bool bDeposited = false;
1114 double dToDepositHere = 0;
1115
1116 // The current elevation is below the Dean elevation, so we have can have beach deposition here
1117 if (dStillToDepositOnProfile > SED_ELEV_TOLERANCE)
1118 {
1119 dToDepositHere = tMin(dElevDiff, dStillToDepositOnProfile, dStillToDepositOnPoly);
1120
1121 // LogStream << " DOWN-COAST nPoly = " << nPoly << " nCoastPoint = " << nCoastPoint << " texture = " << strTexture << " dToDepositHere = " << dToDepositHere << endl;
1122
1123 int const nTopLayer = m_pRasterGrid->m_Cell[nX][nY].nGetNumOfTopLayerAboveBasement();
1124
1125 // Safety check
1126 if (nTopLayer == INT_NODATA)
1127 return RTN_ERR_NO_TOP_LAYER;
1128
1129 if (dToDepositHere > SED_ELEV_TOLERANCE)
1130 {
1131 bDeposited = true;
1132
1133 if (nTexture == TEXTURE_SAND)
1134 {
1135 double const dSandNow = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetSandDepth();
1136
1137 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->SetSandDepth(dSandNow + dToDepositHere);
1138
1139 // Set the changed-this-timestep switch
1140 m_bUnconsChangedThisIter[nTopLayer] = true;
1141
1142 dDepositedOnProfile += dToDepositHere;
1143 dDepositedOnPoly += dToDepositHere;
1144
1145 // LogStream << "XXXXXXX texture = " << strTexture << " dDepositedOnProfile = " << dDepositedOnProfile << " dDepositedOnPoly = " << dDepositedOnPoly << endl;
1146
1147 // Update the cell's beach deposition, and total beach deposition, values
1148 m_pRasterGrid->m_Cell[nX][nY].IncrBeachDeposition(dToDepositHere);
1149
1150 dStillToDepositOnPoly -= dToDepositHere;
1151 dStillToDepositOnProfile -= dToDepositHere;
1152
1153 // if (dDepositedOnPoly > dTargetToDepositOnPoly)
1154 // {
1155 // // LogStream << "££££££ 1 nPoly = " << nPoly << " nCoastPoint = " << nCoastPoint << " dDepositedOnProfile = " << dDepositedOnProfile << " dDepositedOnPoly = " << dDepositedOnPoly << " dTargetToDepositOnPoly = " << dTargetToDepositOnPoly << endl;
1156 // }
1157
1158 // LogStream << m_ulIter << ": nPoly = " << nPoly << " texture = " << strTexture << " dToDepositHere = " << dToDepositHere << " at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "} nCoastPoint = " << nCoastPoint << " nSeawardFromCoast = " << nSeawardFromCoast << endl;
1159 }
1160 else if (nTexture == TEXTURE_COARSE)
1161 {
1162 double const dCoarseNow = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetCoarseDepth();
1163
1164 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->SetCoarseDepth(dCoarseNow + dToDepositHere);
1165
1166 // Set the changed-this-timestep switch
1167 m_bUnconsChangedThisIter[nTopLayer] = true;
1168
1169 // if (dDepositedOnPoly > dTargetToDepositOnPoly)
1170 // {
1171 // // LogStream << "££££££ 2 nPoly = " << nPoly << " dDepositedOnPoly = " << dDepositedOnPoly << " dTargetToDepositOnPoly = " << dTargetToDepositOnPoly << endl;
1172 // }
1173
1174 dDepositedOnProfile += dToDepositHere;
1175 dDepositedOnPoly += dToDepositHere;
1176
1177 // LogStream << "XXXXXXX texture = " << strTexture << " dDepositedOnProfile = " << dDepositedOnProfile << " dDepositedOnPoly = " << dDepositedOnPoly << endl;
1178
1179 // Update the cell's beach deposition, and total beach deposition, values
1180 m_pRasterGrid->m_Cell[nX][nY].IncrBeachDeposition(dToDepositHere);
1181
1182 dStillToDepositOnPoly -= dToDepositHere;
1183 dStillToDepositOnProfile -= dToDepositHere;
1184
1185 // LogStream << m_ulIter << ": nPoly = " << nPoly << " texture = " << strTexture << " dToDepositHere = " << dToDepositHere << " at [" << nX << "][" << nY << "] nCoastPoint = " << nCoastPoint << " nSeawardFromCoast = " << nSeawardFromCoast << endl;
1186 }
1187 }
1188 }
1189
1190 // if (dDepositedOnPoly > dTargetToDepositOnPoly)
1191 // {
1192 // // LogStream << "££££££ 3 nPoly = " << nPoly << " texture = " << strTexture << " dDepositedOnPoly = " << dDepositedOnPoly << " dTargetToDepositOnPoly = " << dTargetToDepositOnPoly << endl;
1193 // }
1194
1195 if (bDeposited)
1196 {
1197 // Update the cell's layer elevations
1198 m_pRasterGrid->m_Cell[nX][nY].CalcAllLayerElevsAndD50();
1199
1200 // Update the cell's sea depth
1201 m_pRasterGrid->m_Cell[nX][nY].SetSeaDepth();
1202
1203 // And set the landform category
1204 int const nCat = pLandform->nGetLFCategory();
1205
1206 if ((nCat != LF_SEDIMENT_INPUT_CONSOLIDATED) && (nCat != LF_SEDIMENT_INPUT_UNCONSOLIDATED))
1207 pLandform->SetLFCategory(LF_DRIFT_BEACH);
1208
1209 // Update this-timestep totals
1210 m_ulThisIterNumBeachDepositionCells++;
1211
1212 if (nTexture == TEXTURE_SAND)
1213 m_dThisIterBeachDepositionSand += dToDepositHere;
1214
1215 else if (nTexture == TEXTURE_COARSE)
1216 m_dThisIterBeachDepositionCoarse += dToDepositHere;
1217 }
1218 }
1219 else if (bElevAboveDeanElev(nX, nY, dElevDiff, pLandform))
1220 {
1221 // The current elevation is higher than the Dean elevation, so we have potential beach erosion (i.e. not constrained by availability of unconsolidated sediment) here
1222 m_ulThisIterNumPotentialBeachErosionCells++;
1223
1224 m_pRasterGrid->m_Cell[nX][nY].SetPotentialBeachErosion(-dElevDiff);
1225
1226 // LogStream << m_ulIter << ": nPoly = " << nPoly << " texture = " << strTexture << " going DOWN-COAST, potential beach erosion = " << -dElevDiff << " at [" << nX << "][" << nY << "] nCoastPoint = " << nCoastPoint << " nSeawardFromCoast = " << nSeawardFromCoast << " dThisElevNow = " << dThisElevNow << " Dean Elev = " << VdParProfileDeanElev[nSeawardFromCoast] << endl;
1227
1228 // Now get the number of the highest layer with non-zero thickness
1229 int const nThisLayer = m_pRasterGrid->m_Cell[nX][nY].nGetTopNonZeroLayerAboveBasement();
1230
1231 // Safety check
1232 if (nThisLayer == INT_NODATA)
1233 return RTN_ERR_NO_TOP_LAYER;
1234
1235 if (nThisLayer != NO_NONZERO_THICKNESS_LAYERS)
1236 {
1237 // We still have at least one layer left with non-zero thickness (i.e. we are not down to basement), and the cell's current elevation is higher than the Dean equilibrium profile elevation. So do some beach erosion on this cell (Note: we are in nDoUnconsDepositionOnPolygon() still)
1238 if (nTexture == TEXTURE_SAND)
1239 {
1240 double dSandRemoved = 0;
1241 ErodeCellBeachSedimentSupplyLimited(nX, nY, nThisLayer, TEXTURE_SAND, -dElevDiff, dSandRemoved);
1242
1243 // Update totals for this parallel profile
1244 dDepositedOnProfile -= dSandRemoved;
1245 dStillToDepositOnProfile += dSandRemoved;
1246
1247 // Update totals for the polygon (note that these may become slightly -ve if erosion occurs early during this routine, but this is not a problem since the values will eventually become +ve again
1248 dDepositedOnPoly -= dSandRemoved;
1249 dStillToDepositOnPoly += dSandRemoved;
1250
1251 // if (dDepositedOnPoly < 0)
1252 // LogStream << m_ulIter << ": BBB LESS THAN ZERO dDepositedOnPoly = " << dDepositedOnPoly << endl;
1253
1254 // Update this-timestep totals
1255 m_ulThisIterNumActualBeachErosionCells++;
1256
1257 // // DEBUG CODE ==================================================================================================
1258 // LogStream << m_ulIter << ": $$$$$$$$$ BEACH 1 Dean profile lower than existing profile, SAND depth eroded on [" << nX << "][" << nY << "] in Poly " << nPoly << " is " << dSandRemoved * m_dCellArea << endl;
1259 // // DEBUG CODE ==================================================================================================
1260
1261 // Store the this-polygon depth of sand sediment eroded during Dean profile deposition of beach unconsolidated sediment
1262 pPolygon->AddBeachSandErodedDeanProfile(dSandRemoved);
1263 }
1264 else if (nTexture == TEXTURE_COARSE)
1265 {
1266 double dCoarseRemoved = 0;
1267 ErodeCellBeachSedimentSupplyLimited(nX, nY, nThisLayer, TEXTURE_COARSE, -dElevDiff, dCoarseRemoved);
1268
1269 // Update totals for this parallel profile
1270 dDepositedOnProfile -= dCoarseRemoved;
1271 dStillToDepositOnProfile += dCoarseRemoved;
1272
1273 // Update totals for the polygon (note that these may become slightly -ve if erosion occurs early during this routine, but this is not a problem since the values will eventually become +ve again
1274 dDepositedOnPoly -= dCoarseRemoved;
1275 dStillToDepositOnPoly += dCoarseRemoved;
1276
1277 // Update this-timestep totals
1278 m_ulThisIterNumActualBeachErosionCells++;
1279
1280 // // DEBUG CODE ==================================================================================================
1281 // LogStream << m_ulIter << ": $$$$$$$$$ BEACH 1 Dean profile lower than existing profile, COARSE depth eroded on [" << nX << "][" << nY << "] in Poly " << nPoly << " is " << dCoarseRemoved * m_dCellArea << endl;
1282 // // DEBUG CODE ==================================================================================================
1283
1284 // Store the this-polygon depth of coarse sediment eroded during Dean profile deposition of beach unconsolidated sediment
1285 pPolygon->AddBeachCoarseErodedDeanProfile(dCoarseRemoved);
1286 }
1287 }
1288
1289 // if (dDepositedOnPoly > dTargetToDepositOnPoly)
1290 // {
1291 // LogStream << "££££££ 4 nPoly = " << nPoly << " texture = " << strTexture << " dDepositedOnPoly = " << dDepositedOnPoly << " dTargetToDepositOnPoly = " << dTargetToDepositOnPoly << endl;
1292 // }
1293
1294 // LogStream << m_ulIter << ": in polygon " << nPoly << " have net deposition for " << strTexture << ", actual beach erosion = " << dSand + dCoarse << " at [" << nX << "][" << nY << "]" << endl;
1295 }
1296 }
1297 }
1298
1299 // LogStream << m_ulIter << ": nPoly = " << nPoly << " texture = " << strTexture << " nCoastPoint = " << nCoastPoint << " nSeawardOffset = " << nSeawardOffset << " dTargetStillToDepositOnProfile = " << dTargetStillToDepositOnProfile << endl;
1300 }
1301
1302 // Have we deposited enough?
1303 if (dTargetToDepositOnPoly > 0)
1304 {
1305 // LogStream << "##### texture = " << strTexture << " dDepositedOnPoly = " << dDepositedOnPoly << " dTargetToDepositOnPoly = " << dTargetToDepositOnPoly << " Not enough deposited, so now going UP-COAST" << endl;
1306 // No, so do the same in an UP-COAST (i.e. decreasing coastpoint indices) direction
1307 // UP-COAST (decreasing coastpoint indices) ==================================================================================
1308
1309 // Start by finding the seaward end point of the down-coast part-profile, as before we are using only part of each profile, seaward as far as the depth of closure
1310 // CGeom2DIPoint PtiDownCoastPartProfileSeawardEnd;
1311 // int nIndex1 = pDownCoastProfile->nGetCellGivenDepth(m_pRasterGrid, m_dDepthOfClosure, &PtiDownCoastPartProfileSeawardEnd);
1312 int const nIndex1 = pDownCoastProfile->nGetCellGivenDepth(m_pRasterGrid, m_dDepthOfClosure);
1313
1314 if (nIndex1 == INT_NODATA)
1315 {
1316 LogStream << m_ulIter << ": " << ERR << "while depositing beach for coast " << nCoast << " polygon " << pPolygon->nGetPolygonCoastID() << ", could not find the seaward end point of the up-coast profile (" << nUpCoastProfile << ") for depth of closure = " << m_dDepthOfClosure << endl;
1317
1318 return RTN_ERR_NO_SEAWARD_END_OF_PROFILE_DOWNCOAST_BEACH_DEPOSITION;
1319 }
1320
1321 // The part-profile length is one greater than nIndex1, since pPtiGetCellGivenDepth() returns the index of the cell at depth of closure. This will be the number of cells in the Dean profile portion of every parallel profile
1322 int const nDownCoastDeanLen = nIndex1 + 1;
1323
1324 // assert(bIsWithinValidGrid(&PtiDownCoastPartProfileSeawardEnd));
1325
1326 // Get the distance between the start and end of the part-profile (the Dean length), in external CRS units
1327 // CGeom2DPoint
1328 // PtDownCoastProfileStart = *pDownCoastProfile->pPtGetPointInProfile(0),
1329 // PtDownCoastProfileEnd = PtGridCentroidToExt(&PtiDownCoastPartProfileSeawardEnd);
1330
1331 // double dDownCoastDeanLen = dGetDistanceBetween(&PtDownCoastProfileStart, &PtDownCoastProfileEnd);
1332
1333 int const nXDownCoastProfileExistingCoastPoint = m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nDownCoastProfileCoastPoint)->nGetX();
1334 int const nYDownCoastProfileExistingCoastPoint = m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nDownCoastProfileCoastPoint)->nGetY();
1335
1336 // int nCoastSegLen;
1337
1338 // Store the coast point numbers for this polygon so that we can shuffle them
1339 nVCoastPoint.resize(0);
1340
1341 if (nUpCoastProfileCoastPoint == 0)
1342 {
1343 // This is the final up-coast polygon, so also include the up-coast polygon boundary
1344 nCoastSegLen = nDownCoastProfileCoastPoint - nUpCoastProfileCoastPoint + 1;
1345
1346 for (int nCoastPoint = nUpCoastProfileCoastPoint; nCoastPoint <= nDownCoastProfileCoastPoint; nCoastPoint++)
1347 nVCoastPoint.push_back(nCoastPoint);
1348 }
1349 else
1350 {
1351 // This is not the final down-coast polygon, so do not include the polygon's down-coast boundary
1352 nCoastSegLen = nDownCoastProfileCoastPoint - nUpCoastProfileCoastPoint;
1353
1354 for (int nCoastPoint = nUpCoastProfileCoastPoint; nCoastPoint < nDownCoastProfileCoastPoint; nCoastPoint++)
1355 nVCoastPoint.push_back(nCoastPoint);
1356 }
1357
1358 // Shuffle the coast points, this is necessary so that leaving the loop does not create sequence-related artefacts
1359 shuffle(nVCoastPoint.begin(), nVCoastPoint.begin() + nCoastSegLen, m_Rand[1]);
1360
1361 // Recalc the targets for deposition per profile
1362 dTargetToDepositOnProfile = dStillToDepositOnPoly / nCoastSegLen;
1363
1364 // Set a switch
1365 bool bEnoughDepositedOnPolygon = false;
1366
1367 // Now traverse the polygon's existing coastline in a random (but broadly UP-COAST, i.e. decreasing coastpoint indices) sequence, fitting a Dean profile at each coast point
1368 for (int n = 0; n < nCoastSegLen; n++)
1369 {
1370 // Check the switch
1371 if (bEnoughDepositedOnPolygon)
1372 break;
1373
1374 // Pick a random coast point
1375 int const nCoastPoint = nVCoastPoint[n];
1376 CGeom2DIPoint const PtiCoastPoint = *m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nCoastPoint);
1377 int const nCoastX = PtiCoastPoint.nGetX();
1378 int const nCoastY = PtiCoastPoint.nGetY();
1379
1380 // Calculate the x-y offset between this coast point, and the coast point of the down-coast normal
1381 int const nXOffset = nCoastX - nXDownCoastProfileExistingCoastPoint;
1382 int const nYOffset = nCoastY - nYDownCoastProfileExistingCoastPoint;
1383 int nSeawardOffset = -1;
1384 // nParProfLen;
1385 vector<CGeom2DIPoint> PtiVParProfile;
1386 vector<double> VdParProfileDeanElev;
1387
1388 // OK, loop until we can deposit sufficient unconsolidated sediment
1389 while (true)
1390 {
1391 // Move seaward by one cell
1392 nSeawardOffset++;
1393
1394 // And lengthen the parallel profile
1395 int nParProfLen = nDownCoastDeanLen + nSeawardOffset;
1396
1397 if (nParProfLen > (pDownCoastProfile->nGetNumCellsInProfile()))
1398 {
1399 // We've reached the seaward end of the down-coast profile, need to quit
1400 // if (m_nLogFileDetail >= LOG_FILE_MIDDLE_DETAIL)
1401 // LogStream << m_ulIter << ": " << WARN << "reached seaward end of down-coast profile during UP-COAST deposition of unconsolidated sediment for coast " << nCoast << " polygon " << nPoly << " (nCoastPoint = " << nCoastPoint << " nSeawardOffset = " << nSeawardOffset << ")" << endl;
1402
1403 break;
1404 }
1405
1406 // Get the x-y coords of a profile starting from this coast point and parallel to the down-coast polygon boundary profile (these are in natural sequence, like the boundary part-profile)
1407 PtiVParProfile.resize(0);
1408
1409 for (int m = 0; m < nParProfLen; m++)
1410 {
1411 CGeom2DIPoint const PtiProf = *pDownCoastProfile->pPtiGetCellInProfile(m);
1412 CGeom2DIPoint const PtiTmp(PtiProf.nGetX() + nXOffset, PtiProf.nGetY() + nYOffset);
1413 PtiVParProfile.push_back(PtiTmp);
1414 }
1415
1416 // Get the existing elevation of the seaward end of the parallel profile
1417 int nSeaEndX = PtiVParProfile.back().nGetX();
1418 int nSeaEndY = PtiVParProfile.back().nGetY();
1419
1420 // Safety check
1421 if (! bIsWithinValidGrid(nSeaEndX, nSeaEndY))
1422 {
1423 // LogStream << WARN << "12 @@@@ while doing UP-COAST deposition on polygon " << nPoly << ", hit edge of grid at [" << nSeaEndX << "][" << nSeaEndY << "] for parallel profile from coast point " << nCoastPoint << " at [" << nCoastX << "][" << nCoastY << "]. Constraining this parallel profile at its seaward end" << endl;
1424
1425 KeepWithinValidGrid(nCoastX, nCoastY, nSeaEndX, nSeaEndY);
1426 PtiVParProfile.back().SetX(nSeaEndX);
1427 PtiVParProfile.back().SetY(nSeaEndY);
1428 }
1429
1430 double const dParProfEndElev = m_pRasterGrid->m_Cell[nSeaEndX][nSeaEndY].dGetAllSedTopElevOmitTalus();
1431
1432 // Set the start elevation for the Dean profile just a bit above mean SWL for this timestep (i.e. so that it is a Bruun profile)
1433 double const dParProfStartElev = m_dThisIterMeanSWL + m_dDeanProfileStartAboveSWL;
1434
1435 // Calculate the total length of the parallel profile, including any seaward offset
1436 double const dParProfLen = dGetDistanceBetween(&PtiVParProfile.front(), &PtiVParProfile.back());
1437
1438 // Now calculate the length of the Dean profile-only part i.e. without any seaward offset. The approach used here is approximate but probably OK
1439 double dParProfDeanLen = dParProfLen - (nSeawardOffset * m_dCellSide);
1440
1441 // Safety check
1442 if (dParProfDeanLen <= 0)
1443 dParProfDeanLen = 1;
1444
1445 // Solve for dA so that the existing elevations at the end of the parallel profile, and at the end of a Dean equilibrium profile on that part-normal, are the same
1446 double const dParProfA = (dParProfStartElev - dParProfEndElev) / pow(dParProfDeanLen, DEAN_POWER);
1447
1448 nParProfLen = static_cast<int>(PtiVParProfile.size());
1449 VdParProfileDeanElev.resize(nParProfLen, 0);
1450
1451 double const dInc = dParProfDeanLen / (nParProfLen - nSeawardOffset - 2);
1452
1453 // The elevation of the coast point in the Dean profile is the same as the elevation of the current coast point TODO 020 Is this correct? Should it be dParProfStartElev?
1454 double const dCoastElev = m_pRasterGrid->m_Cell[nCoastX][nCoastY].dGetAllSedTopElevOmitTalus();
1455
1456 // For this depositing parallel profile, calculate the Dean equilibrium profile of the unconsolidated sediment h(y) = A * y^(2/3) where h(y) is the distance below the highest point in the profile at a distance y from the landward start of the profile
1457 CalcDeanProfile(&VdParProfileDeanElev, dInc, dParProfStartElev, dParProfA, true, nSeawardOffset, dCoastElev);
1458
1459 double dParProfTotDiff = 0;
1460
1461 for (int m = 0; m < nParProfLen; m++)
1462 {
1463 CGeom2DIPoint PtiTmp = PtiVParProfile[m];
1464 int nX = PtiTmp.nGetX();
1465 int nY = PtiTmp.nGetY();
1466
1467 // Safety check
1468 if (! bIsWithinValidGrid(nX, nY))
1469 {
1470 // LogStream << WARN << "13 @@@@ while constructing parallel profile for UP-COAST deposition on polygon " << nPoly << ", hit edge of grid at [" << nX << "][" << nY << "] for parallel profile from coast point " << nCoastPoint << " at [" << nCoastX << "][" << nCoastY << "]. Constraining this parallel profile at its seaward end" << endl;
1471
1472 KeepWithinValidGrid(nCoastX, nCoastY, nX, nY);
1473 PtiTmp.SetX(nX);
1474 PtiTmp.SetY(nY);
1475 }
1476
1477 // Don't do anything to intervention cells
1478 if (bIsInterventionCell(nX, nY))
1479 continue;
1480
1481 // Don't do cells twice
1482 if (! m_pRasterGrid->m_Cell[nX][nY].bBeachErosionOrDepositionThisIter())
1483 {
1484 double const dTmpElev = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus();
1485 double const dDiff = VdParProfileDeanElev[m] - dTmpElev;
1486
1487 dParProfTotDiff += dDiff;
1488 }
1489 }
1490
1491 // // DEBUG CODE -----------------------------------------------------
1492 // LogStream << endl << "\tFor polygon " << nPoly << " doing UP-COAST deposition, nSeawardOffset = " << nSeawardOffset << ", parallel profile from [" << PtiVParProfile[0].nGetX() << "][" << PtiVParProfile[0].nGetY() << "] to [" << PtiVParProfile.back().nGetX() << "][" << PtiVParProfile.back().nGetY() << "], nDownCoastDeanLen = " << nDownCoastDeanLen << " dDownCoastDeanLen = " << dDownCoastDeanLen << " nParProfLen = " << nParProfLen << " dParProfDeanLen = " << dParProfDeanLen << " dInc = " << dInc << " dParProfStartElev = " << dParProfStartElev << " dParProfEndElev = " << dParProfEndElev << " dParProfA = " << dParProfA << endl;
1493 //
1494 // LogStream << "\tExisting profile for deposition = ";
1495 // for (int n = 0; n < nParProfLen; n++)
1496 // {
1497 // CGeom2DIPoint PtiTmp = PtiVParProfile[n];
1498 // int
1499 // nX = PtiTmp.nGetX(),
1500 // nY = PtiTmp.nGetY();
1501 //
1502 // // Safety check
1503 // if (! bIsWithinValidGrid(nX, nY))
1504 // KeepWithinValidGrid(nX, nY);
1505 //
1506 // LogStream << m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus() << " ";
1507 // }
1508 // LogStream << endl;
1509 // LogStream << "\tParallel Dean equilibrium profile for deposition = ";
1510 // for (int n = 0; n < nParProfLen; n++)
1511 // {
1512 // LogStream << VdParProfileDeanElev[n] << " ";
1513 // }
1514 // LogStream << endl;
1515 //
1516 // LogStream << "\tnCoastPoint = " << nCoastPoint << " nSeawardOffset = " << nSeawardOffset << " difference = ";
1517 // for (int n = 0; n < nParProfLen; n++)
1518 // {
1519 // CGeom2DIPoint PtiTmp = PtiVParProfile[n];
1520 // int
1521 // nX = PtiTmp.nGetX(),
1522 // nY = PtiTmp.nGetY();
1523 //
1524 // // Safety check
1525 // if (! bIsWithinValidGrid(nX, nY))
1526 // KeepWithinValidGrid(nX, nY);
1527 //
1528 // double
1529 // dTmpElev = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus(),
1530 // dDiff = VdParProfileDeanElev[n] - dTmpElev;
1531 //
1532 // LogStream << dDiff << " ";
1533 // }
1534 // LogStream << endl;
1535 // // END DEBUG CODE -----------------------------------------------------
1536
1537 // So will we be able to deposit as much as is needed?
1538 if (dParProfTotDiff >= dTargetToDepositOnProfile)
1539 {
1540 // LogStream << " UP-COAST nCoastPoint = " << nCoastPoint << " nSeawardOffset = " << nSeawardOffset << " dParProfTotDiff = " << dParProfTotDiff << " dTargetToDepositOnProfile = " << dTargetToDepositOnProfile << " WILL BE ABLE TO DEPOSIT ENOUGH" << endl;
1541
1542 break;
1543 }
1544 }
1545
1546 // assert(dParProfTotDiff > 0);
1547
1548 // OK, this value of nSeawardOffset gives us enough deposition. So start depositing on the parallel profile from the coastward end
1549 double dDepositedOnProfile = 0; // Total for this parallel profile
1550
1551 dStillToDepositOnProfile = dTargetToDepositOnProfile;
1552
1553 for (unsigned int nSeawardFromCoast = 0; nSeawardFromCoast < PtiVParProfile.size(); nSeawardFromCoast++)
1554 {
1555 // Move along this parallel profile starting from the coast. Leave the loop if we have deposited enough on this polygon
1556 if (dStillToDepositOnPoly < SED_ELEV_TOLERANCE)
1557 {
1558 // LogStream << m_ulIter << ": nPoly = " << nPoly << " texture = " << strTexture << " UP-COAST nCoastPoint = " << nCoastPoint << " nSeawardOffset = " << nSeawardOffset << " leaving loop because enough deposited on polygon, dTargetToDepositOnPoly = " << dTargetToDepositOnPoly << " dDepositedOnPoly = " << dDepositedOnPoly << endl;
1559
1560 bEnoughDepositedOnPolygon = true;
1561
1562 break;
1563 }
1564
1565 // Leave the loop if we have deposited enough on this parallel profile
1566 if (dStillToDepositOnProfile < SED_ELEV_TOLERANCE)
1567 {
1568 // LogStream << m_ulIter << ": nPoly = " << nPoly << " texture = " << strTexture << " UP-COAST nCoastPoint = " << nCoastPoint << " nSeawardOffset = " << nSeawardOffset << " leaving loop because enough deposited on parallel profile, dTargetStillToDepositOnProfile = " << dTargetStillToDepositOnProfile << " dDepositedOnProfile = " << dDepositedOnProfile << endl;
1569
1570 break;
1571 }
1572
1573 // assert(nSeawardFromCoast < PtiVParProfile.size());
1574 CGeom2DIPoint PtiTmp = PtiVParProfile[nSeawardFromCoast];
1575 int nX = PtiTmp.nGetX();
1576 int nY = PtiTmp.nGetY();
1577
1578 // Safety check
1579 if (! bIsWithinValidGrid(nX, nY))
1580 {
1581 // LogStream << WARN << "14 @@@@ while constructing parallel profile for UP-COAST deposition on polygon " << nPoly << ", hit edge of grid at [" << nX << "][" << nY << "] for parallel profile from coast point " << nCoastPoint << " at [" << nCoastX << "][" << nCoastY << "]. Constraining this parallel profile at its seaward end" << endl;
1582
1583 KeepWithinValidGrid(nCoastX, nCoastY, nX, nY);
1584 PtiTmp.SetX(nX);
1585 PtiTmp.SetY(nY);
1586 }
1587
1588 // Don't do anything to intervention cells
1589 if (bIsInterventionCell(nX, nY))
1590 continue;
1591
1592 // Don't do cells twice
1593 if (! m_pRasterGrid->m_Cell[nX][nY].bBeachErosionOrDepositionThisIter())
1594 {
1595 // Get this cell's current elevation
1596 double const dThisElevNow = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus();
1597
1598 // LogStream << "\tnPoly = " << nPoly << " going UP-COAST, [" << nX << "][" << nY << "] nCoastPoint = " << nCoastPoint << " nSeawardFromCoast = " << nSeawardFromCoast << " dThisElevNow = " << dThisElevNow << " Dean Elev = " << VdParProfileDeanElev[nSeawardFromCoast] << endl;
1599
1600 // Subtract the two elevations
1601 double const dElevDiff = VdParProfileDeanElev[nSeawardFromCoast] - dThisElevNow;
1602 CRWCellLandform* pLandform = m_pRasterGrid->m_Cell[nX][nY].pGetLandform();
1603
1604 if (dElevDiff > SED_ELEV_TOLERANCE)
1605 {
1606 bool bDeposited = false;
1607 double dToDepositHere = 0;
1608
1609 // The current elevation is below the Dean elevation, so we have can have beach deposition here
1610 if (dStillToDepositOnProfile > SED_ELEV_TOLERANCE)
1611 {
1612 dToDepositHere = tMin(dElevDiff, dStillToDepositOnProfile, dStillToDepositOnPoly);
1613
1614 // LogStream << " UP-COAST nPoly = " << nPoly << " nCoastPoint = " << nCoastPoint << " texture = " << strTexture << " dToDepositHere = " << dToDepositHere << endl;
1615
1616 int const nTopLayer = m_pRasterGrid->m_Cell[nX][nY].nGetNumOfTopLayerAboveBasement();
1617
1618 // Safety check
1619 if (nTopLayer == INT_NODATA)
1620 return RTN_ERR_NO_TOP_LAYER;
1621
1622 if (dToDepositHere > SED_ELEV_TOLERANCE)
1623 {
1624 bDeposited = true;
1625
1626 if (nTexture == TEXTURE_SAND)
1627 {
1628 double const dSandNow = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetSandDepth();
1629
1630 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->SetSandDepth(dSandNow + dToDepositHere);
1631
1632 // Set the changed-this-timestep switch
1633 m_bUnconsChangedThisIter[nTopLayer] = true;
1634
1635 dDepositedOnProfile += dToDepositHere;
1636 dDepositedOnPoly += dToDepositHere;
1637
1638 // LogStream << "XXXXXXX texture = " << strTexture << " dDepositedOnProfile = " << dDepositedOnProfile << " dDepositedOnPoly = " << dDepositedOnPoly << endl;
1639
1640 // Update the cell's beach deposition, and total beach deposition, values
1641 m_pRasterGrid->m_Cell[nX][nY].IncrBeachDeposition(dToDepositHere);
1642
1643 dStillToDepositOnPoly -= dToDepositHere;
1644 dStillToDepositOnProfile -= dToDepositHere;
1645
1646 // if (dDepositedOnPoly > dTargetToDepositOnPoly)
1647 // {
1648 // // LogStream << "££££££ 5 nPoly = " << nPoly << " nCoastPoint = " << nCoastPoint << " dDepositedOnProfile = " << dDepositedOnProfile << " dDepositedOnPoly = " << dDepositedOnPoly << " dTargetToDepositOnPoly = " << dTargetToDepositOnPoly << endl;
1649 // }
1650
1651 // LogStream << m_ulIter << ": nPoly = " << nPoly << " texture = " << strTexture << " dToDepositHere = " << dToDepositHere << " at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "} nCoastPoint = " << nCoastPoint << " nSeawardFromCoast = " << nSeawardFromCoast << endl;
1652 }
1653 else if (nTexture == TEXTURE_COARSE)
1654 {
1655 double const dCoarseNow = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetCoarseDepth();
1656
1657 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->SetCoarseDepth(dCoarseNow + dToDepositHere);
1658
1659 // Set the changed-this-timestep switch
1660 m_bUnconsChangedThisIter[nTopLayer] = true;
1661
1662 if (dDepositedOnPoly > dTargetToDepositOnPoly)
1663 {
1664 // LogStream << "££££££ 6 nPoly = " << nPoly << " dDepositedOnPoly = " << dDepositedOnPoly << " dTargetToDepositOnPoly = " << dTargetToDepositOnPoly << endl;
1665 }
1666
1667 dDepositedOnProfile += dToDepositHere;
1668 dDepositedOnPoly += dToDepositHere;
1669
1670 // LogStream << "XXXXXXX texture = " << strTexture << " dDepositedOnProfile = " << dDepositedOnProfile << " dDepositedOnPoly = " << dDepositedOnPoly << endl;
1671
1672 // Update the cell's beach deposition, and total beach deposition, values
1673 m_pRasterGrid->m_Cell[nX][nY].IncrBeachDeposition(dToDepositHere);
1674
1675 dStillToDepositOnPoly -= dToDepositHere;
1676 dStillToDepositOnProfile -= dToDepositHere;
1677
1678 // LogStream << m_ulIter << ": nPoly = " << nPoly << " texture = " << strTexture << " dToDepositHere = " << dToDepositHere << " at [" << nX << "][" << nY << "] nCoastPoint = " << nCoastPoint << " nSeawardFromCoast = " << nSeawardFromCoast << endl;
1679 }
1680 }
1681 }
1682
1683 if (bDeposited)
1684 {
1685 // Now update the cell's layer elevations
1686 m_pRasterGrid->m_Cell[nX][nY].CalcAllLayerElevsAndD50();
1687
1688 // Update the cell's sea depth
1689 m_pRasterGrid->m_Cell[nX][nY].SetSeaDepth();
1690
1691 int const nCat = pLandform->nGetLFCategory();
1692
1693 if ((nCat != LF_SEDIMENT_INPUT_CONSOLIDATED) && (nCat != LF_SEDIMENT_INPUT_UNCONSOLIDATED))
1694 pLandform->SetLFCategory(LF_DRIFT_BEACH);
1695
1696 // Update this-timestep totals
1697 m_ulThisIterNumBeachDepositionCells++;
1698
1699 if (nTexture == TEXTURE_SAND)
1700 m_dThisIterBeachDepositionSand += dToDepositHere;
1701
1702 else if (nTexture == TEXTURE_COARSE)
1703 m_dThisIterBeachDepositionCoarse += dToDepositHere;
1704 }
1705 }
1706 else if (bElevAboveDeanElev(nX, nY, dElevDiff, pLandform))
1707 {
1708 // The current elevation is higher than the Dean elevation, so we could have beach erosion here
1709 m_ulThisIterNumPotentialBeachErosionCells++;
1710
1711 m_pRasterGrid->m_Cell[nX][nY].SetPotentialBeachErosion(-dElevDiff);
1712
1713 // LogStream << m_ulIter << ": nPoly = " << nPoly << " texture = " << strTexture << " going UP-COAST, potential beach erosion = " << -dElevDiff << " at [" << nX << "][" << nY << "] nCoastPoint = " << nCoastPoint << " nSeawardFromCoast = " << nSeawardFromCoast << " dThisElevNow = " << dThisElevNow << " Dean Elev = " << VdParProfileDeanElev[nSeawardFromCoast] << endl;
1714
1715 // Now get the number of the highest layer with non-zero thickness
1716 int const nThisLayer = m_pRasterGrid->m_Cell[nX][nY].nGetTopNonZeroLayerAboveBasement();
1717
1718 // Safety check
1719 if (nThisLayer == INT_NODATA)
1720 return RTN_ERR_NO_TOP_LAYER;
1721
1722 if (nThisLayer != NO_NONZERO_THICKNESS_LAYERS)
1723 {
1724 // We still have at least one layer left with non-zero thickness (i.e. we are not down to basement), and the cell's current elevation is higher than the Dean equilibrium profile elevation. So do some beach erosion on this cell (Note: we are in nDoUnconsDepositionOnPolygon() still)
1725 if (nTexture == TEXTURE_SAND)
1726 {
1727 double dSandRemoved = 0;
1728 ErodeCellBeachSedimentSupplyLimited(nX, nY, nThisLayer, TEXTURE_SAND, -dElevDiff, dSandRemoved);
1729
1730 // Update total for this parallel profile
1731 dDepositedOnProfile -= dSandRemoved;
1732 dStillToDepositOnProfile += dSandRemoved;
1733
1734 // Update totals for the polygon (note that these may become slightly -ve if erosion occurs early during this routine, but this is not a problem since the values will eventually become +ve again
1735 dDepositedOnPoly -= dSandRemoved;
1736 dStillToDepositOnPoly += dSandRemoved;
1737
1738 // if (dDepositedOnPoly < 0)
1739 // LogStream << m_ulIter << ": AAA LESS THAN ZERO dDepositedOnPoly = " << dDepositedOnPoly << endl;
1740
1741 // Update this-timestep totals
1742 m_ulThisIterNumActualBeachErosionCells++;
1743
1744 // // DEBUG CODE ==============================================================================================
1745 // LogStream << m_ulIter << ": $$$$$$$$$ BEACH 3 Dean profile lower than existing profile, SAND depth eroded on [" << nX << "][" << nY << "] in Poly " << nPoly << " is " << dSandRemoved * m_dCellArea << endl;
1746 // // DEBUG CODE ==============================================================================================
1747
1748 // Store the this-polygon depth of sand sediment eroded during Dean profile deposition of beach unconsolidated sediment
1749 pPolygon->AddBeachSandErodedDeanProfile(dSandRemoved);
1750 }
1751 else if (nTexture == TEXTURE_COARSE)
1752 {
1753 double dCoarseRemoved = 0;
1754 ErodeCellBeachSedimentSupplyLimited(nX, nY, nThisLayer, TEXTURE_COARSE, -dElevDiff, dCoarseRemoved);
1755
1756 // Update total for this parallel profile
1757 dDepositedOnProfile -= dCoarseRemoved;
1758 dStillToDepositOnProfile += dCoarseRemoved;
1759
1760 // Update totals for the polygon (note that these may become slightly -ve if erosion occurs early during this routine, but this is not a problem since the values will eventually become +ve again
1761 dDepositedOnPoly -= dCoarseRemoved;
1762 dStillToDepositOnPoly += dCoarseRemoved;
1763
1764 // Update this-timestep totals
1765 m_ulThisIterNumActualBeachErosionCells++;
1766
1767 // // DEBUG CODE ==============================================================================================
1768 // LogStream << m_ulIter << ": $$$$$$$$$ BEACH 3 Dean profile lower than existing profile, COARSE depth eroded on [" << nX << "][" << nY << "] in Poly " << nPoly << " is " << dCoarseRemoved * m_dCellArea << endl;
1769 // // DEBUG CODE ==============================================================================================
1770
1771 // Store the this-polygon depth of coarse sediment eroded during Dean profile deposition of beach unconsolidated sediment
1772 pPolygon->AddBeachCoarseErodedDeanProfile(dCoarseRemoved);
1773 }
1774 }
1775
1776 // if (dDepositedOnPoly > dTargetToDepositOnPoly)
1777 // {
1778 // // LogStream << "££££££ 4 nPoly = " << nPoly << " texture = " << strTexture << " dDepositedOnPoly = " << dDepositedOnPoly << " dTargetToDepositOnPoly = " << dTargetToDepositOnPoly << endl;
1779 // }
1780
1781 // LogStream << m_ulIter << ": in polygon " << nPoly << " have net deposition for " << strTexture << ", actual beach erosion = " << dSand + dCoarse << " at [" << nX << "][" << nY << "]" << endl;
1782 }
1783 }
1784 }
1785 }
1786 }
1787
1788 // Store amounts deposited on this polygon
1789 if (nTexture == TEXTURE_SAND)
1790 {
1791 pPolygon->SetBeachDepositionUnconsSand(dDepositedOnPoly);
1792
1793 // Check mass balance for sand deposited
1794 if (m_nLogFileDetail >= LOG_FILE_MIDDLE_DETAIL)
1795 if (! bFPIsEqual(pPolygon->dGetToDoBeachDepositionUnconsSand(), dDepositedOnPoly, MASS_BALANCE_TOLERANCE))
1796 LogStream << m_ulIter << ": \tcoast " << nCoast << " polygon " << pPolygon->nGetPolygonCoastID() << " NOT equal SAND dGetToDoBeachDepositionUnconsSand() = " << pPolygon->dGetToDoBeachDepositionUnconsSand() << " dDepositedOnPoly = " << dDepositedOnPoly << endl;
1797
1798 pPolygon->SetZeroToDoDepositionUnconsSand();
1799 }
1800 else if (nTexture == TEXTURE_COARSE)
1801 {
1802 pPolygon->SetBeachDepositionUnconsCoarse(dDepositedOnPoly);
1803
1804 // Check mass balance for coarse deposited
1805 if (m_nLogFileDetail >= LOG_FILE_MIDDLE_DETAIL)
1806 if (! bFPIsEqual(pPolygon->dGetToDoBeachDepositionUnconsCoarse(), dDepositedOnPoly, MASS_BALANCE_TOLERANCE))
1807 LogStream << m_ulIter << ": \tcoast " << nCoast << " polygon " << pPolygon->nGetPolygonCoastID() << " NOT equal COARSE dGetToDoBeachDepositionUnconsCoarse() = " << pPolygon->dGetToDoBeachDepositionUnconsCoarse() << " dDepositedOnPoly = " << dDepositedOnPoly << endl;
1808
1809 pPolygon->SetZeroToDoDepositionUnconsCoarse();
1810 }
1811
1812 // // DEBUG CODE #####################
1813 // if (m_ulIter == 5)
1814 // {
1815 // LogStream << m_ulIter << ": leaving nDoUnconsDepositionOnPolygon() nCoast = " << nCoast << " nPoly = " << nPoly << " strTexture = " << strTexture << " dTargetToDepositOnPoly = " << dTargetToDepositOnPoly * m_dCellArea << " dDepositedOnPoly = " << dDepositedOnPoly * m_dCellArea << endl;
1816 //
1817 // double dTmpSandUncons = 0;
1818 // for (int nX1 = 0; nX1 < m_nXGridSize; nX1++)
1819 // {
1820 // for (int nY1 = 0; nY1 < m_nYGridSize; nY1++)
1821 // {
1822 // dTmpSandUncons += m_pRasterGrid->m_Cell[nX1][nY1].dGetUnconsSandDepthAllLayers();
1823 // }
1824 // }
1825 //
1826 // LogStream << m_ulIter << ": leaving nDoUnconsDepositionOnPolygon() for nPoly = " << nPoly << " TOTAL UNCONSOLIDATED SAND ON ALL CELLS = " << dTmpSandUncons * m_dCellArea << endl;
1827 // }
1828 // // DEBUG CODE #####################
1829
1830 return RTN_OK;
1831}
1832
1833//===============================================================================================================================
1835//===============================================================================================================================
1836bool CSimulation::bElevAboveDeanElev(int const nX, int const nY, double const dElevDiff, CRWCellLandform const* pLandform)
1837{
1838 // TODO 075 What if it is bedrock that sticks above Dean profile?
1839 if (dElevDiff <= 0)
1840 {
1841 if (m_pRasterGrid->m_Cell[nX][nY].bIsInContiguousSea())
1842 return true;
1843
1844 int const nCat = pLandform->nGetLFCategory();
1845
1846 if ((nCat == LF_DRIFT_BEACH) || (nCat == LF_DRIFT_DUNES) || (nCat == LF_SEDIMENT_INPUT_UNCONSOLIDATED))
1847 return true;
1848 }
1849
1850 return false;
1851}
2d_point.h
Contains CGeom2DPoint definitions.
2di_point.h
Contains CGeom2DIPoint definitions.
CGeom2DIPoint
Geometry class used to represent 2D point objects with integer coordinates.
Definition 2di_point.h:25
CGeom2DIPoint::SetY
void SetY(int const)
The integer parameter sets a value for the CGeom2DIPoint object's Y coordinate.
Definition 2di_point.cpp:75
CGeom2DIPoint::nGetY
int nGetY(void) const
Returns the CGeom2DIPoint object's integer Y coordinate.
Definition 2di_point.cpp:51
CGeom2DIPoint::SetX
void SetX(int const)
The integer parameter sets a value for the CGeom2DIPoint object's X coordinate.
Definition 2di_point.cpp:69
CGeom2DIPoint::nGetX
int nGetX(void) const
Returns the CGeom2DIPoint object's integer X coordinate.
Definition 2di_point.cpp:45
CGeom2DPoint
Geometry class used to represent 2D point objects with floating-point coordinates.
Definition 2d_point.h:25
CGeomCoastPolygon
Geometry class used for coast polygon objects.
Definition coast_polygon.h:30
CGeomCoastPolygon::AddBeachCoarseErodedDeanProfile
void AddBeachCoarseErodedDeanProfile(double const)
Adds a depth (in m) of coarse sediment eroded during beach Dean profile deposition.
Definition coast_polygon.cpp:394
CGeomCoastPolygon::SetBeachDepositionUnconsCoarse
void SetBeachDepositionUnconsCoarse(double const)
Sets the depth of coarse-sized unconsolidated sediment deposited on this polygon during this timestep...
Definition coast_polygon.cpp:334
CGeomCoastPolygon::AddBeachSandErodedDeanProfile
void AddBeachSandErodedDeanProfile(double const)
Adds a depth (in m) of sand sediment eroded during beach Dean profile deposition.
Definition coast_polygon.cpp:382
CGeomCoastPolygon::SetZeroToDoDepositionUnconsSand
void SetZeroToDoDepositionUnconsSand(void)
Re-initialises this timestep's still-to-do deposition of unconsolidated sand sediment (from beach red...
Definition coast_polygon.cpp:346
CGeomCoastPolygon::SetBeachDepositionUnconsSand
void SetBeachDepositionUnconsSand(double const)
Sets the depth of sand-sized unconsolidated sediment deposited on this polygon during this timestep (...
Definition coast_polygon.cpp:322
CGeomCoastPolygon::SetZeroToDoDepositionUnconsCoarse
void SetZeroToDoDepositionUnconsCoarse(void)
Re-initialises this timestep's still-to-do deposition of unconsolidated coarse sediment (from beach r...
Definition coast_polygon.cpp:364
CGeomCoastPolygon::nGetPolygonCoastID
int nGetPolygonCoastID(void) const
Get the coast ID, this is the same as the down-coast sequence of polygons.
Definition coast_polygon.cpp:138
CGeomCoastPolygon::nGetUpCoastProfile
int nGetUpCoastProfile(void) const
Return the number of the up-coast profile.
Definition coast_polygon.cpp:191
CGeomCoastPolygon::dGetToDoBeachDepositionUnconsSand
double dGetToDoBeachDepositionUnconsSand(void) const
Returns this timestep's still-to-do deposition of sand-sized unconsolidated sediment (from beach redi...
Definition coast_polygon.cpp:358
CGeomCoastPolygon::nGetDownCoastProfile
int nGetDownCoastProfile(void) const
Return the number of the down-coast profile.
Definition coast_polygon.cpp:197
CGeomCoastPolygon::dGetToDoBeachDepositionUnconsCoarse
double dGetToDoBeachDepositionUnconsCoarse(void) const
Returns this timestep's still-to-do deposition of coarse unconsolidated sediment (from beach redistri...
Definition coast_polygon.cpp:376
CGeomProfile
Geometry class used to represent coast profile objects.
Definition profile.h:33
CGeomProfile::nGetCellGivenDepth
int nGetCellGivenDepth(CGeomRasterGrid const *, double const)
Returns the index of the cell on this profile which has a sea depth which is just less than a given d...
Definition profile.cpp:545
CGeomProfile::nGetCoastPoint
int nGetCoastPoint(void) const
Returns the coast point at which the profile starts.
Definition profile.cpp:80
CGeomProfile::pPtiGetCellInProfile
CGeom2DIPoint * pPtiGetCellInProfile(int const)
Returns a single cell (grid CRS) in the profile.
Definition profile.cpp:518
CGeomProfile::nGetNumCellsInProfile
int nGetNumCellsInProfile(void) const
Returns the number of cells in the profile.
Definition profile.cpp:539
CRWCellLandform
Real-world class used to represent the landform of a cell.
Definition cell_landform.h:25
CRWCellLandform::nGetLFCategory
int nGetLFCategory(void) const
Get the landform category.
Definition cell_landform.cpp:49
CRWCellLandform::SetLFCategory
void SetLFCategory(int const)
Set the landform category.
Definition cell_landform.cpp:43
CSimulation::m_nLogFileDetail
int m_nLogFileDetail
The level of detail in the log file output. Can be LOG_FILE_LOW_DETAIL, LOG_FILE_MIDDLE_DETAIL,...
Definition simulation.h:588
CSimulation::m_pRasterGrid
CGeomRasterGrid * m_pRasterGrid
Pointer to the raster grid object.
Definition simulation.h:1558
CSimulation::LogStream
ofstream LogStream
Definition simulation.h:1893
CSimulation::m_VCoast
vector< CRWCoast > m_VCoast
The coastline objects.
Definition simulation.h:1561
CSimulation::m_dFineErodibilityNormalized
double m_dFineErodibilityNormalized
Relative erodibility of fine unconsolidated beach sediment, normalized.
Definition simulation.h:810
CSimulation::ErodeCellBeachSedimentSupplyLimited
void ErodeCellBeachSedimentSupplyLimited(int const, int const, int const, int const, double const, double &)
Erodes the unconsolidated beach sediment on a single cell, for a single size class,...
Definition do_beach_within_polygon.cpp:679
CSimulation::dGetDistanceBetween
static double dGetDistanceBetween(CGeom2DPoint const *, CGeom2DPoint const *)
Returns the distance (in external CRS) between two points.
Definition gis_utils.cpp:143
CSimulation::m_dCoarseErodibilityNormalized
double m_dCoarseErodibilityNormalized
Relative erodibility of coarse unconsolidated beach sediment, normalized.
Definition simulation.h:816
CSimulation::KeepWithinValidGrid
void KeepWithinValidGrid(int &, int &) const
Constrains the supplied point (in the grid CRS) to be a valid cell within the raster grid.
Definition gis_utils.cpp:212
CSimulation::dSubtractProfiles
static double dSubtractProfiles(vector< double > const *, vector< double > const *, vector< bool > const *)
Calculate the total elevation difference between every point in two elevation profiles (first profile...
Definition utils.cpp:2717
CSimulation::m_ulThisIterNumBeachDepositionCells
unsigned long m_ulThisIterNumBeachDepositionCells
The number of grid cells on which beach (unconsolidated sediment) deposition occurs,...
Definition simulation.h:639
CSimulation::m_dDeanProfileStartAboveSWL
double m_dDeanProfileStartAboveSWL
Berm height i.e. height above SWL of start of depositional Dean profile.
Definition simulation.h:975
CSimulation::nDoUnconsDepositionOnPolygon
int nDoUnconsDepositionOnPolygon(int const, CGeomCoastPolygon *, int const, double, double &)
Deposits unconsolidated beach sediment (sand or coarse) on the cells within a polygon....
Definition do_beach_within_polygon.cpp:749
CSimulation::m_dSandErodibilityNormalized
double m_dSandErodibilityNormalized
Relative erodibility of sand unconsolidated beach sediment, normalized.
Definition simulation.h:813
CSimulation::m_dThisIterBeachDepositionCoarse
double m_dThisIterBeachDepositionCoarse
Total beach deposition (coarse unconsolidated sediment) for this iteration (depth in m)
Definition simulation.h:876
CSimulation::bElevAboveDeanElev
bool bElevAboveDeanElev(int const, int const, double const, CRWCellLandform const *)
Return true if the given elevation is higher than the Dean elevation (and other land use conditions a...
Definition do_beach_within_polygon.cpp:1836
CSimulation::m_ulThisIterNumActualBeachErosionCells
unsigned long m_ulThisIterNumActualBeachErosionCells
The number of grid cells on which actual beach (unconsolidated sediment) erosion occurs,...
Definition simulation.h:636
CSimulation::m_bUnconsChangedThisIter
vector< bool > m_bUnconsChangedThisIter
One element per layer: has the consolidated sediment of this layer been changed during this iteration...
Definition simulation.h:1387
CSimulation::m_Rand
default_random_engine m_Rand[NUMBER_OF_RNGS]
The c++11 random number generators.
Definition simulation.h:1591
CSimulation::nDoUnconsErosionOnPolygon
int nDoUnconsErosionOnPolygon(int const, CGeomCoastPolygon *, int const, double const, double &)
Erodes unconsolidated beach sediment of one texture class on the cells within a polygon....
Definition do_beach_within_polygon.cpp:39
CSimulation::m_ulThisIterNumPotentialBeachErosionCells
unsigned long m_ulThisIterNumPotentialBeachErosionCells
The number of grid cells on which potential beach (unconsolidated sediment) erosion occurs,...
Definition simulation.h:633
CSimulation::bIsWithinValidGrid
bool bIsWithinValidGrid(int const, int const) const
Checks whether the supplied point (an x-y pair, in the grid CRS) is within the raster grid,...
Definition gis_utils.cpp:184
CSimulation::nDoParallelProfileUnconsErosion
int nDoParallelProfileUnconsErosion(CGeomCoastPolygon *, int const, int const, int const, int const, int const, int const, vector< CGeom2DIPoint > const *, vector< double > const *, double &, double &, double &)
This routine erodes unconsolidated beach sediment (either fine, sand, or coarse) on a parallel profil...
Definition do_beach_within_polygon.cpp:486
CSimulation::m_dThisIterBeachDepositionSand
double m_dThisIterBeachDepositionSand
Total beach deposition (sand unconsolidated sediment) for this iteration (depth in m)
Definition simulation.h:873
CSimulation::PtGridCentroidToExt
CGeom2DPoint PtGridCentroidToExt(CGeom2DIPoint const *) const
Transforms a pointer to a CGeom2DIPoint in the raster grid CRS (assumed to be the centroid of a cell)...
Definition gis_utils.cpp:83
CSimulation::m_ulIter
unsigned long m_ulIter
The number of the current iteration (time step)
Definition simulation.h:609
CSimulation::m_dCellSide
double m_dCellSide
Length of a cell side (in external CRS units)
Definition simulation.h:672
CSimulation::bIsInterventionCell
bool bIsInterventionCell(int const, int const) const
Returns true if the cell is an intervention.
Definition utils.cpp:3066
CSimulation::m_dDepthOfClosure
double m_dDepthOfClosure
Depth of closure (in m) TODO 007 can be calculated using Hallermeier, R.J. (1978) or Birkemeier (1985...
Definition simulation.h:831
CSimulation::CalcDeanProfile
static void CalcDeanProfile(vector< double > *, double const, double const, double const, bool const, int const, double const)
Calculates a Dean equilibrium profile h(y) = A * y^(2/3) where h(y) is the distance below the highest...
Definition utils.cpp:2678
CSimulation::m_dThisIterMeanSWL
double m_dThisIterMeanSWL
The mean still water level for this timestep (does not include tidal changes, but includes any long-t...
Definition simulation.h:729
cme.h
This file contains global definitions for CoastalME.
NO_NONZERO_THICKNESS_LAYERS
int const NO_NONZERO_THICKNESS_LAYERS
Definition cme.h:672
INT_NODATA
int const INT_NODATA
Definition cme.h:380
tMin
T tMin(T a, T b)
Definition cme.h:1175
TOLERANCE
double const TOLERANCE
Definition cme.h:725
RTN_ERR_NO_TOP_LAYER
int const RTN_ERR_NO_TOP_LAYER
Definition cme.h:629
ERR
string const ERR
Definition cme.h:805
LF_SEDIMENT_INPUT_CONSOLIDATED
int const LF_SEDIMENT_INPUT_CONSOLIDATED
Definition cme.h:446
TEXTURE_COARSE
int const TEXTURE_COARSE
Definition cme.h:419
LOG_FILE_MIDDLE_DETAIL
int const LOG_FILE_MIDDLE_DETAIL
Definition cme.h:393
MIN_INLAND_OFFSET_UNCONS_EROSION
int const MIN_INLAND_OFFSET_UNCONS_EROSION
Definition cme.h:384
bFPIsEqual
bool bFPIsEqual(const T d1, const T d2, const T dEpsilon)
Definition cme.h:1213
RTN_ERR_NO_SEAWARD_END_OF_PROFILE_DOWNCOAST_BEACH_DEPOSITION
int const RTN_ERR_NO_SEAWARD_END_OF_PROFILE_DOWNCOAST_BEACH_DEPOSITION
Definition cme.h:627
MASS_BALANCE_TOLERANCE
double const MASS_BALANCE_TOLERANCE
Definition cme.h:727
LOG_FILE_HIGH_DETAIL
int const LOG_FILE_HIGH_DETAIL
Definition cme.h:394
TEXTURE_FINE
int const TEXTURE_FINE
Definition cme.h:417
LF_SEDIMENT_INPUT_UNCONSOLIDATED
int const LF_SEDIMENT_INPUT_UNCONSOLIDATED
Definition cme.h:445
LF_DRIFT_DUNES
int const LF_DRIFT_DUNES
Definition cme.h:441
LOG_FILE_ALL
int const LOG_FILE_ALL
Definition cme.h:395
RTN_OK
int const RTN_OK
Definition cme.h:585
MIN_PARALLEL_PROFILE_SIZE
int const MIN_PARALLEL_PROFILE_SIZE
Definition cme.h:385
DEAN_POWER
double const DEAN_POWER
Definition cme.h:719
LF_DRIFT_BEACH
int const LF_DRIFT_BEACH
Definition cme.h:440
SED_ELEV_TOLERANCE
double const SED_ELEV_TOLERANCE
Definition cme.h:726
RTN_ERR_NO_SEAWARD_END_OF_PROFILE_BEACH_EROSION
int const RTN_ERR_NO_SEAWARD_END_OF_PROFILE_BEACH_EROSION
Definition cme.h:625
TEXTURE_SAND
int const TEXTURE_SAND
Definition cme.h:418
RTN_ERR_NO_SEAWARD_END_OF_PROFILE_UPCOAST_BEACH_DEPOSITION
int const RTN_ERR_NO_SEAWARD_END_OF_PROFILE_UPCOAST_BEACH_DEPOSITION
Definition cme.h:626
coast.h
Contains CRWCoast definitions.
simulation.h
Contains CSimulation definitions.