45using std::reverse_copy;
59 for (
int nCoast = 0; nCoast < static_cast<int>(
m_VCoast.size()); nCoast++)
61 int nDistFromPrevProfile = 0;
62 int nDistToNextProfile = 0;
64 double dPrevProfileDeepWaterWaveHeight = 0;
65 double dPrevProfileDeepWaterWaveAngle = 0;
66 double dPrevProfileDeepWaterWavePeriod = 0;
67 double dNextProfileDeepWaterWaveHeight = 0;
68 double dNextProfileDeepWaterWaveAngle = 0;
69 double dNextProfileDeepWaterWavePeriod = 0;
72 for (
int nPoint = 0; nPoint <
m_VCoast[nCoast].nGetCoastlineSize(); nPoint++)
75 if (
m_VCoast[nCoast].bIsProfileAtCoastPoint(nPoint))
85 m_VCoast[nCoast].SetCoastDeepWaterWaveHeight(nPoint, dThisDeepWaterWaveHeight);
86 m_VCoast[nCoast].SetCoastDeepWaterWaveAngle(nPoint, dThisDeepWaterWaveAngle);
87 m_VCoast[nCoast].SetCoastDeepWaterWavePeriod(nPoint, dThisDeepWaterWavePeriod);
90 nDistFromPrevProfile = 0;
91 dPrevProfileDeepWaterWaveHeight = dThisDeepWaterWaveHeight;
92 dPrevProfileDeepWaterWaveAngle = dThisDeepWaterWaveAngle;
93 dPrevProfileDeepWaterWavePeriod = dThisDeepWaterWavePeriod;
98 if (pNextProfile == NULL)
106 dDist = nDistToNextProfile;
119 nDistFromPrevProfile++;
120 nDistToNextProfile--;
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);
128 m_VCoast[nCoast].SetCoastDeepWaterWaveHeight(nPoint, dThisDeepWaterWaveHeight);
129 m_VCoast[nCoast].SetCoastDeepWaterWaveAngle(nPoint, dThisDeepWaterWaveAngle);
130 m_VCoast[nCoast].SetCoastDeepWaterWavePeriod(nPoint, dThisDeepWaterWavePeriod);
146 vector<bool> VbBreakingAll;
148 vector<double> VdXAll;
149 vector<double> VdYAll;
150 vector<double> VdHeightXAll;
151 vector<double> VdHeightYAll;
154 bool bSomeNonStartOrEndOfCoastProfiles =
false;
156 for (
int nCoast = 0; nCoast < static_cast<int>(
m_VCoast.size()); nCoast++)
158 int const nCoastSize =
m_VCoast[nCoast].nGetCoastlineSize();
159 int const nNumProfiles =
m_VCoast[nCoast].nGetNumProfiles();
161 static bool bDownCoast =
true;
164 for (
int nn = 0; nn < nNumProfiles; nn++)
166 vector<bool> VbBreaking;
169 vector<double> VdHeightX;
170 vector<double> VdHeightY;
175 pProfile =
m_VCoast[nCoast].pGetProfileWithDownCoastSeq(nn);
178 pProfile =
m_VCoast[nCoast].pGetProfileWithUpCoastSeq(nn);
201 if (VbBreaking.empty())
208 bSomeNonStartOrEndOfCoastProfiles =
true;
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());
227 bDownCoast = ! bDownCoast;
231 if (! bSomeNonStartOrEndOfCoastProfiles)
233 LogStream <<
m_ulIter <<
": waves are on-shore only, for start and/or end of coast profiles" << endl;
239 double dDeepWaterWaveX;
240 double dDeepWaterWaveY;
253 int const nPolyID =
m_pRasterGrid->m_Cell[nX][0].nGetPolygonID();
258 VdXAll.push_back(nX);
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);
268 VdHeightXAll.push_back(dDeepWaterWaveX);
269 VdHeightYAll.push_back(dDeepWaterWaveY);
280 VdXAll.push_back(nX);
290 VdHeightXAll.push_back(dDeepWaterWaveX);
291 VdHeightYAll.push_back(dDeepWaterWaveY);
300 int const nPolyID =
m_pRasterGrid->m_Cell[0][nY].nGetPolygonID();
306 VdYAll.push_back(nY);
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);
315 VdHeightXAll.push_back(dDeepWaterWaveX);
316 VdHeightYAll.push_back(dDeepWaterWaveY);
328 VdYAll.push_back(nY);
337 VdHeightXAll.push_back(dDeepWaterWaveX);
338 VdHeightYAll.push_back(dDeepWaterWaveY);
353 if (VbBreakingAll.empty())
504 for (
int nCoast = 0; nCoast < static_cast<int>(
m_VCoast.size()); nCoast++)
506 int const nNumProfiles =
m_VCoast[nCoast].nGetNumProfiles();
508 for (
int nProfile = 0; nProfile < nNumProfiles; nProfile++)
577 for (
int nCoast = 0; nCoast < static_cast<int>(
m_VCoast.size()); nCoast++)
579 int const nCoastSize =
m_VCoast[nCoast].nGetCoastlineSize();
580 int const nNumProfiles =
m_VCoast[nCoast].nGetNumProfiles();
583 for (
int n = 0; n < nNumProfiles - 1; n++)
590 for (
int nCoastPoint = 0; nCoastPoint < nCoastSize; nCoastPoint++)
593 double const dBreakingWaveHeight =
m_VCoast[nCoast].dGetBreakingWaveHeight(nCoastPoint);
594 double const dCoastPointWavePeriod =
m_VCoast[nCoast].dGetCoastDeepWaterWavePeriod(nCoastPoint);
599 m_VCoast[nCoast].SetBreakingWaveHeight(nCoastPoint, 0);
607 double const dWaveEnergy = dErosiveWaveForce *
m_dTimeStep * 3600;
608 m_VCoast[nCoast].SetWaveEnergyAtBreaking(nCoastPoint, dWaveEnergy);
621 double dWaveToNormalAngle = 0;
625 dWaveToNormalAngle = fmod((dWaveAngle - dCoastAngle + 360), 360) - 90;
629 dWaveToNormalAngle = fmod((dWaveAngle - dCoastAngle + 360), 360) - 270;
631 if ((dWaveToNormalAngle >= 90) || (dWaveToNormalAngle <= -90))
634 return dWaveToNormalAngle;
657 int const nSeaHand =
m_VCoast[nCoast].nGetSeaHandedness();
661 double const dFluxOrientationThis =
m_VCoast[nCoast].dGetFluxOrientation(nCoastPoint);
662 double dFluxOrientationPrev = 0;
663 double dFluxOrientationNext = 0;
665 if (nCoastPoint == 0)
667 dFluxOrientationPrev = dFluxOrientationThis;
668 dFluxOrientationNext =
m_VCoast[nCoast].dGetFluxOrientation(1);
671 else if (nCoastPoint == nCoastSize - 1)
673 dFluxOrientationPrev =
m_VCoast[nCoast].dGetFluxOrientation(nCoastPoint - 2);
674 dFluxOrientationNext = dFluxOrientationThis;
679 dFluxOrientationPrev =
m_VCoast[nCoast].dGetFluxOrientation(nCoastPoint - 1);
680 dFluxOrientationNext =
m_VCoast[nCoast].dGetFluxOrientation(nCoastPoint + 1);
701 double dWaveToNormalAnglePrev;
706 double const dPrevDeepWaterWaveAngle =
m_VCoast[nCoast].dGetCoastDeepWaterWaveAngle(nCoastPoint - 1);
713 dWaveToNormalAnglePrev = dWaveToNormalAngle;
722 double dWaveToNormalAngleNext;
724 if (nCoastPoint < nCoastSize - 1)
727 double const dNextDeepWaterWaveAngle =
m_VCoast[nCoast].dGetCoastDeepWaterWaveAngle(nCoastPoint + 1);
734 dWaveToNormalAngleNext = dWaveToNormalAngle;
745 if (dWaveToNormalAngle > 45)
747 if (dWaveToNormalAnglePrev < 45)
749 dWaveToNormalAngle = 45;
755 dWaveToNormalAngle = dWaveToNormalAnglePrev;
763 if (dWaveToNormalAngle < -45)
765 if (dWaveToNormalAngleNext > -45)
767 dWaveToNormalAngle = -45;
773 dWaveToNormalAngle = dWaveToNormalAngleNext;
784 dWaveToNormalAngle = dFluxOrientationPrev;
792 dWaveToNormalAngle = dFluxOrientationNext;
796 bool bBreaking =
false;
799 int nProfileBreakingDist = 0;
801 double dProfileBreakingWaveHeight =
DBL_NODATA;
802 double dProfileBreakingWaveAngle = 0;
803 double dProfileBreakingDepth = 0;
809 vector<bool> VbWaveIsBreaking(nProfileSize, 0);
811 vector<double> VdWaveHeight(nProfileSize, 0);
812 vector<double> VdWaveSetupSurge(nProfileSize, 0);
814 vector<double>
const VdWaveSetupRunUp(nProfileSize, 0);
815 vector<double> VdWaveDirection(nProfileSize, 0);
820 double const dCShoreTimeStep = 3600;
824 vector<double> VdProfileZ;
825 vector<double> VdProfileDistXY;
826 vector<double> VdProfileFrictionFactor;
843 if (VdProfileDistXY.empty())
852 dWaveToNormalAngle =
tMax(dWaveToNormalAngle, -80.0);
853 dWaveToNormalAngle =
tMin(dWaveToNormalAngle, 80.0);
855 int const nProfileDistXYSize =
static_cast<int>(VdProfileDistXY.size());
856 vector<double> VdFreeSurfaceStd(nProfileDistXYSize, 0);
857 vector<double> VdSinWaveAngleRadians(nProfileDistXYSize, 0);
858 vector<double> VdFractionBreakingWaves(nProfileDistXYSize, 0);
861 int const nILine = 1;
862 int const nIProfl = 0;
863 int const nIPerm = 0;
864 int const nIOver = 0;
865 int const nIWCInt = 0;
866 int const nIRoll = 0;
867 int const nIWind = 0;
868 int const nITide = 0;
870 int const nNWave = 1;
871 int const nNSurge = 1;
872 double const dDX = VdProfileDistXY.back() / (
static_cast<double>(VdProfileDistXY.size() - 1));
873 double const dWaveInitTime = 0;
874 double const dSurgeInitTime = 0;
876#if defined CSHORE_FILE_INOUT || CSHORE_BOTH
885#if defined CSHORE_FILE_INOUT
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);
907 strErr = to_string(
m_ulIter) +
": CShore WARNING 1: negative depth at the first node ";
911 strErr = to_string(
m_ulIter) +
": CShore WARNING 2: negative value at end of landward marching computation ";
915 strErr = to_string(
m_ulIter) +
": CShore WARNING 3: large energy gradients at the first node: small waves with short period at sea boundary ";
919 strErr = to_string(
m_ulIter) +
": CShore WARNING 4: zero energy at the first node ";
923 strErr = to_string(
m_ulIter) +
": CShore WARNING 5: at end of landward marching computation, insufficient water depth ";
927 strErr = to_string(
m_ulIter) +
": CShore WARNING 7: did not reach convergence ";
931 strErr +=
"(coast " + to_string(nCoast) +
" profile " + to_string(pProfile->
nGetProfileID()) +
" profile length " + to_string(nOutSize) +
")\n";
940 string strOSETUP =
"OSETUP";
941 string strOYVELO =
"OYVELO";
942 string strOPARAM =
"OPARAM";
944 nRet =
nReadCShoreOutput(nProfile, &strOSETUP, 4, 4, &VdProfileDistXY, &VdFreeSurfaceStd);
949 nRet =
nReadCShoreOutput(nProfile, &strOYVELO, 4, 2, &VdProfileDistXY, &VdSinWaveAngleRadians);
954 nRet =
nReadCShoreOutput(nProfile, &strOPARAM, 4, 3, &VdProfileDistXY, &VdFractionBreakingWaves);
965 nRet = system(
"./clean.bat");
970 string strCommand =
"./save_CShore_output.sh ";
973 strCommand += to_string(nCoast);
975 strCommand += to_string(nProfile);
977 nRet = system(strCommand.c_str());
983 nRet = system(
"./clean.sh");
997#if defined CSHORE_ARG_INOUT || CSHORE_BOTH
1004 vector<double> VdInitTime = {dWaveInitTime, dCShoreTimeStep};
1005 vector<double> VdTPIn = {dDeepWaterWavePeriod, dDeepWaterWavePeriod};
1006 vector<double> VdHrmsIn = {dProfileDeepWaterWaveHeight, dProfileDeepWaterWaveHeight};
1007 vector<double> VdWangIn = {dWaveToNormalAngle, dWaveToNormalAngle};
1008 vector<double> VdTSurg = {dSurgeInitTime, dCShoreTimeStep};
1009 vector<double> VdSWLin = {dSurgeLevel, dSurgeLevel};
1010 vector<double> VdFPInp = VdProfileFrictionFactor;
1041 &nProfileDistXYSize,
1042 &VdProfileDistXY[0],
1047 &VdXYDistFromCShoreOut[0],
1048 &VdFreeSurfaceStdOut[0],
1049 &VdWaveSetupSurgeOut[0],
1050 &VdSinWaveAngleRadiansOut[0],
1051 &VdFractionBreakingWavesOut[0]);
1058 LogStream <<
m_ulIter <<
": " <<
WARN <<
"for coast " << nCoast <<
" profile " << pProfile->
nGetProfileID() <<
", only " << nOutSize <<
" CShore output rows, abandoning this profile" << endl;
1082 strErr = to_string(
m_ulIter) +
": CShore WARNING 1: negative depth at the first node ";
1086 strErr = to_string(
m_ulIter) +
": CShore WARNING 2: negative value at end of landward marching computation ";
1090 strErr = to_string(
m_ulIter) +
": CShore WARNING 3: large energy gradients at the first node: small waves with short period at sea boundary ";
1094 strErr = to_string(
m_ulIter) +
": CShore WARNING 4: zero energy at the first node ";
1098 strErr = to_string(
m_ulIter) +
": CShore WARNING 5: at end of landward marching computation, insufficient water depth ";
1102 strErr = to_string(
m_ulIter) +
": CShore WARNING 7: did not reach convergence ";
1106 strErr +=
"(coast " + to_string(nCoast) +
" profile " + to_string(pProfile->
nGetProfileID()) +
" profile length " + to_string(nOutSize) +
")\n";
1116 InterpolateCShoreOutput(&VdProfileDistXY, nOutSize, nProfileSize, &VdXYDistFromCShoreOut, &VdFreeSurfaceStdOut, &VdWaveSetupSurgeOut, &VdSinWaveAngleRadiansOut, &VdFractionBreakingWavesOut, &VdFreeSurfaceStd, &VdWaveSetupSurge, &VdSinWaveAngleRadians, &VdFractionBreakingWaves);
1119#if defined CSHORE_BOTH
1124 string strCommand =
"./save_CShore_output.sh ";
1127 strCommand += to_string(nCoast);
1129 strCommand += to_string(nProfile);
1131 nRet = system(strCommand.c_str());
1148 for (
int nProfilePoint = (nProfileSize - 1); nProfilePoint >= 0; nProfilePoint--)
1154 if (nProfilePoint >
static_cast<int>(VdFreeSurfaceStd.size()) - 1)
1158 if (isnan(VdFreeSurfaceStd[nProfilePoint]))
1159 VdFreeSurfaceStd[nProfilePoint] = 0;
1161 VdWaveHeight[nProfilePoint] = sqrt(8) * VdFreeSurfaceStd[nProfilePoint];
1164 if (isnan(VdSinWaveAngleRadians[nProfilePoint]))
1166 VdSinWaveAngleRadians[nProfilePoint] = 0;
1167 VdWaveHeight[nProfilePoint] = 0;
1171 if (VdSinWaveAngleRadians[nProfilePoint] < -1)
1172 VdSinWaveAngleRadians[nProfilePoint] = -1;
1174 if (VdSinWaveAngleRadians[nProfilePoint] > 1)
1175 VdSinWaveAngleRadians[nProfilePoint] = 1;
1177 double const dAlpha = asin(VdSinWaveAngleRadians[nProfilePoint]) * (180 /
PI);
1180 VdWaveDirection[nProfilePoint] =
dKeepWithin360(dAlpha + 90 + dFluxOrientationThis);
1183 VdWaveDirection[nProfilePoint] =
dKeepWithin360(dAlpha + 270 + dFluxOrientationThis);
1186 if (isnan(VdFractionBreakingWaves[nProfilePoint]))
1188 VdFractionBreakingWaves[nProfilePoint] = 0;
1189 VdWaveHeight[nProfilePoint] = 0;
1197 dProfileBreakingWaveHeight = VdWaveHeight[nProfilePoint];
1198 dProfileBreakingWaveAngle = VdWaveDirection[nProfilePoint];
1199 dProfileBreakingDepth =
m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth();
1200 nProfileBreakingDist = nProfilePoint + 1;
1205 VbWaveIsBreaking[nProfilePoint] = bBreaking;
1208 if (dProfileBreakingWaveHeight >= dProfileDeepWaterWaveHeight)
1220 for (
int nProfilePoint = (nProfileSize - 1); nProfilePoint >= 0; nProfilePoint--)
1229 double const dSeaDepth =
m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth();
1231 if (dSeaDepth > dDepthLookupMax)
1234 dProfileWaveHeight = dProfileDeepWaterWaveHeight;
1235 dProfileWaveAngle = dProfileDeepWaterWaveAngle;
1243 double const dL =
m_dL_0 * sqrt(tanh((2 *
PI * dSeaDepth) /
m_dL_0));
1244 double const dC =
m_dC_0 * tanh((2 *
PI * dSeaDepth) / dL);
1245 double const dk = 2 *
PI / dL;
1246 double const dn = ((2 * dSeaDepth * dk) / (sinh(2 * dSeaDepth * dk)) + 1) / 2;
1247 double const dKs = sqrt(
m_dC_0 / (dn * dC * 2));
1248 double const dAlpha = (180 /
PI) * asin((dC /
m_dC_0) * sin((
PI / 180) * dWaveToNormalAngle));
1249 double const dKr = sqrt(cos((
PI / 180) * dWaveToNormalAngle) / cos((
PI / 180) * dAlpha));
1250 dProfileWaveHeight = dProfileDeepWaterWaveHeight * dKs * dKr;
1253 dProfileWaveAngle =
dKeepWithin360(dAlpha + 90 + dFluxOrientationThis);
1256 dProfileWaveAngle =
dKeepWithin360(dAlpha + 270 + dFluxOrientationThis);
1261 dProfileBreakingWaveHeight = dProfileWaveHeight;
1262 dProfileBreakingWaveAngle = dProfileWaveAngle;
1272 dProfileWaveAngle = dProfileBreakingWaveAngle;
1274 dProfileWaveHeight = dProfileBreakingWaveHeight * (nProfilePoint / nProfileBreakingDist);
1276 dProfileBreakingDepth = dSeaDepth;
1277 nProfileBreakingDist = nProfilePoint;
1282 VdWaveDirection[nProfilePoint] = dProfileWaveAngle;
1283 VdWaveHeight[nProfilePoint] = dProfileWaveHeight;
1284 VbWaveIsBreaking[nProfilePoint] = bBreaking;
1289 for (
int nProfilePoint = (nProfileSize - 1); nProfilePoint >= 0; nProfilePoint--)
1299 double const dWaveHeight = VdWaveHeight[nProfilePoint];
1300 double const dWaveAngle = VdWaveDirection[nProfilePoint];
1302 bBreaking = VbWaveIsBreaking[nProfilePoint];
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);
1330 double const dDiffProfileDistXY = hypot(dXDist, dYDist);
1333 double const dtanBeta = tan(
tAbs(
m_pRasterGrid->m_Cell[nX1][nY1].dGetSeaDepth() -
m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth()) / dDiffProfileDistXY);
1336 int nValidPointsWaveHeight = 0;
1337 int nValidPointsWaveSetup = 0;
1339 for (
int nPoint = 0; nPoint < static_cast<int>(VdWaveHeight.size()); nPoint++)
1341 if (VdWaveHeight[nPoint] > 1e-4)
1343 nValidPointsWaveHeight += 1;
1352 nValidPointsWaveHeight -= 1;
1354 for (
int nPoint = 0; nPoint < static_cast<int>(VdWaveSetupSurge.size()); nPoint++)
1356 if (
tAbs(VdWaveSetupSurge[nPoint]) < 1)
1358 nValidPointsWaveSetup += 1;
1367 nValidPointsWaveSetup -= 1;
1369 double dWaveHeight = 0;
1374 dWaveHeight = VdWaveHeight[nValidPointsWaveHeight];
1383 dRunUp = 0.36 * pow(9.81, 0.5) * dtanBeta * pow(dWaveHeight, 0.5) * dDeepWaterWavePeriod;
1389 double const dS0 = 2 *
PI * dWaveHeight / (9.81 * dDeepWaterWavePeriod * dDeepWaterWavePeriod);
1390 dRunUp = 1.86 * dWaveHeight * pow(pow(dtanBeta / dS0, 0.5), 0.71);
1396 double const dS0 = 2 *
PI * dWaveHeight / (9.81 * dDeepWaterWavePeriod * dDeepWaterWavePeriod);
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));
1405 if ((
tAbs(dRunUp) < 1e-4) || (isnan(dRunUp)))
1410 double dWaveSetupSurge = 0;
1415 dWaveSetupSurge = VdWaveSetupSurge[nValidPointsWaveSetup];
1418 if ((
tAbs(dWaveSetupSurge) < 1e-4) || (isnan(dWaveSetupSurge)))
1420 dWaveSetupSurge = 0;
1425 m_VCoast[nCoast].SetCoastWaveHeight(nCoastPoint, dWaveHeight);
1426 m_VCoast[nCoast].SetWaveSetupSurge(nCoastPoint, dWaveSetupSurge);
1427 m_VCoast[nCoast].SetRunUp(nCoastPoint, dRunUp);
1429 if (nProfileBreakingDist > 0)
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);
1454#if defined CSHORE_FILE_INOUT
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)
1461 ofstream CShoreOutStream;
1462 CShoreOutStream.open(
CSHORE_INFILE.c_str(), ios::out | ios::app);
1464 if (CShoreOutStream.fail())
1472 CShoreOutStream << 3 << endl;
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;
1487 CShoreOutStream << setw(11) << nILab <<
" -> ILAB" << endl;
1488 CShoreOutStream << setw(11) << nWave <<
" -> NWAVE" << endl;
1489 CShoreOutStream << setw(11) << nSurge <<
" -> NSURGE" << endl;
1492 CShoreOutStream << setw(11) << setprecision(2) << dWaveInitTime;
1493 CShoreOutStream << setw(11) << setprecision(4) << dWavePeriod;
1494 CShoreOutStream << setw(11) << dHrms;
1495 CShoreOutStream << setw(11) << dWaveAngle << endl;
1498 CShoreOutStream << setw(11) << setprecision(2) << dTimestep;
1499 CShoreOutStream << setw(11) << setprecision(4) << dWavePeriod;
1500 CShoreOutStream << setw(11) << dHrms;
1501 CShoreOutStream << setw(11) << dWaveAngle << endl;
1504 CShoreOutStream << setw(11) << setprecision(2) << dSurgeInitTime;
1505 CShoreOutStream << setw(11) << setprecision(4) << dSurgeLevel << endl;
1508 CShoreOutStream << setw(11) << setprecision(2) << dTimestep;
1509 CShoreOutStream << setw(11) << setprecision(4) << dSurgeLevel << endl;
1512 CShoreOutStream << setw(8) << pVdXdist->size() <<
" -> NBINP" << endl;
1514 CShoreOutStream << fixed << setprecision(4);
1516 for (
unsigned int i = 0; i < pVdXdist->size(); i++)
1518 CShoreOutStream << setw(11) << pVdXdist->at(i) << setw(11) << pVdBottomElevation->at(i) << setw(11) << pVdWaveFriction->at(i) << endl;
1520 CShoreOutStream << endl;
1523 CShoreOutStream.close();
1534 bool bIsBehindIntervention =
false;
1541 double dProfileDistXY = 0;
1542 double dProfileFricFact;
1543 double dPrevDist = -1;
1545 for (
int i = nProfSize - 1; i >= 0; i--)
1551 if (i == nProfSize - 1)
1558 dProfileDistXY = dProfileDistXY + hypot(dXDist, dYDist);
1563 dProfileDistXY += 0.1;
1570 int const nTopLayer =
m_pRasterGrid->m_Cell[nX][nY].nGetTopNonZeroLayerAboveBasement();
1581 double const dTopElev =
m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev() +
m_pRasterGrid->m_Cell[nX][nY].dGetInterventionHeight();
1589 VdVZ->push_back(0.1);
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;
1597 VdVZ->push_back(VdProfileZ);
1603 VdVZ->push_back(VdProfileZ);
1607 VdDistXY->push_back(dProfileDistXY);
1610 double const dInterventionHeight =
m_pRasterGrid->m_Cell[nX][nY].dGetInterventionHeight();
1613 if (dInterventionHeight > 0 || bIsBehindIntervention)
1617 bIsBehindIntervention =
true;
1624 VdFricF->push_back(dProfileFricFact);
1627 dPrevDist = dProfileDistXY;
1633#if defined CSHORE_FILE_INOUT
1637int CSimulation::nReadCShoreOutput(
int const nProfile,
string const *strCShoreFilename,
int const nExpectedColumns,
int const nCShorecolumn, vector<double>
const *pVdProfileDistXYCME, vector<double> *pVdInterpolatedValues)
1641 InStream.open(strCShoreFilename->c_str(), ios::in);
1644 if (!InStream.is_open())
1647 LogStream <<
m_ulIter <<
": " <<
ERR <<
"for profile " << nProfile <<
", cannot open " << *strCShoreFilename <<
" for input" << endl;
1653 vector<double> VdXYDistCShore;
1654 vector<double> VdValuesCShore;
1658 int nExpectedRows = 0;
1661 while (getline(InStream, strLineIn))
1672 string strErr =
ERR +
"invalid integer for number of expected rows '" + VstrItems[1] +
"' in " + *strCShoreFilename +
"\n";
1680 nExpectedRows = stoi(VstrItems[1]);
1688 int nCols =
static_cast<int>(VstrItems.size());
1690 if (nCols != nExpectedColumns)
1693 LogStream <<
m_ulIter <<
": " <<
ERR <<
"for profile " << nProfile <<
", expected " << nExpectedColumns <<
" CShore output columns but read " << nCols <<
" columns from header section of file " << *strCShoreFilename << endl;
1699 VdXYDistCShore.push_back(strtod(VstrItems[0].c_str(), NULL));
1700 VdValuesCShore.push_back(strtod(VstrItems[nCShorecolumn - 1].c_str(), NULL));
1705 int nReadRows =
static_cast<int>(VdXYDistCShore.size());
1707 if (nReadRows != nExpectedRows)
1710 LogStream <<
m_ulIter <<
": " <<
ERR <<
"for profile " << nProfile <<
", expected " << nExpectedRows <<
" CShore output rows, but read " << nReadRows <<
" rows from file " << *strCShoreFilename << endl;
1719 LogStream <<
m_ulIter <<
": " <<
WARN <<
"for profile " << nProfile <<
", only " << nReadRows <<
" CShore output rows in file " << *strCShoreFilename << endl;
1722 VdXYDistCShore.push_back(VdXYDistCShore[0]);
1723 VdValuesCShore.push_back(VdValuesCShore[0]);
1730 vector<double> VdXYDistCShoreTmp(nReadRows, 0);
1732 for (
int i = 0; i < nReadRows; i++)
1733 VdXYDistCShoreTmp[i] = VdXYDistCShore[nReadRows - 1] - VdXYDistCShore[i];
1736 reverse(VdXYDistCShoreTmp.begin(), VdXYDistCShoreTmp.end());
1739 reverse(VdValuesCShore.begin(), VdValuesCShore.end());
1742 vector<double> VdDistXYCopy(pVdProfileDistXYCME->begin(), pVdProfileDistXYCME->end());
1751#if defined CSHORE_ARG_INOUT || CSHORE_BOTH
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)
1762 bool bTooShort =
false;
1765 if (nOutSize < nProfileSize)
1768 nTooShort = nProfileSize - nOutSize - 1;
1789 double dLastDiff = pVdXYDistFromCShore->at(nOutSize - 1) - pVdXYDistFromCShore->at(nOutSize - 2);
1791 for (
int n = 0; n < nTooShort; n++)
1792 pVdXYDistFromCShore->at(nOutSize + n) = pVdXYDistFromCShore->at(nOutSize - 1 + n) + dLastDiff;
1794 for (
int n = 0; n < nTooShort; n++)
1795 pVdFreeSurfaceStdCShore->at(nOutSize + n) = pVdFreeSurfaceStdCShore->at(nOutSize - 1);
1797 for (
int n = 0; n < nTooShort; n++)
1798 pVdWaveSetupSurgeCShore->at(nOutSize + n) = pVdWaveSetupSurgeCShore->at(nOutSize - 1);
1802 for (
int n = 0; n < nTooShort; n++)
1803 pVdSinWaveAngleRadiansCShore->at(nOutSize + n) = pVdSinWaveAngleRadiansCShore->at(nOutSize - 1);
1805 dLastDiff = pVdFractionBreakingWavesCShore->at(nOutSize - 1) - pVdFractionBreakingWavesCShore->at(nOutSize - 2);
1807 for (
int n = 0; n < nTooShort; n++)
1808 pVdFractionBreakingWavesCShore->at(nOutSize + n) =
tMin(pVdFractionBreakingWavesCShore->at(nOutSize - 1 + n) + dLastDiff, 1.0);
1819 vector<double> VdXYDistCShoreTmp(nProfileSize);
1820 copy(pVdXYDistFromCShore->begin(), pVdXYDistFromCShore->begin() + nProfileSize, begin(VdXYDistCShoreTmp));
1823 vector<double> VdFreeSurfaceStdCShoreTmp(nProfileSize);
1824 reverse_copy(pVdFreeSurfaceStdCShore->begin(), pVdFreeSurfaceStdCShore->begin() + nProfileSize, begin(VdFreeSurfaceStdCShoreTmp));
1826 vector<double> VdWaveSetupSurgeCShoreTmp(nProfileSize);
1827 reverse_copy(pVdWaveSetupSurgeCShore->begin(), pVdWaveSetupSurgeCShore->begin() + nProfileSize, begin(VdWaveSetupSurgeCShoreTmp));
1829 vector<double> VdSinWaveAngleRadiansCShoreTmp(nProfileSize);
1830 reverse_copy(pVdSinWaveAngleRadiansCShore->begin(), pVdSinWaveAngleRadiansCShore->begin() + nProfileSize, begin(VdSinWaveAngleRadiansCShoreTmp));
1832 vector<double> VdFractionBreakingWavesCShoreTmp(nProfileSize);
1833 reverse_copy(pVdFractionBreakingWavesCShore->begin(), pVdFractionBreakingWavesCShore->begin() + nProfileSize, begin(VdFractionBreakingWavesCShoreTmp));
1870 bool bModfiedWaveHeightisBreaking =
false;
1871 bool bProfileIsinShadowZone =
false;
1874 int nThisBreakingDist =
m_VCoast[nCoast].nGetBreakingDistance(nThisCoastPoint);
1875 double dThisBreakingWaveHeight =
m_VCoast[nCoast].dGetBreakingWaveHeight(nThisCoastPoint);
1876 double dThisBreakingWaveAngle =
m_VCoast[nCoast].dGetBreakingWaveAngle(nThisCoastPoint);
1877 double dThisBreakingDepth =
m_VCoast[nCoast].dGetDepthOfBreaking(nThisCoastPoint);
1880 for (
int nProfilePoint = (nProfileSize - 1); nProfilePoint >= 0; nProfilePoint--)
1888 bProfileIsinShadowZone =
true;
1891 double const dWaveHeight =
m_pRasterGrid->m_Cell[nX][nY].dGetWaveHeight();
1897 bModfiedWaveHeightisBreaking =
true;
1900 dThisBreakingWaveAngle =
m_pRasterGrid->m_Cell[nX][nY].dGetWaveAngle();
1902 nThisBreakingDist = nProfilePoint;
1908 if (bProfileIsinShadowZone && bModfiedWaveHeightisBreaking)
1911 if (dThisBreakingWaveHeight > dThisBreakingDepth * 0.78)
1913 dThisBreakingWaveHeight = dThisBreakingDepth * 0.78;
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);
1925 else if (bProfileIsinShadowZone && (! bModfiedWaveHeightisBreaking))
1953 int const nThisBreakingDist =
m_VCoast[nCoast].nGetBreakingDistance(nThisCoastPoint);
1954 double const dThisBreakingWaveHeight =
m_VCoast[nCoast].dGetBreakingWaveHeight(nThisCoastPoint);
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);
1959 double const dThisRunUp =
m_VCoast[nCoast].dGetRunUp(nThisCoastPoint);
1976 pTmpProfile = pNextProfile;
1987 int const nDistBetween = nNextCoastPoint - nThisCoastPoint;
1990 if (nDistBetween <= 0)
1994 int const nNextBreakingDist =
m_VCoast[nCoast].nGetBreakingDistance(nNextCoastPoint);
1995 double const dNextBreakingWaveHeight =
m_VCoast[nCoast].dGetBreakingWaveHeight(nNextCoastPoint);
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);
2002 for (
int n = nThisCoastPoint; n <= nNextCoastPoint; n++)
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;
2011 dWaveSetupSurge = (dThisWeight * dThisWaveSetupSurge) + (dNextWeight * dNextWaveSetupSurge);
2012 m_VCoast[nCoast].SetWaveSetupSurge(n, dWaveSetupSurge);
2014 dRunUp = (dThisWeight * dThisRunUp) + (dNextWeight * dNextRunUp);
2015 m_VCoast[nCoast].SetRunUp(n, dRunUp);
2027 for (
int n = nThisCoastPoint; n < nNextCoastPoint; n++)
2042 for (
int n = nThisCoastPoint; n < nNextCoastPoint; n++)
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);
2058 for (
int n = nThisCoastPoint + 1; n <= nNextCoastPoint; n++)
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);
2072 for (
int n = nThisCoastPoint + 1; n < nNextCoastPoint; n++)
2074 int const nDist = n - nThisCoastPoint;
2081 if ((dNextBreakingDepth > 0) && (dThisBreakingDepth > 0))
2083 double const dThisWeight = (nDistBetween - nDist) /
static_cast<double>(nDistBetween);
2084 double const dNextWeight = 1 - dThisWeight;
2086 dBreakingWaveHeight = (dThisWeight * dThisBreakingWaveHeight) + (dNextWeight * dNextBreakingWaveHeight),
2087 dBreakingWaveAngle = (dThisWeight * dThisBreakingWaveAngle) + (dNextWeight * dNextBreakingWaveAngle),
2088 dBreakingDepth = (dThisWeight * dThisBreakingDepth) + (dNextWeight * dNextBreakingDepth),
2089 dBreakingDist = (dThisWeight * nThisBreakingDist) + (dNextWeight * nNextBreakingDist);
2092 else if (dNextBreakingDepth > 0)
2094 dBreakingWaveHeight = dNextBreakingWaveHeight,
2095 dBreakingWaveAngle = dNextBreakingWaveAngle,
2096 dBreakingDepth = dNextBreakingDepth,
2097 dBreakingDist = nNextBreakingDist;
2100 else if (dThisBreakingDepth > 0)
2102 dBreakingWaveHeight = dThisBreakingWaveHeight,
2103 dBreakingWaveAngle = dThisBreakingWaveAngle,
2104 dBreakingDepth = dThisBreakingDepth,
2105 dBreakingDist = nThisBreakingDist;
2110 m_VCoast[nCoast].SetBreakingWaveHeight(n, dBreakingWaveHeight);
2111 m_VCoast[nCoast].SetBreakingWaveAngle(n, dBreakingWaveAngle);
2112 m_VCoast[nCoast].SetDepthOfBreaking(n, dBreakingDepth);
2122 int const nCoastPoints =
m_VCoast[nCoast].nGetCoastlineSize();
2125 vector<int> nVCoastWaveHeightX;
2126 vector<double> dVCoastWaveHeightY;
2129 for (
int n = 0; n < nCoastPoints; n++)
2131 double const dCoastWaveHeight =
m_VCoast[nCoast].dGetCoastWaveHeight(n);
2135 nVCoastWaveHeightX.push_back(n);
2136 dVCoastWaveHeightY.push_back(dCoastWaveHeight);
2141 if ((nVCoastWaveHeightX.size() >= 3) && (nVCoastWaveHeightX.size() == dVCoastWaveHeightY.size()))
2143 for (
int n = 0; n < nCoastPoints; n++)
2145 double const dInterpCoastWaveHeight =
dGetInterpolatedValue(&nVCoastWaveHeightX, &dVCoastWaveHeightY, n,
false);
2146 m_VCoast[nCoast].SetCoastWaveHeight(n, dInterpCoastWaveHeight);
2152 for (
int n = 0; n < nCoastPoints; n++)
2154 m_VCoast[nCoast].SetCoastWaveHeight(n, 0);
2164 int const nCoastSize =
m_VCoast[nCoast].nGetCoastlineSize();
2166 for (
int nCoastPoint = 0; nCoastPoint < nCoastSize; nCoastPoint++)
2171 if (nCoastPoint == 0)
2181 else if (nCoastPoint == nCoastSize - 1)
2207 m_VCoast[nCoast].SetFluxOrientation(nCoastPoint, 90);
2211 m_VCoast[nCoast].SetFluxOrientation(nCoastPoint, 270);
2219 m_VCoast[nCoast].SetFluxOrientation(nCoastPoint, 0);
2223 m_VCoast[nCoast].SetFluxOrientation(nCoastPoint, 180);
2229 double dAzimuthAngle;
2232 dAzimuthAngle = (180 /
PI) * (
PI * 0.5 - atan(dYDiff / dXDiff));
2235 dAzimuthAngle = (180 /
PI) * (
PI * 1.5 - atan(dYDiff / dXDiff));
2237 m_VCoast[nCoast].SetFluxOrientation(nCoastPoint, dAzimuthAngle);
2259 int const nPolyID =
m_pRasterGrid->m_Cell[nX][nY].nGetPolygonID();
2262 bool const bActive =
m_pRasterGrid->m_Cell[nX][nY].bIsInActiveZone();
2267 double const dTmpd50 =
m_pRasterGrid->m_Cell[nX][nY].dGetUnconsD50();
2274 VnPolygonD50Count[nPolyID]++;
2275 VdPolygonD50[nPolyID] += dTmpd50;
2300 int nDownDriftNum = 0;
2303 double dWaveHeight = 0;
2304 double dWaveAngle = 0;
2312 if (
m_pRasterGrid->m_Cell[nXTmp][nYTmp].bIsInContiguousSea())
2315 dWaveHeight +=
m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveHeight();
2316 dWaveAngle += (
m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveAngle());
2321 int nTmp =
m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetShadowZoneNumber();
2329 nTmp =
m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetDownDriftZoneNumber();
2334 nDownDriftNum = nTmp;
2348 if (
m_pRasterGrid->m_Cell[nXTmp][nYTmp].bIsInContiguousSea())
2351 dWaveHeight +=
m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveHeight();
2352 dWaveAngle += (
m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveAngle());
2357 int nTmp =
m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetShadowZoneNumber();
2365 nTmp =
m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetDownDriftZoneNumber();
2370 nDownDriftNum = nTmp;
2384 if (
m_pRasterGrid->m_Cell[nXTmp][nYTmp].bIsInContiguousSea())
2387 dWaveHeight +=
m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveHeight();
2388 dWaveAngle += (
m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveAngle());
2393 int nTmp =
m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetShadowZoneNumber();
2401 nTmp =
m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetDownDriftZoneNumber();
2406 nDownDriftNum = nTmp;
2420 if (
m_pRasterGrid->m_Cell[nXTmp][nYTmp].bIsInContiguousSea())
2423 dWaveHeight +=
m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveHeight();
2424 dWaveAngle += (
m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetWaveAngle());
2429 int nTmp =
m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetShadowZoneNumber();
2437 nTmp =
m_pRasterGrid->m_Cell[nXTmp][nYTmp].nGetDownDriftZoneNumber();
2442 nDownDriftNum = nTmp;
2453 dWaveHeight /= nRead;
2454 dWaveAngle /= nRead;
2466 double const dDeepWaterWaveHeight =
m_pRasterGrid->m_Cell[nX][nY].dGetCellDeepWaterWaveHeight();
2474 double const dDeepWaterWaveAngle =
m_pRasterGrid->m_Cell[nX][nY].dGetCellDeepWaterWaveAngle();
2482 int const nShadowZoneCode =
m_pRasterGrid->m_Cell[nX][nY].nGetShadowZoneNumber();
2484 if (nShadowZoneCode <= 0)
2487 if ((nShadow == 4) || (nShadow + nDownDrift == 4) || (nShadow + nCoast == 4))
2489 m_pRasterGrid->m_Cell[nX][nY].SetShadowZoneNumber(nShadowNum);
2496 int const nDownDriftZoneCode =
m_pRasterGrid->m_Cell[nX][nY].nGetDownDriftZoneNumber();
2498 if ((nDownDriftZoneCode == 0) && (nDownDrift == 4))
2500 m_pRasterGrid->m_Cell[nX][nY].SetDownDriftZoneNumber(nDownDriftNum);
2510 for (
int nCoast = 0; nCoast < static_cast<int>(
m_VCoast.size()); nCoast++)
2512 for (
int nPoly = 0; nPoly <
m_VCoast[nCoast].nGetNumPolygons(); nPoly++)
2517 if (VnPolygonD50Count[nPolyID] > 0)
2518 VdPolygonD50[nPolyID] /= VnPolygonD50Count[nPolyID];
Contains CGeom2DPoint definitions.
int nGetY(void) const
Returns the CGeom2DIPoint object's integer Y coordinate.
int nGetX(void) const
Returns the CGeom2DIPoint object's integer X coordinate.
Geometry class used to represent 2D point objects with floating-point coordinates.
double dGetY(void) const
Returns the CGeom2DPoint object's double Y coordinate.
double dGetX(void) const
Returns the CGeom2DPoint object's double X coordinate.
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.
double dGetProfileDeepWaterWaveAngle(void) const
Returns the deep-water wave orientation for this profile.
double dGetProfileDeepWaterWaveHeight(void) const
Returns the deep-water wave height for this profile.
int nGetProfileID(void) const
Returns the profile's this-coast ID.
int nGetCoastPoint(void) const
Returns the coast point at which the profile starts.
CGeomProfile * pGetDownCoastAdjacentProfile(void) const
Returns the down-coast adjacent profile, returns NULL if there is no down-coast adjacent profile.
CGeom2DIPoint * pPtiGetCellInProfile(int const)
Returns a single cell (grid CRS) in the profile.
int nGetNumCellsInProfile(void) const
Returns the number of cells in the profile.
void SetCShoreProblem(bool const)
Sets a switch to indicate whether this profile has a CShore problem.
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...
bool bStartOfCoast(void) const
Returns the switch to indicate whether this is a start-of-coast profile.
vector< CGeom2DIPoint > * pPtiVGetCellsInProfile(void)
Returns all cells (grid CRS) in the profile.
bool bEndOfCoast(void) const
Returns the switch to indicate whether this is an end-of-coast profile.
double dGetProfileDeepWaterWavePeriod(void) const
Returns the deep-water wave period for this profile.
int m_nLogFileDetail
The level of detail in the log file output. Can be LOG_FILE_LOW_DETAIL, LOG_FILE_MIDDLE_DETAIL,...
double m_dAllCellsDeepWaterWaveHeight
Deep water wave height (m) for all sea cells.
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)
double m_dThisIterSWL
The still water level for this timestep (this includes tidal changes and any long-term SWL change)
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.
double m_dL_0
Deep water wave length (m)
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.
int m_nYGridSize
The size of the grid in the y direction.
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.
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.
double m_dBreakingWaveHeightDepthRatio
Breaking wave height-to-depth ratio.
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...
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.
double m_dAllCellsDeepWaterWaveAngle
Deep water wave angle for all sea cells.
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.
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)
unsigned long m_ulIter
The number of the current iteration (time step)
double dGridCentroidXToExtCRSX(int const) const
int m_nNumPolygonGlobal
Number of global (all coasts) polygons.
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...
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...
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
int const WAVE_MODEL_COVE
int const RTN_ERR_NO_TOP_LAYER
int const RTN_ERR_CSHORE_FILE_INPUT
double const CSHORE_FRICTION_FACTOR
int const WAVE_MODEL_CSHORE
int const LOG_FILE_MIDDLE_DETAIL
int const RTN_ERR_CSHORE_ERROR
bool bFPIsEqual(const T d1, const T d2, const T dEpsilon)
double const WALKDEN_HALL_PARAM_2
int const RTN_ERR_READING_CSHORE_FILE_OUTPUT
double const CSHORE_SURGE_LEVEL
bool const SAVE_CSHORE_OUTPUT
int const CSHOREARRAYOUTSIZE
The size of the arrays output by CShore. If you change this, then you must also set the same value on...
double const WALKDEN_HALL_PARAM_1
int const RTN_ERR_CSHORE_EMPTY_PROFILE
string const CSHORE_INFILE
Contains CRWCoast definitions.
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.