Skip to content

Commit ce714bd

Browse files
ZeitsperreZeitsperre
authored
Merge branch 'main' into multiplatform-builds
2 parents fb122a1 + d66c853 commit ce714bd

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

56 files changed

+988
-578
lines changed

‎src/Assimilate.cpp‎

Lines changed: 75 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,14 @@
11
/*----------------------------------------------------------------
22
Raven Library Source Code
3-
Copyright (c) 2008-2021 the Raven Development Team
3+
Copyright (c) 2008-2025 the Raven Development Team
44
----------------------------------------------------------------*/
55
#include "RavenInclude.h"
66
#include "Model.h"
77

88
/*****************************************************************
99
Model Streamflow Assimilation Routines
1010
*****************************************************************/
11-
bool IsContinuousFlowObs2(const CTimeSeriesABC *pObs,long SBID)
11+
bool IsContinuousFlowObs2(const CTimeSeriesABC *pObs,long long SBID)
1212
{
1313
// clears up terribly ugly repeated if statements
1414
if(pObs==NULL) { return false; }
@@ -25,19 +25,25 @@ void CModel::InitializeDataAssimilation(const optStruct &Options)
2525
// Initialize Streamflow assimilation
2626
if(Options.assimilate_flow)
2727
{
28-
_aDAscale =new double[_nSubBasins];
29-
_aDAlength =new double[_nSubBasins];
30-
_aDAtimesince=new double[_nSubBasins];
31-
_aDAoverride =new bool [_nSubBasins];
32-
_aDAobsQ =new double[_nSubBasins];
33-
_aDAlast =new double[_nSubBasins];
28+
_aDAscale =new double[_nSubBasins];
29+
_aDAscale_last=new double[_nSubBasins];
30+
_aDAQadjust =new double[_nSubBasins];
31+
_aDADrainSum =new double[_nSubBasins];
32+
_aDADownSum =new double[_nSubBasins];
33+
_aDAlength =new double[_nSubBasins];
34+
_aDAtimesince =new double[_nSubBasins];
35+
_aDAoverride =new bool [_nSubBasins];
36+
_aDAobsQ =new double[_nSubBasins];
3437
for(int p=0; p<_nSubBasins; p++) {
35-
_aDAscale [p]=1.0;
36-
_aDAlength [p]=0.0;
37-
_aDAtimesince[p]=0.0;
38-
_aDAoverride [p]=false;
39-
_aDAobsQ [p]=0.0;
40-
_aDAlast [p]=1.0;
38+
_aDAscale [p]=1.0;
39+
_aDAscale_last[p]=1.0;
40+
_aDAQadjust [p]=0.0;
41+
_aDADrainSum [p]=0.0;
42+
_aDADownSum [p]=0.0;
43+
_aDAlength [p]=0.0;
44+
_aDAtimesince [p]=0.0;
45+
_aDAoverride [p]=false;
46+
_aDAobsQ [p]=0.0;
4147
}
4248
int count=0;
4349
for(int p=0; p<_nSubBasins; p++) {
@@ -93,17 +99,20 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim
9399
//Option B: instantaneous flow
94100
Qmod = _pSubBasins[p]->GetOutflowRate();
95101
if(Qmod>PRETTY_SMALL) {
96-
_aDAscale[p]=1.0+alpha*((Qobs-Qmod)/Qmod);
102+
_aDAscale [p]=1.0+alpha*((Qobs-Qmod)/Qmod); //if alpha = 1, Q=Qobs in observation basin
103+
_aDAQadjust[p]=alpha*(Qobs-Qmod);
97104
}
98105
else {
99-
_aDAscale[p]=1.0;
106+
_aDAscale [p]=1.0;
107+
_aDAQadjust[p]=0.0;
100108
}
101-
_aDAlast[p]=1.0;//no need to scale previous
109+
_aDAscale_last[p]=1.0;//no need to scale previous
102110
}
103111

104112
// actually scale flows
105113
//---------------------------------------------------------------------
106-
double mass_added=_pSubBasins[p]->ScaleAllFlows(_aDAscale[p]/_aDAlast[p],_aDAoverride[p],Options.timestep,tt.model_time);
114+
double mass_added=_pSubBasins[p]->ScaleAllFlows(_aDAscale[p]/_aDAscale_last[p],_aDAoverride[p],Options.timestep,tt.model_time);
115+
//double mass_added=_pSubBasins[p]->AdjustAllFlows(_aDAadjust[p],_aDAoverride[p],Options.timestep,tt.model_time);
107116
if(mass_added>0.0){_CumulInput +=mass_added/(_WatershedArea*M2_PER_KM2)*MM_PER_METER;}
108117
else {_CumulOutput-=mass_added/(_WatershedArea*M2_PER_KM2)*MM_PER_METER;}
109118

@@ -119,22 +128,24 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim
119128
//
120129
void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
121130
{
122-
if (!Options.assimilate_flow) {return;}
123-
if(tt.model_time<(Options.assimilation_start+Options.timestep/2.0)){return;}//assimilates all data after assimilation date
131+
if (!Options.assimilate_flow) {return;}
132+
if (tt.model_time<(Options.assimilation_start-Options.timestep/2.0)){return;}//assimilates all data after or on assimilation date
124133

125134
int p,pdown;
126135
double Qobs;
127136
double t_observationsOFF=ALMOST_INF;//Only used for debugging - keep as ALMOST_INF otherwise
128137

129-
int nn=(int)((tt.model_time+TIME_CORRECTION)/Options.timestep);//current timestep index
138+
int nn=(int)((tt.model_time+TIME_CORRECTION)/Options.timestep)+1;//end of timestep index
130139

131140
for(p=0; p<_nSubBasins; p++) {
132-
_aDAlast[p]=_aDAscale[p];
141+
_aDAscale_last[p]=_aDAscale[p];
142+
_aDADrainSum[p]=0.0;
133143
}
134144

135145
for(int pp=_nSubBasins-1; pp>=0; pp--)//downstream to upstream
136146
{
137147
p=GetOrderedSubBasinIndex(pp);
148+
pdown = GetSubBasinByID(_pSubBasins[p]->GetDownstreamID())->GetGlobalIndex();
138149

139150
bool ObsExists=false; //observation available in THIS basin
140151
// observations in this basin, determine scaling variables based upon blank/not blank
@@ -145,31 +156,45 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
145156
{
146157
if(IsContinuousFlowObs2(_pObservedTS[i],_pSubBasins[p]->GetID()))//flow observation is available and linked to this subbasin
147158
{
148-
Qobs = _pObservedTS[i]->GetSampledValue(nn+1);//+1: correction for period ending storage of hydrograph
159+
Qobs = _pObservedTS[i]->GetSampledValue(nn); //end of timestep flow
149160

150161
//bool fakeblank=((tt.model_time>30) && (tt.model_time<40)) || ((tt.model_time>45) && (tt.model_time<47));//TMP DEBUG
151162
//if (fakeblank){Qobs=RAV_BLANK_DATA;}
152163

153164
if((Qobs!=RAV_BLANK_DATA) && (tt.model_time<t_observationsOFF))
154165
{
155166
//_aDAscale[p] calculated live in AssimilationOverride when up-to-date modelled flow available
167+
//same with _aDAQadjust[p]
156168
_aDAlength [p]=0.0;
157169
_aDAtimesince[p]=0.0;
158170
_aDAoverride [p]=true;
159171
_aDAobsQ [p]=Qobs;
172+
_aDADrainSum [p]=0.0; //??? maybe doesnt matter
173+
if (pdown != DOESNT_EXIST) {
174+
_aDADrainSum [pdown]+=_pSubBasins[p]->GetDrainageArea(); //DOES THIS HANDLE NESTING RIGHT?
175+
}
160176
}
161177
else
162178
{ //found a blank or zero flow value
163-
_aDAscale [p]=_aDAscale[p];
179+
_aDAscale [p]=_aDAscale[p];//same adjustment as before - scaling persists
180+
_aDAQadjust [p]=_aDAQadjust[p];//same adjustment as before - flow magnitude persists
164181
_aDAtimesince[p]+=Options.timestep;
165182
_aDAlength [p]=0.0;
166183
_aDAoverride [p]=false;
167184
_aDAobsQ [p]=0.0;
185+
if (pdown != DOESNT_EXIST) {
186+
_aDADrainSum[pdown] += _aDADrainSum[p];
187+
}
168188
}
169189
ObsExists=true;
170190
break; //avoids duplicate observations
171191
}
172192
}
193+
}
194+
else {
195+
if (pdown != DOESNT_EXIST) {
196+
_aDADrainSum[pdown] += _aDADrainSum[p];
197+
}
173198
}
174199
// no observations in this basin, get scaling from downstream
175200
//----------------------------------------------------------------
@@ -178,25 +203,50 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
178203
//if ((pdown!=DOESNT_EXIST) && (!_aDAoverride[pdown])){ //alternate - allow information to pass through reservoirs
179204
if((pdown!=DOESNT_EXIST) && (_pSubBasins[p]->GetReservoir()==NULL) && (!_aDAoverride[pdown])) {
180205
_aDAscale [p]= _aDAscale [pdown];
206+
_aDAQadjust [p]= _aDAQadjust [pdown] * (_pSubBasins[p]->GetDrainageArea() / _pSubBasins[pdown]->GetDrainageArea());
181207
_aDAlength [p]+=_pSubBasins [pdown]->GetReachLength();
182208
_aDAtimesince [p]= _aDAtimesince[pdown];
183209
_aDAoverride [p]=false;
184210
}
185211
else{ //Nothing downstream or reservoir present in this basin, no assimilation
186212
_aDAscale [p]=1.0;
213+
_aDAQadjust [p]=0.0;
187214
_aDAlength [p]=0.0;
188215
_aDAtimesince[p]=0.0;
189216
_aDAoverride [p]=false;
190217
}
191218
}
192219
}// end downstream to upstream
193220

221+
for(int pp=0;pp<_nSubBasins; pp++)//upstream to downstream
222+
{
223+
p=GetOrderedSubBasinIndex(pp);
224+
pdown = GetSubBasinByID(_pSubBasins[p]->GetDownstreamID())->GetGlobalIndex();
225+
if (_aDAoverride[p]) {
226+
_aDADrainSum[p]=_pSubBasins[p]->GetDrainageArea();
227+
}
228+
else if (pdown!=DOESNT_EXIST){
229+
_aDADrainSum[p]=_aDADrainSum[pdown];
230+
}
231+
}
232+
233+
194234
// Apply time and space correction factors
195235
//----------------------------------------------------------------
196236
double time_fact = _pGlobalParams->GetParams()->assim_time_decay;
197237
double distfact = _pGlobalParams->GetParams()->assim_upstream_decay/M_PER_KM; //[1/km]->[1/m]
238+
double ECCCwt;
198239
for(p=0; p<_nSubBasins; p++)
199240
{
200-
_aDAscale[p] =1.0+(_aDAscale[p]-1.0)*exp(-distfact*_aDAlength[p])*exp(-time_fact*_aDAtimesince[p]);
241+
_aDAscale [p] =1.0+(_aDAscale[p]-1.0)*exp(-distfact*_aDAlength[p])*exp(-time_fact*_aDAtimesince[p]);
242+
//_aDAQadjust[p] =_aDAQadjust[p]*exp(-distfact*_aDAlength[p])*exp(-time_fact*_aDAtimesince[p]);
243+
244+
if (_aDADrainSum[p]!=0.0){
245+
ECCCwt = (_pSubBasins[p]->GetDrainageArea() - _aDADrainSum[p])/(_aDADownSum[p]-_aDADrainSum[p]);
246+
}
247+
else {
248+
ECCCwt=1.0;
249+
}
250+
_aDAQadjust[p] = _aDAQadjust[p]*ECCCwt;
201251
}
202252
}

‎src/ChannelXSect.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/*----------------------------------------------------------------
22
Raven Library Source Code
3-
Copyright (c) 2008-2022 the Raven Development Team
3+
Copyright (c) 2008-2025 the Raven Development Team
44
----------------------------------------------------------------*/
55
#include "ChannelXSect.h"
66
void TestManningsInfluence(const CChannelXSect *pChan,const double &Qref);
@@ -403,7 +403,7 @@ double CChannelXSect::GetDiffusivity(const double &Q, const double &SB_slope, co
403403
/// \param SB_n [in] subbasin mannings (or AUTO_COMPUTE, if channel mannings to be used)
404404
/// \param SBID [in] subbasin identifier
405405
//
406-
void CChannelXSect::CheckReferenceFlow(const double& Qref,const double& SB_slope,const double& SB_n,const long SBID) const
406+
void CChannelXSect::CheckReferenceFlow(const double& Qref,const double& SB_slope,const double& SB_n,const long long SBID) const
407407
{
408408
string warn;
409409
double junk,Q_mult;

‎src/ChannelXSect.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/*----------------------------------------------------------------
22
Raven Library Source Code
3-
Copyright (c) 2008-2022 the Raven Development Team
3+
Copyright (c) 2008-2025 the Raven Development Team
44
----------------------------------------------------------------
55
ChannelXSect.h
66
------------------------------------------------------------------
@@ -104,6 +104,6 @@ class CChannelXSect
104104
double GetCelerity (const double &Qref, const double &SB_slope,const double &SB_n) const;
105105
double GetDiffusivity(const double &Q, const double &SB_slope, const double &SB_n) const;
106106

107-
void CheckReferenceFlow(const double& Qref,const double& SB_slope,const double& SB_n,const long SBID) const;
107+
void CheckReferenceFlow(const double& Qref,const double& SB_slope,const double& SB_n,const long long SBID) const;
108108
};
109109
#endif

‎src/CommonFunctions.cpp‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -442,7 +442,10 @@ time_struct DateStringToTimeStruct(const string sDate, string sTime, const int c
442442
if (tt.month==12){tt.julian_day+=30;}
443443
if ((tt.leap_yr ) && (tt.month> 2)){tt.julian_day+= 1;}
444444

445-
if (tt.day_of_month > DAYS_PER_MONTH[tt.month - 1]) {
445+
int leap=0;
446+
if (tt.leap_yr){leap=1;}
447+
448+
if (tt.day_of_month > (DAYS_PER_MONTH[tt.month - 1]+leap)) {
446449
ExitGracefully("DateStringToTimeStruct: Invalid time format used - exceeded max day of month",BAD_DATA);
447450
}
448451
if (tt.day_of_month <=0 ) {

‎src/ConstituentModel.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/*----------------------------------------------------------------
22
Raven Library Source Code
3-
Copyright (c) 2008-2022 the Raven Development Team
3+
Copyright (c) 2008-2025 the Raven Development Team
44
----------------------------------------------------------------
55
Constituent Transport/Tracer Model class
66
coordinates information about constituent storage
@@ -12,7 +12,7 @@ coordinates information about constituent storage
1212
#include "Transport.h"
1313
#include "EnergyTransport.h"
1414

15-
bool IsContinuousConcObs(const CTimeSeriesABC *pObs,const long SBID,const int c); //Defined in StandardOutput.cpp
15+
bool IsContinuousConcObs(const CTimeSeriesABC *pObs,const long long SBID,const int c); //Defined in StandardOutput.cpp
1616
void WriteNetCDFGlobalAttributes(const int out_ncid,const optStruct& Options,const string descript);
1717
int NetCDFAddMetadata (const int fileid,const int time_dimid,string shortname,string longname,string units);
1818
int NetCDFAddMetadata2D (const int fileid,const int time_dimid,int nbasins_dimid,string shortname,string longname,string units);

‎src/ControlStructures.cpp‎

Lines changed: 11 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@
1111
/// \param downID [in] downstream subbasin
1212
/// \notes assumes SBID basin has valid reservoir
1313
//
14-
CControlStructure::CControlStructure(string name, const long SBID,const long downID) {
14+
CControlStructure::CControlStructure(string name, const long long SBID,const long long downID) {
1515
_name=name;
1616
_SBID=SBID;
1717
_target_SBID=downID;
@@ -37,13 +37,13 @@ string CControlStructure::GetName () const {return _name;}
3737
/// \brief returns target subbasin ID to which outflow is directed
3838
/// \returns target subbasin ID
3939
//
40-
long CControlStructure::GetTargetBasinID() const {return _target_SBID;}
40+
long long CControlStructure::GetTargetBasinID() const {return _target_SBID;}
4141

4242
//////////////////////////////////////////////////////////////////
4343
/// \brief sets target subbasin ID to which outflow is directed
4444
/// \param SBID [in] - valid subbasin ID (or DOESNT_EXIST, if flow is directed outside model)
4545
//
46-
void CControlStructure::SetTargetBasin(const long SBID)
46+
void CControlStructure::SetTargetBasin(const long long SBID)
4747
{
4848
_target_SBID=SBID;
4949
}
@@ -210,11 +210,11 @@ void COutflowRegime::AddRegimeConstraint(RegimeConstraint* pCond)
210210
//
211211
bool EvaluateCondition(const comparison cond, const double& v, const double& v1, const double& v2)
212212
{
213-
if (cond==COMPARE_IS_EQUAL ){return (v==v1); }
214-
else if (cond==COMPARE_NOT_EQUAL ){return (v!=v1); }
213+
if (cond==COMPARE_IS_EQUAL ){return fabs(v-v1)<REAL_SMALL; }
214+
else if (cond==COMPARE_NOT_EQUAL ){return fabs(v-v1)>REAL_SMALL; }
215215
else if (cond==COMPARE_GREATERTHAN){return (v>v1); }
216216
else if (cond==COMPARE_LESSTHAN ){return (v<v1); }
217-
else if (cond==COMPARE_BETWEEN ){return (v>=v1) && (v<=v2); }
217+
else if (cond==COMPARE_BETWEEN ){return (v>=v1) && (v<=v2); }//inclusive
218218
return false;
219219
}
220220

@@ -228,14 +228,17 @@ bool COutflowRegime::AreConditionsMet(const time_struct& tt) const
228228
for (int i = 0; i < _nConditions; i++)
229229
{
230230
comparison comp=_pConditions[i]->compare;
231-
long SBID=_pConditions[i]->basinID;
231+
long long SBID=_pConditions[i]->basinID;
232232
double v1=_pConditions[i]->compare_val1;
233233
double v2=_pConditions[i]->compare_val2;
234234
string var=_pConditions[i]->variable;
235235

236236
if (var == "DATE") {
237237
if (!EvaluateCondition(comp, tt.model_time, v1, v2)) {return false;}
238238
}
239+
else if (var == "YEAR") {
240+
if (!EvaluateCondition(comp, (double)(tt.year), v1, v2)) {return false;}
241+
}
239242
else if ((var == "DOY") || (var=="DAY_OF_YEAR")) {
240243
if (comp!=COMPARE_BETWEEN){
241244
if (!EvaluateCondition(comp, tt.julian_day, v1, v2)) {return false;}
@@ -286,7 +289,7 @@ bool COutflowRegime::AreConditionsMet(const time_struct& tt) const
286289
/// \param Q_start [in] control structure outflow at start of timestep [m3/s]
287290
/// \returns outflow, in [m3/s]
288291
//
289-
double COutflowRegime::GetOutflow(const double &h, const double &h_start, const double &Q_start, const long &target_SBID, const double &drefelev) const
292+
double COutflowRegime::GetOutflow(const double &h, const double &h_start, const double &Q_start, const long long &target_SBID, const double &drefelev) const
290293
{
291294
double tstep = _pModel->GetOptStruct()->timestep;
292295
double rivdepth = _pModel->GetSubBasinByID(target_SBID)->GetRiverDepth();

‎src/ControlStructures.h‎

Lines changed: 13 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -17,8 +17,8 @@
1717
*****************************************************************/
1818
struct RegimeCondition
1919
{
20-
string variable; ///< one of {DATE, DOY, STAGE, FLOW, ...}
21-
long basinID; ///< defaults to -1 (this basin)
20+
string variable; ///< one of {DATE, DOY, STAGE, FLOW, ...}
21+
long long basinID; ///< defaults to -1 (this basin)
2222
comparison compare; ///< e.g., greater than, between, etc.
2323
double compare_val1; ///< primary comparison value
2424
double compare_val2; ///< secondary comparison value (used for between)
@@ -78,7 +78,7 @@ class COutflowRegime
7878
void AddRegimeConstraint(RegimeConstraint *pCond);
7979

8080
bool AreConditionsMet(const time_struct &tt) const;
81-
double GetOutflow(const double &h, const double &hstart, const double &Qstart, const long &target_SBID, const double &drefelev) const;
81+
double GetOutflow(const double &h, const double &hstart, const double &Qstart, const long long &target_SBID, const double &drefelev) const;
8282
};
8383

8484

@@ -91,28 +91,28 @@ class CControlStructure
9191
{
9292
private:/*-------------------------------------------------------*/
9393
string _name; ///< structure name
94-
long _SBID; ///< subbasin ID
95-
long _target_SBID; ///< outflow directed to this subbasin ID
94+
long long _SBID; ///< subbasin ID
95+
long long _target_SBID; ///< outflow directed to this subbasin ID
9696
double _dRefElev; ///< downstream reference elevation [m] (used in some structures for backwater, limiting flow, etc.)
9797

9898
int _nRegimes; ///< number of flow regimes
9999
COutflowRegime **_aRegimes; ///< array of pointers to flow regimes with unique Q=f(h,Q,...)
100100

101101
public:/*-------------------------------------------------------*/
102102
//Constructors:
103-
CControlStructure(const string name,const long SBID, const long downID);
103+
CControlStructure(const string name,const long long SBID, const long long downID);
104104
~CControlStructure();
105105

106-
void AddRegime (COutflowRegime *pRegime);
107-
void SetTargetBasin (const long SBID);
106+
void AddRegime (COutflowRegime *pRegime);
107+
void SetTargetBasin (const long long SBID);
108108

109-
void SetDownstreamRefElevation (const double dRefElev);
109+
void SetDownstreamRefElevation (const double dRefElev);
110110

111-
string GetName () const;
112-
long GetTargetBasinID() const;
111+
string GetName () const;
112+
long long GetTargetBasinID() const;
113113

114-
double GetOutflow(const double &stage,const double &stage_start,const double &Q_start,const time_struct &tt) const;
114+
double GetOutflow(const double &stage,const double &stage_start,const double &Q_start,const time_struct &tt) const;
115115

116-
string GetCurrentRegimeName(const time_struct &tt) const;
116+
string GetCurrentRegimeName(const time_struct &tt) const;
117117
};
118118
#endif

0 commit comments

Comments
 (0)