CoastalME (Coastal Modelling Environment)
Simulates the long-term behaviour of complex coastlines
Loading...
Searching...
No Matches
calc_waves.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 <cstdio>
27
28#include <cmath>
29using std::sin;
30using std::cos;
31using std::pow;
32using std::fmod;
33
34#include <string>
35using std::to_string;
36
37#include <iostream>
38using std::endl;
39using std::ios;
40
41#include <fstream>
42using std::ifstream;
43
44#include <algorithm>
45using std::reverse_copy;
46
47#include "cme.h"
48#include "coast.h"
49#include "simulation.h"
50#include "cshore/cshore.h"
51#include "2d_point.h"
52
53//===============================================================================================================================
55//===============================================================================================================================
57{
58 // For each coastline, put a value for deep water wave height and direction at each coastline point
59 for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
60 {
61 int nDistFromPrevProfile = 0;
62 int nDistToNextProfile = 0;
63
64 double dPrevProfileDeepWaterWaveHeight = 0;
65 double dPrevProfileDeepWaterWaveAngle = 0;
66 double dPrevProfileDeepWaterWavePeriod = 0;
67 double dNextProfileDeepWaterWaveHeight = 0;
68 double dNextProfileDeepWaterWaveAngle = 0;
69 double dNextProfileDeepWaterWavePeriod = 0;
70 double dDist = 0;
71
72 for (int nPoint = 0; nPoint < m_VCoast[nCoast].nGetCoastlineSize(); nPoint++)
73 {
74 // We are going down-coast (i.e. along the coast in the direction of increasing coastline point numbers)
75 if (m_VCoast[nCoast].bIsProfileAtCoastPoint(nPoint))
76 {
77 // OK, a coastline-normal profile begins at this coastline point, so set the deep water wave values at this coastline point to be the values at the seaward end of the coastline normal
78 CGeomProfile const* pProfile = m_VCoast[nCoast].pGetProfileAtCoastPoint(nPoint);
79 // int nProfile = pProfile->nGetProfileID();
80
81 double const dThisDeepWaterWaveHeight = pProfile->dGetProfileDeepWaterWaveHeight();
82 double const dThisDeepWaterWaveAngle = pProfile->dGetProfileDeepWaterWaveAngle();
83 double const dThisDeepWaterWavePeriod = pProfile->dGetProfileDeepWaterWavePeriod();
84
85 m_VCoast[nCoast].SetCoastDeepWaterWaveHeight(nPoint, dThisDeepWaterWaveHeight);
86 m_VCoast[nCoast].SetCoastDeepWaterWaveAngle(nPoint, dThisDeepWaterWaveAngle);
87 m_VCoast[nCoast].SetCoastDeepWaterWavePeriod(nPoint, dThisDeepWaterWavePeriod);
88
89 // Reset for next time
90 nDistFromPrevProfile = 0;
91 dPrevProfileDeepWaterWaveHeight = dThisDeepWaterWaveHeight;
92 dPrevProfileDeepWaterWaveAngle = dThisDeepWaterWaveAngle;
93 dPrevProfileDeepWaterWavePeriod = dThisDeepWaterWavePeriod;
94
95 // Find the next profile
96 CGeomProfile const* pNextProfile = pProfile->pGetDownCoastAdjacentProfile();
97
98 if (pNextProfile == NULL)
99 {
100 // We are at the end of the coast
101 break;
102 }
103
104 // And the distance (in along-coast points) to the next profile
105 nDistToNextProfile = pNextProfile->nGetCoastPoint() - nPoint;
106 dDist = nDistToNextProfile;
107
108 // And the next profile's deep water wave values
109 dNextProfileDeepWaterWaveHeight = pNextProfile->dGetProfileDeepWaterWaveHeight();
110 dNextProfileDeepWaterWaveAngle = pNextProfile->dGetProfileDeepWaterWaveAngle();
111 dNextProfileDeepWaterWavePeriod = pNextProfile->dGetProfileDeepWaterWavePeriod();
112
113 // LogStream << m_ulIter << ": coast point = " << nPoint << " is start of profile " << pProfile->nGetProfileID() << ", next profile is " << pNextProfile->nGetProfileID() << ", which starts at coast piint " << pNextProfile->nGetCoastPoint() << ", dThisDeepWaterWaveHeight = " << dThisDeepWaterWaveHeight << ", dThisDeepWaterWaveAngle = " << dThisDeepWaterWaveAngle << " nDistToNextProfile = " << nDistToNextProfile << endl;
114 }
115
116 else
117 {
118 // This coast point is not the start of a coastline normal, so set the deep water wave values at this coastline point to be a weighted average of those from the up-coast and down-coast profiles
119 nDistFromPrevProfile++;
120 nDistToNextProfile--;
121
122 double const dPrevWeight = (dDist - nDistFromPrevProfile) / dDist;
123 double const dNextWeight = (dDist - nDistToNextProfile) / dDist;
124 double const dThisDeepWaterWaveHeight = (dPrevWeight * dPrevProfileDeepWaterWaveHeight) + (dNextWeight * dNextProfileDeepWaterWaveHeight);
125 double const dThisDeepWaterWaveAngle = dKeepWithin360((dPrevWeight * dPrevProfileDeepWaterWaveAngle) + (dNextWeight * dNextProfileDeepWaterWaveAngle));
126 double const dThisDeepWaterWavePeriod = (dPrevWeight * dPrevProfileDeepWaterWavePeriod) + (dNextWeight * dNextProfileDeepWaterWavePeriod);
127
128 m_VCoast[nCoast].SetCoastDeepWaterWaveHeight(nPoint, dThisDeepWaterWaveHeight);
129 m_VCoast[nCoast].SetCoastDeepWaterWaveAngle(nPoint, dThisDeepWaterWaveAngle);
130 m_VCoast[nCoast].SetCoastDeepWaterWavePeriod(nPoint, dThisDeepWaterWavePeriod);
131
132 // LogStream << m_ulIter << ": coast point = " << nPoint << " dThisDeepWaterWaveHeight = " << dThisDeepWaterWaveHeight << " dThisDeepWaterWaveAngle = " << dThisDeepWaterWaveAngle << " nDistFromPrevProfile = " << nDistFromPrevProfile << " nDistToNextProfile = " << nDistToNextProfile << endl;
133 }
134 }
135 }
136
137 return RTN_OK;
138}
139
140//===============================================================================================================================
142//===============================================================================================================================
144{
145 // Set up all-profile vectors to hold the wave attribute data at every profile point on all profiles
146 vector<bool> VbBreakingAll;
147
148 vector<double> VdXAll;
149 vector<double> VdYAll;
150 vector<double> VdHeightXAll;
151 vector<double> VdHeightYAll;
152
153 // Calculate wave properties for every coast
154 bool bSomeNonStartOrEndOfCoastProfiles = false;
155
156 for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
157 {
158 int const nCoastSize = m_VCoast[nCoast].nGetCoastlineSize();
159 int const nNumProfiles = m_VCoast[nCoast].nGetNumProfiles();
160
161 static bool bDownCoast = true;
162
163 // Calculate wave properties at every point along each valid profile, and for the cells under the profiles. Do this alternately in up-coast and down-coast sequence
164 for (int nn = 0; nn < nNumProfiles; nn++)
165 {
166 vector<bool> VbBreaking;
167 vector<double> VdX;
168 vector<double> VdY;
169 vector<double> VdHeightX;
170 vector<double> VdHeightY;
171
172 CGeomProfile *pProfile;
173
174 if (bDownCoast)
175 pProfile = m_VCoast[nCoast].pGetProfileWithDownCoastSeq(nn);
176
177 else
178 pProfile = m_VCoast[nCoast].pGetProfileWithUpCoastSeq(nn);
179
180 int const nRet = nCalcWavePropertiesOnProfile(nCoast, nCoastSize, pProfile, &VdX, &VdY, &VdHeightX, &VdHeightY, &VbBreaking);
181
182 if (nRet != RTN_OK)
183 {
184 if (nRet == RTN_ERR_CSHORE_ERROR)
185 {
186 // Abandon calculations on this profile, and flag the profile
187 pProfile->SetCShoreProblem(true);
188
189 // Move on to next profile
190 continue;
191 }
192
193 else
194 {
195 // A serious CShore error, so abort the run
196 return nRet;
197 }
198 }
199
200 // Are the waves off-shore? If so, do nothing more with this profile. The wave values for cells have already been given the off-shore value
201 if (VbBreaking.empty())
202 continue;
203
204 // Is this a start of coast or end of coast profile?
205 if ((! pProfile->bStartOfCoast()) && (! pProfile->bEndOfCoast()))
206 {
207 // It is neither a start of coast or an end of coast profile, so set switch
208 bSomeNonStartOrEndOfCoastProfiles = true;
209 }
210
211 // // DEBUG CODE ===============================================================================================
212 // for (int nn = 0; nn < VdX.size(); nn++)
213 // {
214 // LogStream << "nProfile = " << nProfile << " nn = " << nn << " VdX[nn] = " << VdX[nn] << " VdY[nn] = " << VdY[nn] << " VdHeightX[nn] = " << VdHeightX[nn] << " VdHeightY[nn] = " << VdHeightY[nn] << " VbBreaking[nn] = " << VbBreaking[nn] << endl;
215 // }
216 // LogStream << endl;
217 // // DEBUG CODE ===============================================================================================
218
219 // Append to the all-profile vectors
220 VdXAll.insert(VdXAll.end(), VdX.begin(), VdX.end());
221 VdYAll.insert(VdYAll.end(), VdY.begin(), VdY.end());
222 VdHeightXAll.insert(VdHeightXAll.end(), VdHeightX.begin(), VdHeightX.end());
223 VdHeightYAll.insert(VdHeightYAll.end(), VdHeightY.begin(), VdHeightY.end());
224 VbBreakingAll.insert(VbBreakingAll.end(), VbBreaking.begin(), VbBreaking.end());
225 }
226
227 bDownCoast = ! bDownCoast;
228 }
229
230 // OK, do we have some profiles other than start of coast or end of coast profiles in the all-profile vectors? We need to check this, because GDALGridCreate() in nInterpolateWavePropertiesToWithinPolygonCells() does not work if we give it only a start-of-coast or an end-of-coast profile to work with TODO 006 Is this still true?
231 if (! bSomeNonStartOrEndOfCoastProfiles)
232 {
233 LogStream << m_ulIter << ": waves are on-shore only, for start and/or end of coast profiles" << endl;
234
235 return RTN_OK;
236 }
237
238 // We need to also send the deepwater points from the edge of the grid to nInterpolateWavePropertiesToWithinPolygonCells(), this is necessary to prevent GDALGridCreate() leaving holes in the interpolated grid when the polygons are far from regular
239 double dDeepWaterWaveX;
240 double dDeepWaterWaveY;
241
243 {
244 // Just using the same value of deep water wave hehght and angle for all cells
245 dDeepWaterWaveX = m_dAllCellsDeepWaterWaveHeight * sin(m_dAllCellsDeepWaterWaveAngle * PI / 180),
246 dDeepWaterWaveY = m_dAllCellsDeepWaterWaveHeight * cos(m_dAllCellsDeepWaterWaveAngle * PI / 180);
247 }
248
249 for (int nX = 0; nX < m_nXGridSize; nX++)
250 {
251 if (m_pRasterGrid->m_Cell[nX][0].bIsInContiguousSea())
252 {
253 int const nPolyID = m_pRasterGrid->m_Cell[nX][0].nGetPolygonID();
254
255 if (nPolyID == INT_NODATA)
256 {
257 // Not in a polygon
258 VdXAll.push_back(nX);
259 VdYAll.push_back(0);
260
262 {
263 // Not using the same value of deep water height and angle for all cells, so get this cell's deep water height and angle values
264 dDeepWaterWaveX = m_pRasterGrid->m_Cell[nX][0].dGetCellDeepWaterWaveHeight() * sin(m_pRasterGrid->m_Cell[nX][0].dGetCellDeepWaterWaveAngle() * PI / 180);
265 dDeepWaterWaveY = m_pRasterGrid->m_Cell[nX][0].dGetCellDeepWaterWaveHeight() * cos(m_pRasterGrid->m_Cell[nX][0].dGetCellDeepWaterWaveAngle() * PI / 180);
266 }
267
268 VdHeightXAll.push_back(dDeepWaterWaveX);
269 VdHeightYAll.push_back(dDeepWaterWaveY);
270 }
271 }
272
273 if (m_pRasterGrid->m_Cell[nX][m_nYGridSize - 1].bIsInContiguousSea())
274 {
275 int const nPolyID = m_pRasterGrid->m_Cell[nX][m_nYGridSize - 1].nGetPolygonID();
276
277 if (nPolyID == INT_NODATA)
278 {
279 // Not in a polygon
280 VdXAll.push_back(nX);
281 VdYAll.push_back(m_nYGridSize - 1);
282
284 {
285 // Not using the same value of deep water height and angle for all cells, so get this cell's deep water height and angle values
286 dDeepWaterWaveX = m_pRasterGrid->m_Cell[nX][m_nYGridSize - 1].dGetCellDeepWaterWaveHeight() * sin(m_pRasterGrid->m_Cell[nX][m_nYGridSize - 1].dGetCellDeepWaterWaveAngle() * PI / 180);
287 dDeepWaterWaveY = m_pRasterGrid->m_Cell[nX][m_nYGridSize - 1].dGetCellDeepWaterWaveHeight() * cos(m_pRasterGrid->m_Cell[nX][m_nYGridSize - 1].dGetCellDeepWaterWaveAngle() * PI / 180);
288 }
289
290 VdHeightXAll.push_back(dDeepWaterWaveX);
291 VdHeightYAll.push_back(dDeepWaterWaveY);
292 }
293 }
294 }
295
296 for (int nY = 0; nY < m_nYGridSize; nY++)
297 {
298 if (m_pRasterGrid->m_Cell[0][nY].bIsInContiguousSea())
299 {
300 int const nPolyID = m_pRasterGrid->m_Cell[0][nY].nGetPolygonID();
301
302 if (nPolyID == INT_NODATA)
303 {
304 // Not in a polygon
305 VdXAll.push_back(0);
306 VdYAll.push_back(nY);
307
309 {
310 // Not using the same value of deep water height and angle for all cells, so get this cell's deep water height and angle values
311 dDeepWaterWaveX = m_pRasterGrid->m_Cell[0][nY].dGetCellDeepWaterWaveHeight() * sin(m_pRasterGrid->m_Cell[0][nY].dGetCellDeepWaterWaveAngle() * PI / 180);
312 dDeepWaterWaveY = m_pRasterGrid->m_Cell[0][nY].dGetCellDeepWaterWaveHeight() * cos(m_pRasterGrid->m_Cell[0][nY].dGetCellDeepWaterWaveAngle() * PI / 180);
313 }
314
315 VdHeightXAll.push_back(dDeepWaterWaveX);
316 VdHeightYAll.push_back(dDeepWaterWaveY);
317 }
318 }
319
320 if (m_pRasterGrid->m_Cell[m_nXGridSize - 1][nY].bIsInContiguousSea())
321 {
322 int const nPolyID = m_pRasterGrid->m_Cell[m_nXGridSize - 1][nY].nGetPolygonID();
323
324 if (nPolyID == INT_NODATA)
325 {
326 // Not in a polygon
327 VdXAll.push_back(m_nXGridSize - 1);
328 VdYAll.push_back(nY);
329
331 {
332 // Not using the same value of deep water height and angle for all cells, so get this cell's deep water height and angle values
333 dDeepWaterWaveX = m_pRasterGrid->m_Cell[m_nXGridSize - 1][nY].dGetCellDeepWaterWaveHeight() * sin(m_pRasterGrid->m_Cell[m_nXGridSize - 1][nY].dGetCellDeepWaterWaveAngle() * PI / 180);
334 dDeepWaterWaveY = m_pRasterGrid->m_Cell[m_nXGridSize - 1][nY].dGetCellDeepWaterWaveHeight() * cos(m_pRasterGrid->m_Cell[m_nXGridSize - 1][nY].dGetCellDeepWaterWaveAngle() * PI / 180);
335 }
336
337 VdHeightXAll.push_back(dDeepWaterWaveX);
338 VdHeightYAll.push_back(dDeepWaterWaveY);
339 }
340 }
341 }
342
343 // // DEBUG CODE ============================================================================================================
344 // LogStream << "Out of loop" << endl;
345 // for (int nn = 0; nn < VdXAll.size(); nn++)
346 // {
347 // LogStream << "nn = " << nn << " VdXAll[nn] = " << VdXAll[nn] << " VdYAll[nn] = " << VdYAll[nn] << " VdHeightXAll[nn] = " << VdHeightXAll[nn] << " VdHeightYAll[nn] = " << VdHeightYAll[nn] << " VbBreakingAll[nn] = " << VbBreakingAll[nn] << endl;
348 // }
349 // LogStream << endl;
350 // // DEBUG CODE ============================================================================================================
351
352 // Are the waves off-shore for every profile? If so, do nothing more
353 if (VbBreakingAll.empty())
354 {
355 LogStream << m_ulIter << ": waves off-shore for all profiles" << endl;
356
357 return RTN_OK;
358 }
359
360 // Some waves are on-shore, so interpolate the wave attributes from all profile points to all within-polygon sea cells, also update the active zone status for each cell
361 int nRet = nInterpolateWavesToPolygonCells(&VdXAll, &VdYAll, &VdHeightXAll, &VdHeightYAll);
362
363 if (nRet != RTN_OK)
364 return nRet;
365
366 // // DEBUG CODE ===========================================================================================================
367 // string strOutFile = m_strOutPath;
368 // strOutFile += "sea_wave_height_CHECKPOINT_1_";
369 // strOutFile += to_string(m_ulIter);
370 // strOutFile += ".tif";
371 //
372 // GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
373 // GDALDataset* pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
374 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
375 // pDataSet->SetGeoTransform(m_dGeoTransform);
376 //
377 // int nn = 0;
378 // double* pdRaster = new double[m_nXGridSize * m_nYGridSize];
379 // for (int nY = 0; nY < m_nYGridSize; nY++)
380 // {
381 // for (int nX = 0; nX < m_nXGridSize; nX++)
382 // {
383 // pdRaster[nn++] = m_pRasterGrid->m_Cell[nX][nY].dGetWaveHeight();
384 // }
385 // }
386 //
387 // GDALRasterBand* pBand = pDataSet->GetRasterBand(1);
388 // pBand->SetNoDataValue(m_dMissingValue);
389 // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
390 //
391 // if (nRet == CE_Failure)
392 // return RTN_ERR_GRIDCREATE;
393 //
394 // GDALClose(pDataSet);
395 // delete[] pdRaster;
396 // // DEBUG CODE ===========================================================================================================
397
398 // // DEBUG CODE ===========================================================================================================
399 // strOutFile = m_strOutPath;
400 // strOutFile += "sea_wave_angle_CHECKPOINT_1_";
401 // strOutFile += to_string(m_ulIter);
402 // strOutFile += ".tif";
403 //
404 // pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
405 // pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
406 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
407 // pDataSet->SetGeoTransform(m_dGeoTransform);
408 //
409 // nn = 0;
410 // pdRaster = new double[m_nXGridSize * m_nYGridSize];
411 // for (int nY = 0; nY < m_nYGridSize; nY++)
412 // {
413 // for (int nX = 0; nX < m_nXGridSize; nX++)
414 // {
415 // pdRaster[nn++] = m_pRasterGrid->m_Cell[nX][nY].dGetWaveAngle();
416 // }
417 // }
418 //
419 // pBand = pDataSet->GetRasterBand(1);
420 // pBand->SetNoDataValue(m_dMissingValue);
421 // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
422 //
423 // if (nRet == CE_Failure)
424 // return RTN_ERR_GRIDCREATE;
425 //
426 // GDALClose(pDataSet);
427 // delete[] pdRaster;
428 // // DEBUG CODE ===========================================================================================================
429
430 // Find any shadow zones and then modify waves in and adjacent to them
431 nRet = nDoAllShadowZones();
432
433 if (nRet != RTN_OK)
434 return nRet;
435
436 // Calculate the D50 for each polygon. Also fill in any artefactual 'holes' in active zone and wave property patterns
438
439 // // DEBUG CODE ===========================================================================================================
440 // strOutFile = m_strOutPath;
441 // strOutFile += "sea_wave_height_CHECKPOINT_2_";
442 // strOutFile += to_string(m_ulIter);
443 // strOutFile += ".tif";
444 //
445 // pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
446 // pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
447 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
448 // pDataSet->SetGeoTransform(m_dGeoTransform);
449 //
450 // nn = 0;
451 // pdRaster = new double[m_nXGridSize * m_nYGridSize];
452 // for (int nY = 0; nY < m_nYGridSize; nY++)
453 // {
454 // for (int nX = 0; nX < m_nXGridSize; nX++)
455 // {
456 // pdRaster[nn++] = m_pRasterGrid->m_Cell[nX][nY].dGetWaveHeight();
457 // }
458 // }
459 //
460 // pBand = pDataSet->GetRasterBand(1);
461 // pBand->SetNoDataValue(m_dMissingValue);
462 // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
463 //
464 // if (nRet == CE_Failure)
465 // return RTN_ERR_GRIDCREATE;
466 //
467 // GDALClose(pDataSet);
468 // delete[] pdRaster;
469 // // DEBUG CODE ===========================================================================================================
470
471 // // DEBUG CODE ===========================================================================================================
472 // strOutFile = m_strOutPath;
473 // strOutFile += "sea_wave_angle_CHECKPOINT_2_";
474 // strOutFile += to_string(m_ulIter);
475 // strOutFile += ".tif";
476 //
477 // pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
478 // pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
479 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
480 // pDataSet->SetGeoTransform(m_dGeoTransform);
481 //
482 // nn = 0;
483 // pdRaster = new double[m_nXGridSize * m_nYGridSize];
484 // for (int nY = 0; nY < m_nYGridSize; nY++)
485 // {
486 // for (int nX = 0; nX < m_nXGridSize; nX++)
487 // {
488 // pdRaster[nn++] = m_pRasterGrid->m_Cell[nX][nY].dGetWaveAngle();
489 // }
490 // }
491 //
492 // pBand = pDataSet->GetRasterBand(1);
493 // pBand->SetNoDataValue(m_dMissingValue);
494 // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
495 //
496 // if (nRet == CE_Failure)
497 // return RTN_ERR_GRIDCREATE;
498 //
499 // GDALClose(pDataSet);
500 // delete[] pdRaster;
501 // // DEBUG CODE ===========================================================================================================
502
503 // Modify the wave breaking properties (wave height, wave dir, breaking depth, breaking distance) for coastline points within the shadow zone
504 for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
505 {
506 int const nNumProfiles = m_VCoast[nCoast].nGetNumProfiles();
507
508 for (int nProfile = 0; nProfile < nNumProfiles; nProfile++)
510 }
511
512 // // DEBUG CODE ===========================================================================================================
513 // strOutFile = m_strOutPath;
514 // strOutFile += "sea_wave_height_CHECKPOINT_3_";
515 // strOutFile += to_string(m_ulIter);
516 // strOutFile += ".tif";
517 //
518 // pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
519 // pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
520 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
521 // pDataSet->SetGeoTransform(m_dGeoTransform);
522 //
523 // nn = 0;
524 // pdRaster = new double[m_nXGridSize * m_nYGridSize];
525 // for (int nY = 0; nY < m_nYGridSize; nY++)
526 // {
527 // for (int nX = 0; nX < m_nXGridSize; nX++)
528 // {
529 // pdRaster[nn++] = m_pRasterGrid->m_Cell[nX][nY].dGetWaveHeight();
530 // }
531 // }
532 //
533 // pBand = pDataSet->GetRasterBand(1);
534 // pBand->SetNoDataValue(m_dMissingValue);
535 // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
536 //
537 // if (nRet == CE_Failure)
538 // return RTN_ERR_GRIDCREATE;
539 //
540 // GDALClose(pDataSet);
541 // delete[] pdRaster;
542 // // DEBUG CODE ===========================================================================================================
543
544 // // DEBUG CODE ===========================================================================================================
545 // strOutFile = m_strOutPath;
546 // strOutFile += "sea_wave_angle_CHECKPOINT_3_";
547 // strOutFile += to_string(m_ulIter);
548 // strOutFile += ".tif";
549 //
550 // pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
551 // pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
552 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
553 // pDataSet->SetGeoTransform(m_dGeoTransform);
554 //
555 // nn = 0;
556 // pdRaster = new double[m_nXGridSize * m_nYGridSize];
557 // for (int nY = 0; nY < m_nYGridSize; nY++)
558 // {
559 // for (int nX = 0; nX < m_nXGridSize; nX++)
560 // {
561 // pdRaster[nn++] = m_pRasterGrid->m_Cell[nX][nY].dGetWaveAngle();
562 // }
563 // }
564 //
565 // pBand = pDataSet->GetRasterBand(1);
566 // pBand->SetNoDataValue(m_dMissingValue);
567 // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
568 //
569 // if (nRet == CE_Failure)
570 // return RTN_ERR_GRIDCREATE;
571 //
572 // GDALClose(pDataSet);
573 // delete[] pdRaster;
574 // // DEBUG CODE ===========================================================================================================
575
576 // Interpolate these wave properties for all remaining coastline points. Do this in along-coastline sequence
577 for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
578 {
579 int const nCoastSize = m_VCoast[nCoast].nGetCoastlineSize();
580 int const nNumProfiles = m_VCoast[nCoast].nGetNumProfiles();
581
582 // Interpolate these wave properties for all remaining coastline points. Do this in along-coastline sequence, but do not do this for the end-of-coastline profile (which is the final one)
583 for (int n = 0; n < nNumProfiles - 1; n++)
585
586 // And do the same for the coastline cells
588
589 // Calculate wave energy at breaking for every point on the coastline
590 for (int nCoastPoint = 0; nCoastPoint < nCoastSize; nCoastPoint++)
591 {
592 // Equation 4 from Walkden & Hall, 2005
593 double const dBreakingWaveHeight = m_VCoast[nCoast].dGetBreakingWaveHeight(nCoastPoint);
594 double const dCoastPointWavePeriod = m_VCoast[nCoast].dGetCoastDeepWaterWavePeriod(nCoastPoint);
595
596 // TODO 080 Why do we get -ve dBreakingWaveHeight here?
597 if (bFPIsEqual(dBreakingWaveHeight, DBL_NODATA, TOLERANCE))
598 {
599 m_VCoast[nCoast].SetBreakingWaveHeight(nCoastPoint, 0);
600 }
601
602 else
603 {
604 double const dErosiveWaveForce = pow(dBreakingWaveHeight, WALKDEN_HALL_PARAM_1) * pow(dCoastPointWavePeriod, WALKDEN_HALL_PARAM_2);
605
606 // Calculate total wave energy at this coast point during this timestep
607 double const dWaveEnergy = dErosiveWaveForce * m_dTimeStep * 3600;
608 m_VCoast[nCoast].SetWaveEnergyAtBreaking(nCoastPoint, dWaveEnergy);
609 }
610 }
611 }
612
613 return RTN_OK;
614}
615
616//===============================================================================================================================
618//===============================================================================================================================
619double CSimulation::dCalcWaveAngleToCoastNormal(double const dCoastAngle, double const dWaveAngle, int const nSeaHand)
620{
621 double dWaveToNormalAngle = 0;
622
623 if (nSeaHand == LEFT_HANDED)
624 // Left-handed coast
625 dWaveToNormalAngle = fmod((dWaveAngle - dCoastAngle + 360), 360) - 90;
626
627 else
628 // Right-handed coast
629 dWaveToNormalAngle = fmod((dWaveAngle - dCoastAngle + 360), 360) - 270;
630
631 if ((dWaveToNormalAngle >= 90) || (dWaveToNormalAngle <= -90))
632 dWaveToNormalAngle = DBL_NODATA;
633
634 return dWaveToNormalAngle;
635}
636
637//===============================================================================================================================
639//===============================================================================================================================
640int CSimulation::nCalcWavePropertiesOnProfile(int const nCoast, int const nCoastSize, CGeomProfile *pProfile, vector<double> *pVdX, vector<double> *pVdY, vector<double> *pVdHeightX, vector<double> *pVdHeightY, vector<bool> *pVbBreaking)
641{
642 // Only do this for profiles without problems. Still do start- and end-of-coast profiles however
643 if (! pProfile->bOKIncStartAndEndOfCoast())
644 {
645 // if (m_nLogFileDetail >= LOG_FILE_ALL)
646 // LogStream << m_ulIter << ": Coast " << nCoast << ", profile " << nProfile << " has been marked invalid, will not calc wave properties on this profile" << endl;
647
648 return RTN_OK;
649 }
650
651 // Calculate some wave properties based on the wave period following Airy wave theory
652 double const dDeepWaterWavePeriod = pProfile->dGetProfileDeepWaterWavePeriod();
653
654 m_dC_0 = (m_dG * dDeepWaterWavePeriod) / (2 * PI); // Deep water (offshore) wave celerity (m/s)
655 m_dL_0 = m_dC_0 * dDeepWaterWavePeriod; // Deep water (offshore) wave length (m)
656
657 int const nSeaHand = m_VCoast[nCoast].nGetSeaHandedness();
658 int const nCoastPoint = pProfile->nGetCoastPoint();
659
660 // Get the flux orientation (the orientation of a line which is tangential to the coast) at adjacent coastline points. Note special treatment for the coastline end points
661 double const dFluxOrientationThis = m_VCoast[nCoast].dGetFluxOrientation(nCoastPoint);
662 double dFluxOrientationPrev = 0;
663 double dFluxOrientationNext = 0;
664
665 if (nCoastPoint == 0)
666 {
667 dFluxOrientationPrev = dFluxOrientationThis;
668 dFluxOrientationNext = m_VCoast[nCoast].dGetFluxOrientation(1);
669 }
670
671 else if (nCoastPoint == nCoastSize - 1)
672 {
673 dFluxOrientationPrev = m_VCoast[nCoast].dGetFluxOrientation(nCoastPoint - 2);
674 dFluxOrientationNext = dFluxOrientationThis;
675 }
676
677 else
678 {
679 dFluxOrientationPrev = m_VCoast[nCoast].dGetFluxOrientation(nCoastPoint - 1);
680 dFluxOrientationNext = m_VCoast[nCoast].dGetFluxOrientation(nCoastPoint + 1);
681 }
682
683 // Get the deep water wave orientation for this profile
684 double const dDeepWaterWaveAngle = pProfile->dGetProfileDeepWaterWaveAngle();
685
686 // Calculate the angle between the deep water wave direction and a normal to the coast tangent
687 double dWaveToNormalAngle = dCalcWaveAngleToCoastNormal(dFluxOrientationThis, dDeepWaterWaveAngle, nSeaHand);
688
689 // Are the waves off-shore?
690 if (bFPIsEqual(dWaveToNormalAngle, DBL_NODATA, TOLERANCE))
691 {
692 // They are so, do nothing (each cell under the profile has already been initialised with deep water wave height and wave direction)
693 // LogStream << m_ulIter << ": profile " << nProfile << " has sea to " << (m_VCoast[nCoast].nGetSeaHandedness() == RIGHT_HANDED ? "right" : "left") << " dWaveToNormalAngle = " << dWaveToNormalAngle << " which is off-shore" << endl;
694
695 return RTN_OK;
696 }
697
698 // LogStream << m_ulIter << ": profile = " << nProfile << " has sea to " << (m_VCoast[nCoast].nGetSeaHandedness() == RIGHT_HANDED ? "right" : "left") << " dWaveToNormalAngle = " << dWaveToNormalAngle << " which is " << (dWaveToNormalAngle < 0 ? "DOWN" : "UP") << "-coast" << endl;
699
700 // Calculate the angle between the deep water wave direction and a normal to the coast tangent for the previous coast point
701 double dWaveToNormalAnglePrev;
702
703 if (nCoastPoint > 0)
704 {
705 // Get the deep water wave orientation for the up-coast point
706 double const dPrevDeepWaterWaveAngle = m_VCoast[nCoast].dGetCoastDeepWaterWaveAngle(nCoastPoint - 1);
707
708 dWaveToNormalAnglePrev = dCalcWaveAngleToCoastNormal(dFluxOrientationPrev, dPrevDeepWaterWaveAngle, nSeaHand);
709 }
710
711 else
712 {
713 dWaveToNormalAnglePrev = dWaveToNormalAngle;
714 }
715
716 // if (dWaveToNormalAnglePrev == DBL_NODATA)
717 // LogStream << "\tPrevious profile, dWaveToNormalAnglePrev = " << dWaveToNormalAnglePrev << " which is off-shore" << endl;
718 // else
719 // LogStream << "\tPrevious profile, dWaveToNormalAnglePrev = " << dWaveToNormalAnglePrev << " which is " << (dWaveToNormalAnglePrev < 0 ? "DOWN" : "UP") << "-coast" << endl;
720
721 // Calculate the angle between the deep water wave direction and a normal to the coast tangent for the next coast point
722 double dWaveToNormalAngleNext;
723
724 if (nCoastPoint < nCoastSize - 1)
725 {
726 // Get the deep water wave orientation for the down-coast point
727 double const dNextDeepWaterWaveAngle = m_VCoast[nCoast].dGetCoastDeepWaterWaveAngle(nCoastPoint + 1);
728
729 dWaveToNormalAngleNext = dCalcWaveAngleToCoastNormal(dFluxOrientationNext, dNextDeepWaterWaveAngle, nSeaHand);
730 }
731
732 else
733 {
734 dWaveToNormalAngleNext = dWaveToNormalAngle;
735 }
736
737 // if (dWaveToNormalAngleNext == DBL_NODATA)
738 // LogStream << "\tNext profile, dWaveToNormalAngleNext = " << dWaveToNormalAngleNext << " which is off-shore" << endl;
739 // else
740 // LogStream << "\tNext profile, dWaveToNormalAngleNext = " << dWaveToNormalAngleNext << " which is " << (dWaveToNormalAngleNext < 0 ? "DOWN" : "UP") << "-coast" << endl;
741
742 // Following Ashton and Murray (2006), if we have high-angle waves then use the flux orientation of the previous (up-coast) profile, if transitioning from diffusive to antidiffusive use flux maximizing angle (45 degrees)
743 if ((dWaveToNormalAngle > 0) && (! bFPIsEqual(dWaveToNormalAnglePrev, DBL_NODATA, TOLERANCE)) && (dWaveToNormalAnglePrev > 0))
744 {
745 if (dWaveToNormalAngle > 45)
746 {
747 if (dWaveToNormalAnglePrev < 45)
748 {
749 dWaveToNormalAngle = 45;
750 // LogStream << "\tA1" << endl;
751 }
752
753 else
754 {
755 dWaveToNormalAngle = dWaveToNormalAnglePrev;
756 // LogStream << "\tA2" << endl;
757 }
758 }
759 }
760
761 else if ((dWaveToNormalAngle < 0) && (! bFPIsEqual(dWaveToNormalAngleNext, DBL_NODATA, TOLERANCE)) && (dWaveToNormalAngleNext < 0))
762 {
763 if (dWaveToNormalAngle < -45)
764 {
765 if (dWaveToNormalAngleNext > -45)
766 {
767 dWaveToNormalAngle = -45;
768 // LogStream << "\tB1" << endl;
769 }
770
771 else
772 {
773 dWaveToNormalAngle = dWaveToNormalAngleNext;
774 // LogStream << "\tB2" << endl;
775 }
776 }
777 }
778
779 else if ((dWaveToNormalAngle > 45) && (! bFPIsEqual(dWaveToNormalAnglePrev, DBL_NODATA, TOLERANCE)) && (dWaveToNormalAnglePrev > 0))
780 {
781 // The wave direction here has an up-coast (decreasing indices) component: so for high-angle waves use the orientation from the up-coast (previous) profile
782 // LogStream << "\tCCC" << endl;
783
784 dWaveToNormalAngle = dFluxOrientationPrev;
785 }
786
787 else if ((dWaveToNormalAngle < -45) && (! bFPIsEqual(dWaveToNormalAngleNext, DBL_NODATA, TOLERANCE)) && (dWaveToNormalAngleNext < 0))
788 {
789 // The wave direction here has a down-coast (increasing indices) component: so for high-angle waves use the orientation from the down-coast (next) profile
790 // LogStream << "\tDDD" << endl;
791
792 dWaveToNormalAngle = dFluxOrientationNext;
793 }
794
795 // Initialize the wave properties at breaking for this profile
796 bool bBreaking = false;
797
798 int const nProfileSize = pProfile->nGetNumCellsInProfile();
799 int nProfileBreakingDist = 0;
800
801 double dProfileBreakingWaveHeight = DBL_NODATA;
802 double dProfileBreakingWaveAngle = 0;
803 double dProfileBreakingDepth = 0;
804 double dProfileWaveHeight = DBL_NODATA;
805 double const dProfileDeepWaterWaveHeight = pProfile->dGetProfileDeepWaterWaveHeight();
806 double dProfileWaveAngle = DBL_NODATA;
807 double const dProfileDeepWaterWaveAngle = pProfile->dGetProfileDeepWaterWaveAngle();
808
809 vector<bool> VbWaveIsBreaking(nProfileSize, 0);
810
811 vector<double> VdWaveHeight(nProfileSize, 0);
812 vector<double> VdWaveSetupSurge(nProfileSize, 0);
813 // vector<double> VdStormSurge(nProfileSize, 0);
814 vector<double> const VdWaveSetupRunUp(nProfileSize, 0);
815 vector<double> VdWaveDirection(nProfileSize, 0);
816
818 {
819 // We are using CShore to propagate the waves
820 double const dCShoreTimeStep = 3600; // In seconds, not important because we are not using CShore to erode the profile, just to get the hydrodynamics
821 double const dSurgeLevel = CSHORE_SURGE_LEVEL;
822
823 // Set up vectors for the coastline-normal profile elevations. The length of this vector line is given by the number of cells 'under' the profile. Thus each point on the vector relates to a single cell in the grid. This assumes that all points on the profile vector are equally spaced (not quite true, depends on the orientation of the line segments which comprise the profile)
824 vector<double> VdProfileZ; // Initial (pre-erosion) elevation of both consolidated and unconsolidated sediment for cells 'under' the profile, in CShore units
825 vector<double> VdProfileDistXY; // Along-profile distance measured from the seaward limit, in CShore units
826 vector<double> VdProfileFrictionFactor; // Along-profile friction factor from seaward limit
827
828 // The elevation of each of these profile points is the elevation of the centroid of the cell that is 'under' the point. However we cannot always be confident that this is the 'true' elevation of the point on the vector since (unless the profile runs planview N-S or W-E) the vector does not always run exactly through the centroid of the cell
829 int nRet = nGetThisProfileElevationsForCShore(nCoast, pProfile, nProfileSize, &VdProfileDistXY, &VdProfileZ, &VdProfileFrictionFactor);
830
831 if (nRet != RTN_OK)
832 {
833 // Could not create the profile elevation vectors
834 LogStream << m_ulIter << ": could not create CShore profile elevation vectors for profile " << pProfile->nGetProfileID() << endl;
835
836 return nRet;
837 }
838
839 // assert(static_cast<int>(VdProfileDistXY.size()) == nProfileSize);
840 // assert(static_cast<int>(VdProfileZ.size()) == nProfileSize);
841 // assert(static_cast<int>(VdProfileFrictionFactor.size()) == nProfileSize);
842
843 if (VdProfileDistXY.empty())
844 {
845 // The profile elevation vector was created, but was not populated
846 LogStream << m_ulIter << ": could not populate CShore profile elevation vector for profile " << pProfile->nGetProfileID() << endl;
847
849 }
850
851 // Constrain the wave to normal angle to be between -80 and 80 degrees, this is a requirement of CShore
852 dWaveToNormalAngle = tMax(dWaveToNormalAngle, -80.0);
853 dWaveToNormalAngle = tMin(dWaveToNormalAngle, 80.0);
854
855 int const nProfileDistXYSize = static_cast<int>(VdProfileDistXY.size());
856 vector<double> VdFreeSurfaceStd(nProfileDistXYSize, 0); // This is converted to Hrms by Hrms = sqr(8)*FreeSurfaceStd
857 vector<double> VdSinWaveAngleRadians(nProfileDistXYSize, 0); // This is converted to deg by asin(VdSinWaveAngleRadians)*(180/pi)
858 vector<double> VdFractionBreakingWaves(nProfileDistXYSize, 0); // Is 0 if no wave breaking, and 1 if all waves breaking
859
860 // Now define the other values that CShore requires
861 int const nILine = 1; // This is the number of cross-shore lines i.e. the number of CoastalME profiles. Only one at a time, at present
862 int const nIProfl = 0; // 0 for fixed bottom profile, 1 for profile evolution computation
863 int const nIPerm = 0; // 0 for impermeable bottom, 1 for permeable bottom of stone structure
864 int const nIOver = 0; // 0 for no wave overtopping and overflow on crest, 1 for wave overtopping and overflow
865 int const nIWCInt = 0; // 0 for no wave/current interaction, 1 for wave/current interaction in frequency dispersion, momentum and wave action equations
866 int const nIRoll = 0; // 0 for no roller effects, 1 for roller effects in governing equations
867 int const nIWind = 0; // 0 for no wind effects, 1 for wind shear stresses on momentum equations
868 int const nITide = 0; // 0 for no tidal effect on currents, 1 for longshore and cross-shore tidal currents
869 int const nILab = 0; // 0 for field data set, 1 for laboratory data set
870 int const nNWave = 1; // Number of waves at x = 0 starting from time = 0
871 int const nNSurge = 1; // Number of water levels at x = 0 from time = 0
872 double const dDX = VdProfileDistXY.back() / (static_cast<double>(VdProfileDistXY.size() - 1));
873 double const dWaveInitTime = 0; // CShore wave start time
874 double const dSurgeInitTime = 0; // CShore surge start time
875
876#if defined CSHORE_FILE_INOUT || CSHORE_BOTH
877 // Move to the CShore folder
878 nRet = chdir(CSHORE_DIR.c_str());
879
880 if (nRet != RTN_OK)
881 return nRet;
882
883#endif
884
885#if defined CSHORE_FILE_INOUT
886 // We are communicating with CShore using ASCII files, so create an input file for this profile which will be read by CShore
887 nRet = nCreateCShoreInfile(nCoast, nProfile, nILine, nIProfl, nIPerm, nIOver, nIWCInt, nIRoll, nIWind, nITide, nILab, nNWave, nNSurge, dDX, dCShoreTimeStep, dWaveInitTime, dDeepWaterWavePeriod, dProfileDeepWaterWaveHeight, dWaveToNormalAngle, dSurgeInitTime, dSurgeLevel, &VdProfileDistXY, &VdProfileZ, &VdProfileFrictionFactor);
888
889 if (nRet != RTN_OK)
890 return nRet;
891
892 // Set the error flag: this will be changed to 0 within CShore if CShore returns correctly
893 nRet = -1;
894
895 // Run CShore for this profile
896 CShore(&nRet);
897 // CShore(&nRet, &m_ulIter, &nProfile, &nProfile);
898
899 // Check return code for error
900 if (nRet != 0)
901 {
902 string strErr;
903
904 switch (nRet)
905 {
906 case -1:
907 strErr = to_string(m_ulIter) + ": CShore WARNING 1: negative depth at the first node ";
908 break;
909
910 case 2:
911 strErr = to_string(m_ulIter) + ": CShore WARNING 2: negative value at end of landward marching computation ";
912 break;
913
914 case 3:
915 strErr = to_string(m_ulIter) + ": CShore WARNING 3: large energy gradients at the first node: small waves with short period at sea boundary ";
916 break;
917
918 case 4:
919 strErr = to_string(m_ulIter) + ": CShore WARNING 4: zero energy at the first node ";
920 break;
921
922 case 5:
923 strErr = to_string(m_ulIter) + ": CShore WARNING 5: at end of landward marching computation, insufficient water depth ";
924 break;
925
926 case 7:
927 strErr = to_string(m_ulIter) + ": CShore WARNING 7: did not reach convergence ";
928 break;
929 }
930
931 strErr += "(coast " + to_string(nCoast) + " profile " + to_string(pProfile->nGetProfileID()) + " profile length " + to_string(nOutSize) + ")\n";
932
933 // OK, give up for this profile
934 // LogStream << strErr;
935 //
936 // return RTN_ERR_CSHORE_ERROR;
937 }
938
939 // Fetch the CShore results by reading files written by CShore
940 string strOSETUP = "OSETUP";
941 string strOYVELO = "OYVELO";
942 string strOPARAM = "OPARAM";
943
944 nRet = nReadCShoreOutput(nProfile, &strOSETUP, 4, 4, &VdProfileDistXY, &VdFreeSurfaceStd);
945
946 if (nRet != RTN_OK)
947 return nRet;
948
949 nRet = nReadCShoreOutput(nProfile, &strOYVELO, 4, 2, &VdProfileDistXY, &VdSinWaveAngleRadians);
950
951 if (nRet != RTN_OK)
952 return nRet;
953
954 nRet = nReadCShoreOutput(nProfile, &strOPARAM, 4, 3, &VdProfileDistXY, &VdFractionBreakingWaves);
955
956 if (nRet != RTN_OK)
957 return nRet;
958
959 // Read surge outputs
960 // VdTSurg = {dSurgeInitTime, dCShoreTimeStep}, // Ditto
961 // VdSWLin = {dSurgeLevel, dSurgeLevel}, // Ditto
962
963 // Clean up the CShore outputs
964#ifdef _WIN32
965 nRet = system("./clean.bat");
966#else
967
969 {
970 string strCommand = "./save_CShore_output.sh ";
971 strCommand += to_string(m_ulIter);
972 strCommand += " ";
973 strCommand += to_string(nCoast);
974 strCommand += " ";
975 strCommand += to_string(nProfile);
976
977 nRet = system(strCommand.c_str());
978
979 if (nRet != RTN_OK)
980 return nRet;
981 }
982
983 nRet = system("./clean.sh");
984#endif
985
986 if (nRet != RTN_OK)
987 return nRet;
988
989 // And return to the CoastalME folder
990 nRet = chdir(m_strCMEDir.c_str());
991
992 if (nRet != RTN_OK)
993 return nRet;
994
995#endif
996
997#if defined CSHORE_ARG_INOUT || CSHORE_BOTH
998 // We are communicating with CShore by passing arguments
999 int nOutSize = 0; // CShore will return the size of the output vectors
1000
1001 // Set the error flag: this will be changed within CShore if there is a problem
1002 nRet = 0;
1003
1004 vector<double> VdInitTime = {dWaveInitTime, dCShoreTimeStep}; // Size is nNwave+1, value 1 is for the start of the CShore run, value 2 for end of CShore run
1005 vector<double> VdTPIn = {dDeepWaterWavePeriod, dDeepWaterWavePeriod}; // Ditto
1006 vector<double> VdHrmsIn = {dProfileDeepWaterWaveHeight, dProfileDeepWaterWaveHeight}; // Ditto
1007 vector<double> VdWangIn = {dWaveToNormalAngle, dWaveToNormalAngle}; // Ditto
1008 vector<double> VdTSurg = {dSurgeInitTime, dCShoreTimeStep}; // Ditto
1009 vector<double> VdSWLin = {dSurgeLevel, dSurgeLevel}; // Ditto
1010 vector<double> VdFPInp = VdProfileFrictionFactor; // Set the value for wave friction at every point of the normal profile
1011 vector<double> VdXYDistFromCShoreOut(CSHOREARRAYOUTSIZE, 0); // Output from CShore
1012 vector<double> VdFreeSurfaceStdOut(CSHOREARRAYOUTSIZE, 0); // Ditto
1013 vector<double> VdWaveSetupSurgeOut(CSHOREARRAYOUTSIZE, 0); // Ditto
1014 // vector<double> VdStormSurgeOut(CSHOREARRAYOUTSIZE, 0); // Ditto
1015 vector<double> const VdWaveSetupRunUpOut(CSHOREARRAYOUTSIZE, 0); // Ditto
1016 vector<double> VdSinWaveAngleRadiansOut(CSHOREARRAYOUTSIZE, 0); // Ditto
1017 vector<double> VdFractionBreakingWavesOut(CSHOREARRAYOUTSIZE, 0); // Ditto
1018
1019 // Call CShore using the argument-passing wrapper
1020 // long lIter = static_cast<long>(m_ulIter); // Bodge to get round compiler 'invalid conversion' error
1021
1022 CShoreWrapper(&nILine, /* In_ILINE */
1023 &nIProfl, /* In_IPROFL */
1024 &nIPerm, /* In_IPERM */
1025 &nIOver, /* In_IOVER */
1026 &nIWCInt, /* In_IWCINT */
1027 &nIRoll, /* In_IROLL */
1028 &nIWind, /* In_IWIND */
1029 &nITide, /* In_ITIDE */
1030 &nILab, /* In_ILAB */
1031 &nNWave, /* In_NWAVE */
1032 &nNSurge, /* In_NSURG */
1033 &dDX, /* In_DX */
1034 &m_dBreakingWaveHeightDepthRatio, /* In_GAMMA */
1035 &VdInitTime[0], /* In_TWAVE */
1036 &VdTPIn[0], /* In_TPIN */
1037 &VdHrmsIn[0], /* In_HRMSIN */
1038 &VdWangIn[0], /* In_WANGIN */
1039 &VdTSurg[0], /* In_TSURG */
1040 &VdSWLin[0], /* In_SWLIN */
1041 &nProfileDistXYSize, /* In_NBINP */
1042 &VdProfileDistXY[0], /* In_XBINP */
1043 &VdProfileZ[0], /* In_ZBINP */
1044 &VdFPInp[0], /* In_FBINP */
1045 &nRet, /* Out_IError */
1046 &nOutSize, /* Out_nOutSize */
1047 &VdXYDistFromCShoreOut[0], /* Out_XYDist */
1048 &VdFreeSurfaceStdOut[0], /* Out_FreeSurfaceStd */
1049 &VdWaveSetupSurgeOut[0], /* Out_WaveSetupSurge */
1050 &VdSinWaveAngleRadiansOut[0], /* Out_SinWaveAngleRadians */
1051 &VdFractionBreakingWavesOut[0]); /* Out_FractionBreakingWaves */
1052
1053 // OK, now check for warnings and errors
1054 if (nOutSize < 2)
1055 {
1056 // CShore sometimes returns only one row of results, which contains data only for the seaward point of the profile. This happens when all other (more coastward) points give an invalid result during CShore's calculations. This is a problem. We don't want to abandon the simulation just because of this, so instead we just put some dummy data into the second row, and carry on with these two rows. The profile will get ignored later, since it is too small to be useful
1058 LogStream << m_ulIter << ": " << WARN << "for coast " << nCoast << " profile " << pProfile->nGetProfileID() << ", only " << nOutSize << " CShore output rows, abandoning this profile" << endl;
1059
1060 // // Set dummy data in the second row
1061 // VdXYDistFromCShoreOut[1] = 1e-5; // Dummy data, must not be the same as VdXYDistFromCShoreOut[0] tho', or get crash in linear interpolation routine
1062 // VdFreeSurfaceStdOut[1] = VdFreeSurfaceStdOut[0];
1063 // VdSinWaveAngleRadiansOut[1] = VdSinWaveAngleRadiansOut[0];
1064 // VdFractionBreakingWavesOut[1] = VdFractionBreakingWavesOut[0];
1065 // VdWaveSetupSurgeOut[1] = VdWaveSetupSurgeOut[0];
1066 // // VdStormSurgeOut[1] = VdStormSurgeOut[0];
1067 // // VdWaveSetupRunUpOut[1] = VdWaveSetupRunUpOut[0];
1068 //
1069 // // And increase the expected number of rows
1070 // nOutSize = 2;
1071
1072 return RTN_ERR_CSHORE_ERROR;
1073 }
1074
1075 if (nRet != RTN_OK)
1076 {
1077 string strErr;
1078
1079 switch (nRet)
1080 {
1081 case -1:
1082 strErr = to_string(m_ulIter) + ": CShore WARNING 1: negative depth at the first node ";
1083 break;
1084
1085 case 2:
1086 strErr = to_string(m_ulIter) + ": CShore WARNING 2: negative value at end of landward marching computation ";
1087 break;
1088
1089 case 3:
1090 strErr = to_string(m_ulIter) + ": CShore WARNING 3: large energy gradients at the first node: small waves with short period at sea boundary ";
1091 break;
1092
1093 case 4:
1094 strErr = to_string(m_ulIter) + ": CShore WARNING 4: zero energy at the first node ";
1095 break;
1096
1097 case 5:
1098 strErr = to_string(m_ulIter) + ": CShore WARNING 5: at end of landward marching computation, insufficient water depth ";
1099 break;
1100
1101 case 7:
1102 strErr = to_string(m_ulIter) + ": CShore WARNING 7: did not reach convergence ";
1103 break;
1104 }
1105
1106 strErr += "(coast " + to_string(nCoast) + " profile " + to_string(pProfile->nGetProfileID()) + " profile length " + to_string(nOutSize) + ")\n";
1107 LogStream << strErr;
1108
1109 // OK, give up for this profile
1110 // return RTN_ERR_CSHORE_ERROR;
1111 }
1112
1113 // LogStream << m_ulIter << ": interpolating profile " << nProfile << endl;
1114
1115 // All OK, so interpolate the CShore output, and convert from the CShore convention (cross-shore distance has its origin at the seaward end) to the CoastalME convention (origin at the shoreline)
1116 InterpolateCShoreOutput(&VdProfileDistXY, nOutSize, nProfileSize, &VdXYDistFromCShoreOut, &VdFreeSurfaceStdOut, &VdWaveSetupSurgeOut, &VdSinWaveAngleRadiansOut, &VdFractionBreakingWavesOut, &VdFreeSurfaceStd, &VdWaveSetupSurge, &VdSinWaveAngleRadians, &VdFractionBreakingWaves);
1117#endif
1118
1119#if defined CSHORE_BOTH
1120#if !defined _WIN32
1121
1123 {
1124 string strCommand = "./save_CShore_output.sh ";
1125 strCommand += to_string(m_ulIter);
1126 strCommand += " ";
1127 strCommand += to_string(nCoast);
1128 strCommand += " ";
1129 strCommand += to_string(nProfile);
1130
1131 nRet = system(strCommand.c_str());
1132
1133 if (nRet != RTN_OK)
1134 return nRet;
1135 }
1136
1137#endif
1138
1139 // Return to the CoastalME folder
1140 nRet = chdir(m_strCMEDir.c_str());
1141
1142 if (nRet != RTN_OK)
1143 return nRet;
1144
1145#endif
1146
1147 // Do more safety checks, then convert the CShore output to wave height and wave direction, and update wave profile attributes
1148 for (int nProfilePoint = (nProfileSize - 1); nProfilePoint >= 0; nProfilePoint--)
1149 {
1150 int const nX = pProfile->pPtiGetCellInProfile(nProfilePoint)->nGetX();
1151 int const nY = pProfile->pPtiGetCellInProfile(nProfilePoint)->nGetY();
1152
1153 // Safety check
1154 if (nProfilePoint > static_cast<int>(VdFreeSurfaceStd.size()) - 1)
1155 continue;
1156
1157 // Safety check: deal with NaN values
1158 if (isnan(VdFreeSurfaceStd[nProfilePoint]))
1159 VdFreeSurfaceStd[nProfilePoint] = 0;
1160
1161 VdWaveHeight[nProfilePoint] = sqrt(8) * VdFreeSurfaceStd[nProfilePoint];
1162
1163 // Another safety check: deal with NaN values
1164 if (isnan(VdSinWaveAngleRadians[nProfilePoint]))
1165 {
1166 VdSinWaveAngleRadians[nProfilePoint] = 0;
1167 VdWaveHeight[nProfilePoint] = 0;
1168 }
1169
1170 // More safety checks: constrain to the interval -1 to +1 to keep asin() happy
1171 if (VdSinWaveAngleRadians[nProfilePoint] < -1)
1172 VdSinWaveAngleRadians[nProfilePoint] = -1;
1173
1174 if (VdSinWaveAngleRadians[nProfilePoint] > 1)
1175 VdSinWaveAngleRadians[nProfilePoint] = 1;
1176
1177 double const dAlpha = asin(VdSinWaveAngleRadians[nProfilePoint]) * (180 / PI);
1178
1179 if (nSeaHand == LEFT_HANDED)
1180 VdWaveDirection[nProfilePoint] = dKeepWithin360(dAlpha + 90 + dFluxOrientationThis);
1181
1182 else
1183 VdWaveDirection[nProfilePoint] = dKeepWithin360(dAlpha + 270 + dFluxOrientationThis);
1184
1185 // Yet another safety check: deal with NaN values
1186 if (isnan(VdFractionBreakingWaves[nProfilePoint]))
1187 {
1188 VdFractionBreakingWaves[nProfilePoint] = 0;
1189 VdWaveHeight[nProfilePoint] = 0;
1190 }
1191
1192 // if ((VdFractionBreakingWaves[nProfilePoint] >= 0.10) && (! bBreaking)) // Sometimes is possible that waves break again
1193 if ((VdFractionBreakingWaves[nProfilePoint] >= 0.10) && (m_dDepthOfClosure >= m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth()) && (! bBreaking))
1194 {
1195 bBreaking = true;
1196 // assert(VdWaveHeight[nProfilePoint] >= 0);
1197 dProfileBreakingWaveHeight = VdWaveHeight[nProfilePoint];
1198 dProfileBreakingWaveAngle = VdWaveDirection[nProfilePoint];
1199 dProfileBreakingDepth = m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth(); // Water depth for the cell 'under' this point in the profile
1200 nProfileBreakingDist = nProfilePoint + 1; // At the nearest point nProfilePoint = 0, so, plus one
1201
1202 // LogStream << m_ulIter << ": CShore breaking at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "} nProfile = " << nProfile << ", nProfilePoint = " << nProfilePoint << ", dBreakingWaveHeight = " << dBreakingWaveHeight << ", dBreakingWaveAngle = " << dBreakingWaveAngle << ", dProfileBreakingDepth = " << dProfileBreakingDepth << ", nProfileBreakingDist = " << nProfileBreakingDist << endl;
1203 }
1204
1205 VbWaveIsBreaking[nProfilePoint] = bBreaking;
1206 }
1207
1208 if (dProfileBreakingWaveHeight >= dProfileDeepWaterWaveHeight)
1209 {
1210 dProfileBreakingWaveHeight = DBL_NODATA; // checking poorly conditions profiles problems for cshore
1211 }
1212 }
1213
1215 {
1216 // We are using COVE's linear wave theory to propagate the waves
1217 double const dDepthLookupMax = m_dWaveDepthRatioForWaveCalcs * dProfileDeepWaterWaveHeight;
1218
1219 // Go landwards along the profile, calculating wave height and wave angle for every inundated point on the profile (don't do point zero, this is on the coastline) until the waves start to break after breaking wave height is assumed to decrease linearly to zero at the shoreline and wave angle is equalt to wave angle at breaking
1220 for (int nProfilePoint = (nProfileSize - 1); nProfilePoint >= 0; nProfilePoint--)
1221 {
1222 int const nX = pProfile->pPtiGetCellInProfile(nProfilePoint)->nGetX();
1223 int const nY = pProfile->pPtiGetCellInProfile(nProfilePoint)->nGetY();
1224
1225 // Safety check
1226 if (! m_pRasterGrid->m_Cell[nX][nY].bIsInContiguousSea())
1227 continue;
1228
1229 double const dSeaDepth = m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth(); // Water depth for the cell 'under' this point in the profile
1230
1231 if (dSeaDepth > dDepthLookupMax)
1232 {
1233 // Sea depth is too large relative to wave height to feel the bottom, so do nothing since each cell under the profile has already been initialised with deep water wave height and wave direction
1234 dProfileWaveHeight = dProfileDeepWaterWaveHeight;
1235 dProfileWaveAngle = dProfileDeepWaterWaveAngle;
1236 }
1237
1238 else
1239 {
1240 if (! bBreaking)
1241 {
1242 // Start calculating wave properties using linear wave theory
1243 double const dL = m_dL_0 * sqrt(tanh((2 * PI * dSeaDepth) / m_dL_0)); // Wavelength (m) in intermediate-shallow waters
1244 double const dC = m_dC_0 * tanh((2 * PI * dSeaDepth) / dL); // Wave speed (m/s) set by dSeaDepth, dL and m_dC_0
1245 double const dk = 2 * PI / dL; // Wave number (1/m)
1246 double const dn = ((2 * dSeaDepth * dk) / (sinh(2 * dSeaDepth * dk)) + 1) / 2; // Shoaling factor
1247 double const dKs = sqrt(m_dC_0 / (dn * dC * 2)); // Shoaling coefficient
1248 double const dAlpha = (180 / PI) * asin((dC / m_dC_0) * sin((PI / 180) * dWaveToNormalAngle)); // Calculate angle between wave direction and the normal to the coast tangent
1249 double const dKr = sqrt(cos((PI / 180) * dWaveToNormalAngle) / cos((PI / 180) * dAlpha)); // Refraction coefficient
1250 dProfileWaveHeight = dProfileDeepWaterWaveHeight * dKs * dKr; // Calculate wave height, based on the previous (more seaward) wave height
1251
1252 if (nSeaHand == LEFT_HANDED)
1253 dProfileWaveAngle = dKeepWithin360(dAlpha + 90 + dFluxOrientationThis);
1254
1255 else
1256 dProfileWaveAngle = dKeepWithin360(dAlpha + 270 + dFluxOrientationThis);
1257
1258 // Test to see if the wave breaks at this depth
1259 if (dProfileWaveHeight < (dSeaDepth * m_dBreakingWaveHeightDepthRatio))
1260 {
1261 dProfileBreakingWaveHeight = dProfileWaveHeight;
1262 dProfileBreakingWaveAngle = dProfileWaveAngle;
1263 }
1264 }
1265
1266 else
1267 {
1268 // It does
1269 bBreaking = true;
1270
1271 // Wave has already broken
1272 dProfileWaveAngle = dProfileBreakingWaveAngle; // Wave orientation remains equal to wave orientation at breaking
1273
1274 dProfileWaveHeight = dProfileBreakingWaveHeight * (nProfilePoint / nProfileBreakingDist); // Wave height decreases linearly to zero at shoreline
1275 // dProfileWaveHeight = dSeaDepth * m_dBreakingWaveHeightDepthRatio; // Wave height is limited by depth
1276 dProfileBreakingDepth = dSeaDepth;
1277 nProfileBreakingDist = nProfilePoint;
1278 }
1279 }
1280
1281 // Save current wave attributes
1282 VdWaveDirection[nProfilePoint] = dProfileWaveAngle;
1283 VdWaveHeight[nProfilePoint] = dProfileWaveHeight;
1284 VbWaveIsBreaking[nProfilePoint] = bBreaking;
1285 }
1286 }
1287
1288 // Go landwards along the profile, fetching the calculated wave height and wave angle for every inundated point on this profile
1289 for (int nProfilePoint = (nProfileSize - 1); nProfilePoint >= 0; nProfilePoint--)
1290 {
1291 int nX = pProfile->pPtiGetCellInProfile(nProfilePoint)->nGetX();
1292 int nY = pProfile->pPtiGetCellInProfile(nProfilePoint)->nGetY();
1293
1294 // Safety check
1295 if (! m_pRasterGrid->m_Cell[nX][nY].bIsInContiguousSea())
1296 continue;
1297
1298 // Get the wave attributes calculated for this profile: wave height, wave angle, and whether is in the active zone
1299 double const dWaveHeight = VdWaveHeight[nProfilePoint];
1300 double const dWaveAngle = VdWaveDirection[nProfilePoint];
1301
1302 bBreaking = VbWaveIsBreaking[nProfilePoint];
1303
1304 // And store the wave properties for this point in the all-profiles vectors
1305 pVdX->push_back(nX);
1306 pVdY->push_back(nY);
1307 pVdHeightX->push_back(dWaveHeight * sin(dWaveAngle * PI / 180));
1308 pVdHeightY->push_back(dWaveHeight * cos(dWaveAngle * PI / 180));
1309 pVbBreaking->push_back(bBreaking);
1310 }
1311
1312 // Obtain the profile nodes near the coast
1313 int const nX = pProfile->pPtiGetCellInProfile(nProfileSize - 2)->nGetX();
1314 int const nY = pProfile->pPtiGetCellInProfile(nProfileSize - 2)->nGetY();
1315 int const nX1 = pProfile->pPtiGetCellInProfile(nProfileSize - 1)->nGetX();
1316 int const nY1 = pProfile->pPtiGetCellInProfile(nProfileSize - 1)->nGetY();
1317
1318 // Calculate the horizontal distance between the profile points
1319 double dXDist = tAbs(dGridCentroidXToExtCRSX(nX1) - dGridCentroidXToExtCRSX(nX));
1320 double dYDist = tAbs(dGridCentroidYToExtCRSY(nY1) - dGridCentroidYToExtCRSY(nY));
1321
1322 // Safety check
1323 if (bFPIsEqual(dXDist, 0.0, TOLERANCE))
1324 dXDist = 1e-3;
1325
1326 // Safety check
1327 if (bFPIsEqual(dYDist, 0.0, TOLERANCE))
1328 dYDist = 1e-3;
1329
1330 double const dDiffProfileDistXY = hypot(dXDist, dYDist);
1331
1332 // Compute the beach slope
1333 double const dtanBeta = tan(tAbs(m_pRasterGrid->m_Cell[nX1][nY1].dGetSeaDepth() - m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth()) / dDiffProfileDistXY);
1334
1335 // Compute the wave run-up using NIELSEN & HANSLOW (1991) & DHI (2004)
1336 int nValidPointsWaveHeight = 0;
1337 int nValidPointsWaveSetup = 0;
1338
1339 for (int nPoint = 0; nPoint < static_cast<int>(VdWaveHeight.size()); nPoint++)
1340 {
1341 if (VdWaveHeight[nPoint] > 1e-4)
1342 {
1343 nValidPointsWaveHeight += 1;
1344 }
1345
1346 else
1347 {
1348 break;
1349 }
1350 }
1351
1352 nValidPointsWaveHeight -= 1;
1353
1354 for (int nPoint = 0; nPoint < static_cast<int>(VdWaveSetupSurge.size()); nPoint++)
1355 {
1356 if (tAbs(VdWaveSetupSurge[nPoint]) < 1) // limiting the absolute value of setup + surge if cshore run fails
1357 {
1358 nValidPointsWaveSetup += 1;
1359 }
1360
1361 else
1362 {
1363 break;
1364 }
1365 }
1366
1367 nValidPointsWaveSetup -= 1;
1368
1369 double dWaveHeight = 0;
1370
1371 // Safety checks
1372 if ((nValidPointsWaveHeight >= 0) && (! bFPIsEqual(VdWaveHeight[nValidPointsWaveHeight], DBL_NODATA, TOLERANCE)))
1373 {
1374 dWaveHeight = VdWaveHeight[nValidPointsWaveHeight];
1375 }
1376
1377 // TODO 060 Remove these 'magic numbers'
1378 double dRunUp = 0;
1379
1380 if (m_nRunUpEquation == 0)
1381 {
1382 // Compute the run-up using Nielsen & Hanslow (1991) & DHI (2004)
1383 dRunUp = 0.36 * pow(9.81, 0.5) * dtanBeta * pow(dWaveHeight, 0.5) * dDeepWaterWavePeriod;
1384 }
1385
1386 else if (m_nRunUpEquation == 1)
1387 {
1388 // Compute the run-up using MASE 1989
1389 double const dS0 = 2 * PI * dWaveHeight / (9.81 * dDeepWaterWavePeriod * dDeepWaterWavePeriod);
1390 dRunUp = 1.86 * dWaveHeight * pow(pow(dtanBeta / dS0, 0.5), 0.71);
1391 }
1392
1393 else if (m_nRunUpEquation == 2)
1394 {
1395 // Compute the run-up using STOCKDON (2006)
1396 double const dS0 = 2 * PI * dWaveHeight / (9.81 * dDeepWaterWavePeriod * dDeepWaterWavePeriod);
1397 // dRunUp = 1.1 * ((0.35 * dWaveHeight * (pow((1 / dS0) * dWaveHeight * dWaveHeight, 0.5))) + (((((1 / dS0) * dWaveHeight * dWaveHeight) * (0.563 * dWaveHeight * dWaveHeight + 0.0004)), 0.5)) / 2);
1398
1399 double const dH0OverL0 = (1 / dS0) * dWaveHeight;
1400 double const dTmp1 = 0.35 * dWaveHeight * pow(dH0OverL0, 0.5);
1401 double const dTmp2 = pow(dH0OverL0 * ((0.563 * dWaveHeight * dWaveHeight) + 0.0004), 0.5);
1402 dRunUp = 1.1 * (dTmp1 + (dTmp2 / 2));
1403 }
1404
1405 if ((tAbs(dRunUp) < 1e-4) || (isnan(dRunUp)))
1406 {
1407 dRunUp = 0;
1408 }
1409
1410 double dWaveSetupSurge = 0;
1411
1412 // Safety checks
1413 if ((nValidPointsWaveSetup >= 0) && (! bFPIsEqual(VdWaveSetupSurge[nValidPointsWaveSetup], DBL_NODATA, TOLERANCE)))
1414 {
1415 dWaveSetupSurge = VdWaveSetupSurge[nValidPointsWaveSetup];
1416 }
1417
1418 if ((tAbs(dWaveSetupSurge) < 1e-4) || (isnan(dWaveSetupSurge)))
1419 {
1420 dWaveSetupSurge = 0;
1421 }
1422
1423 // Update wave attributes along the coastline object. Wave height at the coast is always calculated (i.e. whether or not waves are breaking)
1424 // cout << "Wave Height at the coast is " << VdWaveHeight[nProfileSize - 1] << endl;
1425 m_VCoast[nCoast].SetCoastWaveHeight(nCoastPoint, dWaveHeight);
1426 m_VCoast[nCoast].SetWaveSetupSurge(nCoastPoint, dWaveSetupSurge);
1427 m_VCoast[nCoast].SetRunUp(nCoastPoint, dRunUp);
1428
1429 if (nProfileBreakingDist > 0)
1430 {
1431 // This coast point is in the active zone, so set breaking wave height, breaking wave angle, and depth of breaking for the coast point
1432 m_VCoast[nCoast].SetBreakingWaveHeight(nCoastPoint, dProfileBreakingWaveHeight);
1433 m_VCoast[nCoast].SetBreakingWaveAngle(nCoastPoint, dProfileBreakingWaveAngle);
1434 m_VCoast[nCoast].SetDepthOfBreaking(nCoastPoint, dProfileBreakingDepth);
1435 m_VCoast[nCoast].SetBreakingDistance(nCoastPoint, nProfileBreakingDist);
1436
1437 // LogStream << m_ulIter << ": nProfile = " << nProfile << ", nCoastPoint = " << nCoastPoint << " in active zone, dBreakingWaveHeight = " << dBreakingWaveHeight << endl;
1438 }
1439
1440 else
1441 {
1442 // This coast point is not in the active zone, so breaking wave height, breaking wave angle, and depth of breaking are all meaningless
1443 m_VCoast[nCoast].SetBreakingWaveHeight(nCoastPoint, DBL_NODATA);
1444 m_VCoast[nCoast].SetBreakingWaveAngle(nCoastPoint, DBL_NODATA);
1445 m_VCoast[nCoast].SetDepthOfBreaking(nCoastPoint, DBL_NODATA);
1446 m_VCoast[nCoast].SetBreakingDistance(nCoastPoint, INT_NODATA);
1447
1448 // LogStream << m_ulIter << ": nProfile = " << nProfile << ", nCoastPoint = " << nCoastPoint << " NOT in active zone" << endl;
1449 }
1450
1451 return RTN_OK;
1452}
1453
1454#if defined CSHORE_FILE_INOUT
1455//===============================================================================================================================
1457//===============================================================================================================================
1458int CSimulation::nCreateCShoreInfile(int const nCoast, int const nProfile, int const nILine, int const nIProfl, int const nIPerm, int const nIOver, int const nIWcint, int const nIRoll, int const nIWind, int const nITide, int const nILab, int const nWave, int const nSurge, double const dX, double const dTimestep, double const dWaveInitTime, double const dWavePeriod, double const dHrms, double const dWaveAngle, double const dSurgeInitTime, double const dSurgeLevel, vector<double> const *pVdXdist, vector<double> const *pVdBottomElevation, vector<double> const *pVdWaveFriction)
1459{
1460 // Create the CShore input file
1461 ofstream CShoreOutStream;
1462 CShoreOutStream.open(CSHORE_INFILE.c_str(), ios::out | ios::app);
1463
1464 if (CShoreOutStream.fail())
1465 {
1466 // Error, cannot open file for writing
1467 LogStream << m_ulIter << ": " << ERR << "cannot write to CShore input file '" << CSHORE_INFILE << "'" << endl;
1469 }
1470
1471 // And write to the file
1472 CShoreOutStream << 3 << endl; // Number of comment lines
1473 CShoreOutStream << "------------------------------------------------------------" << endl;
1474 CShoreOutStream << "CShore input file created by CoastalME for iteration " << m_ulIter << ", coast " << nCoast << ", profile " << nProfile << endl;
1475 CShoreOutStream << "------------------------------------------------------------" << endl;
1476 CShoreOutStream << nILine << " -> ILINE" << endl;
1477 CShoreOutStream << nIProfl << " -> IPROFL" << endl;
1478 CShoreOutStream << nIPerm << " -> IPERM" << endl;
1479 CShoreOutStream << nIOver << " -> IOVER" << endl;
1480 CShoreOutStream << nIWcint << " -> IWCINT" << endl;
1481 CShoreOutStream << nIRoll << " -> IROLL" << endl;
1482 CShoreOutStream << nIWind << " -> IWIND" << endl;
1483 CShoreOutStream << nITide << " -> ITIDE" << endl;
1484 CShoreOutStream << fixed;
1485 CShoreOutStream << setw(11) << setprecision(4) << dX << " -> DX" << endl;
1486 CShoreOutStream << setw(11) << m_dBreakingWaveHeightDepthRatio << " -> GAMMA" << endl;
1487 CShoreOutStream << setw(11) << nILab << " -> ILAB" << endl;
1488 CShoreOutStream << setw(11) << nWave << " -> NWAVE" << endl;
1489 CShoreOutStream << setw(11) << nSurge << " -> NSURGE" << endl;
1490
1491 // Line 18 of infile
1492 CShoreOutStream << setw(11) << setprecision(2) << dWaveInitTime; // TWAVE(1) in CShore
1493 CShoreOutStream << setw(11) << setprecision(4) << dWavePeriod; // TPIN(1) in CShore
1494 CShoreOutStream << setw(11) << dHrms; // HRMS(1) in CShore
1495 CShoreOutStream << setw(11) << dWaveAngle << endl; // WANGIN(1) in CShore
1496
1497 // Line 19 of infile
1498 CShoreOutStream << setw(11) << setprecision(2) << dTimestep; // TWAVE(2) in CShore
1499 CShoreOutStream << setw(11) << setprecision(4) << dWavePeriod; // TPIN(2) in CShore
1500 CShoreOutStream << setw(11) << dHrms; // HRMS(2) in CShore
1501 CShoreOutStream << setw(11) << dWaveAngle << endl; // WANGIN(2) in CShore
1502
1503 // Line 20 of infile
1504 CShoreOutStream << setw(11) << setprecision(2) << dSurgeInitTime; // TSURG(1) in CShore
1505 CShoreOutStream << setw(11) << setprecision(4) << dSurgeLevel << endl; // SWLIN(1) in CShore
1506
1507 // Line 21 of infile
1508 CShoreOutStream << setw(11) << setprecision(2) << dTimestep; // TSURG(2) in CShore
1509 CShoreOutStream << setw(11) << setprecision(4) << dSurgeLevel << endl; // SWLIN(2) in CShore
1510
1511 // Line 22 of infile
1512 CShoreOutStream << setw(8) << pVdXdist->size() << " -> NBINP" << endl;
1513
1514 CShoreOutStream << fixed << setprecision(4);
1515
1516 for (unsigned int i = 0; i < pVdXdist->size(); i++)
1517 // These are BINP(J,1), ZBINP(J,1), FBINP(J-1,1) in CShore
1518 CShoreOutStream << setw(11) << pVdXdist->at(i) << setw(11) << pVdBottomElevation->at(i) << setw(11) << pVdWaveFriction->at(i) << endl;
1519
1520 CShoreOutStream << endl;
1521
1522 // File written, so close it
1523 CShoreOutStream.close();
1524
1525 return RTN_OK;
1526}
1527#endif
1528
1529//===============================================================================================================================
1531//===============================================================================================================================
1532int CSimulation::nGetThisProfileElevationsForCShore(int const nCoast, CGeomProfile *pProfile, int const nProfSize, vector<double> *VdDistXY, vector<double> *VdVZ, vector<double> *VdFricF)
1533{
1534 bool bIsBehindIntervention = false;
1535
1536 int nX1 = 0;
1537 int nY1 = 0;
1538
1539 double dXDist;
1540 double dYDist;
1541 double dProfileDistXY = 0;
1542 double dProfileFricFact;
1543 double dPrevDist = -1;
1544
1545 for (int i = nProfSize - 1; i >= 0; i--)
1546 {
1547 int const nX = pProfile->pPtiVGetCellsInProfile()->at(i).nGetX();
1548 int const nY = pProfile->pPtiVGetCellsInProfile()->at(i).nGetY();
1549
1550 // Calculate the horizontal distance relative to the most seaward point
1551 if (i == nProfSize - 1)
1552 dProfileDistXY = 0;
1553
1554 else
1555 {
1558 dProfileDistXY = dProfileDistXY + hypot(dXDist, dYDist);
1559 }
1560
1561 // Before we store the X-Y distance, must check that it is not the same as the previously-stored distance (if it is, we get zero-divide errors in CShore). If they are the same, add on a small distance
1562 if (bFPIsEqual(dProfileDistXY, dPrevDist, TOLERANCE))
1563 dProfileDistXY += 0.1; // TODO 084 Improve this
1564
1565 // Update the cell indexes, the initial cell is now the previous one
1566 nX1 = nX;
1567 nY1 = nY;
1568
1569 // Get the number of the highest layer with non-zero thickness
1570 int const nTopLayer = m_pRasterGrid->m_Cell[nX][nY].nGetTopNonZeroLayerAboveBasement();
1571
1572 // Safety checks
1573 if (nTopLayer == INT_NODATA)
1574 return RTN_ERR_NO_TOP_LAYER;
1575
1576 if (nTopLayer == NO_NONZERO_THICKNESS_LAYERS)
1577 // TODO 009 We are down to basement, decide what to do
1578 return RTN_OK;
1579
1580 // Get the elevation for both consolidated and unconsolidated sediment on this cell
1581 double const dTopElev = m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev() + m_pRasterGrid->m_Cell[nX][nY].dGetInterventionHeight();
1582 double const VdProfileZ = dTopElev - m_dThisIterSWL;
1583
1584 // Check that landward elevation is greater than SWL
1585 if (i == 0)
1586 {
1587 if (VdProfileZ < 0)
1588 {
1589 VdVZ->push_back(0.1); // TODO 053 Set it to a small +ve elevation (compared with SWL). However there must be a better way of doing this
1590
1591 // Could not create the profile elevation vectors
1592 LogStream << m_ulIter << ": " << WARN << "for coast " << nCoast << ", profile " << pProfile->nGetProfileID() << ", elevation at the landward end is " << dTopElev << " m. This is lower than this-iteration SWL (" << m_dThisIterSWL << " m). For CShore, changing the landward elevation for profile " << pProfile->nGetProfileID() << " to " << m_dThisIterSWL + 0.1 << "m" << endl;
1593 }
1594
1595 else
1596 {
1597 VdVZ->push_back(VdProfileZ);
1598 }
1599 }
1600
1601 else
1602 {
1603 VdVZ->push_back(VdProfileZ);
1604 }
1605
1606 // Now store the X-Y plane distance from the start of the profile
1607 VdDistXY->push_back(dProfileDistXY);
1608
1609 // Get the landform type at each point along the profile
1610 double const dInterventionHeight = m_pRasterGrid->m_Cell[nX][nY].dGetInterventionHeight();
1611
1612 // Modify default friction factor if a structural intervention is found, otherwise use the default
1613 if (dInterventionHeight > 0 || bIsBehindIntervention)
1614 {
1615 // Use an arbitrarily high value if structure is present
1616 dProfileFricFact = 100 * CSHORE_FRICTION_FACTOR;
1617 bIsBehindIntervention = true;
1618 }
1619
1620 else
1621 dProfileFricFact = CSHORE_FRICTION_FACTOR;
1622
1623 // Store the friction factor
1624 VdFricF->push_back(dProfileFricFact);
1625
1626 // For next time round
1627 dPrevDist = dProfileDistXY;
1628 }
1629
1630 return RTN_OK;
1631}
1632
1633#if defined CSHORE_FILE_INOUT
1634//===============================================================================================================================
1636//===============================================================================================================================
1637int CSimulation::nReadCShoreOutput(int const nProfile, string const *strCShoreFilename, int const nExpectedColumns, int const nCShorecolumn, vector<double> const *pVdProfileDistXYCME, vector<double> *pVdInterpolatedValues)
1638{
1639 // Read in the first column (contains XY distance relative to seaward limit) and CShore column from the CShore output file
1640 ifstream InStream;
1641 InStream.open(strCShoreFilename->c_str(), ios::in);
1642
1643 // Did it open OK?
1644 if (!InStream.is_open())
1645 {
1646 // Error: cannot open CShore file for input
1647 LogStream << m_ulIter << ": " << ERR << "for profile " << nProfile << ", cannot open " << *strCShoreFilename << " for input" << endl;
1648
1650 }
1651
1652 // Opened OK, so set up the vectors to hold the CShore output data
1653 vector<double> VdXYDistCShore;
1654 vector<double> VdValuesCShore;
1655
1656 // And read in the data
1657 int n = -1;
1658 int nExpectedRows = 0;
1659 string strLineIn;
1660
1661 while (getline(InStream, strLineIn))
1662 {
1663 n++;
1664
1665 if (n == 0)
1666 {
1667 // Read in the header line
1668 vector<string> VstrItems = VstrSplit(&strLineIn, SPACE);
1669
1670 if (! bIsStringValidInt(VstrItems[1]))
1671 {
1672 string strErr = ERR + "invalid integer for number of expected rows '" + VstrItems[1] + "' in " + *strCShoreFilename + "\n";
1673 cerr << strErr;
1674 LogStream << strErr;
1675
1677 }
1678
1679 // Get the number of expected rows
1680 nExpectedRows = stoi(VstrItems[1]);
1681 }
1682
1683 else
1684 {
1685 // Read in a data line
1686 vector<string> VstrItems = VstrSplit(&strLineIn, SPACE);
1687
1688 int nCols = static_cast<int>(VstrItems.size());
1689
1690 if (nCols != nExpectedColumns)
1691 {
1692 // Error: did not read the expected number of CShore output columns
1693 LogStream << m_ulIter << ": " << ERR << "for profile " << nProfile << ", expected " << nExpectedColumns << " CShore output columns but read " << nCols << " columns from header section of file " << *strCShoreFilename << endl;
1694
1696 }
1697
1698 // Number of columns is OK
1699 VdXYDistCShore.push_back(strtod(VstrItems[0].c_str(), NULL));
1700 VdValuesCShore.push_back(strtod(VstrItems[nCShorecolumn - 1].c_str(), NULL));
1701 }
1702 }
1703
1704 // Check that we have read nExpectedRows from the file
1705 int nReadRows = static_cast<int>(VdXYDistCShore.size());
1706
1707 if (nReadRows != nExpectedRows)
1708 {
1709 // Error: did not get nExpectedRows CShore output rows
1710 LogStream << m_ulIter << ": " << ERR << "for profile " << nProfile << ", expected " << nExpectedRows << " CShore output rows, but read " << nReadRows << " rows from file " << *strCShoreFilename << endl;
1711
1713 }
1714
1715 if (nReadRows < 2)
1716 {
1717 // CShore sometimes returns only one row, which contains data for the seaward point of the profile. This happens when all other (more coastward) points give an invalid result during CShore's calculations. This is a problem. We don't want to abandon the simulation just because of this, so instead we just duplicate the row, so that the profile will later get marked as invalid
1719 LogStream << m_ulIter << ": " << WARN << "for profile " << nProfile << ", only " << nReadRows << " CShore output rows in file " << *strCShoreFilename << endl;
1720
1721 // Duplicate the data
1722 VdXYDistCShore.push_back(VdXYDistCShore[0]);
1723 VdValuesCShore.push_back(VdValuesCShore[0]);
1724
1725 // And increase the expected number of rows
1726 nReadRows++;
1727 }
1728
1729 // The output is OK, so change the origin of the across-shore distance from the CShore convention to the one used here (i.e. with the origin at the shoreline)
1730 vector<double> VdXYDistCShoreTmp(nReadRows, 0);
1731
1732 for (int i = 0; i < nReadRows; i++)
1733 VdXYDistCShoreTmp[i] = VdXYDistCShore[nReadRows - 1] - VdXYDistCShore[i];
1734
1735 // Reverse the XY-distance and value vectors (i.e. first point is at the shoreline and must be in strictly ascending order)
1736 reverse(VdXYDistCShoreTmp.begin(), VdXYDistCShoreTmp.end());
1737
1738 // Similarly, reverse the CShore output
1739 reverse(VdValuesCShore.begin(), VdValuesCShore.end());
1740
1741 // Using a simple linear interpolation approach
1742 vector<double> VdDistXYCopy(pVdProfileDistXYCME->begin(), pVdProfileDistXYCME->end());
1743
1744 // assertVdXYDistCShoreTmp.size() == VdValuesCShore.size());
1745 *pVdInterpolatedValues = VdInterpolateCShoreProfileOutput(&VdXYDistCShoreTmp, &VdValuesCShore, &VdDistXYCopy);
1746
1747 return RTN_OK;
1748}
1749#endif
1750
1751#if defined CSHORE_ARG_INOUT || CSHORE_BOTH
1752//===============================================================================================================================
1754//===============================================================================================================================
1755void CSimulation::InterpolateCShoreOutput(vector<double> const *pVdProfileDistXYCME, int const nOutSize, int const nProfileSize, vector<double> *pVdXYDistFromCShore, vector<double> *pVdFreeSurfaceStdCShore, vector<double> *pVdWaveSetupSurgeCShore, vector<double> *pVdSinWaveAngleRadiansCShore, vector<double> *pVdFractionBreakingWavesCShore, vector<double> *pVdFreeSurfaceStdCME, vector<double> *pVdWaveSetupSurgeCME, vector<double> *pVdSinWaveAngleRadiansCME, vector<double> *pVdFractionBreakingWavesCME)
1756{
1757 // // DEBUG CODE ========================================================================================================
1758 // LogStream << m_ulIter << ": nOutSize = " << nOutSize << " nProfileSize = " << nProfileSize << " pVdProfileDistXYCME->size() = " << pVdProfileDistXYCME->size() << " pVdXYDistFromCShore->size() = " << pVdXYDistFromCShore->size() << endl;
1759 // // DEBUG CODE ========================================================================================================
1760
1761 // Sometimes the profile length returned from CShore is shorter than the CoastalME profile length
1762 bool bTooShort = false;
1763 int nTooShort = 0;
1764
1765 if (nOutSize < nProfileSize)
1766 {
1767 bTooShort = true;
1768 nTooShort = nProfileSize - nOutSize - 1;
1769
1770 // LogStream << m_ulIter << ": CShore PROFILE IS TOO SHORT" << endl;
1771 }
1772
1773 // // DEBUG CODE ========================================================================================================
1774 // for (int n = 0; n < static_cast<int>(pVdProfileDistXYCME->size()); n++)
1775 // LogStream << "pVdProfileDistXYCME[" << n << "] = " << pVdProfileDistXYCME->at(n) << endl;
1776 //
1777 // LogStream << endl;
1778 //
1779 // LogStream << "ORIGINAL" << endl;
1780 // for (int n = 0; n < nOutSize; n++)
1781 // LogStream << "pVdXYDistFromCShore[" << n << "] = " << pVdXYDistFromCShore->at(n) << " pVdFreeSurfaceStdCShore[" << n << "] = " << pVdFreeSurfaceStdCShore->at(n) << " pVdWaveSetupSurgeCShore[" << n << "] = " << pVdWaveSetupSurgeCShore->at(n) << " pVdSinWaveAngleRadiansCShore[" << n << "] = " << pVdSinWaveAngleRadiansCShore->at(n) << " pVdFractionBreakingWavesCShore[" << n << "] = " << pVdFractionBreakingWavesCShore->at(n) << endl;
1782 //
1783 // LogStream << endl;
1784 // // DEBUG CODE ========================================================================================================
1785
1786 if (bTooShort)
1787 {
1788 // Add extrapolated value(s) to the end of the valid part of each profile vector so that we have nProfileSize valid values
1789 double dLastDiff = pVdXYDistFromCShore->at(nOutSize - 1) - pVdXYDistFromCShore->at(nOutSize - 2);
1790
1791 for (int n = 0; n < nTooShort; n++)
1792 pVdXYDistFromCShore->at(nOutSize + n) = pVdXYDistFromCShore->at(nOutSize - 1 + n) + dLastDiff;
1793
1794 for (int n = 0; n < nTooShort; n++)
1795 pVdFreeSurfaceStdCShore->at(nOutSize + n) = pVdFreeSurfaceStdCShore->at(nOutSize - 1);
1796
1797 for (int n = 0; n < nTooShort; n++)
1798 pVdWaveSetupSurgeCShore->at(nOutSize + n) = pVdWaveSetupSurgeCShore->at(nOutSize - 1);
1799
1800 // TODO 007 Do same for pVdStormSurgeCShore and pVdWaveSetupRunUpCShore ?
1801
1802 for (int n = 0; n < nTooShort; n++)
1803 pVdSinWaveAngleRadiansCShore->at(nOutSize + n) = pVdSinWaveAngleRadiansCShore->at(nOutSize - 1);
1804
1805 dLastDiff = pVdFractionBreakingWavesCShore->at(nOutSize - 1) - pVdFractionBreakingWavesCShore->at(nOutSize - 2);
1806
1807 for (int n = 0; n < nTooShort; n++)
1808 pVdFractionBreakingWavesCShore->at(nOutSize + n) = tMin(pVdFractionBreakingWavesCShore->at(nOutSize - 1 + n) + dLastDiff, 1.0);
1809
1810 // // DEBUG CODE ========================================================================================================
1811 // LogStream << "EXTENDED" << endl;
1812 // for (int n = 0; n < nProfileSize; n++)
1813 // LogStream << "pVdXYDistFromCShore[" << n << "] = " << pVdXYDistFromCShore->at(n) << " pVdFreeSurfaceStdCShore[" << n << "] = " << pVdFreeSurfaceStdCShore->at(n) << " pVdWaveSetupSurgeCShore[" << n << "] = " << pVdWaveSetupSurgeCShore->at(n) << " pVdSinWaveAngleRadiansCShore[" << n << "] = " << pVdSinWaveAngleRadiansCShore->at(n) << " pVdFractionBreakingWavesCShore[" << n << "] = " << pVdFractionBreakingWavesCShore->at(n) << " " << (n > (nOutSize-1) ? "EXTRAPOLATED" : "") << endl;
1814 //
1815 // LogStream << endl;
1816 // // DEBUG CODE ========================================================================================================
1817 }
1818
1819 vector<double> VdXYDistCShoreTmp(nProfileSize);
1820 copy(pVdXYDistFromCShore->begin(), pVdXYDistFromCShore->begin() + nProfileSize, begin(VdXYDistCShoreTmp));
1821
1822 // The CShore cross-shore distance has its origin at the seaward end, but the CoastalME convention has its origin at the shoreline. So we mst reverse the other vectors to conform with the CoastalMEME convention
1823 vector<double> VdFreeSurfaceStdCShoreTmp(nProfileSize);
1824 reverse_copy(pVdFreeSurfaceStdCShore->begin(), pVdFreeSurfaceStdCShore->begin() + nProfileSize, begin(VdFreeSurfaceStdCShoreTmp));
1825
1826 vector<double> VdWaveSetupSurgeCShoreTmp(nProfileSize);
1827 reverse_copy(pVdWaveSetupSurgeCShore->begin(), pVdWaveSetupSurgeCShore->begin() + nProfileSize, begin(VdWaveSetupSurgeCShoreTmp));
1828
1829 vector<double> VdSinWaveAngleRadiansCShoreTmp(nProfileSize);
1830 reverse_copy(pVdSinWaveAngleRadiansCShore->begin(), pVdSinWaveAngleRadiansCShore->begin() + nProfileSize, begin(VdSinWaveAngleRadiansCShoreTmp));
1831
1832 vector<double> VdFractionBreakingWavesCShoreTmp(nProfileSize);
1833 reverse_copy(pVdFractionBreakingWavesCShore->begin(), pVdFractionBreakingWavesCShore->begin() + nProfileSize, begin(VdFractionBreakingWavesCShoreTmp));
1834
1835 // // DEBUG CODE ========================================================================================================
1836 // LogStream << "REVERSED" << endl;
1837 // for (int n = 0; n < nProfileSize; n++)
1838 // LogStream << "VdXYDistCShoreTmp[" << n << "] = " << VdXYDistCShoreTmp.at(n) << " VdFreeSurfaceStdCShoreTmp[" << n << "] = " << VdFreeSurfaceStdCShoreTmp.at(n) << " VdWaveSetupSurgeCShoreTmp[" << n << "] = " << VdWaveSetupSurgeCShoreTmp.at(n) << " VdSinWaveAngleRadiansCShoreTmp[" << n << "] = " << VdSinWaveAngleRadiansCShoreTmp.at(n) << " VdFractionBreakingWavesCShoreTmp[" << n << "] = " << VdFractionBreakingWavesCShoreTmp.at(n) << endl;
1839 //
1840 // LogStream << endl;
1841 // // DEBUG CODE ========================================================================================================
1842
1843 // Finally we linearly interpolate the CShore output vectors to each point on the CoastalME profile. Input parameters x_old, y_old, x_new and the routine returns y_new
1844 *pVdFreeSurfaceStdCME = VdInterpolateCShoreProfileOutput(&VdXYDistCShoreTmp, pVdFreeSurfaceStdCShore, pVdProfileDistXYCME); // was &VdDistXYCopy);
1845 *pVdWaveSetupSurgeCME = VdInterpolateCShoreProfileOutput(&VdXYDistCShoreTmp, pVdWaveSetupSurgeCShore, pVdProfileDistXYCME);
1846 // *pVdStormSurgeCME = VdInterpolateCShoreProfileOutput(&VdXYDistCShoreTmp, pVdStormSurgeCShore, pVdProfileDistXYCME);
1847 // *pVdWaveSetupRunUpCME = VdInterpolateCShoreProfileOutput(&VdXYDistCShoreTmp, pVdWaveSetupRunUpCShore, pVdProfileDistXYCME);
1848 *pVdSinWaveAngleRadiansCME = VdInterpolateCShoreProfileOutput(&VdXYDistCShoreTmp, pVdSinWaveAngleRadiansCShore, pVdProfileDistXYCME);
1849 *pVdFractionBreakingWavesCME = VdInterpolateCShoreProfileOutput(&VdXYDistCShoreTmp, pVdFractionBreakingWavesCShore, pVdProfileDistXYCME);
1850
1851 // LogStream << "INTERPOLATED" << endl;
1852 // for (int n = 0; n < nProfileSize; n++)
1853 // LogStream << "pVdProfileDistXYCME[" << n << "] = " << pVdProfileDistXYCME->at(n) << " pVdFreeSurfaceStdCME[" << n << "] = " << pVdFreeSurfaceStdCME->at(n) << " pVdWaveSetupSurgeCME[" << n << "] = " << pVdWaveSetupSurgeCME->at(n) << " pVdSinWaveAngleRadiansCME[" << n << "] = " << pVdSinWaveAngleRadiansCME->at(n) << " pVdFractionBreakingWavesCME[" << n << "] = " << pVdFractionBreakingWavesCME->at(n) << endl;
1854 //
1855 // LogStream << "================================================ " << endl;
1856}
1857#endif
1858
1859//===============================================================================================================================
1861//===============================================================================================================================
1863{
1864 CGeomProfile *pProfile = m_VCoast[nCoast].pGetProfile(nProfile);
1865
1866 // Only do this for profiles without problems, including the start and end-of-coast profile
1867 if (! pProfile->bOKIncStartAndEndOfCoast())
1868 return;
1869
1870 bool bModfiedWaveHeightisBreaking = false;
1871 bool bProfileIsinShadowZone = false;
1872 int const nThisCoastPoint = pProfile->nGetCoastPoint();
1873 int const nProfileSize = pProfile->nGetNumCellsInProfile();
1874 int nThisBreakingDist = m_VCoast[nCoast].nGetBreakingDistance(nThisCoastPoint);
1875 double dThisBreakingWaveHeight = m_VCoast[nCoast].dGetBreakingWaveHeight(nThisCoastPoint); // This could be DBL_NODATA
1876 double dThisBreakingWaveAngle = m_VCoast[nCoast].dGetBreakingWaveAngle(nThisCoastPoint);
1877 double dThisBreakingDepth = m_VCoast[nCoast].dGetDepthOfBreaking(nThisCoastPoint);
1878
1879 // Traverse the profile landwards, checking if any profile cell is within the shadow zone
1880 for (int nProfilePoint = (nProfileSize - 1); nProfilePoint >= 0; nProfilePoint--)
1881 {
1882 int const nX = pProfile->pPtiGetCellInProfile(nProfilePoint)->nGetX();
1883 int const nY = pProfile->pPtiGetCellInProfile(nProfilePoint)->nGetY();
1884
1885 // If there is any cell profile within the shadow zone and waves are breaking then modify wave breaking properties otherwise continue
1886 if (m_pRasterGrid->m_Cell[nX][nY].bIsinAnyShadowZone())
1887 {
1888 bProfileIsinShadowZone = true;
1889
1890 // Check if the new wave height is breaking
1891 double const dWaveHeight = m_pRasterGrid->m_Cell[nX][nY].dGetWaveHeight();
1892
1893 // Check that wave height at the given point is lower than maximum real wave height. If breaking wave height is expected that no good wave height are obtained, so, do not take it
1894 if (dWaveHeight > (m_dDepthOfClosure * m_dBreakingWaveHeightDepthRatio) && (! bModfiedWaveHeightisBreaking) && (! bFPIsEqual(dThisBreakingWaveHeight, DBL_NODATA, TOLERANCE)))
1895 {
1896 // It is breaking
1897 bModfiedWaveHeightisBreaking = true;
1898
1899 dThisBreakingWaveHeight = m_dDepthOfClosure * m_dBreakingWaveHeightDepthRatio;
1900 dThisBreakingWaveAngle = m_pRasterGrid->m_Cell[nX][nY].dGetWaveAngle();
1901 dThisBreakingDepth = m_dDepthOfClosure;
1902 nThisBreakingDist = nProfilePoint;
1903 }
1904 }
1905 }
1906
1907 // Update breaking wave properties along coastal line object (Wave height, dir, distance). TODO 010 Update the active zone cells
1908 if (bProfileIsinShadowZone && bModfiedWaveHeightisBreaking) // Modified wave height is still breaking
1909 {
1910 // This coast point is in the active zone, so set breaking wave height, breaking wave angle, and depth of breaking for the coast point TODO 007 Where does the 0.78 come from? TODO 011 Should it be an input variable or a named constant?
1911 if (dThisBreakingWaveHeight > dThisBreakingDepth * 0.78)
1912 {
1913 dThisBreakingWaveHeight = dThisBreakingDepth * 0.78; // Likely CShore output wave height is not adequately reproduced due to input profile and wave properties. TODO 007 Info needed. Does something need to be changed then?
1914 }
1915
1916 // assert(dThisBreakingWaveHeight >= 0);
1917 m_VCoast[nCoast].SetBreakingWaveHeight(nThisCoastPoint, dThisBreakingWaveHeight);
1918 m_VCoast[nCoast].SetBreakingWaveAngle(nThisCoastPoint, dThisBreakingWaveAngle);
1919 m_VCoast[nCoast].SetDepthOfBreaking(nThisCoastPoint, dThisBreakingDepth);
1920 m_VCoast[nCoast].SetBreakingDistance(nThisCoastPoint, nThisBreakingDist);
1921
1922 // LogStream << m_ulIter << ": nProfile = " << nProfile << ", nCoastPoint = " << nCoastPoint << " in active zone, dBreakingWaveHeight = " << dBreakingWaveHeight << endl;
1923 }
1924
1925 else if (bProfileIsinShadowZone && (! bModfiedWaveHeightisBreaking))
1926 {
1927 // This coast point is no longer in the active zone
1928 m_VCoast[nCoast].SetBreakingWaveHeight(nThisCoastPoint, DBL_NODATA);
1929 m_VCoast[nCoast].SetBreakingWaveAngle(nThisCoastPoint, DBL_NODATA);
1930 m_VCoast[nCoast].SetDepthOfBreaking(nThisCoastPoint, DBL_NODATA);
1931 m_VCoast[nCoast].SetBreakingDistance(nThisCoastPoint, INT_NODATA);
1932
1933 // LogStream << m_ulIter << ": nProfile = " << nProfile << ", nCoastPoint = " << nCoastPoint << " NOT in active zone" << endl;
1934 }
1935
1936 return;
1937}
1938
1939//===============================================================================================================================
1941//===============================================================================================================================
1942void CSimulation::InterpolateWavePropertiesBetweenProfiles(int const nCoast, int const nCount)
1943{
1944 CGeomProfile *pProfile = m_VCoast[nCoast].pGetProfileWithDownCoastSeq(nCount);
1945
1946 // Only do this for profiles without problems, including the start-of-coast profile (but not the end-of-coast profile)
1947 // if (! pProfile->bOKIncStartOfCoast())
1948 // return;
1949
1950 int const nThisCoastPoint = pProfile->nGetCoastPoint();
1951
1952 // For the breaking wave stuff, to go into the in-between coastline points
1953 int const nThisBreakingDist = m_VCoast[nCoast].nGetBreakingDistance(nThisCoastPoint);
1954 double const dThisBreakingWaveHeight = m_VCoast[nCoast].dGetBreakingWaveHeight(nThisCoastPoint); // This could be DBL_NODATA
1955 double const dThisBreakingWaveAngle = m_VCoast[nCoast].dGetBreakingWaveAngle(nThisCoastPoint);
1956 double const dThisBreakingDepth = m_VCoast[nCoast].dGetDepthOfBreaking(nThisCoastPoint);
1957 double const dThisWaveSetupSurge = m_VCoast[nCoast].dGetWaveSetupSurge(nThisCoastPoint);
1958 // double dThisStormSurge = m_VCoast[nCoast].dGetStormSurge(nThisCoastPoint);
1959 double const dThisRunUp = m_VCoast[nCoast].dGetRunUp(nThisCoastPoint);
1960
1961 // Get the next profile along the coast, in the down-coast direction. If this next profile has a problem, go to the one after that, etc
1962 CGeomProfile const *pTmpProfile = pProfile;
1963 CGeomProfile *pNextProfile;
1964
1965 while (true)
1966 {
1967 pNextProfile = pTmpProfile->pGetDownCoastAdjacentProfile();
1968
1969 if (pNextProfile->bOKIncStartAndEndOfCoast())
1970 {
1971 // The next profile is OK
1972 break;
1973 }
1974
1975 // The next profile is not OK, so prepare to the one after that
1976 pTmpProfile = pNextProfile;
1977
1978 if (pTmpProfile->bEndOfCoast())
1979 {
1980 // Uh-oh, we've reached the down-coast end of the coast without finding an OK down-coast profile. So give up
1981 return;
1982 }
1983 }
1984
1985 // The next profile is OK
1986 int const nNextCoastPoint = pNextProfile->nGetCoastPoint();
1987 int const nDistBetween = nNextCoastPoint - nThisCoastPoint;
1988
1989 // Safety check
1990 if (nDistBetween <= 0)
1991 // Nothing to do
1992 return;
1993
1994 int const nNextBreakingDist = m_VCoast[nCoast].nGetBreakingDistance(nNextCoastPoint);
1995 double const dNextBreakingWaveHeight = m_VCoast[nCoast].dGetBreakingWaveHeight(nNextCoastPoint); // This could be DBL_NODATA
1996 double const dNextBreakingWaveAngle = m_VCoast[nCoast].dGetBreakingWaveAngle(nNextCoastPoint);
1997 double const dNextBreakingDepth = m_VCoast[nCoast].dGetDepthOfBreaking(nNextCoastPoint);
1998 double const dNextWaveSetupSurge = m_VCoast[nCoast].dGetWaveSetupSurge(nNextCoastPoint);
1999 double const dNextRunUp = m_VCoast[nCoast].dGetRunUp(nNextCoastPoint);
2000
2001 // OK, fill coast point between profiles for setupsurge and runup
2002 for (int n = nThisCoastPoint; n <= nNextCoastPoint; n++)
2003 {
2004 // Fill first wave setup and surge
2005 int const nDist = n - nThisCoastPoint;
2006 double const dThisWeight = (nDistBetween - nDist) / static_cast<double>(nDistBetween);
2007 double const dNextWeight = 1 - dThisWeight;
2008 double dWaveSetupSurge = 0;
2009 double dRunUp = 0;
2010
2011 dWaveSetupSurge = (dThisWeight * dThisWaveSetupSurge) + (dNextWeight * dNextWaveSetupSurge);
2012 m_VCoast[nCoast].SetWaveSetupSurge(n, dWaveSetupSurge);
2013
2014 dRunUp = (dThisWeight * dThisRunUp) + (dNextWeight * dNextRunUp);
2015 m_VCoast[nCoast].SetRunUp(n, dRunUp);
2016 }
2017
2018 // int const nNextProfile = pNextProfile->nGetProfileID();
2019
2020 // If both this profile and the next profile are not in the active zone, then do no more
2021 if ((bFPIsEqual(dThisBreakingWaveHeight, DBL_NODATA, TOLERANCE)) && (bFPIsEqual(dNextBreakingWaveHeight, DBL_NODATA, TOLERANCE)))
2022 {
2023 // if (m_nLogFileDetail >= LOG_FILE_HIGH_DETAIL)
2024 // LogStream << m_ulIter << ": both profile " << pProfile->nGetProfileID() << " at coast point " << nThisCoastPoint << ", and profile " << nNextProfile << " at coast point " << nNextCoastPoint << ", are not in the active zone" << endl;
2025
2026 // Set the breaking wave height, breaking wave angle, and depth of breaking to DBL_NODATA
2027 for (int n = nThisCoastPoint; n < nNextCoastPoint; n++)
2028 {
2029 m_VCoast[nCoast].SetBreakingWaveHeight(nThisCoastPoint, DBL_NODATA);
2030 m_VCoast[nCoast].SetBreakingWaveAngle(nThisCoastPoint, DBL_NODATA);
2031 m_VCoast[nCoast].SetDepthOfBreaking(nThisCoastPoint, DBL_NODATA);
2032 m_VCoast[nCoast].SetBreakingDistance(nThisCoastPoint, DBL_NODATA);
2033 }
2034
2035 return;
2036 }
2037
2038 // OK, at least one of the two profiles is in the active zone
2039 if (bFPIsEqual(dThisBreakingWaveHeight, DBL_NODATA, TOLERANCE))
2040 {
2041 // The next profile must be in the active zone, so use values from the next profile
2042 for (int n = nThisCoastPoint; n < nNextCoastPoint; n++)
2043 {
2044 // Set the breaking wave height, breaking wave angle, and depth of breaking for this coast point TODO 056 This assert sometimes fails: why?
2045 // assert(dNextBreakingWaveHeight >= 0);
2046 m_VCoast[nCoast].SetBreakingWaveHeight(n, dNextBreakingWaveHeight);
2047 m_VCoast[nCoast].SetBreakingWaveAngle(n, dNextBreakingWaveAngle);
2048 m_VCoast[nCoast].SetDepthOfBreaking(n, dNextBreakingDepth);
2049 m_VCoast[nCoast].SetBreakingDistance(n, nNextBreakingDist);
2050 }
2051
2052 return;
2053 }
2054
2055 if (bFPIsEqual(dNextBreakingWaveHeight, DBL_NODATA, TOLERANCE))
2056 {
2057 // This profile must be in the active zone, so use values from this profile
2058 for (int n = nThisCoastPoint + 1; n <= nNextCoastPoint; n++)
2059 {
2060 // Set the breaking wave height, breaking wave angle, and depth of breaking for this coast point TODO 056 This assert sometimes fails: why?
2061 // assert(dThisBreakingWaveHeight >= 0);
2062 m_VCoast[nCoast].SetBreakingWaveHeight(n, dThisBreakingWaveHeight);
2063 m_VCoast[nCoast].SetBreakingWaveAngle(n, dThisBreakingWaveAngle);
2064 m_VCoast[nCoast].SetDepthOfBreaking(n, dThisBreakingDepth);
2065 m_VCoast[nCoast].SetBreakingDistance(n, nThisBreakingDist);
2066 }
2067
2068 return;
2069 }
2070
2071 // The values for both this profile point and the next profile point are fine, so do a weighted interpolation between this profile and the next profile
2072 for (int n = nThisCoastPoint + 1; n < nNextCoastPoint; n++)
2073 {
2074 int const nDist = n - nThisCoastPoint;
2075
2076 double dBreakingWaveHeight = DBL_NODATA;
2077 double dBreakingWaveAngle = DBL_NODATA;
2078 double dBreakingDepth = DBL_NODATA;
2079 double dBreakingDist = DBL_NODATA;
2080
2081 if ((dNextBreakingDepth > 0) && (dThisBreakingDepth > 0))
2082 {
2083 double const dThisWeight = (nDistBetween - nDist) / static_cast<double>(nDistBetween);
2084 double const dNextWeight = 1 - dThisWeight;
2085
2086 dBreakingWaveHeight = (dThisWeight * dThisBreakingWaveHeight) + (dNextWeight * dNextBreakingWaveHeight),
2087 dBreakingWaveAngle = (dThisWeight * dThisBreakingWaveAngle) + (dNextWeight * dNextBreakingWaveAngle),
2088 dBreakingDepth = (dThisWeight * dThisBreakingDepth) + (dNextWeight * dNextBreakingDepth),
2089 dBreakingDist = (dThisWeight * nThisBreakingDist) + (dNextWeight * nNextBreakingDist);
2090 }
2091
2092 else if (dNextBreakingDepth > 0)
2093 {
2094 dBreakingWaveHeight = dNextBreakingWaveHeight,
2095 dBreakingWaveAngle = dNextBreakingWaveAngle,
2096 dBreakingDepth = dNextBreakingDepth,
2097 dBreakingDist = nNextBreakingDist;
2098 }
2099
2100 else if (dThisBreakingDepth > 0)
2101 {
2102 dBreakingWaveHeight = dThisBreakingWaveHeight,
2103 dBreakingWaveAngle = dThisBreakingWaveAngle,
2104 dBreakingDepth = dThisBreakingDepth,
2105 dBreakingDist = nThisBreakingDist;
2106 }
2107
2108 // Set the breaking wave height, breaking wave angle, and depth of breaking for this coast point TODO 056 This assert sometimes fails: why?
2109 // assert(dBreakingWaveHeight >= 0);
2110 m_VCoast[nCoast].SetBreakingWaveHeight(n, dBreakingWaveHeight);
2111 m_VCoast[nCoast].SetBreakingWaveAngle(n, dBreakingWaveAngle);
2112 m_VCoast[nCoast].SetDepthOfBreaking(n, dBreakingDepth);
2113 m_VCoast[nCoast].SetBreakingDistance(n, nRound(dBreakingDist));
2114 }
2115}
2116
2117//===============================================================================================================================
2119//===============================================================================================================================
2121{
2122 int const nCoastPoints = m_VCoast[nCoast].nGetCoastlineSize();
2123
2124 // Initialize all vectors pairs (x,y) for each variable
2125 vector<int> nVCoastWaveHeightX;
2126 vector<double> dVCoastWaveHeightY;
2127
2128 // Search all coast points for non NaN values and store them into temporary variables for later interporlation
2129 for (int n = 0; n < nCoastPoints; n++)
2130 {
2131 double const dCoastWaveHeight = m_VCoast[nCoast].dGetCoastWaveHeight(n);
2132
2133 if (! bFPIsEqual(dCoastWaveHeight, DBL_NODATA, TOLERANCE))
2134 {
2135 nVCoastWaveHeightX.push_back(n);
2136 dVCoastWaveHeightY.push_back(dCoastWaveHeight);
2137 }
2138 }
2139
2140 // Interpolate all coast points. Check first that x,y have more than 3 points and are both of equal size
2141 if ((nVCoastWaveHeightX.size() >= 3) && (nVCoastWaveHeightX.size() == dVCoastWaveHeightY.size()))
2142 {
2143 for (int n = 0; n < nCoastPoints; n++)
2144 {
2145 double const dInterpCoastWaveHeight = dGetInterpolatedValue(&nVCoastWaveHeightX, &dVCoastWaveHeightY, n, false);
2146 m_VCoast[nCoast].SetCoastWaveHeight(n, dInterpCoastWaveHeight);
2147 }
2148 }
2149
2150 else
2151 {
2152 for (int n = 0; n < nCoastPoints; n++)
2153 {
2154 m_VCoast[nCoast].SetCoastWaveHeight(n, 0);
2155 }
2156 }
2157}
2158
2159//===============================================================================================================================
2161//===============================================================================================================================
2163{
2164 int const nCoastSize = m_VCoast[nCoast].nGetCoastlineSize();
2165
2166 for (int nCoastPoint = 0; nCoastPoint < nCoastSize; nCoastPoint++)
2167 {
2168 double dXDiff;
2169 double dYDiff;
2170
2171 if (nCoastPoint == 0)
2172 {
2173 // For the point at the start of the coastline: use the straight line from 'this' point to the next point
2174 CGeom2DPoint const PtThis = *m_VCoast[nCoast].pPtGetCoastlinePointExtCRS(nCoastPoint); // In external CRS
2175 CGeom2DPoint const PtAfter = *m_VCoast[nCoast].pPtGetCoastlinePointExtCRS(nCoastPoint + 1); // In external CRS
2176
2177 dXDiff = PtAfter.dGetX() - PtThis.dGetX();
2178 dYDiff = PtAfter.dGetY() - PtThis.dGetY();
2179 }
2180
2181 else if (nCoastPoint == nCoastSize - 1)
2182 {
2183 // For the point at the end of the coastline: use the straight line from the point before to 'this' point
2184 CGeom2DPoint const PtBefore = *m_VCoast[nCoast].pPtGetCoastlinePointExtCRS(nCoastPoint - 1); // In external CRS
2185 CGeom2DPoint const PtThis = *m_VCoast[nCoast].pPtGetCoastlinePointExtCRS(nCoastPoint); // In external CRS
2186
2187 dXDiff = PtThis.dGetX() - PtBefore.dGetX();
2188 dYDiff = PtThis.dGetY() - PtBefore.dGetY();
2189 }
2190
2191 else
2192 {
2193 // For coastline points not at the start or end of the coast: start with a straight line which links the coastline points before and after 'this' coastline point
2194 CGeom2DPoint const PtBefore = *m_VCoast[nCoast].pPtGetCoastlinePointExtCRS(nCoastPoint - 1); // In external CRS
2195 CGeom2DPoint const PtAfter = *m_VCoast[nCoast].pPtGetCoastlinePointExtCRS(nCoastPoint + 1); // In external CRS
2196
2197 dXDiff = PtAfter.dGetX() - PtBefore.dGetX();
2198 dYDiff = PtAfter.dGetY() - PtBefore.dGetY();
2199 }
2200
2201 // Calculate angle between line and north point, measured clockwise (the azimuth)
2202 if (bFPIsEqual(dYDiff, 0.0, TOLERANCE))
2203 {
2204 // The linking line runs either W-E or E-W
2205 if (dXDiff > 0)
2206 // It runs W to E
2207 m_VCoast[nCoast].SetFluxOrientation(nCoastPoint, 90);
2208
2209 else
2210 // It runs E to W
2211 m_VCoast[nCoast].SetFluxOrientation(nCoastPoint, 270);
2212 }
2213
2214 else if (bFPIsEqual(dXDiff, 0.0, TOLERANCE))
2215 {
2216 // The linking line runs N-S or S-N
2217 if (dYDiff > 0)
2218 // It runs S to N
2219 m_VCoast[nCoast].SetFluxOrientation(nCoastPoint, 0);
2220
2221 else
2222 // It runs N to S
2223 m_VCoast[nCoast].SetFluxOrientation(nCoastPoint, 180);
2224 }
2225
2226 else
2227 {
2228 // The linking line runs neither W-E nor N-S so we have to work a bit harder to find the angle between it and the azimuth
2229 double dAzimuthAngle;
2230
2231 if (dXDiff > 0)
2232 dAzimuthAngle = (180 / PI) * (PI * 0.5 - atan(dYDiff / dXDiff));
2233
2234 else
2235 dAzimuthAngle = (180 / PI) * (PI * 1.5 - atan(dYDiff / dXDiff));
2236
2237 m_VCoast[nCoast].SetFluxOrientation(nCoastPoint, dAzimuthAngle);
2238 }
2239
2240 // LogStream << m_ulIter << ": nCoastPoint = " << nCoastPoint << " FluxOrientation = " << m_VCoast[nCoast].dGetFluxOrientation(nCoastPoint) << endl;
2241 }
2242}
2243
2244//===============================================================================================================================
2246//===============================================================================================================================
2248{
2249 vector<int> VnPolygonD50Count(m_nNumPolygonGlobal + 1, 0); // TODO 044
2250 vector<double> VdPolygonD50(m_nNumPolygonGlobal + 1, 0); // TODO 044
2251
2252 for (int nX = 0; nX < m_nXGridSize; nX++)
2253 {
2254 for (int nY = 0; nY < m_nYGridSize; nY++)
2255 {
2256 if (m_pRasterGrid->m_Cell[nX][nY].bIsInContiguousSea())
2257 {
2258 // This is a sea cell, first get polygon ID
2259 int const nPolyID = m_pRasterGrid->m_Cell[nX][nY].nGetPolygonID();
2260
2261 // Is it in the active zone?
2262 bool const bActive = m_pRasterGrid->m_Cell[nX][nY].bIsInActiveZone();
2263
2264 if (bActive)
2265 {
2266 // It is in the active zone. Does it have unconsolidated sediment on it? Test this using the UnconD50 value: if dGetUnconsD50() returns DBL_NODATA, there is no unconsolidated sediment
2267 double const dTmpd50 = m_pRasterGrid->m_Cell[nX][nY].dGetUnconsD50();
2268
2269 if (! bFPIsEqual(dTmpd50, DBL_NODATA, TOLERANCE))
2270 {
2271 // It does have unconsolidated sediment, so which polygon is this cell in?
2272 if (nPolyID != INT_NODATA)
2273 {
2274 VnPolygonD50Count[nPolyID]++;
2275 VdPolygonD50[nPolyID] += dTmpd50;
2276 }
2277 }
2278 }
2279
2280 // // Now fill in wave calc holes
2281 // if (m_pRasterGrid->m_Cell[nX][nY].dGetWaveHeight() == 0)
2282 // {
2283 // if (nID == INT_NODATA)
2284 // m_pRasterGrid->m_Cell[nX][nY].SetWaveHeight(m_dAllCellsDeepWaterWaveHeight);
2285 // }
2286 //
2287 // if (m_pRasterGrid->m_Cell[nX][nY].dGetWaveAngle() == 0)
2288 // {
2289 // if (nID == INT_NODATA)
2290 // m_pRasterGrid->m_Cell[nX][nY].SetWaveAngle(m_dAllCellsDeepWaterWaveAngle);
2291 // }
2292
2293 // Next look at the cell's N-S and W-E neighbours
2294 int nXTmp;
2295 int nYTmp;
2296 int nActive = 0;
2297 int nShadow = 0;
2298 int nShadowNum = 0;
2299 int nDownDrift = 0;
2300 int nDownDriftNum = 0;
2301 int nCoast = 0;
2302 int nRead = 0;
2303 double dWaveHeight = 0;
2304 double dWaveAngle = 0;
2305
2306 // North
2307 nXTmp = nX;
2308 nYTmp = nY - 1;
2309
2310 if (bIsWithinValidGrid(nXTmp, nYTmp))
2311 {
2312 if (m_pRasterGrid->m_Cell[nXTmp][nYTmp].bIsInContiguousSea())
2313 {
2314 nRead++;
2315 dWaveHeight += m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveHeight();
2316 dWaveAngle += (m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveAngle());
2317
2318 if (m_pRasterGrid->m_Cell[nXTmp][nYTmp].bIsInActiveZone())
2319 nActive++;
2320
2321 int nTmp = m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetShadowZoneNumber();
2322
2323 if (nTmp != 0)
2324 {
2325 nShadow++;
2326 nShadowNum = nTmp;
2327 }
2328
2329 nTmp = m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetDownDriftZoneNumber();
2330
2331 if (nTmp > 0)
2332 {
2333 nDownDrift++;
2334 nDownDriftNum = nTmp;
2335 }
2336 }
2337
2338 else
2339 nCoast++;
2340 }
2341
2342 // East
2343 nXTmp = nX + 1;
2344 nYTmp = nY;
2345
2346 if (bIsWithinValidGrid(nXTmp, nYTmp))
2347 {
2348 if (m_pRasterGrid->m_Cell[nXTmp][nYTmp].bIsInContiguousSea())
2349 {
2350 nRead++;
2351 dWaveHeight += m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveHeight();
2352 dWaveAngle += (m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveAngle());
2353
2354 if (m_pRasterGrid->m_Cell[nXTmp][nYTmp].bIsInActiveZone())
2355 nActive++;
2356
2357 int nTmp = m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetShadowZoneNumber();
2358
2359 if (nTmp != 0)
2360 {
2361 nShadow++;
2362 nShadowNum = nTmp;
2363 }
2364
2365 nTmp = m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetDownDriftZoneNumber();
2366
2367 if (nTmp > 0)
2368 {
2369 nDownDrift++;
2370 nDownDriftNum = nTmp;
2371 }
2372 }
2373
2374 else
2375 nCoast++;
2376 }
2377
2378 // South
2379 nXTmp = nX;
2380 nYTmp = nY + 1;
2381
2382 if (bIsWithinValidGrid(nXTmp, nYTmp))
2383 {
2384 if (m_pRasterGrid->m_Cell[nXTmp][nYTmp].bIsInContiguousSea())
2385 {
2386 nRead++;
2387 dWaveHeight += m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveHeight();
2388 dWaveAngle += (m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveAngle());
2389
2390 if (m_pRasterGrid->m_Cell[nXTmp][nYTmp].bIsInActiveZone())
2391 nActive++;
2392
2393 int nTmp = m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetShadowZoneNumber();
2394
2395 if (nTmp != 0)
2396 {
2397 nShadow++;
2398 nShadowNum = nTmp;
2399 }
2400
2401 nTmp = m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetDownDriftZoneNumber();
2402
2403 if (nTmp > 0)
2404 {
2405 nDownDrift++;
2406 nDownDriftNum = nTmp;
2407 }
2408 }
2409
2410 else
2411 nCoast++;
2412 }
2413
2414 // West
2415 nXTmp = nX - 1;
2416 nYTmp = nY;
2417
2418 if (bIsWithinValidGrid(nXTmp, nYTmp))
2419 {
2420 if (m_pRasterGrid->m_Cell[nXTmp][nYTmp].bIsInContiguousSea())
2421 {
2422 nRead++;
2423 dWaveHeight += m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveHeight();
2424 dWaveAngle += (m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveAngle());
2425
2426 if (m_pRasterGrid->m_Cell[nXTmp][nYTmp].bIsInActiveZone())
2427 nActive++;
2428
2429 int nTmp = m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetShadowZoneNumber();
2430
2431 if (nTmp != 0)
2432 {
2433 nShadow++;
2434 nShadowNum = nTmp;
2435 }
2436
2437 nTmp = m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetDownDriftZoneNumber();
2438
2439 if (nTmp > 0)
2440 {
2441 nDownDrift++;
2442 nDownDriftNum = nTmp;
2443 }
2444 }
2445
2446 else
2447 nCoast++;
2448 }
2449
2450 if (nRead > 0)
2451 {
2452 // Calculate the average of neighbours
2453 dWaveHeight /= nRead;
2454 dWaveAngle /= nRead;
2455 dWaveAngle = dKeepWithin360(dWaveAngle);
2456
2457 // If this sea cell has four active-zone neighbours, then it must also be in the active zone: give it wave height and orientation which is the average of its neighbours
2458 if (nActive == 4)
2459 {
2460 m_pRasterGrid->m_Cell[nX][nY].SetInActiveZone(true);
2461 m_pRasterGrid->m_Cell[nX][nY].SetWaveHeight(dWaveHeight);
2462 m_pRasterGrid->m_Cell[nX][nY].SetWaveAngle(dWaveAngle);
2463 }
2464
2465 // If this sea cell has a wave height which is the same as its deep-water wave height, but its neighbours have a different average wave height, then give it the average of its neighbours
2466 double const dDeepWaterWaveHeight = m_pRasterGrid->m_Cell[nX][nY].dGetCellDeepWaterWaveHeight();
2467
2468 if ((bFPIsEqual(m_pRasterGrid->m_Cell[nX][nY].dGetWaveHeight(), dDeepWaterWaveHeight, TOLERANCE)) && (! bFPIsEqual(dDeepWaterWaveHeight, dWaveHeight, TOLERANCE)))
2469 {
2470 m_pRasterGrid->m_Cell[nX][nY].SetWaveHeight(dWaveHeight);
2471 }
2472
2473 // If this sea cell has a wave orientation which is the same as its deep-water wave orientation, but its neighbours have a different average wave orientation, then give it the average of its neighbours
2474 double const dDeepWaterWaveAngle = m_pRasterGrid->m_Cell[nX][nY].dGetCellDeepWaterWaveAngle();
2475
2476 if ((bFPIsEqual(m_pRasterGrid->m_Cell[nX][nY].dGetWaveAngle(), dDeepWaterWaveAngle, TOLERANCE)) && (! bFPIsEqual(dDeepWaterWaveAngle, dWaveAngle, TOLERANCE)))
2477 {
2478 m_pRasterGrid->m_Cell[nX][nY].SetWaveAngle(dWaveAngle);
2479 }
2480
2481 // Is this sea cell is not already marked as in a shadow zone (note could be marked as in a shadow zone but not yet processed: a -ve number)?
2482 int const nShadowZoneCode = m_pRasterGrid->m_Cell[nX][nY].nGetShadowZoneNumber();
2483
2484 if (nShadowZoneCode <= 0)
2485 {
2486 // If the cell has four neighbours which are all in a shadow zone, or four neighbours some of which are shadow zone and the remainder downdrift zone, or four neighbours some of which are shadow zone and the remainder coast; then it should also be in the shadow zone: give it the average of its neighbours
2487 if ((nShadow == 4) || (nShadow + nDownDrift == 4) || (nShadow + nCoast == 4))
2488 {
2489 m_pRasterGrid->m_Cell[nX][nY].SetShadowZoneNumber(nShadowNum);
2490 m_pRasterGrid->m_Cell[nX][nY].SetWaveHeight(dWaveHeight);
2491 m_pRasterGrid->m_Cell[nX][nY].SetWaveAngle(dWaveAngle);
2492 }
2493 }
2494
2495 // If this sea cell is not marked as in a downdrift zone but has four neighbours which are in a downdrift zone, then it should also be in the downdrift zone: give it the average of its neighbours
2496 int const nDownDriftZoneCode = m_pRasterGrid->m_Cell[nX][nY].nGetDownDriftZoneNumber();
2497
2498 if ((nDownDriftZoneCode == 0) && (nDownDrift == 4))
2499 {
2500 m_pRasterGrid->m_Cell[nX][nY].SetDownDriftZoneNumber(nDownDriftNum);
2501 m_pRasterGrid->m_Cell[nX][nY].SetWaveHeight(dWaveHeight);
2502 m_pRasterGrid->m_Cell[nX][nY].SetWaveAngle(dWaveAngle);
2503 }
2504 }
2505 }
2506 }
2507 }
2508
2509 // Calculate the average d50 for every polygon TODO 044
2510 for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
2511 {
2512 for (int nPoly = 0; nPoly < m_VCoast[nCoast].nGetNumPolygons(); nPoly++)
2513 {
2514 CGeomCoastPolygon* pPolygon = m_VCoast[nCoast].pGetPolygon(nPoly);
2515 int const nPolyID = pPolygon->nGetPolygonCoastID();
2516
2517 if (VnPolygonD50Count[nPolyID] > 0)
2518 VdPolygonD50[nPolyID] /= VnPolygonD50Count[nPolyID];
2519
2520 pPolygon->SetAvgUnconsD50(VdPolygonD50[nPolyID]);
2521 }
2522 }
2523}
Contains CGeom2DPoint definitions.
int nGetY(void) const
Returns the CGeom2DIPoint object's integer Y coordinate.
Definition 2di_point.cpp:55
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
double dGetY(void) const
Returns the CGeom2DPoint object's double Y coordinate.
Definition 2d_point.cpp:53
double dGetX(void) const
Returns the CGeom2DPoint object's double X coordinate.
Definition 2d_point.cpp:47
Geometry class used for coast polygon objects.
int nGetPolygonCoastID(void) const
Get the coast ID, this is the same as the down-coast sequence of polygons.
void SetAvgUnconsD50(double const)
Set the average d50 for unconsolidated sediment on this polygon.
Geometry class used to represent coast profile objects.
Definition profile.h:37
double dGetProfileDeepWaterWaveAngle(void) const
Returns the deep-water wave orientation for this profile.
Definition profile.cpp:585
double dGetProfileDeepWaterWaveHeight(void) const
Returns the deep-water wave height for this profile.
Definition profile.cpp:573
int nGetProfileID(void) const
Returns the profile's this-coast ID.
Definition profile.cpp:80
int nGetCoastPoint(void) const
Returns the coast point at which the profile starts.
Definition profile.cpp:86
CGeomProfile * pGetDownCoastAdjacentProfile(void) const
Returns the down-coast adjacent profile, returns NULL if there is no down-coast adjacent profile.
Definition profile.cpp:491
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
void SetCShoreProblem(bool const)
Sets a switch to indicate whether this profile has a CShore problem.
Definition profile.cpp:134
bool bOKIncStartAndEndOfCoast(void) const
Returns true if this is a problem-free profile (however it could be a start-of-coast or an end-of-coa...
Definition profile.cpp:265
bool bStartOfCoast(void) const
Returns the switch to indicate whether this is a start-of-coast profile.
Definition profile.cpp:116
vector< CGeom2DIPoint > * pPtiVGetCellsInProfile(void)
Returns all cells (grid CRS) in the profile.
Definition profile.cpp:515
bool bEndOfCoast(void) const
Returns the switch to indicate whether this is an end-of-coast profile.
Definition profile.cpp:128
double dGetProfileDeepWaterWavePeriod(void) const
Returns the deep-water wave period for this profile.
Definition profile.cpp:597
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
double m_dAllCellsDeepWaterWaveHeight
Deep water wave height (m) for all sea cells.
Definition simulation.h:754
void InterpolateCShoreOutput(vector< double > const *, int const, int const, vector< double > *, vector< double > *, vector< double > *, vector< double > *, vector< double > *, vector< double > *, vector< double > *, vector< double > *, vector< double > *)
double m_dG
Gravitational acceleration (m**2/sec)
Definition simulation.h:811
double m_dThisIterSWL
The still water level for this timestep (this includes tidal changes and any long-term SWL change)
Definition simulation.h:712
CGeomRasterGrid * m_pRasterGrid
Pointer to the raster grid object.
int nCalcWavePropertiesOnProfile(int const, int const, CGeomProfile *, vector< double > *, vector< double > *, vector< double > *, vector< double > *, vector< bool > *)
Calculates wave properties along a coastline-normal profile using either the COVE linear wave theory ...
int m_nXGridSize
The size of the grid in the x direction.
Definition simulation.h:457
double m_dL_0
Deep water wave length (m)
Definition simulation.h:745
ofstream LogStream
vector< CRWCoast > m_VCoast
The coastline objects.
bool m_bSingleDeepWaterWaveValues
Do we have just a point source for (i.e. only a single measurement of) deep water wave values.
Definition simulation.h:385
int m_nYGridSize
The size of the grid in the y direction.
Definition simulation.h:460
int nDoAllShadowZones(void)
Finds wave shadow zones and modifies waves in and near them. Note that where up-coast and down-coast ...
static double dKeepWithin360(double const)
Constrains the supplied angle to be within 0 and 360 degrees.
void CalcCoastTangents(int const)
Calculates tangents to a coastline: the tangent is assumed to be the orientation of energy/sediment f...
double m_dWaveDepthRatioForWaveCalcs
Start depth for wave calculations.
Definition simulation.h:748
int nCreateCShoreInfile(int const, int const, int const, int const, int const, int const, int const, int const, int const, int const, int const, int const, int const, double const, double const, double const, double const, double const, double const, double const, double const, vector< double > const *, vector< double > const *, vector< double > const *)
void InterpolateWavePropertiesBetweenProfiles(int const, int const)
Interpolates wave properties from profiles to the coastline points between two profiles....
int nInterpolateWavesToPolygonCells(vector< double > const *, vector< double > const *, vector< double > const *, vector< double > const *)
Interpolates wave properties from all profiles to all within-polygon sea cells, using GDALGridCreate(...
int m_nRunUpEquation
The run-up equation used TODO 007.
Definition simulation.h:577
double m_dBreakingWaveHeightDepthRatio
Breaking wave height-to-depth ratio.
Definition simulation.h:751
int nReadCShoreOutput(int const, string const *, int const, int const, vector< double > const *, vector< double > *)
double dGridCentroidYToExtCRSY(int const) const
Given the integer Y-axis ordinate of a cell in the raster grid CRS, returns the external CRS Y-axis o...
Definition gis_utils.cpp:78
void CalcD50AndFillWaveCalcHoles(void)
Calculates an average d50 for each polygon. Also fills in 'holes' in active zone and wave calcs i....
string m_strCMEDir
The CME folder.
vector< double > VdInterpolateCShoreProfileOutput(vector< double > const *, vector< double > const *, vector< double > const *)
Returns a linearly interpolated vector of doubles, to make CShore profile output compatible with CME....
int m_nWavePropagationModel
The wave propagation model used. Possible values are WAVE_MODEL_CSHORE and WAVE_MODEL_COVE.
Definition simulation.h:550
double m_dAllCellsDeepWaterWaveAngle
Deep water wave angle for all sea cells.
Definition simulation.h:757
int nSetAllCoastpointDeepWaterWaveValues(void)
Give every coast point a value for deep water wave height and direction TODO 005 This may not be real...
double m_dTimeStep
The length of an iteration (a time step) in hours.
Definition simulation.h:676
int nGetThisProfileElevationsForCShore(int const, CGeomProfile *, int const, vector< double > *, vector< double > *, vector< double > *)
Get profile horizontal distance and bottom elevation vectors in CShore units.
bool bIsWithinValidGrid(int const, int const) const
void ModifyBreakingWavePropertiesWithinShadowZoneToCoastline(int const, int const)
Modifies the wave breaking properties at coastline points of profiles within the shadow zone.
double m_dC_0
Deep water wave speed (m/s)
Definition simulation.h:742
unsigned long m_ulIter
The number of the current iteration (time step)
Definition simulation.h:598
double dGridCentroidXToExtCRSX(int const) const
Definition gis_utils.cpp:68
int m_nNumPolygonGlobal
Number of global (all coasts) polygons.
Definition simulation.h:523
double dGetInterpolatedValue(vector< double > const *, vector< double > const *, double, bool)
int nDoAllPropagateWaves(void)
Simulates wave propagation along all coastline-normal profiles, on all coasts.
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
void InterpolateWaveHeightToCoastPoints(int const)
Linearly interpolates wave properties from profiles to the coastline cells between two profiles.
static vector< string > * VstrSplit(string const *, char const, vector< string > *)
From http://stackoverflow.com/questions/236129/split-a-string-in-c They implement (approximately) Pyt...
Definition utils.cpp:2466
static double dCalcWaveAngleToCoastNormal(double const, double const, int const)
Calculates the angle between the wave direction and a normal to the coastline tangent....
This file contains global definitions for CoastalME.
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
int const WAVE_MODEL_COVE
Definition cme.h:791
double const TOLERANCE
Definition cme.h:823
int const RTN_ERR_NO_TOP_LAYER
Definition cme.h:738
string const CSHORE_DIR
Definition cme.h:899
int const RTN_ERR_CSHORE_FILE_INPUT
Definition cme.h:746
string const ERR
Definition cme.h:902
double const CSHORE_FRICTION_FACTOR
Definition cme.h:820
int const WAVE_MODEL_CSHORE
Definition cme.h:792
int const LOG_FILE_MIDDLE_DETAIL
Definition cme.h:491
double const PI
Definition cme.h:805
int const RTN_ERR_CSHORE_ERROR
Definition cme.h:751
bool bFPIsEqual(const T d1, const T d2, const T dEpsilon)
Definition cme.h:1303
double const WALKDEN_HALL_PARAM_2
Definition cme.h:813
T tMax(T a, T b)
Definition cme.h:1252
int const RTN_ERR_READING_CSHORE_FILE_OUTPUT
Definition cme.h:747
double const CSHORE_SURGE_LEVEL
Definition cme.h:821
string const WARN
Definition cme.h:903
int const RTN_OK
Definition cme.h:694
bool const SAVE_CSHORE_OUTPUT
Definition cme.h:459
double const DBL_NODATA
Definition cme.h:834
int const CSHOREARRAYOUTSIZE
The size of the arrays output by CShore. If you change this, then you must also set the same value on...
Definition cme.h:472
double const WALKDEN_HALL_PARAM_1
Definition cme.h:812
T tAbs(T a)
Definition cme.h:1277
int const RTN_ERR_CSHORE_EMPTY_PROFILE
Definition cme.h:745
string const CSHORE_INFILE
Definition cme.h:900
int const LEFT_HANDED
Definition cme.h:512
char const SPACE
Definition cme.h:453
Contains CRWCoast definitions.
void CShore(int *)
void CShoreWrapper(int const *, int const *, int const *, int const *, int const *, int const *, int const *, int const *, int const *, int const *, int const *, double const *, double const *, double[], double[], double[], double[], double[], double[], int const *, double[], double[], double[], int *, int *, double[], double[], double[], double[], double[])
Contains CSimulation definitions.
bool bIsStringValidInt(string &str)
Checks to see if a string can be read as a valid integer, from https://stackoverflow....
int nRound(double const d)
Version of the above that returns an int.