From b085a47f90760c1f4b07f8d0bb7ecd3b3d4f80b9 Mon Sep 17 00:00:00 2001 From: lawfordp2017 Date: Mon, 14 Feb 2022 13:36:04 -0600 Subject: [PATCH 1/5] Initial add of corrected Muskingum implementation, corrections to Makefile --- crhmcode/Makefile | 12 +- crhmcode/src/core/ClassMuskingum2.cpp | 180 ++++++ crhmcode/src/core/ClassMuskingum2.h | 42 ++ crhmcode/src/modules/ClassNetroute_M2_D.cpp | 598 ++++++++++++++++++++ crhmcode/src/modules/ClassNetroute_M2_D.h | 119 ++++ 5 files changed, 945 insertions(+), 6 deletions(-) create mode 100644 crhmcode/src/core/ClassMuskingum2.cpp create mode 100644 crhmcode/src/core/ClassMuskingum2.h create mode 100644 crhmcode/src/modules/ClassNetroute_M2_D.cpp create mode 100644 crhmcode/src/modules/ClassNetroute_M2_D.h diff --git a/crhmcode/Makefile b/crhmcode/Makefile index ba09182e9..c4caf51be 100644 --- a/crhmcode/Makefile +++ b/crhmcode/Makefile @@ -13,7 +13,7 @@ CFLAGS = -Wno-switch -Wno-c++11-extensions -Wno-logical-op-parentheses \ CFLAGS = -I. -I$(BOOST) -I$(SDIR)/core -I$(SDIR)/macro \ -I$(SDIR)/modules/waterquality -I$(SDIR)/modules/newmodules \ --I$(SDIR)/modules/coremodules -I$(SDIR)/gcc -std=c++14 -static -static-libgcc -static-libstdc++ +-I$(SDIR)/modules/coremodules -I$(SDIR)/gcc -std=c++14 -static-libgcc -static-libstdc++ -lfmt CORE_OBJS = GlobalDll.o StandardConverterUtility.o @@ -91,7 +91,7 @@ OBJ = $(patsubst %,$(ODIR)/%,$(ALL_OBS)) crhm: $(OBJ) $(CC) -o $@ $^ $(CFLAGS) $(DFLAGS) $(LIBS) -$(ODIR)/main.o: $(SDIR)/gcc/main.cpp | $(ODIR) +$(ODIR)/main.o: $(SDIR)/main.cpp | $(ODIR) $(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS) $(ODIR)/%.o: $(SDIR)/core/%.cpp $(SDIR)/core/%.h | $(ODIR) @@ -103,7 +103,7 @@ $(ODIR)/%.o: $(SDIR)/core/ClassModule/%.cpp $(SDIR)/core/ClassModule/%.h | $(ODI $(ODIR)/%.o: $(SDIR)/core/ClassCRHM/%.cpp $(SDIR)/core/ClassCRHM/%.h | $(ODIR) $(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS) -$(ODIR)/%.o: $(SDIR)/core/ClassCRHM/Filters/%.cpp $(SDIR)/core/ClassCRHM/Filters/%.h | $(ODIR) +$(ODIR)/%.o: $(SDIR)/core/filter/%.cpp $(SDIR)/core/filter/%.h | $(ODIR) $(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS) $(ODIR)/%.o: $(SDIR)/core/TStringList/%.cpp $(SDIR)/core/TStringList/%.h | $(ODIR) @@ -112,10 +112,10 @@ $(ODIR)/%.o: $(SDIR)/core/TStringList/%.cpp $(SDIR)/core/TStringList/%.h | $(ODI $(ODIR)/%.o: $(SDIR)/core/Common/%.cpp $(SDIR)/core/Common/%.h | $(ODIR) $(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS) -$(ODIR)/%.o: $(SDIR)/core/CRHM_parse/%.cpp $(SDIR)/core/CRHM_parse/%.h | $(ODIR) +$(ODIR)/%.o: $(SDIR)/core/def/%.cpp $(SDIR)/core/def/%.h | $(ODIR) $(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS) -$(ODIR)/%.o: $(SDIR)/core/CRHM_parse/exec/%.cpp $(SDIR)/core/CRHM_parse/exec/%.h | $(ODIR) +$(ODIR)/%.o: $(SDIR)/core/exec/%.cpp $(SDIR)/core/exec/%.h | $(ODIR) $(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS) $(ODIR)/%.o: $(SDIR)/core/CRHMmain/%.cpp $(SDIR)/core/CRHMmain/%.h | $(ODIR) @@ -136,7 +136,7 @@ $(ODIR)/%.o: $(SDIR)/modules/waterquality/%.cpp $(SDIR)/modules/waterquality/%.h $(ODIR)/%.o: $(SDIR)/modules/newmodules/%.cpp $(SDIR)/modules/newmodules/%.h | $(ODIR) $(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS) -$(ODIR)/%.o: $(SDIR)/modules/coremodules/%.cpp $(SDIR)/modules/coremodules/%.h | $(ODIR) +$(ODIR)/%.o: $(SDIR)/modules/%.cpp $(SDIR)/modules/%.h | $(ODIR) $(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS) clean: diff --git a/crhmcode/src/core/ClassMuskingum2.cpp b/crhmcode/src/core/ClassMuskingum2.cpp new file mode 100644 index 000000000..f56a84463 --- /dev/null +++ b/crhmcode/src/core/ClassMuskingum2.cpp @@ -0,0 +1,180 @@ +#include "ClassMuskingum2.h" + + +//--------------------------------------------------------------------------- +ClassMuskingum2::ClassMuskingum2(const double* inVar, double* outVar, const double* k, const double* X_M, const double* lag, const long nhru, const long setlag) + : inVar(inVar), outVar(outVar), nhru(nhru) { + + // !!! UNITS !!! + + // kstorage (days) + // lag (hours) + + LastIn = new double[nhru]; + LastOut = new double[nhru]; + + c0 = new double[nhru]; + c1 = new double[nhru]; + c2 = new double[nhru]; + + /* + * ilag for routing water along a stream, it is used to simulate the delay caused by the travel time of water. + * ilag is the length of the stream (the number of elements or letters). + * If on every cycle the water moves one step then the number of elements determines how many steps + * are required to get to the end of the stream. + */ + ilag = new long[nhru]; + maxlag = new long[nhru]; + ulag = new long[nhru]; + long Biggest = 0; + + for (long hh = 0; hh < nhru; hh++) { + + c0[hh] = (Global::Interval - 2.0 * k[hh] * X_M[hh]) / + (2.0 * k[hh] * (1.0 - X_M[hh]) + Global::Interval); // units of Global::Interval (days) + + c1[hh] = (Global::Interval + 2.0 * k[hh] * X_M[hh]) / + (2.0 * k[hh] * (1.0 - X_M[hh]) + Global::Interval); // units of kstorage (days) + + c2[hh] = (2.0 * k[hh] * (1.0 - X_M[hh]) - Global::Interval) / + (2.0 * k[hh] * (1.0 - X_M[hh]) + Global::Interval); // units of kstorage (days) + + ilag[hh] = (long)(max(lag[hh], 0.0) / 24.0 * Global::Freq + 1.1); // =1 for lag of zero + + if (setlag == -1 || ilag[hh] > setlag) + maxlag[hh] = ilag[hh]; + else + maxlag[hh] = setlag; + + ulag[hh] = 0; + + LastIn[hh] = 0.0; // zero initial conditions + + LastOut[hh] = 0.0; // zero initial conditions + + if (maxlag[hh] > Biggest) Biggest = maxlag[hh]; + } + + LagArray = new double* [nhru]; // create lag array + + for (long hh = 0; hh < nhru; hh++) { + LagArray[hh] = new double[maxlag[hh]]; + for (long jj = 0; jj < maxlag[hh]; jj++) + LagArray[hh][jj] = 0.0; + } +} + +ClassMuskingum2::~ClassMuskingum2() { + delete[] LastIn; + delete[] LastOut; + delete[] c0; + delete[] c1; + delete[] c2; + delete[] ilag; + delete[] maxlag; + delete[] ulag; + + for (long hh = 0; hh < nhru; hh++) + delete[] LagArray[hh]; + delete[] LagArray; +} + +void ClassMuskingum2::ChangeLag(const double* newlag, const long hh) +{ + + long newilag = (long)(max(newlag[hh], 0.0) / 24.0 * Global::Freq + 1.1); // =1 for lag of zero + + if (newilag == ilag[hh]) { + return; + } + + double* AccArray = new double[ilag[hh]]; // work area for ChangeLag + + AccArray[0] = 0.0; + + for (int ii = 1; ii < ilag[hh]; ++ii) + { + AccArray[ii] = AccArray[ii - 1] + LagArray[hh][(ulag[hh] + ii) % ilag[hh]]; // accumulate storage + } + + delete[] LagArray[hh]; // delete previous length + + LagArray[hh] = new double[newilag]; // create new length + + ulag[hh] = 0; // next input value save here. + LagArray[hh][0] = 0.0; // looks better + + double LastValue = 0.0; + + for (int mm = 1; mm < newilag - 1; ++mm) + { + double Y = double(mm) / ((long long)newilag - 1ll) * ((long long)ilag[hh] - 1ll); + int Yint = (int)(Y + 0.0001); + if ((Yint + 1) > ilag[hh] - 1) + { + CRHMException Except("Attempting to read out of bounds array address", TExcept::TERMINATE); + LogError(Except); + throw(Except); + } + double Ydif = Y - Yint; + + + double NewValue = AccArray[Yint] + Ydif * (AccArray[Yint + 1] - AccArray[Yint]); + + + + + LagArray[hh][(ulag[hh] + mm) % newilag] = NewValue - LastValue; + + LastValue = NewValue; + } + + LagArray[hh][(ulag[hh] + newilag - 1) % newilag] = AccArray[ilag[hh] - 1] - LastValue; // final values + + delete[] AccArray; // work area for ChangeLag + + ilag[hh] = newilag; // assign new lag +} + +void ClassMuskingum2::DoMuskingum() { + + for (long hh = 0; hh < nhru; hh++) { + + LagArray[hh][ulag[hh]] = inVar[hh]; + + ulag[hh] = ++ulag[hh] % ilag[hh]; + + outVar[hh] = c0[hh] * LagArray[hh][ulag[hh]] + c1[hh] * LastIn[hh] + c2[hh] * LastOut[hh]; + + LastIn[hh] = LagArray[hh][ulag[hh]]; + + LastOut[hh] = outVar[hh]; + } +} + +void ClassMuskingum2::DoMuskingum(const long hh) { + + LagArray[hh][ulag[hh]] = inVar[hh]; + + ulag[hh] = ++ulag[hh] % ilag[hh]; // now points to fully delayed value + + outVar[hh] = c0[hh] * LagArray[hh][ulag[hh]] + c1[hh] * LastIn[hh] + c2[hh] * LastOut[hh]; + + LastIn[hh] = LagArray[hh][ulag[hh]]; + + LastOut[hh] = outVar[hh]; +} + +double ClassMuskingum2::Left(int hh) { + + double Slag = 0; + + for (int ii = 1; ii < ilag[hh]; ++ii) + Slag += LagArray[hh][(ulag[hh] + ii) % ilag[hh]]; + + double Sstorage = (1.0 / (1.0 - c2[hh])) * (c1[hh] * LastIn[hh] + c2[hh] * outVar[hh]); + + return Slag + Sstorage; +} + +//--------------------------------------------------------------------------- diff --git a/crhmcode/src/core/ClassMuskingum2.h b/crhmcode/src/core/ClassMuskingum2.h new file mode 100644 index 000000000..947df6da2 --- /dev/null +++ b/crhmcode/src/core/ClassMuskingum2.h @@ -0,0 +1,42 @@ +#pragma once + +#ifndef ClassMuskingum2 +#define ClassMuskingum2 + +#include "ClassModule.h" + + +class ClassMuskingum2 { + +public: + ClassMuskingum2(const double* inVar, double* outVar, const double* kstorage, const double* route_X_M, const double* lag, const long nhru, const long setlag = -1); + ~ClassMuskingum2(); + void DoMuskingum(); + void DoMuskingum(const long hh); + void ChangeLag(const double* newlag, const long hh); + double Left(int hh); + + + double* c0; // storage constant from K + double* c1; // storage constant from K + double* c2; // storage constant from K + double prevdate{ 0.0 }; + +private: + const double* kstorage{ NULL }; + const double* inVar; + double* outVar; + + double** LagArray; + + double* LastIn; // + double* LastOut; // + + long nhru; + long* maxlag; // maximum lag - i.e. storage + + long* ilag; // lag interval (hours) + long* ulag; // lag interval (#intervals) +}; + +#endif // !ClassMuskingum2 diff --git a/crhmcode/src/modules/ClassNetroute_M2_D.cpp b/crhmcode/src/modules/ClassNetroute_M2_D.cpp new file mode 100644 index 000000000..67975fb54 --- /dev/null +++ b/crhmcode/src/modules/ClassNetroute_M2_D.cpp @@ -0,0 +1,598 @@ +//created by Manishankar Mondal + +#include +#include +#include +#include +#include +#include + +#include "ClassNetroute_M2_D.h" +#include "../core/GlobalDll.h" +#include "../core/ClassCRHM.h" +#include "newmodules/SnobalDefines.h" + + +using namespace CRHM; + +ClassNetroute_M2_D* ClassNetroute_M2_D::klone(string name) const{ + return new ClassNetroute_M2_D(name); +} + +double ClassNetroute_M2_D::Function1(double *I, long hh) { + return runDelay->ChangeLag(I, hh); +} + +double ClassNetroute_M2_D::Function2(double *X, long hh) { + return runDelay->ChangeStorage(X, hh); +} + +void ClassNetroute_M2_D::decl(void) { + + Description = "'Handles the routing of surface runoff, subsurface runoff and HRU routing using the lag and route and Muskingum method for \"outflow\". Distributed flow.'"; + + declvar("inflow", TDim::NHRU, "inflow from other HRUs", "(mm*km^2/int)", &inflow); + + declstatdiag("cuminflow", TDim::NHRU, "cumulative inflow from other HRUs", "(mm*km^2)", &cuminflow); + + declvar("outflow", TDim::NHRU, "HRU outflow", "(mm*km^2/int)", &outflow); + + declstatdiag("cumoutflow", TDim::NHRU, "cumulative HRU outflow", "(mm*km^2)", &cumoutflow); + + decldiag("outflow_diverted", TDim::NHRU, "HRU outflow diverted to another HRU", "(mm*km^2/int)", &outflow_diverted); + + declstatdiag("cumoutflow_diverted", TDim::NHRU, "cumulative HRU outflow diverted to another HRU", "(mm*km^2/int)", &cumoutflow_diverted); + + declstatdiag("cum_to_Sd", TDim::NHRU, "cumulative other HRU to depressional storage (Sd) of this HRU", "(mm)", &cum_to_Sd); + + declstatdiag("cum_to_soil_rechr", TDim::NHRU, "cumulative other HRU recycled tosoil_rechr of this HRU", "(mm)", &cum_to_soil_rechr); + + declvar("gwinflow", TDim::NHRU, "ground water inflow", "(mm*km^2/int)", &gwinflow); + + declstatdiag("gwcuminflow", TDim::NHRU, "cumulative gw inflow", "(mm*km^2)", &gwcuminflow); + + declvar("gwoutflow", TDim::NHRU, "HRU gw outflow", "(mm*km^2/int)", &gwoutflow); + + declstatdiag("gwcumoutflow", TDim::NHRU, "cumulative HRU gw outflow", "(mm*km^2)", &gwcumoutflow); + + decldiag("gwoutflow_diverted", TDim::NHRU, "HRU gw outflow diverted to another HRU", "(mm*km^2/int)", &gwoutflow_diverted); + + declstatdiag("gwcumoutflow_diverted", TDim::NHRU, "cumulative HRU gw outflow diverted to another HRU", "(mm*km^2/int)", &gwcumoutflow_diverted); + + declvar("ssrinflow", TDim::NHRU, "inflow from other HRUs", "(mm*km^2/int)", &ssrinflow); + + declstatdiag("ssrcuminflow", TDim::NHRU, "cumulative inflow from other HRUs", "(mm*km^2)", &ssrcuminflow); + + declvar("ssroutflow", TDim::NHRU, "HRU outflow", "(mm*km^2/int)", &ssroutflow); + + declstatdiag("ssrcumoutflow", TDim::NHRU, "cumulative HRU outflow", "(mm*km^2)", &ssrcumoutflow); + + declstatdiag("HRU_cumbasinflow", TDim::NHRU, "cumulative HRU to basinflow", "(mm*km^2)", &HRU_cumbasinflow); + + declstatdiag("cum_preferential_flow_to_gw", TDim::NHRU, "cumulative other HRU's runoff to this HRU's gw via preferential flow path", "(mm)", &cum_preferential_flow_to_gw); + + declparam("preferential_flow", TDim::NHRU, "[0]", "0", "1","0 - no preferential and remain as runoff routing to other HRU, 1 - preferential flow and route runoff to other HRU's gw.", "()", &preferential_flow); + + declvar("runinflow", TDim::NHRU, "inflow from other HRUs", "(mm*km^2/int)", &runinflow); + + declstatdiag("runcuminflow", TDim::NHRU, "cumulative inflow from other HRUs", "(mm*km^2)", &runcuminflow); + + declvar("runoutflow", TDim::NHRU, "HRU outflow", "(mm*km^2/int)", &runoutflow); + + declstatdiag("runcumoutflow", TDim::NHRU, "cumulative HRU outflow", "(mm*km^2)", &runcumoutflow); + + declstatdiag("cumscaling_boost", TDim::NHRU, "cumulative amout inflow boosted", "(mm*km^2)", &cumscaling_boost); + + declvar("basinflow", TDim::BASIN, "basin surface and sub-surface outflow", "(m^3/int)", &basinflow); + + decldiag("basinflow_s", TDim::BASIN, "basin surface and sub-surface outflow", "(m^3/s)", &basinflow_s); + + declstatdiag("cumbasinflow", TDim::BASIN, "cumulative basin surface and sub-surface outflow", "(m^3)", &cumbasinflow); + + declvar("basingw", TDim::BASIN, "cumulative basin groundwater outflow", "(m^3/int)", &basingw); + + decldiag("basingw_s", TDim::BASIN, "cumulative basin groundwater outflow", "(m^3/s)", &basingw_s); + + declstatdiag("cumbasingw", TDim::BASIN, "cumulative basin groundwater outflow", "(m^3)", &cumbasingw); + + decllocal("soil_ssr_Buf", TDim::NHRU, "buffer subsurface runoff", "(mm/d)", &soil_ssr_Buf); + + decllocal("soil_runoff_Buf", TDim::NHRU, "buffer rain runoff", "(mm/d)", &soil_runoff_Buf); + + decllocal("soil_gw_Buf", TDim::NHRU, "buffer rain runoff", "(mm/d)", &soil_gw_Buf); + + decllocal("Ktravel", TDim::NHRU, "travel time", "(d)", &Ktravel); + + decllocal("distrib_sum", TDim::NHRU, "HRU distribution sum", "()", &distrib_sum); + + + declparam("basin_area", TDim::BASIN, "3", "1e-6", "1e09", "Total basin area", "(km^2)", &basin_area); + + declparam("hru_area", TDim::NHRU, "[1]", "1e-6", "1e09", "HRU area", "(km^2)", &hru_area); + + declparam("Lag", TDim::NHRU, "[0.0]", "0.0","1.0E4.0", "lag delay", "(h)", &Lag); + + declparam("route_n", TDim::NHRU, "[0.025]", "0.016","0.2", "Manning roughness coefficient", "()", &route_n); + + declparam("route_R", TDim::NHRU, "[0.5]", "0.01","1.0E4", "hydraulic radius", "(m)", &route_R); + + declparam("route_S0", TDim::NHRU, "[1e-3]", "1e-6","1.0", "longitudinal channel slope", "()", &route_S0); + + declparam("route_L", TDim::NHRU, "[200.0]", "0.01","1.0E10", "routing length", "(m)", &route_L); + + declparam("route_X_M", TDim::NHRU, "[0.25]", "0.0","0.5", "dimensionless weighting factor", "()", &route_X_M); + + declparam("ssrKstorage", TDim::NHRU, "[0.0]", "0.0","200.0", "subsurface runoff storage constant", "(d)", &ssrKstorage); + + declparam("ssrLag", TDim::NHRU, "[0.0]", "0.0","1.0E4.0", "subsurface runoff lag delay", "(h)", &ssrLag); + + declparam("runKstorage", TDim::NHRU, "[0.0]", "0.0","200.0", "runoff storage constant", "(d)", &runKstorage); + + declparam("runLag", TDim::NHRU, "[0.0]", "0.0","1.0E4", "runoff lag delay", "(h)", &runLag); + + declparam("gwKstorage", TDim::NHRU, "[0.0]", "0.0","200.0", "gw storage constant", "(d)", &gwKstorage); + + declparam("gwLag", TDim::NHRU, "[0.0]", "0.0","1.0E4", "gw lag delay", "(h)", &gwLag); + + declparam("distrib_Route", TDim::NDEFN, "[0.0]", "-1.0E6.0", "1.0E6.0", "route this HRU to these HRUs", "()", &distrib, &distrib_hru, nhru); + + declparam("distrib_Basin", TDim::NHRU, "[1.0]", "0.0", "100.0", "route this HRU to basin (and other HRU(s) determined by 'distrib_Route')", "()", &distrib_Basin); + + declparam("gwwhereto", TDim::NHRU, "[0]", "-1000", "1000", "send to: 0 - basingw, >0 - other HRU surface input <0 - other abs(-HRU) gw input, or (< -HRUmax or > +HRUmax) - surface basinflow", "()", &gwwhereto); + + declparam("Sdmax", TDim::NHRU, "[0]", "0.0", "1000.0","Maximum depression storage", "(mm)", &Sdmax); + + declparam("soil_rechr_max", TDim::NHRU, "[60.0]", "0.0", "350.0", + "Maximum available water holding capacity for soil recharge zone (upper portion of soil_moist where losses occur as both evaporation "// + "and transpiration). Must be less than or equal to soil_moist.", "(mm)", &soil_rechr_max); + + decldiagparam("Sd_ByPass", TDim::NHRU, "[0]", "0", "1","0 - normal, 1 - Bypass Pond/Depressional storage (i.e. Sd).", "()", &Sd_ByPass); + + decldiagparam("soil_rechr_ByPass", TDim::NHRU, "[1]", "0", "1","0 - normal, 1 - Bypass recharge layer (i.e. soil_rechr).", "()", &soil_rechr_ByPass); + + declparam("Channel_shp", TDim::NHRU, "[0]", "0", "2", "rectangular - 0/parabolic - 1/triangular - 2", "()", &route_Cshp); + + decldiagparam("scaling_factor", TDim::NHRU, "[1.0]", "0.0", "1.0E6","multiplies the input to Muskingum by this scaling factor.", "()", &scaling_factor); // temporary + + declparam("order", TDim::NHRU, "[1,2,3,4,5!]", "1","1000", "HRU routing process order", "()", &order); + + + soil_gwDiv = declgetvar("*", "gw_flow", "(mm/int)", &soil_gw); + + soil_ssrDiv = declgetvar("*", "soil_ssr", "(mm/int)", &soil_ssr); + + soil_runoffDiv = declgetvar("*", "soil_runoff", "(mm/int)", &soil_runoff); + + + declputvar("*", "Sd", "(mm)", &Sd); + + declputvar("*", "soil_moist", "(mm)", &soil_moist); + + declputvar("*", "soil_rechr", "(mm)", &soil_rechr); + + declputvar("*", "redirected_residual", "(mm*km^2/int)", &redirected_residual); + + declputvar("*", "gw", "(mm)", &gw); +} + +void ClassNetroute_M2_D::init(void) { + + nhru = getdim(TDim::NHRU); + + if(soil_ssrDiv > 1){ + string S = "Netroute_M2_D: \"soil_ssr\". Converting to mm/int"; + CRHMException TExcept(S.c_str(), TExcept::WARNING); + LogError(TExcept); + } + + if(soil_runoffDiv > 1){ + string S = "Netroute_M2_D: \"soil_runoff\". Converting to mm/int"; + CRHMException TExcept(S.c_str(), TExcept::WARNING); + LogError(TExcept); + } + + if(soil_gwDiv > 1){ + string S = "Netroute_M2_D: \"gw_flow\". Converting to mm/int"; + CRHMException TExcept(S.c_str(), TExcept::WARNING); + LogError(TExcept); + } + +// Refer to EM 1110-2-1417 for these channel Vw/V definitions + const double Vw[3] = {1.67, 1.44, 1.33, 1.5}; // rectangular - 0/parabolic - 1/triangular - 2/natural - 3 + + for(hh = 0; hh < nhru; ++hh){ + double Vavg = (1.0/route_n[hh])*pow(route_R[hh], 2.0/3.0)*pow(route_S0[hh], 0.5f); // (m/s) + Ktravel[hh] = route_L[hh]/(Vw[route_Cshp[hh]]*Vavg)/86400.0; // (d) + } + + hruDelay = new ClassMuskingum2(inflow, outflow, Ktravel, route_X_M, Lag, nhru); + ssrDelay = new ClassClark(ssrinflow, ssroutflow, ssrKstorage, ssrLag, nhru); + runDelay = new ClassClark(runinflow, runoutflow, runKstorage, runLag, nhru); + gwDelay = new ClassClark(gwinflow, gwoutflow, gwKstorage, gwLag, nhru); + + for(hh = 0; hh < nhru; ++hh){ + if(Ktravel[hh] >= (Global::Interval/(2.0*route_X_M[hh]))){ + string S = string("'" + Name + " (Netroute_M2_D) Muskingum coefficient negative in HRU ").c_str() + to_string(hh+1); + CRHMException TExcept(S.c_str(), TExcept::WARNING); + LogError(TExcept); + } + + if(Ktravel[hh] < (Global::Interval/(2.0*(1.0-route_X_M[hh])))){ // if(hruDelay->c0[hh] < 0.0) + hruDelay->c0[hh] = 0.0; + hruDelay->c1[hh] = 1.0; + hruDelay->c2[hh] = 0.0; + } + } + + basinflow[0] = 0.0; + basinflow_s[0] = 0.0; + cumbasinflow[0] = 0.0; + basingw[0] = 0.0; + basingw_s[0] = 0.0; + cumbasingw[0] = 0.0; + + for(hh = 0; hh < nhru; ++hh) { + inflow[hh] = 0.0; + cuminflow[hh] = 0.0; + outflow[hh] = 0.0; + outflow_diverted[hh] = 0.0; + cumoutflow_diverted[hh] = 0.0; + cumoutflow[hh] = 0.0; + cumscaling_boost[hh] = 0.0; + cum_to_Sd[hh] = 0.0; + cum_to_soil_rechr[hh] = 0.0; + cum_preferential_flow_to_gw[hh] = 0.0; + + gwinflow[hh] = 0.0; + gwcuminflow[hh] = 0.0; + gwcumoutflow[hh] = 0.0; + gwcumoutflow_diverted[hh] = 0.0; + HRU_cumbasinflow[hh] = 0.0; + + ssrinflow[hh] = 0.0; + ssrcuminflow[hh] = 0.0; + ssroutflow[hh] = 0.0; + ssrcumoutflow[hh] = 0.0; + + runinflow[hh] = 0.0; + runcuminflow[hh] = 0.0; + runoutflow[hh] = 0.0; + runcumoutflow[hh] = 0.0; + + soil_ssr_Buf[hh] = 0.0; + soil_runoff_Buf[hh] = 0.0; + soil_gw_Buf[hh] = 0.0; + + bool OK = false; + for(long jj = 0; chkStruct(jj); ++jj) + if(order[jj] - 1 == hh){ + OK = true; + break; + } + + if(!OK){ + string SS = string("'" + Name + " (Netroute_M2_D)' the 'order' parameter does not have a unique value for each HRU"); + CRHMException Except(SS.c_str() , TExcept::ERR); + LogError(Except); + throw Except; + } + } +} + +void ClassNetroute_M2_D::run(void) { + + long nstep = getstep(); + nstep = nstep%Global::Freq; + + basinflow[0] = 0.0; + basingw[0] = 0.0; + + if(getstep() == 1){ + for(hh = 0; chkStruct(hh); ++hh) { // do HRUs in sequence. + distrib_sum[hh] = 0.0; + + for(long hhh = 0; chkStruct(hhh); ++hhh) { // do HRUs in sequence + if(distrib_hru[hh][hhh] < 0.0) + const_cast (distrib_hru) [hh][hhh] = -distrib_hru[hh][hhh]*hru_area[hh]; + distrib_sum[hh] += distrib_hru[hh][hhh]; + } + + if(distrib_sum[hh] <= 0 && distrib_Basin[hh] <= 0.0){ + const_cast (distrib_Basin) [hh] = 1; + } + + distrib_sum[hh] += distrib_Basin[hh]; + } + } + + double gw_amount; + + for(long jj = 0; chkStruct(jj); ++jj){ // do HRUs in sequence + + for(hh = 0; chkStruct(hh); ++hh){ + if(order[hh] - 1 == jj) + break; + } + + if(soil_gwDiv == 1) // interval value + soil_gw_Buf[hh] = soil_gw[hh]; + + gwinflow[hh] = soil_gw_Buf[hh]*hru_area[hh]; + + gwoutflow_diverted[hh] = 0.0; + + gw_amount = 0.0; + + for(long hhh = 0; chkStruct(hhh); ++hhh) { + if(gwoutflow[hhh] > 0.0 && gwwhereto[hhh] && (abs(gwwhereto[hhh])-1 == hh || abs(gwwhereto[hhh]) > nhru)){ // handles "gwwhereto" <> 0 + gwoutflow_diverted[hhh] = gwoutflow[hhh]; + gw_amount = gwoutflow_diverted[hhh]; + gwoutflow[hhh] = 0.0; + gwcumoutflow_diverted[hhh] += gwoutflow_diverted[hhh]; + + if(abs(gwwhereto[hhh]) <= nhru){ + if(gwwhereto[hhh] > 0){ // direct to HRU surface + double free = soil_rechr_max[hh] - soil_rechr[hh]; + if(free > 0.0 && !soil_rechr_ByPass[hh]){ + if(free > gw_amount/hru_area[hh]){ // outflow (mm*km^2/int) + soil_rechr[hh] += gw_amount/hru_area[hh]; + soil_moist[hh] += gw_amount/hru_area[hh]; + cum_to_soil_rechr[hh] += gw_amount/hru_area[hh]; + gw_amount = 0.0; + } + else { + gw_amount = (gw_amount/hru_area[hh] - free)*hru_area[hh]; + cum_to_soil_rechr[hh] += free; + soil_moist[hh] += free; + soil_rechr[hh] = soil_rechr_max[hh]; + } + } + + free = Sdmax[hh] - Sd[hh]; + if(free > 0.0 && gw_amount > 0.0 && !Sd_ByPass[hh]){ + if(free > gw_amount/hru_area[hh]){ // outflow (mm*km^2/int) + Sd[hh] += gw_amount/hru_area[hh]; + cum_to_Sd[hh] += gw_amount/hru_area[hh]; + gw_amount = 0.0; + } + else { + gw_amount = (gw_amount/hru_area[hh] - free)*hru_area[hh]; + cum_to_Sd[hh] += free; + Sd[hh] = Sdmax[hh]; + } + } + } // hh > 0 + else{ // hh < 0 + gw[hh] += gw_amount/hru_area[hh]; + gw_amount = 0.0; + } + } // is HRU + else{ // > nhru then put in basinflow + basinflow[0] += gw_amount*1000; // (m3) + HRU_cumbasinflow[hh] = gw_amount; + cumoutflow[hh] += gw_amount; + gw_amount = 0.0; + } + } // match hh + } // for hhh + + gwcuminflow[hh] += gwinflow[hh]; + //gwcumoutflow[hh] += gwoutflow[hh]; Removing to match borland version of code - jhs507 + + inflow[hh] = gw_amount; + } // for jj + + for(hh = 0; chkStruct(hh); ++hh) { // do HRUs in sequence + + if(soil_ssrDiv == 1) // interval value + soil_ssr_Buf[hh] = soil_ssr[hh]; + + if(soil_runoffDiv == 1) // interval value + soil_runoff_Buf[hh] = soil_runoff[hh]; + + runinflow[hh] = soil_runoff_Buf[hh]*hru_area[hh]; // includes melt and rain runoff + runcuminflow[hh] += runinflow[hh]; + + ssrinflow[hh] = soil_ssr_Buf[hh]*hru_area[hh]; // subsurface runoff + ssrcuminflow[hh] += ssrinflow[hh]; + + ssrcumoutflow[hh] += ssroutflow[hh]; + runcumoutflow[hh] += runoutflow[hh]; + + inflow[hh] += scaling_factor[hh]*(runoutflow[hh] + ssroutflow[hh]); // add this HRU runoff and subsurface flow temporary change 10/07/11 + + if(outflow[hh] > 0.0){ + + double Used = outflow[hh]*distrib_Basin[hh]/distrib_sum[hh]; + if(distrib_Basin[hh] > 0.0){ // direct to basin + + basinflow[0] += Used*1000; // (m3) + HRU_cumbasinflow[hh] += Used; + cumoutflow[hh] += Used; + } + + for(long To = 0; chkStruct(To); ++To) { // distribute outflow of HRUs + + if(hh != To && distrib_hru[hh][To] > 0.0){ + double Amount = (outflow[hh]-Used)/hru_area[To]*distrib_hru[hh][To]/(distrib_sum[hh]-distrib_Basin[hh]); // outflow (mm*km^2/int) + + if(preferential_flow[hh]) { + gw[To] += Amount; + cum_preferential_flow_to_gw[To] += Amount; + Amount = 0.0; + } + else { + double free = soil_rechr_max[To] - soil_rechr[To]; + if(free > 0.0 && !soil_rechr_ByPass[To]){ + if(free > Amount){ // outflow (mm*km^2/int) + soil_rechr[To] += Amount; + soil_moist[To] += Amount; + cum_to_soil_rechr[To] += Amount; + Amount = 0.0; + } + else { + Amount -= free; + cum_to_soil_rechr[To] += free; + soil_moist[To] += free; + soil_rechr[To] = soil_rechr_max[To]; + } + } // if + + free = Sdmax[To] - Sd[To]; + if(free > 0.0 && !Sd_ByPass[To] && Amount > 0.0){ + if(free > Amount){ // outflow (mm*km^2/int) + Sd[To] += Amount; + cum_to_Sd[To] += Amount; + Amount = 0.0; + } + else { + Amount -= free; + cum_to_Sd[To] += free; + Sd[To] = Sdmax[To]; + } + } // if + } // else + + if(Amount > 0.0) + redirected_residual[To] += Amount*hru_area[To]; // Return to module soil. Previously handled by Netroute. (mm*km^2/int) + + } // contribute to this HRU + } // distribute outflow over HRUs + } // outflow > 0.0 + + if(gwwhereto[hh] == 0) { // move to basin gw + basingw[0] += gwoutflow[hh]*1000; // (m3) end of every day + gwcumoutflow[hh] += gwoutflow[hh]; + } + + if(nstep == 0){ // end of every day + if(soil_ssrDiv > 1) // daily value - ready for next day + soil_ssr_Buf[hh] = soil_ssr[hh]/soil_ssrDiv; + + if(soil_runoffDiv > 1) // daily value - ready for next day + soil_runoff_Buf[hh] = soil_runoff[hh]/soil_runoffDiv; + + if(soil_gwDiv > 1) // daily value - ready for next day + soil_gw_Buf[hh] = soil_gw[hh]/soil_gwDiv; + } // end if + + } // for hh + + for(hh = 0; chkStruct(hh); ++hh) { + cuminflow[hh] += inflow[hh]; + cumscaling_boost[hh] += inflow[hh]*(scaling_factor[hh] - 1.0); + + outflow_diverted[hh] = 0.0; + if(distrib_sum[hh] > 0.0){ // does not apply to last HRU + for(long hhh = 0; chkStruct(hhh); ++hhh){ + outflow_diverted[hh] += outflow[hh]*distrib_hru[hh][hhh]/distrib_sum[hh]; + } + } + cumoutflow_diverted[hh] += outflow_diverted[hh]; + } // end for + + hruDelay->DoMuskingum(); + runDelay->DoClark(); + ssrDelay->DoClark(); + gwDelay->DoClark(); + + basinflow_s[0] = basinflow[0]*Global::Freq/86400.0; + basingw_s[0] = basingw[0]*Global::Freq/86400.0; + + cumbasinflow[0] += basinflow[0]; + cumbasingw[0] += basingw[0]; +} + +void ClassNetroute_M2_D::finish(bool good) { + + double Allcuminflow = 0.0; + double Allcumoutflow = 0.0; + double Allcumoutflowdiverted = 0.0; + + double Allgwcuminflow = 0.0; + double Allgwcumoutflow = 0.0; + double Allgwcumoutflowdiverted = 0.0; + + double Allssrcuminflow = 0.0; + double Allssrcumoutflow = 0.0; + double Allruncuminflow = 0.0; + double Allruncumoutflow = 0.0; + + double AllSdcuminflow = 0.0; + double Allrechrcuminflow = 0.0; + + for(hh = 0; chkStruct(); ++hh) { + + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' cuminflow (mm) (mm*km^2) (mm*basin): ").c_str(), cuminflow[hh]/hru_area[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' cumoutflow (mm) (mm*km^2) (mm*basin): ").c_str(), cumoutflow[hh]/hru_area[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' cumoutflow_diverted (mm) (mm*km^2) (mm*basin): ").c_str(), cumoutflow_diverted[hh]/hru_area[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' hruDelay_in_storage (mm) (mm*km^2) (mm*basin): ").c_str(), hruDelay->Left(hh)/hru_area[hh], hru_area[hh], basin_area[0]); + + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' scaling_boost (mm) (mm*km^2) (mm*basin): ").c_str(), cumscaling_boost[hh]/hru_area[hh], hru_area[hh], basin_area[0]); + + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' ssrcuminflow (mm) (mm*km^2) (mm*basin): ").c_str(), ssrcuminflow[hh]/hru_area[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' ssrcumoutflow (mm) (mm*km^2) (mm*basin): ").c_str(), ssrcumoutflow[hh]/hru_area[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' ssrDelay_in_storage (mm) (mm*km^2) (mm*basin): ").c_str(), ssrDelay->Left(hh)/hru_area[hh], hru_area[hh], basin_area[0]); + + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' runoffcuminflow (mm) (mm*km^2) (mm*basin): ").c_str(), runcuminflow[hh]/hru_area[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' runoffcumoutflow (mm) (mm*km^2) (mm*basin): ").c_str(), runcumoutflow[hh]/hru_area[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' runDelay_in_storage (mm) (mm*km^2) (mm*basin): ").c_str(), runDelay->Left(hh)/hru_area[hh], hru_area[hh], basin_area[0]); + + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' gwcuminflow (mm) (mm*km^2) (mm*basin): ").c_str(), gwcuminflow[hh]/hru_area[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' gwcumoutflow (mm) (mm*km^2) (mm*basin): ").c_str(), gwcumoutflow[hh]/hru_area[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' gwcumdiverted (mm) (mm*km^2) (mm*basin): ").c_str(), gwcumoutflow_diverted[hh]/hru_area[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' gwDelay_in_storage (mm) (mm*km^2) (mm*basin): ").c_str(), gwDelay->Left(hh)/hru_area[hh], hru_area[hh], basin_area[0]); + + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' cum_to_Sd (mm) (mm*km^2) (mm*basin): ").c_str(), cum_to_Sd[hh], hru_area[hh], basin_area[0], " *** Added to this HRU Sd"); + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' cum_to_soil_rechr (mm) (mm*km^2) (mm*basin): ").c_str(), cum_to_soil_rechr[hh], hru_area[hh], basin_area[0], " *** Added to this HRU recharge"); + LogMessageA(hh, string("'" + Name + " (Netroute_M2_D)' HRU_cumbasinflow (mm) (mm*km^2) (mm*basin): ").c_str(), HRU_cumbasinflow[hh]/hru_area[hh], hru_area[hh], basin_area[0]); + LogDebug(" "); + + Allcuminflow += cuminflow[hh]; + Allcumoutflow += cumoutflow[hh]; + Allcumoutflowdiverted += cumoutflow_diverted[hh]; + + Allgwcuminflow += gwcuminflow[hh]; + Allgwcumoutflow += gwcumoutflow[hh]; + Allgwcumoutflowdiverted += gwcumoutflow_diverted[hh]; + + Allssrcumoutflow += ssrcumoutflow[hh]; + Allssrcuminflow += ssrcuminflow[hh]; + Allruncuminflow += runcuminflow[hh]; + Allruncumoutflow += runcumoutflow[hh]; + + AllSdcuminflow += cum_to_Sd[hh]*hru_area[hh]; + Allrechrcuminflow += cum_to_soil_rechr[hh]*hru_area[hh]; + } + + LogMessage(string("'" + Name + " (Netroute_M2_D)' cumulative surface & subsurface basinflow(m^3): ").c_str(), cumbasinflow[0]); + LogMessage(string("'" + Name + " (Netroute_M2_D)' cumulative basingw (m^3): ").c_str(), cumbasingw[0]); + LogDebug(" "); + + LogMessage(string("'" + Name + " (Netroute_M2_D)' Allgwcuminflow (mm*basin): ").c_str(), Allgwcuminflow); + LogMessage(string("'" + Name + " (Netroute_M2_D)' Allgwcumoutflow (mm*basin): ").c_str(), Allgwcumoutflow); + LogMessage(string("'" + Name + " (Netroute_M2_D)' Allgwcumoutflowdiverted (mm*basin): ").c_str(), Allgwcumoutflowdiverted); + LogDebug(" "); + + LogMessage(string("'" + Name + " (Netroute_M2_D)' Allcuminflow (mm*basin): ").c_str(), Allcuminflow); + LogMessage(string("'" + Name + " (Netroute_M2_D)' Allcumoutflow (mm*basin): ").c_str(), Allcumoutflow); + LogMessage(string("'" + Name + " (Netroute_M2_D)' Allcumoutflowdiverted (mm*basin): ").c_str(), Allcumoutflowdiverted); + LogDebug(" "); + + LogMessage(string("'" + Name + " (Netroute_M2_D)' Allssrcuminflow (mm*basin): ").c_str(), Allssrcuminflow); + LogMessage(string("'" + Name + " (Netroute_M2_D)' Allssrcumoutflow (mm*basin): ").c_str(), Allssrcumoutflow); + LogDebug(" "); + + LogMessage(string("'" + Name + " (Netroute_M2_D)' Allruncuminflow (mm*basin): ").c_str(), Allruncuminflow); + LogMessage(string("'" + Name + " (Netroute_M2_D)' Allruncumoutflow (mm*basin): ").c_str(), Allruncumoutflow); + LogDebug(" "); + + LogMessage(string("'" + Name + " (Netroute_M2_D)' AllSdcuminflow (mm*basin): ").c_str(), AllSdcuminflow); + LogMessage(string("'" + Name + " (Netroute_M2_D)' Allrechrcuminflow (mm*basin): ").c_str(), Allrechrcuminflow); + LogDebug(" "); + + delete hruDelay; + delete ssrDelay; + delete runDelay; + delete gwDelay; +} diff --git a/crhmcode/src/modules/ClassNetroute_M2_D.h b/crhmcode/src/modules/ClassNetroute_M2_D.h new file mode 100644 index 000000000..b683bbed0 --- /dev/null +++ b/crhmcode/src/modules/ClassNetroute_M2_D.h @@ -0,0 +1,119 @@ +//created by Manishankar Mondal + +#include "../core/ClassModule.h" + +class ClassNetroute_M2_D : public ClassModule { +public: + +ClassNetroute_M2_D(string Name, string Version = "undefined", LMODULE Lvl = LMODULE::PROTO) : ClassModule(Name, Version, Lvl) {}; + +long meltrunoffDiv{0}; +long soil_ssrDiv{0}; +long soil_runoffDiv{0}; +long soil_gwDiv{0}; + +// declared variables +double *inflow{ NULL }; // [nhru] +double *cuminflow{ NULL }; // [nhru] +double *outflow{ NULL }; // [nhru] +double *outflow_diverted{ NULL }; // [nhru] +double *cumoutflow_diverted{ NULL }; // [nhru] +double *cumoutflow{ NULL }; // [nhru] +double *cumscaling_boost{ NULL }; // [nhru] +double *gwinflow{ NULL }; // [nhru] +double *gwoutflow_diverted{ NULL }; // [nhru] +double *gwcumoutflow_diverted{ NULL }; // [nhru] +double *HRU_cumbasinflow{ NULL }; // [nhru] + +double *ssrinflow{ NULL }; // [nhru] +double *ssrcuminflow{ NULL }; // [nhru] +double *ssroutflow{ NULL }; // [nhru] +double *ssrcumoutflow{ NULL }; // [nhru] + +double *runinflow{ NULL }; // [nhru] +double *runcuminflow{ NULL }; // [nhru] +double *runoutflow{ NULL }; // [nhru] +double *runcumoutflow{ NULL }; // [nhru] + +double *gwoutflow{ NULL }; // [nhru] +double *gwcuminflow{ NULL }; // [nhru] +double *gwcumoutflow{ NULL }; // [nhru] + +double *basinflow{ NULL }; // [BASIN] all HRUs +double *basinflow_s{ NULL }; // [BASIN] all HRUs +double *cumbasinflow{ NULL }; // [BASIN] all HRUs +double *basingw{ NULL }; // [BASIN} all HRUs +double *basingw_s{ NULL }; // [BASIN} all HRUs +double *cumbasingw{ NULL }; // [BASIN} all HRUs + +double *soil_ssr_Buf{ NULL }; // buffered +double *soil_runoff_Buf{ NULL }; // buffered +double *soil_gw_Buf{ NULL }; // buffered + +double *cum_to_Sd{ NULL }; // [nhru] +double *cum_to_soil_rechr{ NULL }; // [nhru] + +double *Ktravel{ NULL }; // [nhru] Muskingum + +double *distrib_sum{ NULL }; +double *cum_preferential_flow_to_gw{ NULL }; + +ClassMuskingum *hruDelay{ NULL }; +ClassClark *ssrDelay{ NULL }; +ClassClark *runDelay{ NULL }; +ClassClark *gwDelay{ NULL }; + +// declared parameters +const double *route_n{ NULL }; // [nhru] +const double *route_R{ NULL }; // [nhru] +const double *route_S0{ NULL }; // [nhru] +const double *route_L{ NULL }; // [nhru] +const double *route_X_M{ NULL }; // [nhru] +const long *route_Cshp{ NULL }; // [nhru] + +const double *Lag{ NULL }; // [nhru] +const double *ssrKstorage{ NULL }; // [nhru] +const double *ssrLag{ NULL }; // [nhru] +const double *runKstorage{ NULL }; // [nhru] +const double *runLag{ NULL }; // [nhru] +const double *gwKstorage{ NULL }; // [nhru] +const double *gwLag{ NULL }; // [nhru] +const long *gwwhereto{ NULL }; // [nhru] + +const double *basin_area{ NULL }; // [BASIN] +const double *hru_area{ NULL }; // [nhru] +const double *distrib{ NULL }; +const double *distrib_Basin{ NULL }; +const double **distrib_hru{ NULL }; +const double *Sdmax{ NULL }; // [nhru] +const double *soil_rechr_max{ NULL }; // [nhru] +const long *Sd_ByPass{ NULL }; // [nhru] +const long *soil_rechr_ByPass{ NULL }; // [nhru] +const long *order{ NULL }; +const long *preferential_flow{ NULL }; // [nhru] + +const double *scaling_factor{ NULL }; // temporary modification + +// variable inputs +const double *soil_gw{ NULL }; // [nhru] +const double *soil_ssr{ NULL }; // [nhru] +const double *soil_runoff{ NULL }; // [nhru] + +// variable puts +double *Sd{ NULL }; +double *soil_moist{ NULL }; +double *soil_rechr{ NULL }; +double *redirected_residual{ NULL }; +double *gw{ NULL }; + +// local allocated arrays + +void decl(void); +void init(void); +void run(void); +void finish(bool good); +virtual double Function1(double *I, long hh); +virtual double Function2(double *X, long hh); + +ClassNetroute_M2_D* klone(string name) const; +}; From b579f2ec01ac9d36a3a78d6552ceddeab539c29c Mon Sep 17 00:00:00 2001 From: lawfordp2017 Date: Mon, 14 Feb 2022 21:15:41 -0600 Subject: [PATCH 2/5] New muskingum now compiles --- crhmcode/Makefile | 4 +- crhmcode/src/core/ClassModule.h | 1 + crhmcode/src/core/ClassMuskingum2.cpp | 161 +++++------------- crhmcode/src/core/ClassMuskingum2.h | 18 +- crhmcode/src/modules/ClassNetroute_M2_D.cpp | 6 +- crhmcode/src/modules/ClassNetroute_M2_D.h | 4 +- .../src/modules/newmodules/NewModules.cpp | 2 + 7 files changed, 62 insertions(+), 134 deletions(-) diff --git a/crhmcode/Makefile b/crhmcode/Makefile index c4caf51be..3f329d52a 100644 --- a/crhmcode/Makefile +++ b/crhmcode/Makefile @@ -30,7 +30,7 @@ Classramp.o Classrandom.o Classrefwind.o ClassReplace.o Classrh.o \ ClassRH_WtoI.o ClassSim.o Classsin.o ClassSmear.o Classsquare.o Classsub.o \ ClasssubV.o Classtime.o ClassTimeshift.o -ClassModule_OBJS = Administer.o ClassClark.o ClassModule.o ClassMuskingum.o Myparser.o +ClassModule_OBJS = Administer.o ClassClark.o ClassModule.o ClassMuskingum.o ClassMuskingum2.o Myparser.o Common_OBJS = Common.o snowcover.o @@ -61,7 +61,7 @@ ClassGrow_Crop.o ClassHMSA.o ClassHtobs.o ClassIceBulb.o ClassICEflow.o \ ClassIntcp.o Classinterception.o ClassK_Estimate.o ClassKevin.o Classlake.o \ ClassLongVt.o ClassMeltRunoff_Kstorage.o ClassMeltRunoff_Lag.o ClassMod_Exec.o \ ClassNeedle.o Classnetall.o ClassNetroute.o ClassNetroute_D.o ClassNetroute_M.o \ -ClassNetroute_M_D.o ClassNO_pbsm.o ClassNOP.o ClassObs.o ClassObstoPar.o \ +ClassNetroute_M_D.o ClassNetroute_M2_D.o ClassNO_pbsm.o ClassNOP.o ClassObs.o ClassObstoPar.o \ Classpbsm.o Classpbsm_M.o ClasspbsmSnobal.o ClassPrairieInfil.o ClassPSPnew.o \ Classqdrift.o Classqmelt.o Classquinton.o ClassREWroute.o ClassREWroute_stream.o \ ClassREWroute2.o Classsbsm.o ClassSetSoil.o Classshared.o ClassShutWall.o \ diff --git a/crhmcode/src/core/ClassModule.h b/crhmcode/src/core/ClassModule.h index 6bfe73c72..630449df1 100644 --- a/crhmcode/src/core/ClassModule.h +++ b/crhmcode/src/core/ClassModule.h @@ -17,6 +17,7 @@ #include "Myparser.h" #include "ClassClark.h" #include "ClassMuskingum.h" +#include "ClassMuskingum2.h" diff --git a/crhmcode/src/core/ClassMuskingum2.cpp b/crhmcode/src/core/ClassMuskingum2.cpp index f56a84463..2f0c2028f 100644 --- a/crhmcode/src/core/ClassMuskingum2.cpp +++ b/crhmcode/src/core/ClassMuskingum2.cpp @@ -1,8 +1,11 @@ #include "ClassMuskingum2.h" +#include // std::max +#include // std::exception //--------------------------------------------------------------------------- -ClassMuskingum2::ClassMuskingum2(const double* inVar, double* outVar, const double* k, const double* X_M, const double* lag, const long nhru, const long setlag) +ClassMuskingum2::ClassMuskingum2(const double* inVar, double* outVar, + const double* k, const double* X_M, const long nhru) : inVar(inVar), outVar(outVar), nhru(nhru) { // !!! UNITS !!! @@ -10,171 +13,97 @@ ClassMuskingum2::ClassMuskingum2(const double* inVar, double* outVar, const doub // kstorage (days) // lag (hours) - LastIn = new double[nhru]; - LastOut = new double[nhru]; - c0 = new double[nhru]; c1 = new double[nhru]; c2 = new double[nhru]; /* - * ilag for routing water along a stream, it is used to simulate the delay caused by the travel time of water. - * ilag is the length of the stream (the number of elements or letters). * If on every cycle the water moves one step then the number of elements determines how many steps * are required to get to the end of the stream. */ - ilag = new long[nhru]; - maxlag = new long[nhru]; - ulag = new long[nhru]; - long Biggest = 0; - - for (long hh = 0; hh < nhru; hh++) { - c0[hh] = (Global::Interval - 2.0 * k[hh] * X_M[hh]) / - (2.0 * k[hh] * (1.0 - X_M[hh]) + Global::Interval); // units of Global::Interval (days) + buff_sz = new long[nhru]; + K_sub = new double[nhru]; + x_atten = new double[nhru]; - c1[hh] = (Global::Interval + 2.0 * k[hh] * X_M[hh]) / - (2.0 * k[hh] * (1.0 - X_M[hh]) + Global::Interval); // units of kstorage (days) - c2[hh] = (2.0 * k[hh] * (1.0 - X_M[hh]) - Global::Interval) / - (2.0 * k[hh] * (1.0 - X_M[hh]) + Global::Interval); // units of kstorage (days) + for (long hh = 0; hh < nhru; hh++) { - ilag[hh] = (long)(max(lag[hh], 0.0) / 24.0 * Global::Freq + 1.1); // =1 for lag of zero +// Add 1 to buff_sz because the first buffer location just holds the inflow values +// actual outflow values start from the second buffer location + buff_sz[hh] = 1 + std::max(1l, lround(k[hh]/Global::Interval)); + K_sub[hh] = k[hh] / buff_sz[hh]; + x_atten[hh] = X_M[hh]; - if (setlag == -1 || ilag[hh] > setlag) - maxlag[hh] = ilag[hh]; - else - maxlag[hh] = setlag; - ulag[hh] = 0; + c0[hh] = (Global::Interval - 2.0 * K_sub[hh] * X_M[hh]) / + (2.0 * K_sub[hh] * (1.0 - X_M[hh]) + Global::Interval); // units of Global::Interval (days) - LastIn[hh] = 0.0; // zero initial conditions + c1[hh] = (Global::Interval + 2.0 * K_sub[hh] * X_M[hh]) / + (2.0 * K_sub[hh] * (1.0 - X_M[hh]) + Global::Interval); // units of kstorage (days) - LastOut[hh] = 0.0; // zero initial conditions + c2[hh] = (2.0 * K_sub[hh] * (1.0 - X_M[hh]) - Global::Interval) / + (2.0 * K_sub[hh] * (1.0 - X_M[hh]) + Global::Interval); // units of kstorage (days) - if (maxlag[hh] > Biggest) Biggest = maxlag[hh]; } - LagArray = new double* [nhru]; // create lag array + buff_q = new double* [nhru]; // create mem storage for segment discharge for (long hh = 0; hh < nhru; hh++) { - LagArray[hh] = new double[maxlag[hh]]; - for (long jj = 0; jj < maxlag[hh]; jj++) - LagArray[hh][jj] = 0.0; + buff_q[hh] = new double[buff_sz[hh]]; + for (long jj = 0; jj < buff_sz[hh]; jj++) + buff_q[hh][jj] = 0.0; } } ClassMuskingum2::~ClassMuskingum2() { - delete[] LastIn; - delete[] LastOut; delete[] c0; delete[] c1; delete[] c2; - delete[] ilag; - delete[] maxlag; - delete[] ulag; + delete[] buff_sz; + delete[] K_sub; + delete[] x_atten; for (long hh = 0; hh < nhru; hh++) - delete[] LagArray[hh]; - delete[] LagArray; + delete[] buff_q[hh]; + delete[] buff_q; } void ClassMuskingum2::ChangeLag(const double* newlag, const long hh) { - - long newilag = (long)(max(newlag[hh], 0.0) / 24.0 * Global::Freq + 1.1); // =1 for lag of zero - - if (newilag == ilag[hh]) { - return; - } - - double* AccArray = new double[ilag[hh]]; // work area for ChangeLag - - AccArray[0] = 0.0; - - for (int ii = 1; ii < ilag[hh]; ++ii) - { - AccArray[ii] = AccArray[ii - 1] + LagArray[hh][(ulag[hh] + ii) % ilag[hh]]; // accumulate storage - } - - delete[] LagArray[hh]; // delete previous length - - LagArray[hh] = new double[newilag]; // create new length - - ulag[hh] = 0; // next input value save here. - LagArray[hh][0] = 0.0; // looks better - - double LastValue = 0.0; - - for (int mm = 1; mm < newilag - 1; ++mm) - { - double Y = double(mm) / ((long long)newilag - 1ll) * ((long long)ilag[hh] - 1ll); - int Yint = (int)(Y + 0.0001); - if ((Yint + 1) > ilag[hh] - 1) - { - CRHMException Except("Attempting to read out of bounds array address", TExcept::TERMINATE); - LogError(Except); - throw(Except); - } - double Ydif = Y - Yint; - - - double NewValue = AccArray[Yint] + Ydif * (AccArray[Yint + 1] - AccArray[Yint]); - - - - - LagArray[hh][(ulag[hh] + mm) % newilag] = NewValue - LastValue; - - LastValue = NewValue; - } - - LagArray[hh][(ulag[hh] + newilag - 1) % newilag] = AccArray[ilag[hh] - 1] - LastValue; // final values - - delete[] AccArray; // work area for ChangeLag - - ilag[hh] = newilag; // assign new lag + // This is only ever used by ClassMeltRunoff_Lag + throw std::runtime_error("Not Implemented"); } void ClassMuskingum2::DoMuskingum() { for (long hh = 0; hh < nhru; hh++) { - - LagArray[hh][ulag[hh]] = inVar[hh]; - - ulag[hh] = ++ulag[hh] % ilag[hh]; - - outVar[hh] = c0[hh] * LagArray[hh][ulag[hh]] + c1[hh] * LastIn[hh] + c2[hh] * LastOut[hh]; - - LastIn[hh] = LagArray[hh][ulag[hh]]; - - LastOut[hh] = outVar[hh]; + DoMuskingum(hh); } } void ClassMuskingum2::DoMuskingum(const long hh) { - LagArray[hh][ulag[hh]] = inVar[hh]; - - ulag[hh] = ++ulag[hh] % ilag[hh]; // now points to fully delayed value - - outVar[hh] = c0[hh] * LagArray[hh][ulag[hh]] + c1[hh] * LastIn[hh] + c2[hh] * LastOut[hh]; - - LastIn[hh] = LagArray[hh][ulag[hh]]; - - LastOut[hh] = outVar[hh]; + double q_prev = buff_q[hh][0]; + buff_q[hh][0] = inVar[hh]; + for (int i=1; i Date: Tue, 15 Feb 2022 01:59:18 -0600 Subject: [PATCH 3/5] Validated new Muskingum implementation --- crhmcode/src/core/ClassMuskingum2.cpp | 2 +- crhmcode/src/modules/ClassNetroute_M2_D.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/crhmcode/src/core/ClassMuskingum2.cpp b/crhmcode/src/core/ClassMuskingum2.cpp index 2f0c2028f..3ffb2d8f4 100644 --- a/crhmcode/src/core/ClassMuskingum2.cpp +++ b/crhmcode/src/core/ClassMuskingum2.cpp @@ -32,7 +32,7 @@ ClassMuskingum2::ClassMuskingum2(const double* inVar, double* outVar, // Add 1 to buff_sz because the first buffer location just holds the inflow values // actual outflow values start from the second buffer location buff_sz[hh] = 1 + std::max(1l, lround(k[hh]/Global::Interval)); - K_sub[hh] = k[hh] / buff_sz[hh]; + K_sub[hh] = k[hh] / ( buff_sz[hh] - 1); // The 1 additional should not be included here x_atten[hh] = X_M[hh]; diff --git a/crhmcode/src/modules/ClassNetroute_M2_D.cpp b/crhmcode/src/modules/ClassNetroute_M2_D.cpp index e45773c5d..736be467a 100644 --- a/crhmcode/src/modules/ClassNetroute_M2_D.cpp +++ b/crhmcode/src/modules/ClassNetroute_M2_D.cpp @@ -148,7 +148,7 @@ void ClassNetroute_M2_D::decl(void) { decldiagparam("soil_rechr_ByPass", TDim::NHRU, "[1]", "0", "1","0 - normal, 1 - Bypass recharge layer (i.e. soil_rechr).", "()", &soil_rechr_ByPass); - declparam("Channel_shp", TDim::NHRU, "[0]", "0", "2", "rectangular - 0/parabolic - 1/triangular - 2", "()", &route_Cshp); + declparam("Channel_shp", TDim::NHRU, "[0]", "0", "3", "rectangular - 0/parabolic - 1/triangular - 2/natural - 3", "()", &route_Cshp); decldiagparam("scaling_factor", TDim::NHRU, "[1.0]", "0.0", "1.0E6","multiplies the input to Muskingum by this scaling factor.", "()", &scaling_factor); // temporary From 16a23cd2e40b8255448efbabba4bdb7a4a8657b8 Mon Sep 17 00:00:00 2001 From: jhs507 Date: Wed, 16 Feb 2022 16:21:54 -0600 Subject: [PATCH 4/5] Updated CMakeLists.txt to compile objects for ClassMuskingum2.cpp and ClassNetroute_M2_D.cpp --- crhmcode/src/core/CMakeLists.txt | 1 + crhmcode/src/modules/CMakeLists.txt | 1 + 2 files changed, 2 insertions(+) diff --git a/crhmcode/src/core/CMakeLists.txt b/crhmcode/src/core/CMakeLists.txt index bcc40cc4b..faac0a9ac 100644 --- a/crhmcode/src/core/CMakeLists.txt +++ b/crhmcode/src/core/CMakeLists.txt @@ -27,6 +27,7 @@ add_library(core OBJECT ClassModel.cpp ClassModule.cpp ClassMuskingum.cpp + ClassMuskingum2.cpp ClassPar.cpp ClassVar.cpp Common.cpp diff --git a/crhmcode/src/modules/CMakeLists.txt b/crhmcode/src/modules/CMakeLists.txt index 9fb57ee58..e5e97415d 100644 --- a/crhmcode/src/modules/CMakeLists.txt +++ b/crhmcode/src/modules/CMakeLists.txt @@ -55,6 +55,7 @@ add_library(modules OBJECT Classnetall.cpp ClassNetroute_D.cpp ClassNetroute_M_D.cpp + ClassNetroute_M2_D.cpp ClassNetroute_M.cpp ClassNetroute.cpp ClassNO_pbsm.cpp From badbca0cac6103fb07fe3e39880402504c599b1f Mon Sep 17 00:00:00 2001 From: jhs507 Date: Wed, 16 Feb 2022 16:22:40 -0600 Subject: [PATCH 5/5] Updated visual studio solution to include ClassNetroute_M2_D.cpp, ClassNetroute_M2_D.h, ClassMuskingum2.cpp, and ClassMuskingum2.h --- crhmcode/vcc/CRHM_GUI.vcxproj | 6 +++++- crhmcode/vcc/CRHM_GUI.vcxproj.filters | 12 ++++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/crhmcode/vcc/CRHM_GUI.vcxproj b/crhmcode/vcc/CRHM_GUI.vcxproj index 82ce638e6..860fd8dac 100644 --- a/crhmcode/vcc/CRHM_GUI.vcxproj +++ b/crhmcode/vcc/CRHM_GUI.vcxproj @@ -205,6 +205,7 @@ + @@ -332,6 +333,7 @@ + @@ -424,6 +426,7 @@ + @@ -591,6 +594,7 @@ + @@ -1133,4 +1137,4 @@ - + \ No newline at end of file diff --git a/crhmcode/vcc/CRHM_GUI.vcxproj.filters b/crhmcode/vcc/CRHM_GUI.vcxproj.filters index 705dd71f5..f8e6da513 100644 --- a/crhmcode/vcc/CRHM_GUI.vcxproj.filters +++ b/crhmcode/vcc/CRHM_GUI.vcxproj.filters @@ -694,6 +694,12 @@ Source Files + + Source Files + + + Source Files + @@ -2829,5 +2835,11 @@ Resource Files\Header Files + + Resource Files\Header Files + + + Resource Files\Header Files + \ No newline at end of file