v506 [BENCHMARKED] - ECCC assimilation, :AnnualCycleStep, speedup, NOT_UPSTREAM_OF

Fully tested ECCC propagation/assimilation. (Assimilate.cpp;Model.h)
-new routine AssimilationBackPropagate() (Solvers.cpp)
-proper read :AssimilationMethod command (ParseInput.cpp)
-write assimilation corrections (StandardOutput.cpp)

Support :AnnualCycleStep input (TimeSeries.cpp)
-Support provision of interpolation method for monthly data (Radiation.cpp;Reservoir.cpp;VegetationParams.cpp)
-changed InterpolateMo (CommonFunctions.cpp;RavenInclude.h)

Significant speedup for partially disabled models (100x speedup! for Athabasca) :
-ignore non-zero gauge weights for disabled HRUs (ForcingGrid.cpp/.h;ParseTimeSeriesFile.cpp)

new 'NOT_UPSTREAM_OF' Subbasin populate command (ParseHRUFile.cpp)

cleanup
-remove 'DEPTHRESHHOLD' spelling error support
-dont include MULTIPLIERS in template file (ParsePropertyFile.cpp)
-don't write first line of ForcingFunctions.csv (StandardOutput.cpp)

QA/QC:
-better zero time interval message (NetCDFReading.cpp)
-check for valid subbasin property in multiplier (ParseHRUFile.cpp)
-improved handling of bad subbasin in :OverrideStageDischargeCurve (ParseManagementFile.cpp)
-check for WETLAND in :SoilProfiles warnings (ParsePropertyFile.cpp)

Author: James Craig

src/Assimilate.cpp

@@ -84,7 +84,8 @@ void CModel::InitializeDataAssimilation(const optStruct &Options)
 //
 void CModel::AssimilationOverride(const int p,const optStruct& Options,const time_struct& tt)
 {
-  if(!Options.assimilate_flow)                                        { return; }
+  if (!Options.assimilate_flow)    { return; }
+  if (!_pSubBasins[p]->IsEnabled()){ return; }
 
   //Get current, up-to-date flow and calculate scaling factor
   //---------------------------------------------------------------------
@@ -93,13 +94,14 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim
     double Qobs,Qmod,Qmodlast;
     double alpha = _pGlobalParams->GetParams()->assimilation_fact;
 
-    Qobs = _aDAobsQ[p];
-    Qmod = _pSubBasins[p]->GetOutflowRate();
+    Qobs    = _aDAobsQ[p];
+    Qmod    = _pSubBasins[p]->GetOutflowRate();
     Qmodlast= _pSubBasins[p]->GetLastOutflowRate();
     if(Qmod>PRETTY_SMALL) {
       _aDAscale  [p]=1.0+alpha*((Qobs-Qmod)/Qmod); //if alpha = 1, Q=Qobs in observation basin
-      //_aDAQadjust[p]=alpha*(Qobs-Qmod);//Option A: instantaneous flow
-      _aDAQadjust[p]=alpha*(2.0*Qobs-Qmodlast-Qmod);//Option B: mean flow
+      //_aDAQadjust[p]=alpha*(Qobs-Qmod);//Option A: instantaneous flow (should set second argument to AdjustAllFlows() to true)
+      _aDAQadjust[p]=0.5*alpha*(2.0*Qobs-Qmodlast-Qmod);//Option B: mean flow
+      //_aDAQadjust[p]=1.0;//TMP DEBUG - TESTING
     }
     else {
       _aDAscale  [p]=1.0;
@@ -115,7 +117,9 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim
     mass_added=_pSubBasins[p]->ScaleAllFlows(_aDAscale[p]/_aDAscale_last[p],_aDAoverride[p],Options.timestep,tt.model_time);
   }
   else if (Options.assim_method==DA_ECCC) {
-    mass_added=_pSubBasins[p]->AdjustAllFlows(_aDAQadjust[p],_aDAoverride[p],Options.timestep,tt.model_time);
+    if(_aDAoverride[p]){
+      mass_added=_pSubBasins[p]->AdjustAllFlows(_aDAQadjust[p],true,Options.timestep,tt.model_time);
+    }
   }
 
   // 
@@ -141,127 +145,188 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
   double Qobs;
   double t_observationsOFF=ALMOST_INF;//Only used for debugging - keep as ALMOST_INF otherwise
 
-  int nn=(int)((tt.model_time+TIME_CORRECTION)/Options.timestep)+1;//end of timestep index
-
   for(p=0; p<_nSubBasins; p++) {
     _aDAscale_last[p]=_aDAscale[p];
     _aDADrainSum[p]=0.0;
   }
 
+  int nn=(int)((tt.model_time+TIME_CORRECTION)/Options.timestep)+1;//end of timestep index
+
   for(int pp=_nSubBasins-1; pp>=0; pp--)//downstream to upstream
   {
     p=GetOrderedSubBasinIndex(pp);
 
-    pdown=DOESNT_EXIST;
-    if (_pSubBasins[p]->GetDownstreamID()!=DOESNT_EXIST){
-      pdown = GetSubBasinByID(_pSubBasins[p]->GetDownstreamID())->GetGlobalIndex();
-    }
+    pdown=GetDownstreamBasin(p);
+
     bool ObsExists=false; //observation available in THIS basin
-    // observations in this basin, determine scaling variables based upon blank/not blank
-    //----------------------------------------------------------------
+
     if(_pSubBasins[p]->UseInFlowAssimilation())
     {
       for(int i=0; i<_nObservedTS; i++) //determine whether flow observation is available
       {
         if(IsContinuousFlowObs2(_pObservedTS[i],_pSubBasins[p]->GetID()))//flow observation is available and linked to this subbasin
         {
           Qobs = _pObservedTS[i]->GetSampledValue(nn); //mean timestep flow
-
-          //bool fakeblank=((tt.model_time>30) && (tt.model_time<40)) || ((tt.model_time>45) && (tt.model_time<47));//TMP DEBUG
-          //if (fakeblank){Qobs=RAV_BLANK_DATA;}
-
-          if((Qobs!=RAV_BLANK_DATA) && (tt.model_time<t_observationsOFF))
-          {
-            //_aDAscale[p] calculated live in AssimilationOverride when up-to-date modelled flow available
-            //same with _aDAQadjust[p]
-            _aDAlength   [p]=0.0;
-            _aDAtimesince[p]=0.0;
-            _aDAoverride [p]=true;
-            _aDAobsQ     [p]=Qobs;
-            _aDADrainSum [p]=0.0; //??? maybe doesnt matter
-            if (pdown != DOESNT_EXIST) {      
-              _aDADrainSum [pdown]+=_pSubBasins[p]->GetDrainageArea(); //DOES THIS HANDLE NESTING RIGHT?
-            }
-          }
-          else
-          { //found a blank or zero flow value
-            _aDAscale    [p]=_aDAscale[p];//same adjustment as before - scaling persists
-            _aDAQadjust  [p]=_aDAQadjust[p];//same adjustment as before - flow magnitude persists
-            _aDAtimesince[p]+=Options.timestep;
-            _aDAlength   [p]=0.0;
-            _aDAoverride [p]=false;
-            _aDAobsQ     [p]=0.0;
-            if (pdown != DOESNT_EXIST) {        
-              _aDADrainSum[pdown] += _aDADrainSum[p];
-            }
-          }
           ObsExists=true;
           break; //avoids duplicate observations
         }
       }
     } 
-    else {
-      if (pdown != DOESNT_EXIST) {        
-        _aDADrainSum[pdown] += _aDADrainSum[p];
+
+    pdown=GetDownstreamBasin(p);
+
+    // observations in this basin, determine scaling variables based upon blank/not blank
+    //----------------------------------------------------------------
+    if (ObsExists) {
+      if((Qobs!=RAV_BLANK_DATA) && (tt.model_time<t_observationsOFF))
+      {
+        //_aDAscale[p] calculated live in AssimilationOverride when up-to-date modelled flow available
+        _aDAlength   [p]=0.0;
+        _aDAtimesince[p]=0.0;
+        _aDAoverride [p]=true;
+        _aDAobsQ     [p]=Qobs;
+      }
+      else
+      { //found a blank or zero flow value
+        _aDAscale    [p]=_aDAscale[p];//same adjustment as before - scaling persists
+        _aDAlength   [p]=0.0;
+        _aDAtimesince[p]+=Options.timestep;
+        _aDAoverride [p]=false;
+        _aDAobsQ     [p]=0.0;
       }
     }
     // no observations in this basin, get scaling from downstream
     //----------------------------------------------------------------
-    pdown=GetDownstreamBasin(p);
-    if(ObsExists==false) { //observations may be downstream, propagate scaling upstream
+    else if(!ObsExists)  //observations may be downstream, propagate scaling upstream
+    {
       //if ((pdown!=DOESNT_EXIST) && (!_aDAoverride[pdown])){ //alternate - allow information to pass through reservoirs
-      if((pdown!=DOESNT_EXIST) && (_pSubBasins[p]->GetReservoir()==NULL) && (!_aDAoverride[pdown])) {
+      if( (pdown!=DOESNT_EXIST) && (_pSubBasins[p]->GetReservoir()==NULL) && (!_aDAoverride[p]) && (_pSubBasins[p]->IsEnabled()) && (_pSubBasins[pdown]->IsEnabled())) {
         _aDAscale      [p]= _aDAscale    [pdown];
-        _aDAQadjust    [p]= _aDAQadjust  [pdown] * (_pSubBasins[p]->GetDrainageArea() / _pSubBasins[pdown]->GetDrainageArea());
         _aDAlength     [p]+=_pSubBasins  [pdown]->GetReachLength();
         _aDAtimesince  [p]= _aDAtimesince[pdown];
-        _aDAoverride   [p]=false;
+        _aDAoverride   [p]=false; 
       }
       else{ //Nothing downstream or reservoir present in this basin, no assimilation
         _aDAscale    [p]=1.0;
-        _aDAQadjust  [p]=0.0;
         _aDAlength   [p]=0.0;
         _aDAtimesince[p]=0.0;
         _aDAoverride [p]=false;
       }
     }
   }// end downstream to upstream
-
+  
+  //Calculate _aDADrainSum, sum of assimilated drainage areas upstream of a subbasin outlet
+  // and _aDADownSum, drainage of nearest assimilated flow observation
+  // dynamic because data can disappear mid simulation
+  //-------------------------------------------------------------------
+  for(int pp=0;pp<_nSubBasins; pp++)
+  {
+    _aDADrainSum[p]=0.0;
+    _aDADownSum [p]=0.0;
+  }
   for(int pp=0;pp<_nSubBasins; pp++)//upstream to downstream
   {
     p=GetOrderedSubBasinIndex(pp);
+    pdown=GetDownstreamBasin(p);
 
-    pdown=DOESNT_EXIST;
-    if (_pSubBasins[p]->GetDownstreamID()!=DOESNT_EXIST){
-      pdown = GetSubBasinByID(_pSubBasins[p]->GetDownstreamID())->GetGlobalIndex();
-    }
     if (_aDAoverride[p]) {
       _aDADrainSum[p]=_pSubBasins[p]->GetDrainageArea();
     }
+    /*
+    else if (_pSubBasins[p]->GetReservoir()!=NULL){ //??
+      _aDADrainSum[p]=0.0;
+    }
+    */
     else if (pdown!=DOESNT_EXIST){
-      _aDADrainSum[p]=_aDADrainSum[pdown];
+      _aDADrainSum[pdown]+=_aDADrainSum[p];
     }
   }
+  for(int pp=_nSubBasins-1;pp>=0; pp--)// downstream to upstream 
+  {
+    p=GetOrderedSubBasinIndex(pp);
+    pdown=GetDownstreamBasin(p);
 
+    if (_aDAoverride[p]) {
+      _aDADownSum[p]=_pSubBasins[p]->GetDrainageArea();
+    }
+    else if (pdown!=DOESNT_EXIST){
+      _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]
-  double ECCCwt;
   for(p=0; p<_nSubBasins; p++)
   {
-    _aDAscale  [p] =1.0+(_aDAscale[p]-1.0)*exp(-distfact*_aDAlength[p])*exp(-time_fact*_aDAtimesince[p]);
-    //_aDAQadjust[p] =_aDAQadjust[p]*exp(-distfact*_aDAlength[p])*exp(-time_fact*_aDAtimesince[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
+/// \note most scale coefficients will be 1.0  except:
+/// (1) gauges with valid observation data, which scale to override OR
+/// (2) basins at or upstream of missing observation data
+/// actual scaling performed in CModel::AssimilationOverride() during routing
+/// \param tt [in] current model time structure
+/// \param Options [in] current model options structure
+//
+void CModel::AssimilationBackPropagate(const optStruct &Options,const time_struct &tt)
+{
+  int p,pdown;
 
+  if (!Options.assimilate_flow)  {return;}
+
+  for(int pp=_nSubBasins-1; pp>=0; pp--)//downstream to upstream
+  {
+    p=GetOrderedSubBasinIndex(pp);
+
+    pdown=DOESNT_EXIST;
+    if (_pSubBasins[p]->GetDownstreamID()!=DOESNT_EXIST){
+      pdown = GetSubBasinByID(_pSubBasins[p]->GetDownstreamID())->GetGlobalIndex();
+    }
+    
+    // no observations in this basin, get scaling from downstream
+    //----------------------------------------------------------------
+    pdown=GetDownstreamBasin(p);
+    long long SBID=_pSubBasins[p]->GetID();
+    if(!_aDAoverride[p]) { //observations may be downstream, propagate scaling upstream
+      if( (pdown!=DOESNT_EXIST) && (_pSubBasins[p]->GetReservoir()==NULL) && (!_aDAoverride[p]) && (_pSubBasins[p]->IsEnabled()) && (_pSubBasins[pdown]->IsEnabled())) {
+        _aDAQadjust    [p]= _aDAQadjust  [pdown] * (_pSubBasins[p]->GetDrainageArea() / _pSubBasins[pdown]->GetDrainageArea());
+        //cout<<"PROPAGATING "<<SBID<<": " <<setprecision(3)<< _aDAQadjust[p] << " from " << _aDAQadjust[pdown] << endl;
+      }
+      else{ //Nothing downstream or reservoir present in this basin, no assimilation
+        _aDAQadjust  [p]=0.0;
+      }
+    }
+    else if (_pSubBasins[p]->IsEnabled()) {
+      //cout<<"ASSIMILATING AT "<<SBID<<": " << _aDAQadjust[p] << endl;
+    }
+  } // end downstream to upstream
+
+  // 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]
+  double ECCCwt;
+  for(p=0; p<_nSubBasins; p++)
+  {
+    //_aDAQadjust[p] *=exp(-distfact*_aDAlength[p])*exp(-time_fact*_aDAtimesince[p]);
+    ECCCwt=1.0;
     if (_aDADrainSum[p]!=0.0){
       ECCCwt = (_pSubBasins[p]->GetDrainageArea() - _aDADrainSum[p])/(_aDADownSum[p]-_aDADrainSum[p]);
     }
-    else {
-      ECCCwt=1.0;
-    }
+
     if (!_aDAoverride[p]) { //no scaling for stations being overridden, only upstream
       _aDAQadjust[p] = _aDAQadjust[p]*ECCCwt;
     }
   }
-}
+
+  // Actually update flows 
+  //----------------------------------------------------------------
+  for(p=0; p<_nSubBasins; p++)
+  {
+    _pSubBasins[p]->AdjustAllFlows(_aDAQadjust[p],false,Options.timestep,tt.model_time);
+  }
+}

src/CommonFunctions.cpp

@@ -856,11 +856,12 @@ string GetCurrentMachineTime(void)
 ///
 /// \param aVal     [in] Array of doubles representing monthly data
 /// \param &tt      [in] Time structure which specifies interpolation
-/// \param &Options [in] Global model options information
+/// \param &method [in] Interpolation method
 /// \return Interpolated value at time denoted by &tt
 double InterpolateMo(const double       aVal[12],
                      const time_struct &tt,
-                     const optStruct   &Options)
+                     const monthly_interp &method,
+                     const optStruct &Options)
 {
   double wt;
   int leap(0),mo,nextmo;
@@ -870,11 +871,11 @@ double InterpolateMo(const double       aVal[12],
   month  =tt.month;
   year   =tt.year;
 
-  if      (Options.month_interp==MONTHINT_UNIFORM)//uniform over month
+  if      (method==MONTHINT_UNIFORM)//uniform over month
   {
     return aVal[month-1];
   }
-  else if (Options.month_interp==MONTHINT_LINEAR_FOM)//linear from first of month
+  else if (method==MONTHINT_LINEAR_FOM)//linear from first of month
   {
     mo=month-1;
     nextmo=mo+1;
@@ -883,13 +884,13 @@ double InterpolateMo(const double       aVal[12],
     wt=1.0-(double)(day)/(DAYS_PER_MONTH[mo]+leap);
     return wt*aVal[mo]+(1-wt)*aVal[nextmo];
   }
-  else if ((Options.month_interp==MONTHINT_LINEAR_21) ||
-           (Options.month_interp==MONTHINT_LINEAR_MID))
+  else if ((method==MONTHINT_LINEAR_21) ||
+           (method==MONTHINT_LINEAR_MID))
     //linear from 21st of month to 21st of next month (e.g., UBC_WM) or other day
   {
     double pivot=0.0;
-    if      (Options.month_interp==MONTHINT_LINEAR_21){pivot=21;}
-    else if (Options.month_interp==MONTHINT_LINEAR_MID){
+    if      (method==MONTHINT_LINEAR_21){pivot=21;}
+    else if (method==MONTHINT_LINEAR_MID){
       pivot=0.5*DAYS_PER_MONTH[month-1];
       if ((IsLeapYear(year,Options.calendar)) && (month==2)){pivot+=0.5;}
     }

src/ForcingGrid.cpp

@@ -1232,7 +1232,7 @@ void CForcingGrid::SetGridDims(const int GridDims[3])
 /// \param nHydroUnits number of HRUs
 /// \param nGridCells number of grid cells
 //
-void CForcingGrid::SetIdxNonZeroGridCells(const int nHydroUnits, const int nGridCells, const optStruct &Options)
+void CForcingGrid::SetIdxNonZeroGridCells(const int nHydroUnits, const int nGridCells, const bool  *disabledHRUs,const optStruct &Options)
 {
   int row, col;
   int minrow=_GridDims[1];
@@ -1249,7 +1249,7 @@ void CForcingGrid::SetIdxNonZeroGridCells(const int nHydroUnits, const int nGrid
   if (_GridWeight != NULL){
     for(int k=0; k<nHydroUnits; k++) {  // loop over HRUs
       for(int i=0; i<_nWeights[k]; i++) { // loop over all cells of NetCDF
-        if(_GridWeight[k][i] > 0.00001) {
+        if ((_GridWeight[k][i] > 0.00001) && (!disabledHRUs[k])){ // AND HRU not disabled!!!
           nonzero[_GridWtCellIDs[k][i]] = true;
           CellIdxToRowCol(_GridWtCellIDs[k][i],row,col);
           if(row>maxrow) { maxrow=row; }

src/ForcingGrid.h

@@ -1,6 +1,6 @@
 /*----------------------------------------------------------------
   Raven Library Source Code
-  Copyright (c) 2008-2024 the Raven Development Team
+  Copyright (c) 2008-2025 the Raven Development Team
   ----------------------------------------------------------------*/
 #ifndef FORCINGGRID_H
 #define FORCINGGRID_H
@@ -188,7 +188,9 @@ class CForcingGrid //: public CForcingGridABC
   void         SetAsPeriodEnding           ();                                      ///< set _period_ending
   void         SetIs3D(                    const bool   is3D);                      ///< set _is3D                      of class
   void         SetIdxNonZeroGridCells(     const int    nHydroUnits,
-                                           const int    nGridCells, const optStruct &Options);
+                                           const int    nGridCells, 
+                                           const bool  *disabledHRUs,
+                                           const optStruct &Options);
                                                                                     ///< set _IdxNonZeroGridCells       of class
   void         CalculateChunkSize         (const optStruct &Options);
 

src/LandUseClass.cpp

@@ -481,7 +481,6 @@ void  CLandUseClass::SetSurfaceProperty(surface_struct &S,
   else if (!name.compare("DEP_MAX"                )){S.dep_max =value;}
   else if (!name.compare("DEP_MAX_FLOW"           )){S.dep_max_flow =value;}
   else if (!name.compare("DEP_N"                  )){S.dep_n =value;}
-  else if (!name.compare("DEP_THRESHHOLD"         )){S.dep_threshold =value;}/*old typo-backward compat*/
   else if (!name.compare("DEP_THRESHOLD"          )){S.dep_threshold =value;}
   else if (!name.compare("DEP_CRESTRATIO"         )){S.dep_crestratio =value;}
   else if (!name.compare("PDMROF_B"               )){S.PDMROF_b =value; }
@@ -588,7 +587,6 @@ double CLandUseClass::GetSurfaceProperty(const surface_struct &S, string param_n
   else if (!name.compare("DEP_MAX"                )){return S.dep_max ;}
   else if (!name.compare("DEP_MAX_FLOW"           )){return S.dep_max_flow;}
   else if (!name.compare("DEP_N"                  )){return S.dep_n;}
-  else if (!name.compare("DEP_THRESHHOLD"         )){return S.dep_threshold;}/*old typo*/
   else if (!name.compare("DEP_THRESHOLD"          )){return S.dep_threshold;}
   else if (!name.compare("DEP_K"                  )){return S.dep_k;}
   else if (!name.compare("DEP_SEEP_K"             )){return S.dep_seep_k;}

src/Model.h

@@ -544,6 +544,8 @@ class CModel: public CModelABC
   void         AssimilationOverride      (const int p,
                                           const optStruct &Options, const time_struct &tt);
   void         PrepareAssimilation       (const optStruct &Options, const time_struct &tt);
+  void         AssimilationBackPropagate (const optStruct &Options, const time_struct &tt);
+
   void         PrepareForcingPerturbation(const optStruct &Options, const time_struct &tt);
   void         ApplyForcingPerturbation  (const forcing_type f, force_struct &F, const int k, const optStruct& Options, const time_struct& tt);