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