From 0479810a8442c1759d1655f6e5fb986429d0d4af Mon Sep 17 00:00:00 2001 From: David Huard Date: Thu, 2 Apr 2026 14:21:06 -0400 Subject: [PATCH 01/13] Added deflate instructions and global NETCDF_DEFLATE_LEVEL option --- src/CustomOutput.cpp | 8 +++++++- src/RavenInclude.h | 1 + 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/src/CustomOutput.cpp b/src/CustomOutput.cpp index a2c7d60..921c8a7 100644 --- a/src/CustomOutput.cpp +++ b/src/CustomOutput.cpp @@ -579,6 +579,8 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) // (b) Define the time variable. dimids1[0] = time_dimid; retval = nc_def_var(_netcdf_ID, "time", NC_DOUBLE, ndims1,dimids1, &varid_time); HandleNetCDFErrors(retval); + // Enable deflate compression for time variable (shuffle, zlib, deflate_level) + retval = nc_def_var_deflate(_netcdf_ID, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); // (c) Assign units attributes to the netCDF VARIABLES. // --> converts start day into "hours since YYYY-MM-DD HH:MM:SS" @@ -616,6 +618,8 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) dimids1[0] = ndata_dimid; retval = nc_def_var(_netcdf_ID, group_name.c_str(), NC_STRING, ndims1, dimids1, &varid_grps); HandleNetCDFErrors(retval); + // Enable deflate compression for group variable + retval = nc_def_var_deflate(_netcdf_ID, varid_grps, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); //(c) set some attributes to variable "HRUID" or "SBID" tmp =long_name; @@ -632,6 +636,8 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) dimids2[0] = time_dimid; dimids2[1] = ndata_dimid; retval = nc_def_var(_netcdf_ID, netCDFtag.c_str(), NC_DOUBLE, ndims2, dimids2, &varid_data); HandleNetCDFErrors(retval); + // Enable deflate compression for data variable + retval = nc_def_var_deflate(_netcdf_ID, varid_data, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); //(f) set some attributes to variable _netCDFtag tmp=_timeAggStr+" "+_statStr+" "+_varName+" "+_spaceAggStr; @@ -1267,4 +1273,4 @@ CCustomOutput *CCustomOutput::ParseCustomOutputCommand(char *s[MAXINPUTITEMS], c } return pCustom; -} +} diff --git a/src/RavenInclude.h b/src/RavenInclude.h index 22bd49e..f099892 100644 --- a/src/RavenInclude.h +++ b/src/RavenInclude.h @@ -251,6 +251,7 @@ const double USE_TEMPLATE_VALUE =-55555.5; const double NOT_NEEDED =-66666.6; ///< arbitrary value indicating that a non-auto parameter is not needed for the current model configuration const double NOT_NEEDED_AUTO =-77777.7; ///< arbitrary value indicating that a autogeneratable parameter is not needed for the current model configuration const double NETCDF_BLANK_VALUE =-9999.0; ///< NetCDF flag for blank value +const int NETCDF_DEFLATE_LEVEL = 6; ///< NetCDF deflate level (compression 1-9) const double RAV_BLANK_DATA =-1.2345; ///< double corresponding to blank/void data item (also used in input files) const double DIRICHLET_TEMP =-9999.0; ///< dirichlet concentration flag corresponding to air temperature const int FROM_STATION_VAR =-55; ///< special flag indicating that NetCDF indices should be looked up from station attribute table From 3fd23eb1aa509cac43e3e871349d31015ade6be5 Mon Sep 17 00:00:00 2001 From: David Huard Date: Thu, 2 Apr 2026 16:02:19 -0400 Subject: [PATCH 02/13] set compression for standard outputs --- src/CustomOutput.cpp | 3 +++ src/StandardOutput.cpp | 26 ++++++++++++++++++++++++-- 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/src/CustomOutput.cpp b/src/CustomOutput.cpp index 921c8a7..ab849e0 100644 --- a/src/CustomOutput.cpp +++ b/src/CustomOutput.cpp @@ -581,6 +581,9 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) retval = nc_def_var(_netcdf_ID, "time", NC_DOUBLE, ndims1,dimids1, &varid_time); HandleNetCDFErrors(retval); // Enable deflate compression for time variable (shuffle, zlib, deflate_level) retval = nc_def_var_deflate(_netcdf_ID, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); + // Set chunksize for time variable + + // (c) Assign units attributes to the netCDF VARIABLES. // --> converts start day into "hours since YYYY-MM-DD HH:MM:SS" diff --git a/src/StandardOutput.cpp b/src/StandardOutput.cpp index 6660c84..ab54143 100644 --- a/src/StandardOutput.cpp +++ b/src/StandardOutput.cpp @@ -2008,6 +2008,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) int dimids1[ndims1]; // array which will contain all dimension ids for a variable int ncid, varid_pre; // When we create netCDF variables and dimensions, we get back an ID for each one. int time_dimid, varid_time; // dimension ID (holds number of time steps) and variable ID (holds time values) for time + int ntime; // Number of time steps int nSim, nbasins_dimid, varid_bsim,varid_bsim2; // // # of sub-basins with simulated outflow, dimension ID, and // // variable to write basin IDs for simulated outflows @@ -2021,6 +2022,9 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) string tmpFilename; int p; // loop over all sub-basins string tmp,tmp2,tmp3; + size_t time_chunksize[1]; + + // initialize all potential file IDs with -9 == "not existing and hence not opened" _HYDRO_ncid = -9; // output file ID for Hydrographs.nc (-9 --> not opened) @@ -2055,7 +2059,8 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // time // ---------------------------------------------------------- // (a) Define the DIMENSIONS. NetCDF will hand back an ID - retval = nc_def_dim(_HYDRO_ncid, "time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); + ntime = ceil(Options.duration / Options.timestep); + retval = nc_def_dim(_HYDRO_ncid, "time", ntime, &time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. Assign units attributes to the netCDF VARIABLES. dimids1[0] = time_dimid; @@ -2064,6 +2069,13 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) retval = nc_put_att_text(_HYDRO_ncid, varid_time, "calendar", strlen("gregorian"), "gregorian"); HandleNetCDFErrors(retval); retval = nc_put_att_text(_HYDRO_ncid, varid_time, "standard_name", strlen("time"), "time"); HandleNetCDFErrors(retval); + // Enable deflate compression for time variable (shuffle, zlib, deflate_level) + retval = nc_def_var_deflate(_HYDRO_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); + + // Set chunksize to the number of time steps + size_t chunksize_time = ntime; + retval = nc_def_var_chunking(_HYDRO_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval); + // define precipitation variable varid_pre= NetCDFAddMetadata(_HYDRO_ncid, time_dimid,"precip","Precipitation","mm d**-1"); @@ -2081,7 +2093,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // (b) create dimension "nbasins" retval = nc_def_dim(_HYDRO_ncid, "nbasins", nSim, &nbasins_dimid); HandleNetCDFErrors(retval); - // (c) create variable and set attributes for"basin_name" + // (c) create variable and set attributes for "basin_name" dimids1[0] = nbasins_dimid; retval = nc_def_var(_HYDRO_ncid, "basin_name", NC_STRING, ndims1, dimids1, &varid_bsim); HandleNetCDFErrors(retval); tmp ="ID of sub-basins with simulated outflows"; @@ -3049,6 +3061,7 @@ int NetCDFAddMetadata2D(const int fileid,const int time_dimid,int nbasins_dimid, int retval; int dimids2[2]; string tmp; + size_t chunksize2[2]; static double fill_val[] = {NETCDF_BLANK_VALUE}; static double miss_val[] = {NETCDF_BLANK_VALUE}; @@ -3058,6 +3071,15 @@ int NetCDFAddMetadata2D(const int fileid,const int time_dimid,int nbasins_dimid, // (a) create variable retval = nc_def_var(fileid,shortname.c_str(),NC_DOUBLE,2,dimids2,&varid); HandleNetCDFErrors(retval); + // Set compression + retval = nc_def_var_deflate(fileid, varid, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); + // Set time chunksize to number of time steps + retval = nc_inq_dimlen(fileid, time_dimid, &chunksize2[0]); HandleNetCDFErrors(retval); + // Set nbasins chunksize to number ensuring that chunks have approximately 10 MB of data (assuming double precision) + retval = nc_inq_dimlen(fileid, nbasins_dimid, &chunksize2[1]); HandleNetCDFErrors(retval); + chunksize2[1] = max((size_t)1, min(chunksize2[1], (size_t)(10 * 1024 * 1024 / sizeof(double) / chunksize2[0]))); // Ensure at least one basin per chunk + // Set chunksize + retval = nc_def_var_chunking(fileid, varid, NC_CHUNKED, chunksize2); HandleNetCDFErrors(retval); tmp = "basin_name"; From 7f18bde99763367076aa0680e6acc6c31cb086d6 Mon Sep 17 00:00:00 2001 From: David Huard Date: Thu, 2 Apr 2026 17:58:05 -0400 Subject: [PATCH 03/13] parameterized chunksize_mb --- src/CustomOutput.cpp | 14 +++++++++++++- src/RavenInclude.h | 1 + src/StandardOutput.cpp | 4 ++-- 3 files changed, 16 insertions(+), 3 deletions(-) diff --git a/src/CustomOutput.cpp b/src/CustomOutput.cpp index ab849e0..aa5bda2 100644 --- a/src/CustomOutput.cpp +++ b/src/CustomOutput.cpp @@ -554,6 +554,7 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) int retval; // error value for NetCDF routines size_t start[1], count[1]; // determines where and how much will be written to NetCDF + size_t chunksize2[2]; string tmp,tmp2,tmp3,tmp4; bool cant_support=(_aggstat==AGG_RANGE || _aggstat==AGG_95CI || _aggstat==AGG_QUARTILES || _aggstat==AGG_HISTOGRAM); @@ -574,15 +575,19 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) // time // ---------------------------------------------------------- // (a) Define the DIMENSIONS. NetCDF will hand back an ID for each. + + // TODO: Set dimension size instead of unlimited. retval = nc_def_dim(_netcdf_ID, "time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); // (b) Define the time variable. dimids1[0] = time_dimid; retval = nc_def_var(_netcdf_ID, "time", NC_DOUBLE, ndims1,dimids1, &varid_time); HandleNetCDFErrors(retval); + // Enable deflate compression for time variable (shuffle, zlib, deflate_level) retval = nc_def_var_deflate(_netcdf_ID, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); - // Set chunksize for time variable + // TODO: Set chunksize to len(time) + // retval = nc_def_var_chunking(_netcdf_ID, varid_time, NC_CHUNKED, -1); HandleNetCDFErrors(retval); // (c) Assign units attributes to the netCDF VARIABLES. @@ -621,6 +626,7 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) dimids1[0] = ndata_dimid; retval = nc_def_var(_netcdf_ID, group_name.c_str(), NC_STRING, ndims1, dimids1, &varid_grps); HandleNetCDFErrors(retval); + // Enable deflate compression for group variable retval = nc_def_var_deflate(_netcdf_ID, varid_grps, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); @@ -639,9 +645,15 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) dimids2[0] = time_dimid; dimids2[1] = ndata_dimid; retval = nc_def_var(_netcdf_ID, netCDFtag.c_str(), NC_DOUBLE, ndims2, dimids2, &varid_data); HandleNetCDFErrors(retval); + // Enable deflate compression for data variable retval = nc_def_var_deflate(_netcdf_ID, varid_data, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); + // TODO: Set chunksizes for data variable (time, ndata) + //chunksize2[0]=_nDataItems; //chunk along time dimension + //chunksize2[1] = max((size_t)1, min(_nData, (size_t)(NETCDF_CHUNKSIZE_MB * 1024 * 1024 / sizeof(double) / chunksize2[0]))); // Ensure at least one basin per chunk + //retval = nc_def_var_chunking(_netcdf_ID, varid_data, NC_CHUNKED, chunksize2); HandleNetCDFErrors(retval); + //(f) set some attributes to variable _netCDFtag tmp=_timeAggStr+" "+_statStr+" "+_varName+" "+_spaceAggStr; tmp2=_varUnits; diff --git a/src/RavenInclude.h b/src/RavenInclude.h index f099892..21867bb 100644 --- a/src/RavenInclude.h +++ b/src/RavenInclude.h @@ -252,6 +252,7 @@ const double NOT_NEEDED =-66666.6; const double NOT_NEEDED_AUTO =-77777.7; ///< arbitrary value indicating that a autogeneratable parameter is not needed for the current model configuration const double NETCDF_BLANK_VALUE =-9999.0; ///< NetCDF flag for blank value const int NETCDF_DEFLATE_LEVEL = 6; ///< NetCDF deflate level (compression 1-9) +const int NETCDF_CHUNKSIZE_MB = 10; ///< NetCDF chunk size in MB (must be >1 for compression to work) const double RAV_BLANK_DATA =-1.2345; ///< double corresponding to blank/void data item (also used in input files) const double DIRICHLET_TEMP =-9999.0; ///< dirichlet concentration flag corresponding to air temperature const int FROM_STATION_VAR =-55; ///< special flag indicating that NetCDF indices should be looked up from station attribute table diff --git a/src/StandardOutput.cpp b/src/StandardOutput.cpp index ab54143..90cb05a 100644 --- a/src/StandardOutput.cpp +++ b/src/StandardOutput.cpp @@ -2059,7 +2059,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // time // ---------------------------------------------------------- // (a) Define the DIMENSIONS. NetCDF will hand back an ID - ntime = ceil(Options.duration / Options.timestep); + ntime = (int)(ceil((Options.duration+TIME_CORRECTION)/Options.timestep)); retval = nc_def_dim(_HYDRO_ncid, "time", ntime, &time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. Assign units attributes to the netCDF VARIABLES. @@ -3077,7 +3077,7 @@ int NetCDFAddMetadata2D(const int fileid,const int time_dimid,int nbasins_dimid, retval = nc_inq_dimlen(fileid, time_dimid, &chunksize2[0]); HandleNetCDFErrors(retval); // Set nbasins chunksize to number ensuring that chunks have approximately 10 MB of data (assuming double precision) retval = nc_inq_dimlen(fileid, nbasins_dimid, &chunksize2[1]); HandleNetCDFErrors(retval); - chunksize2[1] = max((size_t)1, min(chunksize2[1], (size_t)(10 * 1024 * 1024 / sizeof(double) / chunksize2[0]))); // Ensure at least one basin per chunk + chunksize2[1] = max((size_t)1, min(chunksize2[1], (size_t)(NETCDF_CHUNKSIZE_MB * 1024 * 1024 / sizeof(double) / chunksize2[0]))); // Ensure at least one basin per chunk // Set chunksize retval = nc_def_var_chunking(fileid, varid, NC_CHUNKED, chunksize2); HandleNetCDFErrors(retval); From 1abceabc4ff2b8ea29e595426697c654137125d7 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 7 Apr 2026 14:12:26 +0000 Subject: [PATCH 04/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/Assimilate.cpp | 2 +- src/ChannelXSect.cpp | 2 +- src/ParsePropertyFile.cpp | 18 +++++++++--------- src/Reservoir.cpp | 8 ++++---- src/Reservoir.h | 2 +- 5 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/Assimilate.cpp b/src/Assimilate.cpp index c080e6a..73fedc8 100644 --- a/src/Assimilate.cpp +++ b/src/Assimilate.cpp @@ -157,7 +157,7 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt) { Qobs = _pObservedTS[i]->GetSampledValue(nn); //mean timestep flow Qobs2 = _pObservedTS[i]->GetSampledValue(nn+1); //mean timestep flow - + //override initial conditions directly if ((nn==1) && (Qobs!=RAV_BLANK_DATA)){ _pSubBasins[p]->SetQout(Qobs); diff --git a/src/ChannelXSect.cpp b/src/ChannelXSect.cpp index 0ad00a6..d950a88 100644 --- a/src/ChannelXSect.cpp +++ b/src/ChannelXSect.cpp @@ -389,7 +389,7 @@ double CChannelXSect::GetDiffusivity(const double &Q, const double &SB_slope, co { ExitGracefullyIf(Q<=0,"CChannelXSect::GetDiffusivity: Invalid channel flowrate",BAD_DATA); - ///< diffusivity from Roberson et al. 1995, Hydraulic Engineering + ///< diffusivity from Roberson et al. 1995, Hydraulic Engineering double slope_mult=1.0; double Q_mult =1.0; GetFlowCorrections(SB_slope,SB_n,slope_mult,Q_mult); diff --git a/src/ParsePropertyFile.cpp b/src/ParsePropertyFile.cpp index ddaf48c..3911732 100644 --- a/src/ParsePropertyFile.cpp +++ b/src/ParsePropertyFile.cpp @@ -38,10 +38,10 @@ void AddNewSoilClass (CSoilClass **&pSoilClasses, soil_struct **&parsed_s string *&soiltags, int &num_parsed_soils, int nConstits, const string name, bool isdefault, CModel *pModel); void AddNewLULTClass (CLandUseClass **&pLandUseClasses, surface_struct **&parsed_lults, - string *&lulttags, int &num_parsed_lults, + string *&lulttags, int &num_parsed_lults, const string name, bool isdefault, CModel *pModel); void AddNewVegClass (CVegetationClass **&pVegClasses, veg_struct **&parsed_vegs, - string *&vegtags, int &num_parsed_vegs, + string *&vegtags, int &num_parsed_vegs, const string name, bool isdefault, CModel *pModel); ////////////////////////////////////////////////////////////////// /// \brief This method parses the class properties .rvp file @@ -86,7 +86,7 @@ bool ParseClassPropertiesFile(CModel *&pModel, CSoilProfile **pProfiles=NULL; int num_parsed_aqstacks=0; - + bool invalid_index; int num_read; int *indices=NULL; @@ -518,7 +518,7 @@ bool ParseClassPropertiesFile(CModel *&pModel, invalid_index=ParsePropArray(p,indices,properties,num_read,soiltags,nParamStrings,num_parsed_soils,aAliases,nAliases); ExitGracefullyIf(invalid_index, "ParseClassPropertiesFile: Invalid soiltype code in SoilParameterList command",BAD_DATA); - + if (Options.noisy){ for (int j=0;jAutoCalculateSoilProps (*parsed_soils[c],*parsed_soils[0],pModel->GetTransportModel()->GetNumConstituents()); } for (int c=1;cAutoCalculateLandUseProps (*parsed_surf[c], *parsed_surf[0]); + pLUClasses [c]->AutoCalculateLandUseProps (*parsed_surf[c], *parsed_surf[0]); } for (int c=1;cAutoCalculateTerrainProps ( parsed_terrs[c], parsed_terrs[0]); @@ -1577,7 +1577,7 @@ bool ParseClassPropertiesFile(CModel *&pModel, pModel->CheckForChannelXSectsDuplicates(Options); - + delete p; return true; @@ -1648,7 +1648,7 @@ bool ParsePropArray(CParser *p, //parser tmpproperties[num_read][j-1]=AutoOrDoubleOrAlias(s[j],pAliases,nAliases); } DeletePropArray(indices,properties,num_read); - + properties=tmpproperties; indices=tmpind; num_read++; @@ -1808,7 +1808,7 @@ void AddNewSoilClass(CSoilClass **&pSoilClasses, /// \brief dynamically adds a new LULT class to the array of parsed LULT classes. // void AddNewLULTClass (CLandUseClass **&pLandUseClasses, surface_struct **&parsed_lults, - string *&lulttags, int &num_parsed_lults, + string *&lulttags, int &num_parsed_lults, const string name, bool isdefault, CModel *pModel){ CLandUseClass *pLC; pLC = new CLandUseClass(name, pModel); @@ -1838,7 +1838,7 @@ void AddNewLULTClass (CLandUseClass **&pLandUseClasses, surface_struct ** /// \brief dynamically adds a new veg class to the array of parsed vegetation classes. // void AddNewVegClass (CVegetationClass **&pVegClasses, veg_struct **&parsed_vegs, - string *&vegtags, int &num_parsed_vegs, + string *&vegtags, int &num_parsed_vegs, const string name, bool isdefault, CModel *pModel){ CVegetationClass *pVC; pVC = new CVegetationClass(name, pModel); diff --git a/src/Reservoir.cpp b/src/Reservoir.cpp index 186b0f4..aa36466 100644 --- a/src/Reservoir.cpp +++ b/src/Reservoir.cpp @@ -511,17 +511,17 @@ double CReservoir::GetLakeConvectionCoeff() const { return _lake_convcoeff; } ////////////////////////////////////////////////////////////////// /// \returns current outflow rate [m3/s] // -double CReservoir::GetOutflowRate (const bool adjusted) const { +double CReservoir::GetOutflowRate (const bool adjusted) const { if (adjusted){return max(_Qout+_DAadjust,0.0);} - else {return _Qout;} + else {return _Qout;} } ////////////////////////////////////////////////////////////////// /// \returns previous outflow rate [m3/s] // -double CReservoir::GetOldOutflowRate (const bool adjusted) const { +double CReservoir::GetOldOutflowRate (const bool adjusted) const { if (adjusted){return max(_Qout_last+_DAadjust_last,0.0);} - else {return _Qout_last;} + else {return _Qout_last;} } ////////////////////////////////////////////////////////////////// diff --git a/src/Reservoir.h b/src/Reservoir.h index 4578b0b..26bc365 100644 --- a/src/Reservoir.h +++ b/src/Reservoir.h @@ -99,7 +99,7 @@ class CReservoir double _DAadjust; //< outflow adjustment - used for reporting overriden flows [m3/s] double _DAadjust_last; //< outflow adjustment [m3/s] from previous time step - + int _dry_timesteps; //< number of time steps this reservoir dried out during simulation //state variables : From e6dd68205a533d097427060c68c02a6a7eced7ec Mon Sep 17 00:00:00 2001 From: David Huard Date: Tue, 7 Apr 2026 17:58:57 -0400 Subject: [PATCH 05/13] add function to count the number of custom output time steps --- src/CustomOutput.cpp | 54 +++++++++++++++++++++++++++++++++++++++++++- src/CustomOutput.h | 3 +++ src/Makefile | 4 ++-- 3 files changed, 58 insertions(+), 3 deletions(-) diff --git a/src/CustomOutput.cpp b/src/CustomOutput.cpp index aa5bda2..c211ca0 100644 --- a/src/CustomOutput.cpp +++ b/src/CustomOutput.cpp @@ -217,7 +217,57 @@ void CCustomOutput::SetHistogramParams(const double minv,const double maxv, cons _hist_max=maxv; _nBins=numBins; } - +////////////////////////////////////////////////////////////////// +/// \brief Compute the approximate number of time steps for the output array +/// \param Options [in] Global model options information +/// \return Number of time steps (int) +int CCustomOutput::ComputeNumTimeSteps(const optStruct &Options) const +{ + // Example assumes Options has julian_start_day, timestep, and custom_interval + + double julian_end_day; + int julian_end_year; + time_struct start, end; + + JulianConvert(0, Options.julian_start_day, Options.julian_start_year, Options.calendar, start); + JulianConvert(Options.duration, Options.julian_start_day, Options.julian_start_year, Options.calendar, end); + + // Compute the end date in Julian day and year format + // AddTime(Options.julian_start_day, Options.julian_start_year, Options.duration, Options.calendar, julian_end_day, julian_end_year); + // ntime = (int)(ceil((Options.duration+TIME_CORRECTION)/Options.timestep)); + double dt = Options.timestep; + int nsteps = 0; + switch(_timeAgg) { + case YEARLY: + nsteps = int(end.year - start.year); + break; + case WATER_YEARLY: + // wateryr_mo: Starting month of the water year + nsteps = int(end.year - start.year); + if ((start.month < end.month) & (start.month < Options.wateryr_mo) & (end.month >= Options.wateryr_mo)) nsteps++; + if ((start.month > end.month) & (start.month >= Options.wateryr_mo) & (end.month < Options.wateryr_mo)) nsteps--; + break; + case MONTHLY: + nsteps = int((end.year - start.year) * 12 + (end.month - start.month)); + break; + case DAILY: + nsteps = int(ceil(Options.duration)); + break; + case EVERY_NDAYS: + nsteps = int(ceil(Options.duration / Options.custom_interval)); + break; + case EVERY_TSTEP: + nsteps = int(ceil(Options.duration / dt)); + break; + case ENTIRE_SIM: + nsteps = 1; + break; + default: + nsteps = 0; + break; + } + return max(0, nsteps); +} /////////////////////////////////////////////////////////////////// /// \brief Allocates memory and initialize data storage of a CCustomOutput object /// \remarks Called prior to simulation. Determines size of and allocates memory for (member) data[][] array needed in statistical calculations @@ -577,6 +627,8 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) // (a) Define the DIMENSIONS. NetCDF will hand back an ID for each. // TODO: Set dimension size instead of unlimited. + cout << "Would create time dimension with size " << ComputeNumTimeSteps(Options) << endl; + retval = nc_def_dim(_netcdf_ID, "time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); // (b) Define the time variable. diff --git a/src/CustomOutput.h b/src/CustomOutput.h index d9fb177..5fea67c 100644 --- a/src/CustomOutput.h +++ b/src/CustomOutput.h @@ -126,7 +126,10 @@ class CCustomOutput void SetHistogramParams(const double min,const double max, const int numBins); + int ComputeNumTimeSteps(const optStruct &Options) const; + spatial_agg GetSpatialAgg() const; + void InitializeCustomOutput(const optStruct &Options); void WriteFileHeader (const optStruct &Options); diff --git a/src/Makefile b/src/Makefile index d897009..f135074 100644 --- a/src/Makefile +++ b/src/Makefile @@ -26,8 +26,8 @@ LDFLAGS := CXXFLAGS += -std=c++11 -fPIC # OPTION 1) include netcdf - uncomment following two commands (assumes netCDF path = /usr/local): -#CXXFLAGS += -Dnetcdf -#LDLIBS += -L/usr/local -lnetcdf +CXXFLAGS += -Dnetcdf +LDLIBS += -L/usr/local -lnetcdf # OPTION 1b) include netcdf - for newer MacOS with Apple Silicon (use with option 1 also uncommented) #CXXFLAGS += -I/opt/homebrew/include From 10ad2076daaa0c6c4488531a565c79db96df90c8 Mon Sep 17 00:00:00 2001 From: David Huard Date: Tue, 7 Apr 2026 18:12:04 -0400 Subject: [PATCH 06/13] Add chunking to CustomOutput --- src/CustomOutput.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/CustomOutput.cpp b/src/CustomOutput.cpp index c211ca0..1a1f91f 100644 --- a/src/CustomOutput.cpp +++ b/src/CustomOutput.cpp @@ -604,6 +604,7 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) int retval; // error value for NetCDF routines size_t start[1], count[1]; // determines where and how much will be written to NetCDF + size_t ntime; // Number of time steps size_t chunksize2[2]; string tmp,tmp2,tmp3,tmp4; @@ -627,9 +628,9 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) // (a) Define the DIMENSIONS. NetCDF will hand back an ID for each. // TODO: Set dimension size instead of unlimited. - cout << "Would create time dimension with size " << ComputeNumTimeSteps(Options) << endl; + ntime = ComputeNumTimeSteps(Options); - retval = nc_def_dim(_netcdf_ID, "time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_netcdf_ID, "time", ntime, &time_dimid); HandleNetCDFErrors(retval); // (b) Define the time variable. dimids1[0] = time_dimid; @@ -638,8 +639,8 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) // Enable deflate compression for time variable (shuffle, zlib, deflate_level) retval = nc_def_var_deflate(_netcdf_ID, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); - // TODO: Set chunksize to len(time) - // retval = nc_def_var_chunking(_netcdf_ID, varid_time, NC_CHUNKED, -1); HandleNetCDFErrors(retval); + // Set chunksize to len(time) + retval = nc_def_var_chunking(_netcdf_ID, varid_time, NC_CHUNKED, &ntime); HandleNetCDFErrors(retval); // (c) Assign units attributes to the netCDF VARIABLES. @@ -701,10 +702,10 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) // Enable deflate compression for data variable retval = nc_def_var_deflate(_netcdf_ID, varid_data, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); - // TODO: Set chunksizes for data variable (time, ndata) - //chunksize2[0]=_nDataItems; //chunk along time dimension - //chunksize2[1] = max((size_t)1, min(_nData, (size_t)(NETCDF_CHUNKSIZE_MB * 1024 * 1024 / sizeof(double) / chunksize2[0]))); // Ensure at least one basin per chunk - //retval = nc_def_var_chunking(_netcdf_ID, varid_data, NC_CHUNKED, chunksize2); HandleNetCDFErrors(retval); + // Set chunksizes for data variable (time, ndata) + chunksize2[0] = ntime; //chunk along time dimension + chunksize2[1] = max((size_t)1, min(_nData, (size_t)(NETCDF_CHUNKSIZE_MB * 1024 * 1024 / sizeof(double) / chunksize2[0]))); // Ensure at least one basin per chunk + retval = nc_def_var_chunking(_netcdf_ID, varid_data, NC_CHUNKED, chunksize2); HandleNetCDFErrors(retval); //(f) set some attributes to variable _netCDFtag tmp=_timeAggStr+" "+_statStr+" "+_varName+" "+_spaceAggStr; From f2b63166bb28d384f06333ecd63b2cd54bd138b1 Mon Sep 17 00:00:00 2001 From: David Huard Date: Wed, 8 Apr 2026 09:05:50 -0400 Subject: [PATCH 07/13] try to extend compression logic to other output variables (constituent, storage, etc) --- src/ConstituentModel.cpp | 8 +++++--- src/CustomOutput.cpp | 21 +++++---------------- src/ModelInitialize.cpp | 1 + src/ParseInput.cpp | 6 ++++++ src/RavenInclude.h | 1 + src/StandardOutput.cpp | 22 ++++++++++++---------- 6 files changed, 30 insertions(+), 29 deletions(-) diff --git a/src/ConstituentModel.cpp b/src/ConstituentModel.cpp index 7bac925..b1ae654 100644 --- a/src/ConstituentModel.cpp +++ b/src/ConstituentModel.cpp @@ -1004,6 +1004,7 @@ void CConstituentModel::WriteNetCDFOutputFileHeaders(const optStruct& Options) const int ndims2 = 2; int dimids1[ndims1]; // array which will contain all dimension ids for a variable int retval,ncid; + size_t ntime; // Size of the time dimension //converts start day into "hours since YYYY-MM-DD HH:MM:SS" (model start time) char starttime[200]; // start time string in format 'hours since YYY-MM-DD HH:MM:SS' @@ -1042,7 +1043,8 @@ void CConstituentModel::WriteNetCDFOutputFileHeaders(const optStruct& Options) // time vector // ---------------------------------------------------------- // Define the DIMENSIONS. NetCDF will hand back an ID - retval = nc_def_dim(_CONC_ncid,"time",NC_UNLIMITED,&time_dimid); HandleNetCDFErrors(retval); + size_t ntime = min(1, Options.n_out_time); + retval = nc_def_dim(_CONC_ncid,"time",ntime,&time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. dimids1[0] = time_dimid; @@ -1094,7 +1096,7 @@ void CConstituentModel::WriteNetCDFOutputFileHeaders(const optStruct& Options) // time vector // ---------------------------------------------------------- // Define the DIMENSIONS. NetCDF will hand back an ID - retval = nc_def_dim(_POLLUT_ncid,"time",NC_UNLIMITED,&time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_POLLUT_ncid,"time",ntime,&time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. dimids1[0] = time_dimid; @@ -1182,7 +1184,7 @@ void CConstituentModel::WriteNetCDFOutputFileHeaders(const optStruct& Options) // time vector // ---------------------------------------------------------- // Define the DIMENSIONS. NetCDF will hand back an ID - retval = nc_def_dim(_LOADING_ncid,"time",NC_UNLIMITED,&time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_LOADING_ncid,"time",ntime,&time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. dimids1[0] = time_dimid; diff --git a/src/CustomOutput.cpp b/src/CustomOutput.cpp index 1a1f91f..a2de976 100644 --- a/src/CustomOutput.cpp +++ b/src/CustomOutput.cpp @@ -223,20 +223,12 @@ void CCustomOutput::SetHistogramParams(const double minv,const double maxv, cons /// \return Number of time steps (int) int CCustomOutput::ComputeNumTimeSteps(const optStruct &Options) const { - // Example assumes Options has julian_start_day, timestep, and custom_interval - - double julian_end_day; - int julian_end_year; + int nsteps = 0; // return value time_struct start, end; JulianConvert(0, Options.julian_start_day, Options.julian_start_year, Options.calendar, start); JulianConvert(Options.duration, Options.julian_start_day, Options.julian_start_year, Options.calendar, end); - // Compute the end date in Julian day and year format - // AddTime(Options.julian_start_day, Options.julian_start_year, Options.duration, Options.calendar, julian_end_day, julian_end_year); - // ntime = (int)(ceil((Options.duration+TIME_CORRECTION)/Options.timestep)); - double dt = Options.timestep; - int nsteps = 0; switch(_timeAgg) { case YEARLY: nsteps = int(end.year - start.year); @@ -244,8 +236,8 @@ int CCustomOutput::ComputeNumTimeSteps(const optStruct &Options) const case WATER_YEARLY: // wateryr_mo: Starting month of the water year nsteps = int(end.year - start.year); - if ((start.month < end.month) & (start.month < Options.wateryr_mo) & (end.month >= Options.wateryr_mo)) nsteps++; - if ((start.month > end.month) & (start.month >= Options.wateryr_mo) & (end.month < Options.wateryr_mo)) nsteps--; + if ((start.month < end.month) && (start.month < Options.wateryr_mo) && (end.month >= Options.wateryr_mo)) nsteps++; + if ((start.month > end.month) && (start.month >= Options.wateryr_mo) && (end.month < Options.wateryr_mo)) nsteps--; break; case MONTHLY: nsteps = int((end.year - start.year) * 12 + (end.month - start.month)); @@ -257,7 +249,7 @@ int CCustomOutput::ComputeNumTimeSteps(const optStruct &Options) const nsteps = int(ceil(Options.duration / Options.custom_interval)); break; case EVERY_TSTEP: - nsteps = int(ceil(Options.duration / dt)); + nsteps = int(ceil(Options.duration / Options.timestep)); break; case ENTIRE_SIM: nsteps = 1; @@ -628,7 +620,7 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) // (a) Define the DIMENSIONS. NetCDF will hand back an ID for each. // TODO: Set dimension size instead of unlimited. - ntime = ComputeNumTimeSteps(Options); + ntime = max(1, ComputeNumTimeSteps(Options)); retval = nc_def_dim(_netcdf_ID, "time", ntime, &time_dimid); HandleNetCDFErrors(retval); @@ -680,9 +672,6 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) dimids1[0] = ndata_dimid; retval = nc_def_var(_netcdf_ID, group_name.c_str(), NC_STRING, ndims1, dimids1, &varid_grps); HandleNetCDFErrors(retval); - // Enable deflate compression for group variable - retval = nc_def_var_deflate(_netcdf_ID, varid_grps, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); - //(c) set some attributes to variable "HRUID" or "SBID" tmp =long_name; tmp2="timeseries_id"; diff --git a/src/ModelInitialize.cpp b/src/ModelInitialize.cpp index 6a69e75..93fff82 100644 --- a/src/ModelInitialize.cpp +++ b/src/ModelInitialize.cpp @@ -70,6 +70,7 @@ void CModel::Initialize(const optStruct &Options) } } } + // reserve memory for mass balance arrays //-------------------------------------------------------------- _aCumulativeBal = new double * [_nHydroUnits]; diff --git a/src/ParseInput.cpp b/src/ParseInput.cpp index 0983819..71ff525 100644 --- a/src/ParseInput.cpp +++ b/src/ParseInput.cpp @@ -3729,6 +3729,12 @@ bool ParseMainInputFile (CModel *&pModel, ExitGracefullyIf(Options.duration<0, "ParseMainInputFile::Model duration less than zero. Make sure :EndDate is after :StartDate.",BAD_DATA_WARN); + // Compute the size of the output's time dimension + Options.n_out_time = (int)(floor(ceil(Options.duration/Options.timestep))/Options.output_interval); + if (Options.n_out_time==0) { + WriteAdvisory("ParseMainInputFile::Number of output time steps is zero. Check :Duration, :Timestep, and :OutputInterval commands.",BAD_DATA); + } + if((Options.nNetCDFattribs>0) && (Options.output_format!=OUTPUT_NETCDF)){ WriteAdvisory("ParseMainInputFile: NetCDF attributes were specified but output format is not NetCDF.",Options.noisy); } diff --git a/src/RavenInclude.h b/src/RavenInclude.h index 6c90350..56780ba 100644 --- a/src/RavenInclude.h +++ b/src/RavenInclude.h @@ -1065,6 +1065,7 @@ struct optStruct double convergence_crit; ///< convergence criteria double max_iterations; ///< maximum number of iterations for iterative solver method double timestep; ///< numerical method timestep (in days) + int n_out_time; ///< size of output time dimension double output_interval; ///< write to output file every x number of timesteps ensemble_type ensemble; ///< ensemble type (or ENSEMBLE_NONE if single model) string external_script; ///< call to external script/.exe once per timestep (or "" if none) diff --git a/src/StandardOutput.cpp b/src/StandardOutput.cpp index cf105db..d8a0b1a 100644 --- a/src/StandardOutput.cpp +++ b/src/StandardOutput.cpp @@ -2007,7 +2007,6 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) int dimids1[ndims1]; // array which will contain all dimension ids for a variable int ncid, varid_pre; // When we create netCDF variables and dimensions, we get back an ID for each one. int time_dimid, varid_time; // dimension ID (holds number of time steps) and variable ID (holds time values) for time - int ntime; // Number of time steps int nSim, nbasins_dimid, varid_bsim,varid_bsim2; // // # of sub-basins with simulated outflow, dimension ID, and // // variable to write basin IDs for simulated outflows @@ -2021,8 +2020,6 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) string tmpFilename; int p; // loop over all sub-basins string tmp,tmp2,tmp3; - size_t time_chunksize[1]; - // initialize all potential file IDs with -9 == "not existing and hence not opened" @@ -2058,8 +2055,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // time // ---------------------------------------------------------- // (a) Define the DIMENSIONS. NetCDF will hand back an ID - ntime = (int)(ceil((Options.duration+TIME_CORRECTION)/Options.timestep)); - retval = nc_def_dim(_HYDRO_ncid, "time", ntime, &time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_HYDRO_ncid, "time", Options.n_out_time, &time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. Assign units attributes to the netCDF VARIABLES. dimids1[0] = time_dimid; @@ -2072,7 +2068,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) retval = nc_def_var_deflate(_HYDRO_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); // Set chunksize to the number of time steps - size_t chunksize_time = ntime; + size_t chunksize_time = min(1, Options.n_out_time); retval = nc_def_var_chunking(_HYDRO_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval); // define precipitation variable @@ -2152,7 +2148,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // time // ---------------------------------------------------------- // (a) Define the DIMENSIONS. NetCDF will hand back an ID - retval = nc_def_dim(_RESSTAGE_ncid,"time",NC_UNLIMITED,&time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_RESSTAGE_ncid,"time",ntime,&time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. Assign units attributes to the netCDF VARIABLES. dimids1[0] = time_dimid; @@ -2234,7 +2230,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // time vector // ---------------------------------------------------------- // Define the DIMENSIONS. NetCDF will hand back an ID - retval = nc_def_dim(_STORAGE_ncid, "time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_STORAGE_ncid, "time", ntime, &time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. dimids1[0] = time_dimid; @@ -2284,7 +2280,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // ---------------------------------------------------------- // time vector // ---------------------------------------------------------- - retval = nc_def_dim(_FORCINGS_ncid,"time",NC_UNLIMITED,&time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_FORCINGS_ncid,"time",ntime,&time_dimid); HandleNetCDFErrors(retval); dimids1[0] = time_dimid; retval = nc_def_var(_FORCINGS_ncid,"time",NC_DOUBLE,ndims1,dimids1,&varid_time); HandleNetCDFErrors(retval); retval = nc_put_att_text(_FORCINGS_ncid,varid_time,"units",strlen(starttime),starttime); HandleNetCDFErrors(retval); @@ -2332,7 +2328,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // ---------------------------------------------------------- // time vector // ---------------------------------------------------------- - retval = nc_def_dim(_RESMB_ncid,"time",NC_UNLIMITED,&time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_RESMB_ncid,"time",ntime,&time_dimid); HandleNetCDFErrors(retval); dimids1[0] = time_dimid; retval = nc_def_var(_RESMB_ncid,"time",NC_DOUBLE,ndims1,dimids1,&varid_time); HandleNetCDFErrors(retval); retval = nc_put_att_text(_RESMB_ncid,varid_time,"units",strlen(starttime),starttime); HandleNetCDFErrors(retval); @@ -3076,6 +3072,12 @@ int NetCDFAddMetadata2D(const int fileid,const int time_dimid,int nbasins_dimid, retval = nc_inq_dimlen(fileid, time_dimid, &chunksize2[0]); HandleNetCDFErrors(retval); // Set nbasins chunksize to number ensuring that chunks have approximately 10 MB of data (assuming double precision) retval = nc_inq_dimlen(fileid, nbasins_dimid, &chunksize2[1]); HandleNetCDFErrors(retval); + + // Failsafe if time dimension is unlimited. + if (chunksize2[0] == 0) { + chunksize2[0] = 1; // Avoid division by zero if time dimension is unlimited and has no records yet + } + chunksize2[1] = max((size_t)1, min(chunksize2[1], (size_t)(NETCDF_CHUNKSIZE_MB * 1024 * 1024 / sizeof(double) / chunksize2[0]))); // Ensure at least one basin per chunk // Set chunksize retval = nc_def_var_chunking(fileid, varid, NC_CHUNKED, chunksize2); HandleNetCDFErrors(retval); From f63733c0069df41102f8f2d03c646ae35ec035f7 Mon Sep 17 00:00:00 2001 From: David Huard Date: Wed, 8 Apr 2026 09:07:06 -0400 Subject: [PATCH 08/13] Revert "Merge branch 'deflate' of github.com:CSHS-CWRA/RavenHydroFramework into deflate" This reverts commit 6b0fca9d43575a7a4259d24b678500c9a1fb7a7e, reversing changes made to 10ad2076daaa0c6c4488531a565c79db96df90c8. --- src/Assimilate.cpp | 57 ++++---- src/ChannelXSect.cpp | 28 ++-- src/Model.cpp | 4 + src/Model.h | 4 +- src/ParseInitialConditionFile.cpp | 2 +- src/ParseInput.cpp | 5 +- src/ParseLib.h | 2 +- src/ParseManagementFile.cpp | 1 - src/ParsePropertyFile.cpp | 230 +++++++++++------------------- src/RavenInclude.h | 20 +-- src/Reservoir.cpp | 59 +++++--- src/Reservoir.h | 28 ++-- src/StandardOutput.cpp | 1 + src/SubBasin.cpp | 52 ++++++- src/SubBasin.h | 1 + 15 files changed, 245 insertions(+), 249 deletions(-) diff --git a/src/Assimilate.cpp b/src/Assimilate.cpp index 73fedc8..383b65e 100644 --- a/src/Assimilate.cpp +++ b/src/Assimilate.cpp @@ -1,6 +1,6 @@ /*---------------------------------------------------------------- Raven Library Source Code - Copyright (c) 2008-2026 the Raven Development Team + Copyright (c) 2008-2025 the Raven Development Team ----------------------------------------------------------------*/ #include "RavenInclude.h" #include "Model.h" @@ -25,6 +25,8 @@ void CModel::InitializeDataAssimilation(const optStruct &Options) // Initialize Streamflow assimilation if(Options.assimilate_flow) { + _aDAscale =new double[_nSubBasins]; + _aDAscale_last=new double[_nSubBasins]; _aDAQadjust =new double[_nSubBasins]; _aDADrainSum =new double[_nSubBasins]; _aDADownSum =new double[_nSubBasins]; @@ -35,6 +37,8 @@ void CModel::InitializeDataAssimilation(const optStruct &Options) _aDAobsQ2 =new double[_nSubBasins]; _aDASinceLastBlank=new double [_nSubBasins]; for(int p=0; p<_nSubBasins; p++) { + _aDAscale [p]=1.0; + _aDAscale_last[p]=1.0; _aDAQadjust [p]=0.0; _aDADrainSum [p]=0.0; _aDADownSum [p]=0.0; @@ -77,15 +81,14 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim double Qobs2,Qmod,Qmodlast; double alpha = _pGlobalParams->GetParams()->assimilation_fact; + //alpha*=(1.0-exp(-0.06*_aDASinceLastBlank[p])); //shock/oscillation prevention + Qobs = _aDAobsQ[p]; Qobs2 = _aDAobsQ2[p]; Qmod = _pSubBasins[p]->GetOutflowRate(); Qmodlast= _pSubBasins[p]->GetLastOutflowRate(); - if (_pSubBasins[p]->GetReservoir()!=NULL){ - Qmod = _pSubBasins[p]->GetReservoir()->GetOutflowRate (false); //unadjusted flows - Qmodlast= _pSubBasins[p]->GetReservoir()->GetOldOutflowRate(false); - } if((Qmod>PRETTY_SMALL) && (Qobs!=RAV_BLANK_DATA)){ + _aDAscale [p]=1.0+alpha*((Qobs-Qmod)/Qmod); //if alpha = 1, Q=Qobs in observation basin //_aDAQadjust[p]=alpha*(Qobs-Qmod); //Option A: end of time step flow //_aDAQadjust[p]=0.5*alpha*(2.0*Qobs-Qmodlast-Qmod);//Option B: mean flow - rapidly oscillatory Q - Qmean is perfect if (Qobs2 == RAV_BLANK_DATA) { //Option C: aim for midpoint of two observation flows (this is the way) @@ -96,14 +99,19 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim } } else { + _aDAscale [p]=1.0; _aDAQadjust[p]=0.0; } + _aDAscale_last[p]=1.0;//no need to scale previous } - // actually adjust flows + // actually scale flows //--------------------------------------------------------------------- double mass_added=0.0; - if (Options.assim_method==DA_ECCC) { + if (Options.assim_method==DA_RAVEN_DEFAULT){ + mass_added=_pSubBasins[p]->ScaleAllFlows(_aDAscale[p]/_aDAscale_last[p],_aDAoverride[p],Options.timestep,tt.model_time); + } + else if (Options.assim_method==DA_ECCC) { if(_aDAoverride[p]){ mass_added=_pSubBasins[p]->AdjustAllFlows(_aDAQadjust[p],true,true,Options.timestep,tt.model_time); } @@ -133,6 +141,7 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt) bool ObsExists; //observation available in THIS basin for(p=0; p<_nSubBasins; p++) { + _aDAscale_last[p]=_aDAscale[p]; _aDADrainSum [p]=0.0; _aDAlength [p]=0.0;//reboot every timestep } @@ -155,27 +164,8 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt) { if(IsContinuousFlowObs2(_pObservedTS[i],_pSubBasins[p]->GetID()))//flow observation is available and linked to this subbasin { - Qobs = _pObservedTS[i]->GetSampledValue(nn); //mean timestep flow + Qobs = _pObservedTS[i]->GetSampledValue(nn); //mean timestep flow Qobs2 = _pObservedTS[i]->GetSampledValue(nn+1); //mean timestep flow - - //override initial conditions directly - if ((nn==1) && (Qobs!=RAV_BLANK_DATA)){ - _pSubBasins[p]->SetQout(Qobs); - if (_pSubBasins[p]->GetReservoir()!=NULL) - { - _pSubBasins[p]->GetReservoir()->SetInitialFlow(Qobs,Qobs,true,tt,Options); - } - else{ - _pSubBasins[p]->SetQout(Qobs); - - int N=_pSubBasins[p]->GetInflowHistorySize(); - double *aQobs=new double [N]; - for (int jj=0;jjSetQinHist(N,aQobs); - delete[] aQobs; - } - } - ObsExists=true; break; //avoids duplicate observations } @@ -189,6 +179,7 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt) if (ObsExists) { if((Qobs!=RAV_BLANK_DATA) && (tt.model_timeGetReachLength(); //length propagates from below } @@ -220,11 +212,13 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt) (_pSubBasins[p]->GetReservoir()==NULL) && (!_aDAoverride[p]) && (_pSubBasins[p]->IsEnabled()) && (_pSubBasins[pdown]->IsEnabled())) { + _aDAscale [p]= _aDAscale [pdown]; _aDAlength [p]+=_pSubBasins [pdown]->GetReachLength(); _aDAtimesince [p]= _aDAtimesince[pdown]; _aDAoverride [p]=false; } else{ //Nothing downstream or reservoir present in this basin, no assimilation + _aDAscale [p]=1.0; _aDAlength [p]=0.0; _aDAtimesince[p]=0.0; _aDAoverride [p]=false; @@ -270,6 +264,15 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt) _aDADownSum[p]=_aDADownSum[pdown]; } } + + // Apply time and space correction factors + //---------------------------------------------------------------- + double time_fact = _pGlobalParams->GetParams()->assim_time_decay; + double distfact = _pGlobalParams->GetParams()->assim_upstream_decay/M_PER_KM; //[1/km]->[1/m] + for(p=0; p<_nSubBasins; p++) + { + _aDAscale[p] =1.0+(_aDAscale[p]-1.0)*exp(-distfact*_aDAlength[p])*exp(-time_fact*_aDAtimesince[p]); + } } ///////////////////////////////////////////////////////////////// /// \brief Calculates updated DA scaling coefficients for all subbasins for forthcoming timestep diff --git a/src/ChannelXSect.cpp b/src/ChannelXSect.cpp index d950a88..8c9f8b5 100644 --- a/src/ChannelXSect.cpp +++ b/src/ChannelXSect.cpp @@ -260,17 +260,17 @@ double CChannelXSect::GetNPoints() const { return _nPoints; } /// \brief Returns Riverbed slope /// \return Riverbed slope [m/m] // -double CChannelXSect::GetBedslope () const { return _bedslope;} +double CChannelXSect::GetBedslope() const { return _bedslope;} -double CChannelXSect::GetAQAt (const int i) const { return _aQ[i];} +double CChannelXSect::GetAQAt(const int i) const { return _aQ[i];} -double CChannelXSect::GetAStageAt (const int i) const { return _aStage[i];} +double CChannelXSect::GetAStageAt(const int i) const { return _aStage[i];} double CChannelXSect::GetATopWidthAt(const int i) const { return _aTopWidth[i];} -double CChannelXSect::GetAXAreaAt (const int i) const { return _aXArea[i];} +double CChannelXSect::GetAXAreaAt(const int i) const { return _aXArea[i];} -double CChannelXSect::GetAPerimAt (const int i) const { return _aPerim[i];} +double CChannelXSect::GetAPerimAt(const int i) const { return _aPerim[i];} /***************************************************************** Rating Curve Interpolation functions @@ -381,15 +381,15 @@ double CChannelXSect::GetCelerity(const double &Q, const double &SB_slope,const ////////////////////////////////////////////////////////////////// /// \brief Returns diffusivity of channel [m/s] -/// \param Q [in] Channel flowrate [m3/s] -/// \param bedslope [in] Slope of riverbed [m/m] +/// \param &Q [in] Channel flowrate [m3/s] +/// \param &bedslope [in] Slope of riverbed [m/m] /// \return Diffusivity of channel at specified flowrate and slope [m3/s] // double CChannelXSect::GetDiffusivity(const double &Q, const double &SB_slope, const double &SB_n) const { ExitGracefullyIf(Q<=0,"CChannelXSect::GetDiffusivity: Invalid channel flowrate",BAD_DATA); - ///< diffusivity from Roberson et al. 1995, Hydraulic Engineering + ///< diffusivity from Roberson et al. 1995, Hydraulic Engineering \cite Roberson1998 double slope_mult=1.0; double Q_mult =1.0; GetFlowCorrections(SB_slope,SB_n,slope_mult,Q_mult); @@ -398,7 +398,7 @@ double CChannelXSect::GetDiffusivity(const double &Q, const double &SB_slope, co ////////////////////////////////////////////////////////////////// /// \brief checks to see if reference flow is in excess of channel maximum; throws warning if it is -/// \param Qref [in] Reference flowrate [m3/s] +/// \param &Qref [in] Reference flowrate [m3/s] /// \param SB_slope [in] subbasin slope (or AUTO_COMPUTE, if channel slope to be used) /// \param SB_n [in] subbasin mannings (or AUTO_COMPUTE, if channel mannings to be used) /// \param SBID [in] subbasin identifier @@ -454,15 +454,15 @@ void CChannelXSect::GetPropsFromProfile(const double &elev, Pi=sqrt(dx*dx+(zu-zl)*(zu-zl)); if ((i>0) && (_aElev[i-1]>_aElev[i]) && (_aX[i-1]==_aX[i])) //handles straight adjacent sides (left) { - if(elev<=_aElev[i-1]){Pi+=(elev -_aElev[i]);} - else {Pi+=(_aElev[i-1]-_aElev[i]);} + if(elev<=_aElev[i-1]){Pi+=(elev -_aElev[i]);} + else {Pi+=(_aElev[i-1]-_aElev[i]);} } if ((i<(_nSurveyPts-2)) && (_aElev[i+2]>_aElev[i+1]) && (_aX[i+2]==_aX[i+1])) //handles straight adjacent sides (right) { - if(elev<=_aElev[i+2]){Pi+=(elev -_aElev[i+1]);} - else {Pi+=(_aElev[i+2]-_aElev[i+1]);} + if(elev<=_aElev[i+2]){Pi+=(elev -_aElev[i+1]);} + else {Pi+=(_aElev[i+2]-_aElev[i+1]);} } - if (i==0) {Pi+=(elev-_aElev[i] );} + if (i==0) {Pi+=(elev-_aElev[i] );} if (i==_nSurveyPts-2) {Pi+=(elev-_aElev[i+1]);} } else //partially wet part of profile (includes riverbank) diff --git a/src/Model.cpp b/src/Model.cpp index 0c4fe51..ca92b36 100644 --- a/src/Model.cpp +++ b/src/Model.cpp @@ -58,6 +58,8 @@ CModel::CModel(const int nsoillayers, _aOrderedSBind =NULL; _aDownstreamInds=NULL; + _aDAscale =NULL; //Initialized in InitializeDataAssimilation + _aDAscale_last =NULL; _aDAQadjust =NULL; _aDADrainSum =NULL; _aDADownSum =NULL; @@ -214,6 +216,8 @@ CModel::~CModel() delete [] _aOutputTimes; _aOutputTimes=NULL; delete [] _aObsIndex; _aObsIndex=NULL; + delete [] _aDAscale; _aDAscale=NULL; + delete [] _aDAscale_last; _aDAscale_last=NULL; delete [] _aDAQadjust; _aDAQadjust=NULL; delete [] _aDADrainSum; _aDADrainSum=NULL; delete [] _aDADownSum; _aDADownSum=NULL; diff --git a/src/Model.h b/src/Model.h index b55659b..a6ea20a 100644 --- a/src/Model.h +++ b/src/Model.h @@ -1,6 +1,6 @@ /*---------------------------------------------------------------- Raven Library Source Code - Copyright (c) 2008-2026 the Raven Development Team + Copyright (c) 2008-2025 the Raven Development Team ----------------------------------------------------------------*/ #ifndef MODEL_H #define MODEL_H @@ -143,6 +143,7 @@ class CModel: public CModelABC int _nAggDiagnostics; ///< number of aggregated diagnostics //Data Assimilation + double *_aDAscale; ///< array of data assimilation flow scaling parameters [size: _nSubBasins] (NULL w/o DA) double *_aDAlength; ///< array of downstream distance to nearest DA observation [m] [size: _nSubBasins] (NULL w/o DA) double *_aDAQadjust; ///< array of flow adjustments [m3/s] [size: _nSubBasins] double *_aDADrainSum; ///< sum of assimilated drainage areas upstream of a subbasin outlet [km2] [size: _nSubBasins] @@ -152,6 +153,7 @@ class CModel: public CModelABC double *_aDAobsQ; ///< array of observed flow values in basins [size: _nSubBasins] (NULL w/o DA) double *_aDAobsQ2; ///< array of observed flow values in next time step [size: _nSubBasins] (NULL w/o DA) double *_aDASinceLastBlank; ///< array of time steps since last blank value occurred [size: _nSubBasins] (NULL w/o DA) + double *_aDAscale_last; ///< array of scale factors from previous time step [size: _nSubBasins] (NULL w/o DA) force_perturb**_pPerturbations; ///< array of pointers to perturbation data; defines which forcing functions to perturb and how [size: _nPerturbations] int _nPerturbations; ///< number of forcing functions to perturb diff --git a/src/ParseInitialConditionFile.cpp b/src/ParseInitialConditionFile.cpp index db06b12..511d418 100644 --- a/src/ParseInitialConditionFile.cpp +++ b/src/ParseInitialConditionFile.cpp @@ -647,7 +647,7 @@ bool ParseInitialConditionsFile(CModel *&pModel, const optStruct &Options) pBasin->GetReservoir()->SetControlFlow(s_to_i(s[1]),s_to_d(s[2]),s_to_d(s[3])); } } - else if(!strcmp(s[0],":ResDAadj")) + else if(!strcmp(s[0],":ResDAscale")) { if(Len>=3) { pBasin->GetReservoir()->SetDataAssimFactors(s_to_d(s[1]),s_to_d(s[2])); diff --git a/src/ParseInput.cpp b/src/ParseInput.cpp index 71ff525..ea5d76e 100644 --- a/src/ParseInput.cpp +++ b/src/ParseInput.cpp @@ -332,7 +332,7 @@ bool ParseMainInputFile (CModel *&pModel, Options.aNetCDFattribs =NULL; Options.assimilate_flow =false; Options.assimilate_stage =false; - Options.assim_method =DA_ECCC; + Options.assim_method =DA_RAVEN_DEFAULT; Options.assimilation_start =-1.0; //start before simulation Options.sv_override_endtime =ALMOST_INF; Options.time_zone =0; @@ -1812,7 +1812,8 @@ bool ParseMainInputFile (CModel *&pModel, {/*:AssimilationMethod [method]*/ if(Options.noisy) { cout << "Assimilation Method" << endl; } if (Len<2){ImproperFormatWarning(":AssimilationMethod",p,Options.noisy); break;} - if (!strcmp(s[1],"DA_ECCC" )){Options.assim_method=DA_ECCC;} + if (!strcmp(s[1],"DA_RAVEN_DEFAULT" )){Options.assim_method=DA_RAVEN_DEFAULT;} + else if (!strcmp(s[1],"DA_ECCC" )){Options.assim_method=DA_ECCC;} break; } case(97): //-------------------------------------------- diff --git a/src/ParseLib.h b/src/ParseLib.h index 38d0375..5ffedce 100644 --- a/src/ParseLib.h +++ b/src/ParseLib.h @@ -28,7 +28,7 @@ typedef const double * const * const Ironclad2DArray; ///< unmodifiable 2d a //************************************************************************ const int MAXINPUTITEMS = 1000; ///< maximum delimited input items per line const int MAXCHARINLINE = 10000; ///< maximum characters in line -const bool parserdebug=false; ///< turn to true for debugging of parser +const bool parserdebug=false; ///< turn to true for debugging of parser /////////////////////////////////////////////////////// /// \brief Types of parsing errors diff --git a/src/ParseManagementFile.cpp b/src/ParseManagementFile.cpp index 1c0077a..93835b7 100644 --- a/src/ParseManagementFile.cpp +++ b/src/ParseManagementFile.cpp @@ -189,7 +189,6 @@ void ParseManagementFile(CModel *&pModel,const optStruct &Options) else if(!strcmp(s[0],":UserTimeSeries")) { code=70; } else if(!strcmp(s[0],":NonlinearVariable")) { code=80; } - else if(!strcmp(s[0],":NonLinearVariable")) { code=80; } else if(!strcmp(s[0],":MaxNLSolverIterations")) { code=81; } else if(!strcmp(s[0],":NLSolverTolerance")) { code=82; } else if(!strcmp(s[0],":NLRelaxationCoeff")) { code=83; } diff --git a/src/ParsePropertyFile.cpp b/src/ParsePropertyFile.cpp index 3911732..7d4597f 100644 --- a/src/ParsePropertyFile.cpp +++ b/src/ParsePropertyFile.cpp @@ -9,8 +9,6 @@ #include "GlobalParams.h" #include "SoilAndLandClasses.h" #include - - struct val_alias { string tag; @@ -25,10 +23,9 @@ double AutoOrDoubleOrAlias(const string s, val_alias **aAl,const int nAl) } return AutoOrDouble(s); } -bool ParsePropArray (CParser *p,int *&indices,double **&properties, - int &num_read,string *tags,const int line_length,const int max_classes, - val_alias **pAliases, const int nAliases); -void DeletePropArray (int *indices,double **properties,int num_read); +bool ParsePropArray(CParser *p,int *indices,double **properties, + int &num_read,string *tags,const int line_length,const int max_classes, + val_alias **pAliases, const int nAliases); void RVPParameterWarning (string *aP, class_type *aPC, int &nP, CModel* pModel); void CreateRVPTemplate (string *aP, class_type *aPC, int &nP, const optStruct &Options); void ImproperFormatWarning(string command,CParser *p,bool noisy); @@ -37,12 +34,7 @@ void AddToMasterParamList (string *&aPm, class_type *&aPCm, int void AddNewSoilClass (CSoilClass **&pSoilClasses, soil_struct **&parsed_soils, string *&soiltags, int &num_parsed_soils, int nConstits, const string name, bool isdefault, CModel *pModel); -void AddNewLULTClass (CLandUseClass **&pLandUseClasses, surface_struct **&parsed_lults, - string *&lulttags, int &num_parsed_lults, - const string name, bool isdefault, CModel *pModel); -void AddNewVegClass (CVegetationClass **&pVegClasses, veg_struct **&parsed_vegs, - string *&vegtags, int &num_parsed_vegs, - const string name, bool isdefault, CModel *pModel); + ////////////////////////////////////////////////////////////////// /// \brief This method parses the class properties .rvp file /// @@ -63,14 +55,14 @@ bool ParseClassPropertiesFile(CModel *&pModel, global_struct parsed_globals; int num_parsed_veg=0; - CVegetationClass **pVegClasses=NULL; - veg_struct **parsed_veg =NULL; - string *vegtags =NULL; + CVegetationClass *pVegClasses[MAX_VEG_CLASSES]; + veg_struct parsed_veg [MAX_VEG_CLASSES]; + string vegtags [MAX_VEG_CLASSES]; int num_parsed_lult = 0; // counter for the number of LULT classes parsed - surface_struct **parsed_surf=NULL; - CLandUseClass **pLUClasses =NULL; - string *lulttags =NULL; + surface_struct parsed_surf [MAX_LULT_CLASSES]; + CLandUseClass *pLUClasses [MAX_LULT_CLASSES]; + string lulttags [MAX_LULT_CLASSES]; int num_parsed_soils=0; CSoilClass **pSoilClasses=NULL; @@ -83,14 +75,17 @@ bool ParseClassPropertiesFile(CModel *&pModel, string terraintags [MAX_TERRAIN_CLASSES]; int num_parsed_profiles=0; - CSoilProfile **pProfiles=NULL; + CSoilProfile *pProfiles [MAX_SOIL_PROFILES]; int num_parsed_aqstacks=0; + const int MAX_PROPERTIES_PER_LINE=20; + + int MAX_NUM_IN_CLASS=max(max(MAX_VEG_CLASSES,MAX_LULT_CLASSES),2000); bool invalid_index; int num_read; - int *indices=NULL; - double **properties=NULL; + int *indices; + double **properties; string aParamStrings[MAXINPUTITEMS]; int nParamStrings=0; @@ -100,15 +95,26 @@ bool ParseClassPropertiesFile(CModel *&pModel, pModel->GetGlobalParams()->InitializeGlobalParameters(global_template,true); pModel->GetGlobalParams()->InitializeGlobalParameters(parsed_globals,false); + CVegetationClass::InitializeVegetationProps ("[DEFAULT]",parsed_veg[0],true);//zero-index canopy is template + vegtags [0]="[DEFAULT]"; + num_parsed_veg++; - AddNewVegClass (pVegClasses, parsed_veg, vegtags, num_parsed_veg,"[DEFAULT]",true,pModel); AddNewSoilClass(pSoilClasses, parsed_soils, soiltags, num_parsed_soils, pModel->GetTransportModel()->GetNumConstituents(), "[DEFAULT]", true, pModel);//zero-index soil is template - AddNewLULTClass(pLUClasses, parsed_surf, lulttags, num_parsed_lult,"[DEFAULT]",true,pModel); + + CLandUseClass::InitializeSurfaceProperties("[DEFAULT]", parsed_surf[0], true); // zero-index LULT is template + lulttags [0]="[DEFAULT]"; + num_parsed_lult++; CTerrainClass::InitializeTerrainProperties(parsed_terrs[0],true);//zero-index terrain is template terraintags [0]="[DEFAULT]"; num_parsed_terrs++; + indices =new int [MAX_NUM_IN_CLASS]; + properties=new double *[MAX_NUM_IN_CLASS]; + for (int i=0;i=2) { - //dynamically add to stack - int tmp=num_parsed_profiles; - CSoilProfile *pNewSoil= new CSoilProfile(s[0],pModel); - DynArrayAppend((void**&)(pProfiles),(void*)(pNewSoil),tmp); + if (num_parsed_profiles>=MAX_SOIL_PROFILES-1){ + ExitGracefully("ParseClassPropertiesFile: exceeded maximum # of soil profiles",BAD_DATA);} + + pProfiles[num_parsed_profiles] = new CSoilProfile(s[0], pModel); int nhoriz=s_to_i(s[1]); ExitGracefullyIf(nhoriz<0, @@ -518,7 +524,6 @@ bool ParseClassPropertiesFile(CModel *&pModel, invalid_index=ParsePropArray(p,indices,properties,num_read,soiltags,nParamStrings,num_parsed_soils,aAliases,nAliases); ExitGracefullyIf(invalid_index, "ParseClassPropertiesFile: Invalid soiltype code in SoilParameterList command",BAD_DATA); - if (Options.noisy){ for (int j=0;j=MAX_LULT_CLASSES-1){ + ExitGracefully("ParseClassPropertiesFile: exceeded maximum # of LU/LT classes",BAD_DATA);} - AddNewLULTClass(pLUClasses,parsed_surf,lulttags,num_parsed_lult,s[0],false,pModel); - - parsed_surf[num_parsed_lult-1]->impermeable_frac = s_to_d(s[1]); + CLandUseClass::InitializeSurfaceProperties(s[0], parsed_surf[num_parsed_lult], false); + lulttags [num_parsed_lult] = s[0]; + pLUClasses[num_parsed_lult-1] = new CLandUseClass(s[0], pModel); + parsed_surf[num_parsed_lult].impermeable_frac = s_to_d(s[1]); if (Len>=3){ - parsed_surf[num_parsed_lult-1]->forest_coverage = s_to_d(s[2]); + parsed_surf[num_parsed_lult].forest_coverage = s_to_d(s[2]); } + num_parsed_lult++; } else{ ImproperFormatWarning(":LandUseClasses",p,Options.noisy); @@ -618,12 +626,11 @@ bool ParseClassPropertiesFile(CModel *&pModel, for (int i=0; i=MAX_VEG_CLASSES-1){ + ExitGracefully("ParseClassPropertiesFile: exceeded maximum # of vegetation classes",BAD_DATA);} - parsed_veg [num_parsed_veg-1]->max_height = s_to_d(s[1]); - parsed_veg [num_parsed_veg-1]->max_LAI = s_to_d(s[2]); - parsed_veg [num_parsed_veg-1]->max_leaf_cond = s_to_d(s[3]); + CVegetationClass::InitializeVegetationProps(s[0], parsed_veg [num_parsed_veg], false); + pVegClasses[num_parsed_veg-1] = new CVegetationClass(s[0], pModel); + vegtags [num_parsed_veg] = s[0]; + parsed_veg [num_parsed_veg].max_height = s_to_d(s[1]); + parsed_veg [num_parsed_veg].max_LAI = s_to_d(s[2]); + parsed_veg [num_parsed_veg].max_leaf_cond = s_to_d(s[3]); + + num_parsed_veg++; } else{ ImproperFormatWarning(":VegetationClasses",p,Options.noisy); break; @@ -691,10 +704,9 @@ bool ParseClassPropertiesFile(CModel *&pModel, "ParseClassPropertiesFile: Invalid vegetation code in SeasonalCanopyLAI command",BAD_DATA); for (int i=0;irelative_LAI[mon]=properties[i][mon]; + parsed_veg[indices[i]].relative_LAI[mon]=properties[i][mon]; } } - DeletePropArray(indices,properties,num_read); break; } case(202): //---------------------------------------------- @@ -708,10 +720,9 @@ bool ParseClassPropertiesFile(CModel *&pModel, "ParseClassPropertiesFile: Invalid vegetation code in SeasonalCanopyHeight command",BAD_DATA); for (int i=0;irelative_ht[mon]=properties[i][mon]; + parsed_veg[indices[i]].relative_ht[mon]=properties[i][mon]; } } - DeletePropArray(indices,properties,num_read); break; } case(203): //---------------------------------------------- @@ -761,13 +772,11 @@ bool ParseClassPropertiesFile(CModel *&pModel, { for (int j=0;jGetGlobalParams()->AutoCalculateGlobalParams(parsed_globals,global_template); for (int c=1;cAutoCalculateVegetationProps (*parsed_veg[c], *parsed_veg[0]); + pVegClasses [c-1]->AutoCalculateVegetationProps (parsed_veg[c],parsed_veg[0]); } for (int c=1;cAutoCalculateSoilProps (*parsed_soils[c],*parsed_soils[0],pModel->GetTransportModel()->GetNumConstituents()); + pSoilClasses[c]->AutoCalculateSoilProps (*parsed_soils[c],*parsed_soils[0],pModel->GetTransportModel()->GetNumConstituents()); } for (int c=1;cAutoCalculateLandUseProps (*parsed_surf[c], *parsed_surf[0]); + pModel->GetLanduseClass(c-1)->AutoCalculateLandUseProps(parsed_surf[c], parsed_surf[0]); } for (int c=1;cAutoCalculateTerrainProps ( parsed_terrs[c], parsed_terrs[0]); + pTerrClasses[c-1]->AutoCalculateTerrainProps (parsed_terrs[c],parsed_terrs[0]); } if (!Options.silent){ @@ -1577,7 +1588,8 @@ bool ParseClassPropertiesFile(CModel *&pModel, pModel->CheckForChannelXSectsDuplicates(Options); - + delete [] indices; + for (int i=0;i -void DynAppend(T*& arr, const T& item, const int &size) -{ - // Allocate the new array; guard it so we don't leak if something throws - std::unique_ptr new_arr(new T[size + 1]); - - for (int i = 0; i < size; ++i) { - new_arr[i] = arr[i]; - } - new_arr[size+1] = item; - - delete[] arr; - arr = new_arr.release(); //required for new arr to not be deleted when out of scope - size++; -} //Threshold Correction Functions----------------------------------- double threshPositive(const double &val); diff --git a/src/Reservoir.cpp b/src/Reservoir.cpp index aa36466..3e8f3aa 100644 --- a/src/Reservoir.cpp +++ b/src/Reservoir.cpp @@ -93,8 +93,8 @@ void CReservoir::BaseConstructor(const string Name,const long long SBID) _assimilate_stage=false; _assim_blank=true; _pObsStage=NULL; - _DAadjust=0.0; - _DAadjust_last=0.0; + _DAscale=1.0; + _DAscale_last=1.0; _dry_timesteps=0; } @@ -511,18 +511,12 @@ double CReservoir::GetLakeConvectionCoeff() const { return _lake_convcoeff; } ////////////////////////////////////////////////////////////////// /// \returns current outflow rate [m3/s] // -double CReservoir::GetOutflowRate (const bool adjusted) const { - if (adjusted){return max(_Qout+_DAadjust,0.0);} - else {return _Qout;} -} +double CReservoir::GetOutflowRate () const { return _DAscale*_Qout; } ////////////////////////////////////////////////////////////////// /// \returns previous outflow rate [m3/s] // -double CReservoir::GetOldOutflowRate (const bool adjusted) const { - if (adjusted){return max(_Qout_last+_DAadjust_last,0.0);} - else {return _Qout_last;} -} +double CReservoir::GetOldOutflowRate () const { return _DAscale_last*_Qout_last; } ////////////////////////////////////////////////////////////////// /// \returns current outflow rate from control structure i[m3/s] @@ -534,7 +528,7 @@ double CReservoir::GetControlOutflow (const int i) const { return _aQstruct[ // double CReservoir::GetIntegratedOutflow(const double& tstep) const { - return 0.5*(max(_Qout+_DAadjust,0.0)+max(_Qout_last+_DAadjust_last,0.0))*(tstep*SEC_PER_DAY); //integrated + return 0.5*(_DAscale*_Qout+_DAscale_last*_Qout_last)*(tstep*SEC_PER_DAY); //integrated } ////////////////////////////////////////////////////////////////// @@ -1149,10 +1143,10 @@ void CReservoir::SetDemandMultiplier(const double &value) ////////////////////////////////////////////////////////////////// /// \brief sets data assimilation scale factors (read from .rvc file) // -void CReservoir::SetDataAssimFactors(const double& da_adj,const double& da_adj_last) +void CReservoir::SetDataAssimFactors(const double& da_scale,const double& da_scale_last) { - _DAadjust =da_adj; - _DAadjust_last=da_adj_last; + _DAscale =da_scale; + _DAscale_last=da_scale_last; } ////////////////////////////////////////////////////////////////// @@ -1229,19 +1223,39 @@ void CReservoir::DisableOutflow() for (int i=0;i<_Np;i++){_aQ[i]=0.0;_aQunder[i]=0.0;} } ////////////////////////////////////////////////////////////////// -/// \brief adjust all internal flows by adjustment amount (for assimilation/nudging) +/// \brief scales all internal flows by scale factor (for assimilation/nudging) /// \remark Messes with mass balance something fierce! /// /// \return mass added to system [m3] // -double CReservoir::AdjustFlow(const double& Qadjust, const bool overriding, const double& tstep, const double& t) +double CReservoir::ScaleFlow(const double& scale,const bool overriding, const double& tstep, const double &t) { - _DAadjust_last=_DAadjust; - _DAadjust=Qadjust; + double va=0.0; //volume added + double sf=(scale-1.0)/scale; + + _DAscale=scale; //Estimate volume added through scaling + va+=0.5*(_Qout_last+_Qout)*sf*tstep*SEC_PER_DAY; + + return va; +} +double CReservoir::AdjustFlow(const double& Qadjust, const bool overriding, const double& tstep, const double& t) +{ + + _DAscale=1.0; + _DAscale_last=1.0; + + double scale=(_Qout+Qadjust)/_Qout; + if (_Qout==0.0){scale=1.0;} + double va=0.0; //volume added - va+=0.5*(_DAadjust_last+_DAadjust)*tstep*SEC_PER_DAY; + double sf=(scale-1.0)/scale; + + _DAscale=scale; + + //Estimate volume added through scaling + va+=0.5*(_Qout_last+_Qout)*sf*tstep*SEC_PER_DAY; return va; } @@ -1276,6 +1290,9 @@ void CReservoir::UpdateStage(const double &new_stage,const double &res_outflow, _aQstruct_last[i]=_aQstruct[i]; _aQstruct [i]=aQstruct_new[i]; } + + _DAscale_last=_DAscale; + _DAscale =1.0; } ////////////////////////////////////////////////////////////////// /// \brief returns AET [mm/d], meant to be called at end of time step @@ -1856,8 +1873,8 @@ void CReservoir::WriteToSolutionFile (ofstream &RVC) const { RVC<<" :ResFlow, "<<_Qout<<","<<_Qout_last<IsGauged()) && (_pSubBasins[p]->IsEnabled())){ if (Options.assim_method==DA_ECCC){ASSIM<<","<<_aDAQadjust[p]<<" ";} + else {ASSIM<<","<<_aDAscale [p]<<" ";} } } ASSIM<ScaleFlow(scale,overriding,tstep,t); + } + + //Estivate volume added through scaling + _channel_storage*=scale; va+=_channel_storage*sf; + _rivulet_storage*=scale; va+=_rivulet_storage*sf; + return va; +} + ////////////////////////////////////////////////////////////////// /// \brief adjusts all internal flows by corresponding magnitude (for assimilation/nudging) /// \remark Messes with mass balance something fierce! @@ -1540,7 +1588,7 @@ double CSubBasin::AdjustAllFlows(const double &adjust, const bool overriding,con va+=adjust*tstep*SEC_PER_DAY; } } - else if((overriding) && (assimsite)) //just update Qout + else if((overriding) && (assimsite)) //just update Qin { for (int i=0;i<_nSegments;i++) { diff --git a/src/SubBasin.h b/src/SubBasin.h index d51e81b..48115fd 100644 --- a/src/SubBasin.h +++ b/src/SubBasin.h @@ -317,6 +317,7 @@ class CSubBasin void SetGauged (const bool isgauged); void Disable (); void Enable (); + double ScaleAllFlows (const double &scale_factor, const bool scale_last, const double &tstep, const double &t); double AdjustAllFlows (const double &adjustment, const bool overriding, const bool assimsite, const double &tstep, const double &t); void SetUnusableFlowPercentage(const double &val); void IncludeInAssimilation (); From f72a9e077032d82b050de175e0f954dfb65348a8 Mon Sep 17 00:00:00 2001 From: David Huard Date: Wed, 8 Apr 2026 13:45:47 -0400 Subject: [PATCH 09/13] add n_out_time in Options struct --- src/RavenInclude.h | 2 +- src/StandardOutput.cpp | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/RavenInclude.h b/src/RavenInclude.h index 850dc45..346c9c6 100644 --- a/src/RavenInclude.h +++ b/src/RavenInclude.h @@ -1069,7 +1069,7 @@ struct optStruct double convergence_crit; ///< convergence criteria double max_iterations; ///< maximum number of iterations for iterative solver method double timestep; ///< numerical method timestep (in days) - int n_out_time; ///< size of output time dimension + size_t n_out_time; ///< size of output time dimension double output_interval; ///< write to output file every x number of timesteps ensemble_type ensemble; ///< ensemble type (or ENSEMBLE_NONE if single model) string external_script; ///< call to external script/.exe once per timestep (or "" if none) diff --git a/src/StandardOutput.cpp b/src/StandardOutput.cpp index 6b4b28d..13e833a 100644 --- a/src/StandardOutput.cpp +++ b/src/StandardOutput.cpp @@ -2021,7 +2021,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) string tmpFilename; int p; // loop over all sub-basins string tmp,tmp2,tmp3; - + size_t chunksize_time; // Size of output time dimension // initialize all potential file IDs with -9 == "not existing and hence not opened" _HYDRO_ncid = -9; // output file ID for Hydrographs.nc (-9 --> not opened) @@ -2069,7 +2069,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) retval = nc_def_var_deflate(_HYDRO_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); // Set chunksize to the number of time steps - size_t chunksize_time = min(1, Options.n_out_time); + chunksize_time = min(1, Options.n_out_time); retval = nc_def_var_chunking(_HYDRO_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval); // define precipitation variable @@ -2149,7 +2149,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // time // ---------------------------------------------------------- // (a) Define the DIMENSIONS. NetCDF will hand back an ID - retval = nc_def_dim(_RESSTAGE_ncid,"time",ntime,&time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_RESSTAGE_ncid,"time",Options.n_out_time,&time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. Assign units attributes to the netCDF VARIABLES. dimids1[0] = time_dimid; @@ -2231,7 +2231,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // time vector // ---------------------------------------------------------- // Define the DIMENSIONS. NetCDF will hand back an ID - retval = nc_def_dim(_STORAGE_ncid, "time", ntime, &time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_STORAGE_ncid, "time", Options.n_out_time, &time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. dimids1[0] = time_dimid; @@ -2281,7 +2281,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // ---------------------------------------------------------- // time vector // ---------------------------------------------------------- - retval = nc_def_dim(_FORCINGS_ncid,"time",ntime,&time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_FORCINGS_ncid,"time",Options.n_out_time,&time_dimid); HandleNetCDFErrors(retval); dimids1[0] = time_dimid; retval = nc_def_var(_FORCINGS_ncid,"time",NC_DOUBLE,ndims1,dimids1,&varid_time); HandleNetCDFErrors(retval); retval = nc_put_att_text(_FORCINGS_ncid,varid_time,"units",strlen(starttime),starttime); HandleNetCDFErrors(retval); @@ -2329,7 +2329,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // ---------------------------------------------------------- // time vector // ---------------------------------------------------------- - retval = nc_def_dim(_RESMB_ncid,"time",ntime,&time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_RESMB_ncid,"time",Options.n_out_time,&time_dimid); HandleNetCDFErrors(retval); dimids1[0] = time_dimid; retval = nc_def_var(_RESMB_ncid,"time",NC_DOUBLE,ndims1,dimids1,&varid_time); HandleNetCDFErrors(retval); retval = nc_put_att_text(_RESMB_ncid,varid_time,"units",strlen(starttime),starttime); HandleNetCDFErrors(retval); From c1d0585a337e59796aa78f548ceb4ed109e7ac48 Mon Sep 17 00:00:00 2001 From: David Huard Date: Thu, 9 Apr 2026 13:48:18 -0400 Subject: [PATCH 10/13] Set time dimension back to illimited. No need to define it for chunking. Add compression to other outputs (storage, rainfall) --- src/ConstituentModel.cpp | 8 ++---- src/CustomOutput.cpp | 20 ++++++------- src/CustomOutput.h | 2 +- src/Makefile | 4 +-- src/ParseInput.cpp | 2 +- src/StandardOutput.cpp | 62 ++++++++++++++++++++++++++++------------ 6 files changed, 58 insertions(+), 40 deletions(-) diff --git a/src/ConstituentModel.cpp b/src/ConstituentModel.cpp index b1ae654..52ab642 100644 --- a/src/ConstituentModel.cpp +++ b/src/ConstituentModel.cpp @@ -1004,7 +1004,6 @@ void CConstituentModel::WriteNetCDFOutputFileHeaders(const optStruct& Options) const int ndims2 = 2; int dimids1[ndims1]; // array which will contain all dimension ids for a variable int retval,ncid; - size_t ntime; // Size of the time dimension //converts start day into "hours since YYYY-MM-DD HH:MM:SS" (model start time) char starttime[200]; // start time string in format 'hours since YYY-MM-DD HH:MM:SS' @@ -1043,8 +1042,7 @@ void CConstituentModel::WriteNetCDFOutputFileHeaders(const optStruct& Options) // time vector // ---------------------------------------------------------- // Define the DIMENSIONS. NetCDF will hand back an ID - size_t ntime = min(1, Options.n_out_time); - retval = nc_def_dim(_CONC_ncid,"time",ntime,&time_dimid); HandleNetCDFErrors(retval); + retval = (_CONC_ncid,"time",NC_UNLIMITED,&time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. dimids1[0] = time_dimid; @@ -1096,7 +1094,7 @@ void CConstituentModel::WriteNetCDFOutputFileHeaders(const optStruct& Options) // time vector // ---------------------------------------------------------- // Define the DIMENSIONS. NetCDF will hand back an ID - retval = nc_def_dim(_POLLUT_ncid,"time",ntime,&time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_POLLUT_ncid,"time",NC_UNLIMITED,&time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. dimids1[0] = time_dimid; @@ -1184,7 +1182,7 @@ void CConstituentModel::WriteNetCDFOutputFileHeaders(const optStruct& Options) // time vector // ---------------------------------------------------------- // Define the DIMENSIONS. NetCDF will hand back an ID - retval = nc_def_dim(_LOADING_ncid,"time",ntime,&time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_LOADING_ncid,"time",NC_UNLIMITED,&time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. dimids1[0] = time_dimid; diff --git a/src/CustomOutput.cpp b/src/CustomOutput.cpp index a2de976..60d14ca 100644 --- a/src/CustomOutput.cpp +++ b/src/CustomOutput.cpp @@ -218,10 +218,10 @@ void CCustomOutput::SetHistogramParams(const double minv,const double maxv, cons _nBins=numBins; } ////////////////////////////////////////////////////////////////// -/// \brief Compute the approximate number of time steps for the output array +/// \brief Compute the approximate number of time steps for the output array. This is not meant to be exact. Do not allocate memory based on this value. /// \param Options [in] Global model options information /// \return Number of time steps (int) -int CCustomOutput::ComputeNumTimeSteps(const optStruct &Options) const +int CCustomOutput::ApproximateNumTimeSteps(const optStruct &Options) const { int nsteps = 0; // return value time_struct start, end; @@ -258,7 +258,7 @@ int CCustomOutput::ComputeNumTimeSteps(const optStruct &Options) const nsteps = 0; break; } - return max(0, nsteps); + return max(0, nsteps); // Ensure non-negative and account for partial time step at end of simulation } /////////////////////////////////////////////////////////////////// /// \brief Allocates memory and initialize data storage of a CCustomOutput object @@ -595,9 +595,8 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) int ndata_dimid,varid_data,varid_grps(0); // dimension ID, variable ID for simulated data and group (HRU/SB) info int retval; // error value for NetCDF routines - size_t start[1], count[1]; // determines where and how much will be written to NetCDF - size_t ntime; // Number of time steps - size_t chunksize2[2]; + size_t start[1], count[1]; // determines where and how much will be written to NetCDF // Number of time steps + size_t chunksize2[2]; // Chunksize (time, data) string tmp,tmp2,tmp3,tmp4; bool cant_support=(_aggstat==AGG_RANGE || _aggstat==AGG_95CI || _aggstat==AGG_QUARTILES || _aggstat==AGG_HISTOGRAM); @@ -619,10 +618,8 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) // ---------------------------------------------------------- // (a) Define the DIMENSIONS. NetCDF will hand back an ID for each. - // TODO: Set dimension size instead of unlimited. - ntime = max(1, ComputeNumTimeSteps(Options)); - - retval = nc_def_dim(_netcdf_ID, "time", ntime, &time_dimid); HandleNetCDFErrors(retval); + chunksize2[0] = ApproximateNumTimeSteps(Options) + 1; + retval = nc_def_dim(_netcdf_ID, "time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); // (b) Define the time variable. dimids1[0] = time_dimid; @@ -632,7 +629,7 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) retval = nc_def_var_deflate(_netcdf_ID, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); // Set chunksize to len(time) - retval = nc_def_var_chunking(_netcdf_ID, varid_time, NC_CHUNKED, &ntime); HandleNetCDFErrors(retval); + retval = nc_def_var_chunking(_netcdf_ID, varid_time, NC_CHUNKED, &chunksize2[0]); HandleNetCDFErrors(retval); // (c) Assign units attributes to the netCDF VARIABLES. @@ -692,7 +689,6 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) retval = nc_def_var_deflate(_netcdf_ID, varid_data, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); // Set chunksizes for data variable (time, ndata) - chunksize2[0] = ntime; //chunk along time dimension chunksize2[1] = max((size_t)1, min(_nData, (size_t)(NETCDF_CHUNKSIZE_MB * 1024 * 1024 / sizeof(double) / chunksize2[0]))); // Ensure at least one basin per chunk retval = nc_def_var_chunking(_netcdf_ID, varid_data, NC_CHUNKED, chunksize2); HandleNetCDFErrors(retval); diff --git a/src/CustomOutput.h b/src/CustomOutput.h index 5fea67c..461a98f 100644 --- a/src/CustomOutput.h +++ b/src/CustomOutput.h @@ -126,7 +126,7 @@ class CCustomOutput void SetHistogramParams(const double min,const double max, const int numBins); - int ComputeNumTimeSteps(const optStruct &Options) const; + int ApproximateNumTimeSteps(const optStruct &Options) const; spatial_agg GetSpatialAgg() const; diff --git a/src/Makefile b/src/Makefile index f135074..2539849 100644 --- a/src/Makefile +++ b/src/Makefile @@ -18,9 +18,9 @@ appname := Raven.exe CXX := g++ -CXXFLAGS := -Wno-deprecated +CXXFLAGS := -Wno-deprecated -g LDLIBS := -LDFLAGS := +LDFLAGS := -g # OPTION 0) some compilers require the c++11 flag, some may not CXXFLAGS += -std=c++11 -fPIC diff --git a/src/ParseInput.cpp b/src/ParseInput.cpp index ea5d76e..9087a5f 100644 --- a/src/ParseInput.cpp +++ b/src/ParseInput.cpp @@ -3731,7 +3731,7 @@ bool ParseMainInputFile (CModel *&pModel, "ParseMainInputFile::Model duration less than zero. Make sure :EndDate is after :StartDate.",BAD_DATA_WARN); // Compute the size of the output's time dimension - Options.n_out_time = (int)(floor(ceil(Options.duration/Options.timestep))/Options.output_interval); + Options.n_out_time = (int)(floor(ceil((Options.duration + TIME_CORRECTION)/Options.timestep))/Options.output_interval); if (Options.n_out_time==0) { WriteAdvisory("ParseMainInputFile::Number of output time steps is zero. Check :Duration, :Timestep, and :OutputInterval commands.",BAD_DATA); } diff --git a/src/StandardOutput.cpp b/src/StandardOutput.cpp index 13e833a..77bf9f5 100644 --- a/src/StandardOutput.cpp +++ b/src/StandardOutput.cpp @@ -2021,7 +2021,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) string tmpFilename; int p; // loop over all sub-basins string tmp,tmp2,tmp3; - size_t chunksize_time; // Size of output time dimension + size_t chunksize_time; // Size of output time dimension // initialize all potential file IDs with -9 == "not existing and hence not opened" _HYDRO_ncid = -9; // output file ID for Hydrographs.nc (-9 --> not opened) @@ -2056,22 +2056,22 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // time // ---------------------------------------------------------- // (a) Define the DIMENSIONS. NetCDF will hand back an ID - retval = nc_def_dim(_HYDRO_ncid, "time", Options.n_out_time, &time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_HYDRO_ncid, "time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. Assign units attributes to the netCDF VARIABLES. dimids1[0] = time_dimid; + // Set chunksize to the number of time steps + chunksize_time = Options.n_out_time + 1; + retval = nc_def_var(_HYDRO_ncid, "time", NC_DOUBLE, ndims1,dimids1, &varid_time); HandleNetCDFErrors(retval); + // Enable compression and chunking + retval = nc_def_var_deflate(_HYDRO_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); + retval = nc_def_var_chunking(_HYDRO_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval); + retval = nc_put_att_text(_HYDRO_ncid, varid_time, "units" , strlen(starttime) , starttime); HandleNetCDFErrors(retval); retval = nc_put_att_text(_HYDRO_ncid, varid_time, "calendar", strlen("gregorian"), "gregorian"); HandleNetCDFErrors(retval); retval = nc_put_att_text(_HYDRO_ncid, varid_time, "standard_name", strlen("time"), "time"); HandleNetCDFErrors(retval); - // Enable deflate compression for time variable (shuffle, zlib, deflate_level) - retval = nc_def_var_deflate(_HYDRO_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); - - // Set chunksize to the number of time steps - chunksize_time = min(1, Options.n_out_time); - retval = nc_def_var_chunking(_HYDRO_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval); - // define precipitation variable varid_pre= NetCDFAddMetadata(_HYDRO_ncid, time_dimid,"precip","Precipitation","mm d**-1"); @@ -2149,11 +2149,15 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // time // ---------------------------------------------------------- // (a) Define the DIMENSIONS. NetCDF will hand back an ID - retval = nc_def_dim(_RESSTAGE_ncid,"time",Options.n_out_time,&time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_RESSTAGE_ncid,"time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. Assign units attributes to the netCDF VARIABLES. dimids1[0] = time_dimid; retval = nc_def_var (_RESSTAGE_ncid,"time",NC_DOUBLE,ndims1,dimids1,&varid_time); HandleNetCDFErrors(retval); + // Enable compression and chunking + retval = nc_def_var_deflate(_RESSTAGE_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); + retval = nc_def_var_chunking(_RESSTAGE_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval); + retval = nc_put_att_text(_RESSTAGE_ncid,varid_time,"units",strlen(starttime),starttime); HandleNetCDFErrors(retval); retval = nc_put_att_text(_RESSTAGE_ncid,varid_time,"calendar",strlen("gregorian"),"gregorian"); HandleNetCDFErrors(retval); retval = nc_put_att_text(_RESSTAGE_ncid,varid_time,"standard_name",strlen("time"),"time"); HandleNetCDFErrors(retval); @@ -2231,11 +2235,15 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // time vector // ---------------------------------------------------------- // Define the DIMENSIONS. NetCDF will hand back an ID - retval = nc_def_dim(_STORAGE_ncid, "time", Options.n_out_time, &time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_STORAGE_ncid, "time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); - /// Define the time variable. + /// Define the time variable. dimids1[0] = time_dimid; retval = nc_def_var(_STORAGE_ncid, "time", NC_DOUBLE, ndims1,dimids1, &varid_time); HandleNetCDFErrors(retval); + // Enable compression and chunking + retval = nc_def_var_deflate(_STORAGE_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); + retval = nc_def_var_chunking(_STORAGE_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval); + retval = nc_put_att_text(_STORAGE_ncid, varid_time, "units" , strlen(starttime) , starttime); HandleNetCDFErrors(retval); retval = nc_put_att_text(_STORAGE_ncid, varid_time, "calendar", strlen("gregorian"), "gregorian"); HandleNetCDFErrors(retval); @@ -2281,9 +2289,14 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // ---------------------------------------------------------- // time vector // ---------------------------------------------------------- - retval = nc_def_dim(_FORCINGS_ncid,"time",Options.n_out_time,&time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_FORCINGS_ncid,"time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); + dimids1[0] = time_dimid; retval = nc_def_var(_FORCINGS_ncid,"time",NC_DOUBLE,ndims1,dimids1,&varid_time); HandleNetCDFErrors(retval); + // Enable compression and chunking + retval = nc_def_var_deflate(_FORCINGS_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); + retval = nc_def_var_chunking(_FORCINGS_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval); + retval = nc_put_att_text(_FORCINGS_ncid,varid_time,"units",strlen(starttime),starttime); HandleNetCDFErrors(retval); retval = nc_put_att_text(_FORCINGS_ncid,varid_time,"calendar",strlen("gregorian"),"gregorian"); HandleNetCDFErrors(retval); @@ -2329,9 +2342,16 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // ---------------------------------------------------------- // time vector // ---------------------------------------------------------- - retval = nc_def_dim(_RESMB_ncid,"time",Options.n_out_time,&time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_RESMB_ncid,"time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); + + dimids1[0] = time_dimid; retval = nc_def_var(_RESMB_ncid,"time",NC_DOUBLE,ndims1,dimids1,&varid_time); HandleNetCDFErrors(retval); + + // Enable compression and chunking + retval = nc_def_var_deflate(_RESMB_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); + retval = nc_def_var_chunking(_RESMB_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval); + retval = nc_put_att_text(_RESMB_ncid,varid_time,"units",strlen(starttime),starttime); HandleNetCDFErrors(retval); retval = nc_put_att_text(_RESMB_ncid,varid_time,"calendar",strlen("gregorian"),"gregorian"); HandleNetCDFErrors(retval); @@ -3026,6 +3046,7 @@ int NetCDFAddMetadata(const int fileid,const int time_dimid,string shortname,str int retval; int dimids[1]; dimids[0] = time_dimid; + size_t chunksize[1]; static double fill_val[] = {NETCDF_BLANK_VALUE}; static double miss_val[] = {NETCDF_BLANK_VALUE}; @@ -3033,6 +3054,13 @@ int NetCDFAddMetadata(const int fileid,const int time_dimid,string shortname,str // (a) create variable precipitation retval = nc_def_var(fileid,shortname.c_str(),NC_DOUBLE,1,dimids,&varid); HandleNetCDFErrors(retval); + // Set compression + retval = nc_def_var_deflate(fileid, varid, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); + // Get the chunksize of the time dimension + retval = nc_inq_var_chunking(fileid, time_dimid, NULL, &chunksize[0]); HandleNetCDFErrors(retval); + // Set chunksize to number of time steps, or 1 if time dimension is unlimited + retval = nc_def_var_chunking(fileid, varid, NC_CHUNKED, chunksize); HandleNetCDFErrors(retval); + // (b) add attributes to variable retval = nc_put_att_text (fileid,varid,"units",units.length(),units.c_str()); HandleNetCDFErrors(retval); retval = nc_put_att_text (fileid,varid,"long_name",longname.length(),longname.c_str()); HandleNetCDFErrors(retval); @@ -3070,14 +3098,10 @@ int NetCDFAddMetadata2D(const int fileid,const int time_dimid,int nbasins_dimid, // Set compression retval = nc_def_var_deflate(fileid, varid, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); // Set time chunksize to number of time steps - retval = nc_inq_dimlen(fileid, time_dimid, &chunksize2[0]); HandleNetCDFErrors(retval); + retval = nc_inq_var_chunking(fileid, time_dimid, NULL, &chunksize2[0]); HandleNetCDFErrors(retval); // Set nbasins chunksize to number ensuring that chunks have approximately 10 MB of data (assuming double precision) retval = nc_inq_dimlen(fileid, nbasins_dimid, &chunksize2[1]); HandleNetCDFErrors(retval); - // Failsafe if time dimension is unlimited. - if (chunksize2[0] == 0) { - chunksize2[0] = 1; // Avoid division by zero if time dimension is unlimited and has no records yet - } chunksize2[1] = max((size_t)1, min(chunksize2[1], (size_t)(NETCDF_CHUNKSIZE_MB * 1024 * 1024 / sizeof(double) / chunksize2[0]))); // Ensure at least one basin per chunk // Set chunksize From 9e49de0c2d5b5f74a4b9bdd5f8b4f80e49ec959e Mon Sep 17 00:00:00 2001 From: David Huard Date: Thu, 9 Apr 2026 13:49:31 -0400 Subject: [PATCH 11/13] oups --- src/ConstituentModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ConstituentModel.cpp b/src/ConstituentModel.cpp index 52ab642..7bac925 100644 --- a/src/ConstituentModel.cpp +++ b/src/ConstituentModel.cpp @@ -1042,7 +1042,7 @@ void CConstituentModel::WriteNetCDFOutputFileHeaders(const optStruct& Options) // time vector // ---------------------------------------------------------- // Define the DIMENSIONS. NetCDF will hand back an ID - retval = (_CONC_ncid,"time",NC_UNLIMITED,&time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_CONC_ncid,"time",NC_UNLIMITED,&time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. dimids1[0] = time_dimid; From 54667b6236c1ea7a8f31a8162cc6a66aa447becb Mon Sep 17 00:00:00 2001 From: David Huard Date: Thu, 9 Apr 2026 13:56:13 -0400 Subject: [PATCH 12/13] clean-up --- src/CustomOutput.cpp | 4 ++-- src/StandardOutput.cpp | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/CustomOutput.cpp b/src/CustomOutput.cpp index 60d14ca..6cad104 100644 --- a/src/CustomOutput.cpp +++ b/src/CustomOutput.cpp @@ -595,7 +595,7 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) int ndata_dimid,varid_data,varid_grps(0); // dimension ID, variable ID for simulated data and group (HRU/SB) info int retval; // error value for NetCDF routines - size_t start[1], count[1]; // determines where and how much will be written to NetCDF // Number of time steps + size_t start[1], count[1]; // determines where and how much will be written to NetCDF size_t chunksize2[2]; // Chunksize (time, data) string tmp,tmp2,tmp3,tmp4; @@ -618,7 +618,6 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) // ---------------------------------------------------------- // (a) Define the DIMENSIONS. NetCDF will hand back an ID for each. - chunksize2[0] = ApproximateNumTimeSteps(Options) + 1; retval = nc_def_dim(_netcdf_ID, "time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); // (b) Define the time variable. @@ -629,6 +628,7 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options) retval = nc_def_var_deflate(_netcdf_ID, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval); // Set chunksize to len(time) + chunksize2[0] = ApproximateNumTimeSteps(Options) + 1; retval = nc_def_var_chunking(_netcdf_ID, varid_time, NC_CHUNKED, &chunksize2[0]); HandleNetCDFErrors(retval); diff --git a/src/StandardOutput.cpp b/src/StandardOutput.cpp index 77bf9f5..c0580b4 100644 --- a/src/StandardOutput.cpp +++ b/src/StandardOutput.cpp @@ -2061,7 +2061,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) /// Define the time variable. Assign units attributes to the netCDF VARIABLES. dimids1[0] = time_dimid; // Set chunksize to the number of time steps - chunksize_time = Options.n_out_time + 1; + chunksize_time = Options.n_out_time; retval = nc_def_var(_HYDRO_ncid, "time", NC_DOUBLE, ndims1,dimids1, &varid_time); HandleNetCDFErrors(retval); // Enable compression and chunking @@ -2149,7 +2149,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // time // ---------------------------------------------------------- // (a) Define the DIMENSIONS. NetCDF will hand back an ID - retval = nc_def_dim(_RESSTAGE_ncid,"time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_RESSTAGE_ncid,"time",NC_UNLIMITED,&time_dimid); HandleNetCDFErrors(retval); /// Define the time variable. Assign units attributes to the netCDF VARIABLES. dimids1[0] = time_dimid; @@ -2289,7 +2289,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // ---------------------------------------------------------- // time vector // ---------------------------------------------------------- - retval = nc_def_dim(_FORCINGS_ncid,"time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_FORCINGS_ncid,"time",NC_UNLIMITED,&time_dimid); HandleNetCDFErrors(retval); dimids1[0] = time_dimid; retval = nc_def_var(_FORCINGS_ncid,"time",NC_DOUBLE,ndims1,dimids1,&varid_time); HandleNetCDFErrors(retval); @@ -2342,7 +2342,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options) // ---------------------------------------------------------- // time vector // ---------------------------------------------------------- - retval = nc_def_dim(_RESMB_ncid,"time", NC_UNLIMITED, &time_dimid); HandleNetCDFErrors(retval); + retval = nc_def_dim(_RESMB_ncid,"time",NC_UNLIMITED,&time_dimid); HandleNetCDFErrors(retval); dimids1[0] = time_dimid; From 1fd66f160dc81bd3ab0406019c456a1f2a50b1a8 Mon Sep 17 00:00:00 2001 From: David Huard Date: Thu, 9 Apr 2026 13:59:27 -0400 Subject: [PATCH 13/13] reset the makefile --- src/Makefile | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/Makefile b/src/Makefile index 2539849..738277a 100644 --- a/src/Makefile +++ b/src/Makefile @@ -15,19 +15,20 @@ # Feb 2025 James Craig and Bohao Tang - support for new Mac (Applie Silicon) ### # appname := raven_rev$(shell svnversion -n . | sed -e 's/:/-/g')_$(shell uname -o).exe +# To debug, add -g to CXXFLAGS and LDFLAGS appname := Raven.exe CXX := g++ -CXXFLAGS := -Wno-deprecated -g +CXXFLAGS := -Wno-deprecated LDLIBS := -LDFLAGS := -g +LDFLAGS := # OPTION 0) some compilers require the c++11 flag, some may not CXXFLAGS += -std=c++11 -fPIC # OPTION 1) include netcdf - uncomment following two commands (assumes netCDF path = /usr/local): -CXXFLAGS += -Dnetcdf -LDLIBS += -L/usr/local -lnetcdf +#CXXFLAGS += -Dnetcdf +#LDLIBS += -L/usr/local -lnetcdf # OPTION 1b) include netcdf - for newer MacOS with Apple Silicon (use with option 1 also uncommented) #CXXFLAGS += -I/opt/homebrew/include