Skip to content

fix: NetCDF forcing month-range off-by-one + CSV quantization matching #78

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
DankerMu merged 1 commit into from
Feb 14, 2026
Merged

fix: NetCDF forcing month-range off-by-one + CSV quantization matching #78

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/ModelData/MD_readin.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -846,7 +846,7 @@ void Model_Data::validateTimeStamps()
} else {
for (int i = 0; i < NumForc; i++) {
const double forcMin = tsd_weather[i].getMinTime();
const double forcMax = tsd_weather[i].getMaxTime();
const double forcMax = tsd_weather[i].getMaxTimeCovered();

const bool startCovered = (forcMin - simStart) <= 1e-6;
const bool endCovered = (simEnd - forcMax) <= 1e-6;
Expand Down
78 changes: 69 additions & 9 deletions src/classes/NetcdfForcingProvider.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -778,7 +778,14 @@ struct NetcdfForcingProvider::Impl {
const long long base_days = daysFromCivil(base_y, base_m, base_d);

const long long start_days = base_days + (long long)std::floor(sim_start_min / 1440.0);
const long long end_days = base_days + (long long)std::floor(sim_end_min / 1440.0);
// Simulation END is an exclusive bound for forcing file discovery: the model never
// queries forcing at t == sim_end_min, only for times in [sim_start_min, sim_end_min).
// Avoid requesting the month of the END boundary when it lands exactly on a month/day boundary.
double end_min_excl = sim_end_min;
if (sim_end_min > sim_start_min + 1e-12) {
end_min_excl = std::nextafter(sim_end_min, -std::numeric_limits<double>::infinity());
}
const long long end_days = base_days + (long long)std::floor(end_min_excl / 1440.0);
int y0 = 0, y1 = 0;
unsigned m0 = 0, m1 = 0;
unsigned dtmp = 0;
Expand Down Expand Up @@ -1253,15 +1260,22 @@ struct NetcdfForcingProvider::Impl {
if (!std::isfinite(prec_mm_day)) {
prec_mm_day = 0.0;
}
const double min_mm_day = 0.0001; // match AutoSHUD behavior
if (prec_mm_day < min_mm_day) {
if (prec_mm_day < 0.0) {
prec_mm_day = 0.0;
}
if (prec_mm_day < 0.0) {
// Match AutoSHUD baseline forcing CSV quantization (4 decimals) before thresholding.
prec_mm_day = std::nearbyint(prec_mm_day * 10000.0) / 10000.0;
const double min_mm_day = 0.0001; // match AutoSHUD behavior
if (prec_mm_day < min_mm_day) {
prec_mm_day = 0.0;
}

const double temp_c_val = temp_k - 273.15;
double temp_c_val = temp_k - 273.15;
if (!std::isfinite(temp_c_val)) {
temp_c_val = 0.0;
}
// Match AutoSHUD baseline forcing CSV quantization (2 decimals).
temp_c_val = std::nearbyint(temp_c_val * 100.0) / 100.0;

double rh_percent = 0.263 * pres * shum /
std::exp(17.67 * (temp_k - 273.15) / (temp_k - 29.65));
Expand All @@ -1274,13 +1288,42 @@ struct NetcdfForcingProvider::Impl {
if (rh_percent < 0.0) {
rh_percent = 0.0;
}
const double rh_val = rh_percent / 100.0;
double rh_val = rh_percent / 100.0;
// Match AutoSHUD baseline forcing CSV quantization (4 decimals).
rh_val = std::nearbyint(rh_val * 10000.0) / 10000.0;
if (rh_val > 1.0) {
rh_val = 1.0;
}
if (rh_val < 0.0) {
rh_val = 0.0;
}

double wind_val = std::fabs(wind);
if (!std::isfinite(wind_val)) {
wind_val = 0.0;
}
// Match AutoSHUD baseline forcing CSV quantization (2 decimals) + min clamp.
wind_val = std::nearbyint(wind_val * 100.0) / 100.0;
const double min_wind_ms = 0.05;
if (wind_val < min_wind_ms) {
wind_val = min_wind_ms;
}

double rn_val = srad;
if (!std::isfinite(rn_val)) {
rn_val = 0.0;
}
if (rn_val < 0.0) {
rn_val = 0.0;
}
// Match AutoSHUD baseline forcing CSV quantization (integer).
rn_val = std::nearbyint(rn_val);

prcp_mm_day[i] = prec_mm_day;
temp_c[i] = temp_c_val;
rh_1[i] = rh_val;
wind_ms[i] = std::fabs(wind);
rn_wm2[i] = srad;
wind_ms[i] = wind_val;
rn_wm2[i] = rn_val;
}
loaded_time_idx = t_idx;
}
Expand Down Expand Up @@ -1716,7 +1759,24 @@ struct NetcdfForcingProvider::Impl {
}

double minTimeMin() const { return time_min.empty() ? NA_VALUE : time_min.front(); }
double maxTimeMin() const { return time_min.empty() ? NA_VALUE : time_min.back(); }
double maxTimeMin() const
{
if (time_min.empty()) {
return NA_VALUE;
}
const double t_max = time_min.back();
// Step-function forcing: the last record covers one more interval beyond its timestamp.
// Estimate that final interval width from the last positive dt in the time axis.
double dt_last = 0.0;
for (size_t i = time_min.size(); i > 1; i--) {
const double dt = time_min[i - 1] - time_min[i - 2];
if (dt > 1e-9) {
dt_last = dt;
break;
}
}
return t_max + dt_last;
}
};

NetcdfForcingProvider::NetcdfForcingProvider(const char *forcing_cfg_path,
Expand Down
22 changes: 22 additions & 0 deletions src/classes/TimeSeriesData.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <cmath>
#include <limits>
#include "TimeSeriesData.hpp"

Expand All @@ -16,6 +17,7 @@ _TimeSeriesData::_TimeSeriesData(){
timeRangeCached = 0;
minTime = NA_VALUE;
maxTime = NA_VALUE;
lastDtMin = 0.0;
}

_TimeSeriesData::~_TimeSeriesData(){
Expand All @@ -38,6 +40,7 @@ void _TimeSeriesData::initialize(int n)
timeRangeCached = 0;
minTime = NA_VALUE;
maxTime = NA_VALUE;
lastDtMin = 0.0;
for (int i = 0; i < MAXQUE + 1; i++) {
pRing[i] = i + 1;
}
Expand Down Expand Up @@ -79,6 +82,7 @@ void _TimeSeriesData::computeTimeRange() const

double localMin = std::numeric_limits<double>::infinity();
double localMax = -std::numeric_limits<double>::infinity();
double lastPositiveDtMin = 0.0;

while (std::getline(file, line)) {
lineNo++;
Expand Down Expand Up @@ -118,6 +122,11 @@ void _TimeSeriesData::computeTimeRange() const
myexit(ERRDATAIN);
}

const double dtMin = timeMin - prevTimeMin;
if (dtMin > 1e-12) {
lastPositiveDtMin = dtMin;
}

if (timeMin < localMin) localMin = timeMin;
if (timeMin > localMax) localMax = timeMin;
prevTimeMin = timeMin;
Expand All @@ -133,6 +142,7 @@ void _TimeSeriesData::computeTimeRange() const

minTime = localMin;
maxTime = localMax;
lastDtMin = lastPositiveDtMin;
timeRangeCached = 1;
}

Expand All @@ -148,6 +158,18 @@ double _TimeSeriesData::getMaxTime() const
return maxTime;
}

double _TimeSeriesData::getMaxTimeCovered() const
{
computeTimeRange();
if (!std::isfinite(maxTime)) {
return maxTime;
}
if (!(lastDtMin > 0.0) || !std::isfinite(lastDtMin)) {
return maxTime;
}
return maxTime + lastDtMin;
}

void _TimeSeriesData::read_csv()
{
if (!eof) {
Expand Down
2 changes: 2 additions & 0 deletions src/classes/TimeSeriesData.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ class _TimeSeriesData {
long getStartTime() const;
double getMinTime() const;
double getMaxTime() const;
double getMaxTimeCovered() const;
double xyz[3]={NA_VALUE, NA_VALUE, NA_VALUE};
double lon_deg = NA_VALUE;
double lat_deg = NA_VALUE;
Expand All @@ -44,6 +45,7 @@ class _TimeSeriesData {
mutable int timeRangeCached = 0;
mutable double minTime = NA_VALUE;
mutable double maxTime = NA_VALUE;
mutable double lastDtMin = 0.0;
int ncol = 0;
int Length;
int eof;
Expand Down