diff --git a/src/bmi/ControlFile.cxx b/src/bmi/ControlFile.cxx index 6e4e6d1..dede173 100644 --- a/src/bmi/ControlFile.cxx +++ b/src/bmi/ControlFile.cxx @@ -28,6 +28,8 @@ ueb::ControlFile::~ControlFile() { } void ueb::ControlFile::loadControlFile(std::string const& contrl_file) { + _conFilename = contrl_file; + std::stringstream std_ss(""); std_ss << "Control File is " << contrl_file << std::endl; @@ -310,6 +312,24 @@ int ueb::ControlFile::getStepsInADay() const { return nstepinaDay; } +void ueb::ControlFile::overrideModelTiming( + int startYear, int startMonth, int startDay, double startHour, + int endYear, int endMonth, int endDay, double endHour, + double dt_hours +) { + if (dt_hours <= 0.0) { + throw std::invalid_argument("UEB ControlFile overrideModelTiming received non-positive dt_hours"); + } + + _ModelStartDate = {startYear, startMonth, startDay}; + _ModelStartHour = startHour; + + _ModelEndDate = {endYear, endMonth, endDay}; + _ModelEndHour = endHour; + + _ModelDt = dt_hours; +} + std::ostream& operator<<(std::ostream& os, ueb::ControlFile f) { // operator<< os << "Control file: " << f._conFilename << std::endl; os << "Parameter file: " << f._paramFilename << std::endl; diff --git a/src/bmi/ControlFile.hxx b/src/bmi/ControlFile.hxx index a84f865..cbd6db1 100644 --- a/src/bmi/ControlFile.hxx +++ b/src/bmi/ControlFile.hxx @@ -106,6 +106,9 @@ class ControlFile { int getStepsInADay() const; + void overrideModelTiming(int startYear, int startMonth, int startDay, double startHour, + int endYear, int endMonth, int endDay, double endHour, + double dt_hours); friend std::ostream& ::operator<<(std::ostream& os, ControlFile f); }; diff --git a/src/bmi/bmi_ueb.cxx b/src/bmi/bmi_ueb.cxx index 9f95ccd..16cf0aa 100644 --- a/src/bmi/bmi_ueb.cxx +++ b/src/bmi/bmi_ueb.cxx @@ -20,6 +20,9 @@ namespace { const auto SERIALIZATION_STATE = "serialization_state"; const auto SERIALIZATION_FREE = "serialization_free"; const auto RESET_TIME = "reset_time"; + const auto NGEN_REALIZATION_START_TIME = "ngen_realization_start_time"; + const auto NGEN_REALIZATION_END_TIME = "ngen_realization_end_time"; + const auto NGEN_REALIZATION_DT = "ngen_realization_dt"; } std::stringstream bmi_ueb_ss(""); @@ -36,6 +39,25 @@ void ueb::BmiUEB::Initialize(std::string config_file) { if (config_file.compare("") != 0) { _confile = ControlFile(config_file); + + { + std::stringstream ts_ss; + const auto start_date = _confile.getModelStartDate(); + const auto end_date = _confile.getModelEndDate(); + + ts_ss << "UEB effective control timing after Initialize load:" << std::endl; + ts_ss << " control file: " << _confile.getControlFile() << std::endl; + ts_ss << " start: " + << start_date[0] << "-" << start_date[1] << "-" << start_date[2] + << " " << _confile.getModelStartHour() << std::endl; + ts_ss << " end: " + << end_date[0] << "-" << end_date[1] << "-" << end_date[2] + << " " << _confile.getModelEndHour() << std::endl; + ts_ss << " dt_hours: " << _confile.getModelDt() << std::endl; + ts_ss << " total_timesteps: " << _confile.getModelTotalTimeSteps() << std::endl; + LOG(ts_ss.str(), LogLevel::INFO); + } + _ws = Watershed( _confile.getWatershedFile(), _confile.getWsvarName(), @@ -46,7 +68,6 @@ void ueb::BmiUEB::Initialize(std::string config_file) { _sitevars = SiteVariables(_confile.getSitevarFile()); _forcings = ForcingVariables( _confile.getInputconFile(), - //_forcings = new ForcingVariables( _confile.getInputconFile(), _ws.getActiveCells(_sitevars.getSiteVars().data()), _confile.getWsycorName(), _confile.getWsxcorName() @@ -73,59 +94,56 @@ void ueb::BmiUEB::Initialize(std::string config_file) { int numTotalTs = _confile.getModelTotalTimeSteps(); int nstepinaDay = _confile.getStepsInADay(); + if (numTotalTs <= 0 || nstepinaDay <= 0) { + LOG("UEB Initialize: static config timing is incomplete; using temporary allocation until ngen realization time is supplied through BMI.", LogLevel::INFO); + numTotalTs = 1; + nstepinaDay = 1; + } + + { + std::stringstream ss_info; + ss_info << "UEB allocation sizes:" << std::endl; + ss_info << " numActiveCells=" << numOfActiveCells << std::endl; + ss_info << " numTotalTs=" << numTotalTs << std::endl; + ss_info << " nstepinaDay=" << nstepinaDay << std::endl; + ss_info << " numOut=" << numOut << std::endl; + LOG(ss_info.str(), LogLevel::INFO); + } + _tsprevday.resize(numOfActiveCells); _taveprevday.resize(numOfActiveCells); float Us, Ws, Wc, Apr, cg, rhog, de, tave, WGT; - // WGT=WATER EQUIVALENT GLACIER THICKNESS auto Param = _parms.getParams(); _outvarArray.resize(numOfActiveCells); - for (int cell = 0; cell < numOfActiveCells; ++cell) { - auto sitev = this->getSitevForCell(cell); - auto SiteState = this->getSiteState(cell); - float ts_last = SiteState[30]; + for (int cell = 0; cell < numOfActiveCells; cell++) { + auto sitev = getSitevForCell(cell); - _Ws1[cell] = _statev[cell][1]; - _Wc1[cell] = _statev[cell][3]; + _tsprevday[cell].resize(nstepinaDay, 0.f); + _taveprevday[cell].resize(nstepinaDay, 0.f); + + Us = _statev[cell][0]; + Ws = _statev[cell][1]; + Wc = _statev[cell][3]; + Apr = sitev[1]; + cg = Param[3]; + rhog = Param[7]; + de = Param[10]; if (sitev[9] == 0 || sitev[9] == 3) - WGT = 0.0; + WGT = 0.0f; else - WGT = 1.0; - - if (sitev[9] != 3) // Only do this work for non accumulation cells where model is run - { - _tsprevday[cell].resize(nstepinaDay, -9999.f); - _taveprevday[cell].resize(nstepinaDay, -9999.f); - - // Take surface temperature as 0 where it is unknown the previous time step - // This is for first day of the model to get the force restore going - // #$#$#$#$#_is this all the time steps or the last time? - if (ts_last <= -9999) - _tsprevday[cell][nstepinaDay - 1] = 0; - else - _tsprevday[cell][nstepinaDay - 1] = ts_last; - - // compute Ave.Temp for previous day - Us = _statev[cell][0]; // Ub in UEB - Ws = _statev[cell][1]; // W in UEB - Wc = _statev[cell][3]; // Canopy SWE - Apr = sitev[1]; // Atm. Pressure [PR in UEB] - cg = Param[3]; // Ground heat capacity [nominally 2.09 KJ/kg/C] - rhog = Param[7]; // Soil Density [nominally 1700 kg/m^3] - de = Param[10]; // Thermally active depth of soil (0.1 m) - - tave = TAVG(Us, Ws + WGT, Rho_w, C_s, T_0, rhog, de, cg, H_f); - _taveprevday[cell][nstepinaDay - 1] = tave; - } + WGT = 1.0f; + + tave = TAVG(Us, Ws + WGT, Rho_w, C_s, T_0, rhog, de, cg, H_f); + _taveprevday[cell][nstepinaDay - 1] = tave; _outvarArray[cell].resize(numOut); for (auto& v : _outvarArray[cell]) { #ifdef UEB_SUPPRESS_OUTPUTS - // only need one space when not recording all results v.resize(1); #else v.resize(numTotalTs); @@ -347,6 +365,10 @@ void ueb::BmiUEB::Finalize() { int ueb::BmiUEB::GetVarGrid(std::string name) { int lastgrid = 0; std::array strsvArray = _sitevars.getSiteVars(); + + if (name == NGEN_REALIZATION_START_TIME || name == NGEN_REALIZATION_END_TIME || name == NGEN_REALIZATION_DT) { + return 0; + } for (int i = 0; i < ueb::SiteVariables::nsitevars; ++i) { if (strsvArray[i].svType == 1) { if (name.compare(strsvArray[i].svName) == 0) { @@ -369,7 +391,10 @@ int ueb::BmiUEB::GetVarGrid(std::string name) { } std::string ueb::BmiUEB::GetVarType(std::string name) { - if (name.compare(SERIALIZATION_CREATE) == 0) { + if (name == NGEN_REALIZATION_START_TIME || name == NGEN_REALIZATION_END_TIME || name == NGEN_REALIZATION_DT) { + return "double"; + } + else if (name.compare(SERIALIZATION_CREATE) == 0) { return "uint64_t"; } else if (name.compare(SERIALIZATION_SIZE) == 0) { return "uint64_t"; @@ -380,7 +405,6 @@ std::string ueb::BmiUEB::GetVarType(std::string name) { } else if (name.compare(RESET_TIME) == 0) { return "double"; } - auto it_site = std::find( ueb::SiteVariables::site_var_names.begin(), ueb::SiteVariables::site_var_names.end(), @@ -431,6 +455,12 @@ std::string ueb::BmiUEB::GetVarType(std::string name) { } int ueb::BmiUEB::GetVarItemsize(std::string name) { + if (name.compare(NGEN_REALIZATION_START_TIME) == 0 || + name.compare(NGEN_REALIZATION_END_TIME) == 0 || + name.compare(NGEN_REALIZATION_DT) == 0) { + return sizeof(double); + } + if (name.compare(SERIALIZATION_CREATE) == 0) { return sizeof(uint64_t); } else if (name.compare(SERIALIZATION_SIZE) == 0) { @@ -494,6 +524,12 @@ int ueb::BmiUEB::GetVarItemsize(std::string name) { std::string ueb::BmiUEB::GetVarUnits(std::string name) { auto it_site = ueb::SiteVariables::site_var_units.find(name); + if (name.compare(NGEN_REALIZATION_START_TIME) == 0 || + name.compare(NGEN_REALIZATION_END_TIME) == 0 || + name.compare(NGEN_REALIZATION_DT) == 0) { + return "s"; + } + if (it_site != ueb::SiteVariables::site_var_units.end()) { return it_site->second; } @@ -514,7 +550,11 @@ std::string ueb::BmiUEB::GetVarUnits(std::string name) { } int ueb::BmiUEB::GetVarNbytes(std::string name) { - if (name.compare(SERIALIZATION_CREATE) == 0) { + if (name.compare(NGEN_REALIZATION_START_TIME) == 0 || + name.compare(NGEN_REALIZATION_END_TIME) == 0 || + name.compare(NGEN_REALIZATION_DT) == 0) { + return sizeof(double); + } else if (name.compare(SERIALIZATION_CREATE) == 0) { return sizeof(uint64_t); } else if (name.compare(SERIALIZATION_SIZE) == 0) { return sizeof(uint64_t); @@ -698,6 +738,14 @@ void* ueb::BmiUEB::GetValuePtr(std::string name) { return (void*)this->m_serialized.data(); } + if (name.compare(NGEN_REALIZATION_START_TIME) == 0) { + return (void*)&this->_ngen_realization_start_time; + } else if (name.compare(NGEN_REALIZATION_END_TIME) == 0) { + return (void*)&this->_ngen_realization_end_time; + } else if (name.compare(NGEN_REALIZATION_DT) == 0) { + return (void*)&this->_ngen_realization_dt; + } + auto it_par = std::find( ueb::Parameters::parameter_names.begin(), ueb::Parameters::parameter_names.end(), @@ -862,6 +910,21 @@ void ueb::BmiUEB::GetValueAtIndices(std::string name, void* dest, int* inds, int void ueb::BmiUEB::SetValue(std::string name, void* src) { // special cases for serialized state + + if (name.compare(NGEN_REALIZATION_START_TIME) == 0) { + _ngen_realization_start_time = *static_cast(src); + apply_ngen_realization_time(); + return; + } else if (name.compare(NGEN_REALIZATION_END_TIME) == 0) { + _ngen_realization_end_time = *static_cast(src); + apply_ngen_realization_time(); + return; + } else if (name.compare(NGEN_REALIZATION_DT) == 0) { + _ngen_realization_dt = *static_cast(src); + apply_ngen_realization_time(); + return; + } + if (name.compare(SERIALIZATION_STATE) == 0) { this->load_serialized((char*)src); return; @@ -875,6 +938,7 @@ void ueb::BmiUEB::SetValue(std::string name, void* src) { this->reset_time(); return; } + void* dest = NULL; dest = this->GetValuePtr(name); @@ -915,7 +979,7 @@ int ueb::BmiUEB::GetInputItemCount() { // return Parameters::npar + // NSITEVARS + // NFORCS; - return 8; + return 11; } int ueb::BmiUEB::GetOutputItemCount() { @@ -940,6 +1004,10 @@ std::vector ueb::BmiUEB::GetInputVarNames() { names.push_back("uebv2d"); names.push_back("uebu2d"); + names.push_back(NGEN_REALIZATION_START_TIME); + names.push_back(NGEN_REALIZATION_END_TIME); + names.push_back(NGEN_REALIZATION_DT); + // ngen check all input variables to see if there is // a provider. So here we can't list all input variables. // We will need to a way to get all input values into @@ -1934,6 +2002,129 @@ void ueb::BmiUEB::reset_time() { // indexing into values seems to derive from _currentModelDateTime, so this should be the only thing that needs to be reset } +void ueb::BmiUEB::apply_ngen_realization_time() +{ + if (this->realization_time_applied) { + return; + } + + if (_ngen_realization_start_time <= 0.0 || + _ngen_realization_end_time <= 0.0 || + _ngen_realization_dt <= 0.0) { + return; + } + + time_t start_t = static_cast(_ngen_realization_start_time); + time_t end_t = static_cast(_ngen_realization_end_time); + + struct tm start_tm_struct; + struct tm end_tm_struct; + + gmtime_r(&start_t, &start_tm_struct); + gmtime_r(&end_t, &end_tm_struct); + + int startYear = start_tm_struct.tm_year + 1900; + int startMonth = start_tm_struct.tm_mon + 1; + int startDay = start_tm_struct.tm_mday; + double startHour = start_tm_struct.tm_hour + + start_tm_struct.tm_min / 60.0 + + start_tm_struct.tm_sec / 3600.0; + + int endYear = end_tm_struct.tm_year + 1900; + int endMonth = end_tm_struct.tm_mon + 1; + int endDay = end_tm_struct.tm_mday; + double endHour = end_tm_struct.tm_hour + + end_tm_struct.tm_min / 60.0 + + end_tm_struct.tm_sec / 3600.0; + + double dt_hours = _ngen_realization_dt / 3600.0; + + _confile.overrideModelTiming( + startYear, startMonth, startDay, startHour, + endYear, endMonth, endDay, endHour, + dt_hours + ); + + this->reset_time(); + + int numTotalTs = _confile.getModelTotalTimeSteps(); + int nstepinaDay = _confile.getStepsInADay(); + + if (numTotalTs <= 0) { + throw std::runtime_error("UEB realization time produced numTotalTs <= 0"); + } + if (nstepinaDay <= 0) { + throw std::runtime_error("UEB realization time produced nstepinaDay <= 0"); + } + + auto numOfActiveCells = _ws.getActiveCells(_sitevars.getSiteVars().data()).size(); + + _tsprevday.resize(numOfActiveCells); + _taveprevday.resize(numOfActiveCells); + _outvarArray.resize(numOfActiveCells); + + float Us, Ws, Wc, Apr, cg, rhog, de, tave, WGT; + + auto Param = _parms.getParams(); + + for (int cell = 0; cell < numOfActiveCells; cell++) { + auto sitev = this->getSitevForCell(cell); + auto SiteState = this->getSiteState(cell); + float ts_last = SiteState[30]; + + _Ws1[cell] = _statev[cell][1]; + _Wc1[cell] = _statev[cell][3]; + + if (sitev[9] == 0 || sitev[9] == 3) + WGT = 0.0f; + else + WGT = 1.0f; + + if (sitev[9] != 3) { + _tsprevday[cell].assign(nstepinaDay, -9999.f); + _taveprevday[cell].assign(nstepinaDay, -9999.f); + + if (ts_last <= -9999) + _tsprevday[cell][nstepinaDay - 1] = 0; + else + _tsprevday[cell][nstepinaDay - 1] = ts_last; + + Us = _statev[cell][0]; + Ws = _statev[cell][1]; + Wc = _statev[cell][3]; + Apr = sitev[1]; + cg = Param[3]; + rhog = Param[7]; + de = Param[10]; + + tave = TAVG(Us, Ws + WGT, Rho_w, C_s, T_0, rhog, de, cg, H_f); + _taveprevday[cell][nstepinaDay - 1] = tave; + } + + _outvarArray[cell].resize(numOut); + for (auto& v : _outvarArray[cell]) { + #ifdef UEB_SUPPRESS_OUTPUTS + v.assign(1, 0.0f); + #else + v.assign(numTotalTs, 0.0f); + #endif + } + } + + this->realization_time_applied = true; + + { + std::stringstream ss; + ss << "UEB realization time applied:" << std::endl; + ss << " startdate=" << startYear << "-" << startMonth << "-" << startDay << " " << startHour << std::endl; + ss << " enddate=" << endYear << "-" << endMonth << "-" << endDay << " " << endHour << std::endl; + ss << " dt_hours=" << dt_hours << std::endl; + ss << " total_timesteps=" << numTotalTs << std::endl; + ss << " steps_in_day=" << nstepinaDay << std::endl; + LOG(ss.str(), LogLevel::INFO); + } +} + template void ueb::BmiUEB::serialize(Archive& ar, const unsigned int version) { ar & this->_currentModelDateTime; diff --git a/src/bmi/bmi_ueb.hxx b/src/bmi/bmi_ueb.hxx index 0d9bbef..5ae7e1b 100644 --- a/src/bmi/bmi_ueb.hxx +++ b/src/bmi/bmi_ueb.hxx @@ -108,6 +108,15 @@ class BmiUEB : public bmi::Bmi { OutControl _outcontrol; double _currentModelDateTime; // this is Julian date in days + double realization_start_time = -1; + double realization_end_time = -1; + double realization_dt = -1; + bool realization_time_set = false; + double _ngen_realization_start_time = -1.0; + double _ngen_realization_end_time = -1.0; + double _ngen_realization_dt = -1.0; + + bool realization_time_applied = false; std::vector> _statev; @@ -137,6 +146,8 @@ class BmiUEB : public bmi::Bmi { void setStatev(); + void apply_ngen_realization_time(); + void prepareInputForPoint( double const& UTCHour, // input double const& dHour, // input