v525 - state variable overriding, divert_fract [BENCHMARKED]

Support for state variable overriding
-new members _nStateVarOverrides/_pStateVarOverrides (Model.h/cpp)
-new routine CModel::ApplyStateOverrides() (Model.cpp)
-new struct sv_over (Model.h)
-new options item sv_override_endtime (RavenInclude.h; ParseInput.cpp)
-new commands :StateOverrideEndTime & :OverrideState (ParseInput.cpp;ParseTimeSeriesFile.cpp)
-adjust model initialization (ModelInitialize.cpp)
-incorporate in main loop (RavenMain.cpp; Raven_BMI.cpp)

divert_fract as subbasin property (SubBasin.cpp/.h)
bug fix
- correct handling of temperature perturbations in AdjustDailyHRUForcings() (HydroUnits.cpp)

Author: James Craig

src/HydroUnits.cpp

@@ -506,7 +506,7 @@ void CHydroUnit::AdjustDailyHRUForcings(const forcing_type Ftyp, force_struct &F
     {
       double adjust=0;
       for (int n=0;n<nStepsPerDay;n++){
-        adjust+=1.0*epsilon[n]; //mean multiplicative perturbation over the day -exact for single time step
+        adjust+=epsilon[n]/nStepsPerDay; //mean multiplicative perturbation over the day -exact for single time step
       }
       F.temp_daily_min*=adjust; //approximate - uses mean perturbation to adjust extremes
       F.temp_daily_max*=adjust; //approximate - uses mean perturbation to adjust extremes
@@ -516,7 +516,7 @@ void CHydroUnit::AdjustDailyHRUForcings(const forcing_type Ftyp, force_struct &F
     {
       double adjust=0;
       for (int n=0;n<nStepsPerDay;n++){
-        adjust+=epsilon[n];
+        adjust+=epsilon[n]/nStepsPerDay;
       }
       F.temp_daily_ave+=adjust; //exact
       F.temp_daily_min+=adjust; //approximate - uses mean perturbation to adjust extremes

src/LatFlush.cpp

@@ -223,6 +223,7 @@ void CmvLatFlush::GetLateralExchange( const double * const     *state_vars, //ar
   double stor,Afrom,Ato;
   double to_stor,max_to_stor,max_rate;
   double mult=1.0; //default: flush 100%
+  int    p;
 
   for(int q=0; q<_nLatConnections; q++)
   {
@@ -235,6 +236,8 @@ void CmvLatFlush::GetLateralExchange( const double * const     *state_vars, //ar
 
     if (_divert) {
       mult=pHRUs[_kFrom[q]]->GetSurfaceProps()->divert_fract;
+      p=pHRUs[_kFrom[q]]->GetSubBasinIndex();
+      mult*=_pModel->GetSubBasin(p)->GetDivertFract();
       //cout<<" Lat Flush "<<q<<" kFrom: "<<_kFrom[q]<<" "<<_kTo[q] <<" iFrom: "<<_iFromLat[q]<<" "<< mult << " " << _aFrac[q] << " " << stor << endl;
     }
 

src/Model.cpp

@@ -33,6 +33,7 @@ CModel::CModel(const int        nsoillayers,
   _nTransParams=0;    _pTransParams=NULL;
   _nClassChanges=0;   _pClassChanges=NULL;
   _nParamOverrides=0; _pParamOverrides=NULL;
+  _nStateVarOverrides=0;_pStateVarOverrides=NULL; 
   _nObservedTS=0;     _pObservedTS=NULL; _pModeledTS=NULL; _aObsIndex=NULL;
   _nObsWeightTS =0;   _pObsWeightTS=NULL;
   _nDiagnostics=0;    _pDiagnostics=NULL;
@@ -192,12 +193,12 @@ CModel::~CModel()
   if (_aShouldApplyProcess!=NULL){
     for (k=0;k<_nProcesses;   k++){delete [] _aShouldApplyProcess[k]; } delete [] _aShouldApplyProcess;  _aShouldApplyProcess=NULL;
   }
-  for (kk=0;kk<_nHRUGroups;kk++)  {delete _pHRUGroups[kk];    } delete [] _pHRUGroups;      _pHRUGroups  =NULL;
-  for (kk=0;kk<_nSBGroups;kk++ )  {delete _pSBGroups[kk];     } delete [] _pSBGroups;       _pSBGroups  =NULL;
-  for (j=0;j<_nTransParams;j++)   {delete _pTransParams[j];   } delete [] _pTransParams;    _pTransParams=NULL;
-  for (j=0;j<_nClassChanges;j++)  {delete _pClassChanges[j];  } delete [] _pClassChanges;   _pClassChanges=NULL;
-  for (j=0;j<_nParamOverrides;j++){delete _pParamOverrides[j];} delete [] _pParamOverrides; _pParamOverrides=NULL;
-
+  for (kk=0;kk<_nHRUGroups;kk++)     {delete _pHRUGroups[kk];       } delete [] _pHRUGroups;      _pHRUGroups  =NULL;
+  for (kk=0;kk<_nSBGroups;kk++ )     {delete _pSBGroups[kk];        } delete [] _pSBGroups;       _pSBGroups  =NULL;
+  for (j=0;j<_nTransParams;j++)      {delete _pTransParams[j];      } delete [] _pTransParams;    _pTransParams=NULL;
+  for (j=0;j<_nClassChanges;j++)     {delete _pClassChanges[j];     } delete [] _pClassChanges;   _pClassChanges=NULL;
+  for (j=0;j<_nParamOverrides;j++)   {delete _pParamOverrides[j];   } delete [] _pParamOverrides; _pParamOverrides=NULL;
+  for (j=0;j<_nStateVarOverrides;j++){delete _pStateVarOverrides[j];} delete [] _pStateVarOverrides; _pStateVarOverrides=NULL;
   for (i=0;i<_nPerturbations;   i++)
   {
     delete [] _pPerturbations[i]->eps;
@@ -1308,6 +1309,16 @@ void CModel::AddParameterOverride(param_override   *pPO)
     ExitGracefully("CModel::AddParameterOverride: adding NULL override",BAD_DATA);}
 }
 //////////////////////////////////////////////////////////////////
+/// \brief Adds state override to model
+///
+/// \param *pSO [in] (valid) pointer to state override structure to be added to model
+//
+void CModel::AddStateVarOverride(sv_over   *pSO)
+{
+  if (!DynArrayAppend((void**&)(_pStateVarOverrides),(void*)(pSO),_nStateVarOverrides)){
+    ExitGracefully("CModel::AddStateVarOverride: adding NULL override",BAD_DATA);}
+}
+//////////////////////////////////////////////////////////////////
 /// \brief Adds class change to model
 ///
 /// \param *pTP [in] (valid) pointer to transient parameter to be added to model
@@ -2699,7 +2710,32 @@ void CModel::ApplyLocalParamOverrrides(const int k, const bool revert)
     }
   }
 }
-
+//////////////////////////////////////////////////////////////////
+/// \brief overrides state variables with time series values
+/// \notes called at start of timestep prior to mass energy balance
+///
+/// \param &Options [in] Global model options information
+/// \param &tt [in] Current time structure
+//
+void CModel::ApplyStateOverrrides(const optStruct   &Options,const time_struct &tt)
+{
+  if (tt.model_time>Options.sv_override_endtime-0.5*Options.timestep){return;}
+  int i,k;
+  double val;
+  int nn=(int)((tt.model_time+TIME_CORRECTION)/Options.timestep);
+  for (int j=0;j<_nStateVarOverrides;j++)
+  {
+    CHRUGroup *pGrp=GetHRUGroup(_pStateVarOverrides[j]->kk);
+    for (int k_loc=0;k_loc<pGrp->GetNumHRUs();k_loc++){
+        k=pGrp->GetHRU(k_loc)->GetGlobalIndex();
+        i= _pStateVarOverrides[j]->sv_ind;
+        val=_pStateVarOverrides[j]->pTS->GetSampledValue(nn);
+        if (val!=RAV_BLANK_DATA){
+          _pHydroUnits[k]->SetStateVarValue(i,val); //DOESNT CHECK FOR MAX OR MIN!
+        }
+    }
+  }
+}
 //////////////////////////////////////////////////////////////////
 /// \brief Recalculates HRU derived parameters
 /// \details Recalculate HRU derived parameters that are based upon time-of-year/day and SVs (storage, temp)

src/Model.h

@@ -48,6 +48,15 @@ struct class_change;
 class CTransientParam;
 class CDemandOptimizer;
 
+struct sv_over { // state variable override structure
+  CTimeSeries *pTS;    ///< time series of overridden state
+  int          kk;     ///< HRU Group
+  int          sv_ind; ///< state variable index i
+  sv_over(CTimeSeries* pT, const int kk_in, const int i) {
+    pTS=pT; kk=kk_in; sv_ind=i;
+  }
+  ~sv_over() {delete pTS;}
+};
 ////////////////////////////////////////////////////////////////////
 /// \brief Data abstraction for water surface model
 /// \details Stores and organizes HRUs and basins, provides access to all
@@ -79,6 +88,9 @@ class CModel: public CModelABC
   sv_type       *_aStateVarType;  ///< type of state variable in unit i  [size:_nStateVars]
   int          *_aStateVarLayer;  ///< index of state variable for multilayer variables (e.g., SOIL); [size:_nStateVars] value=DOESNT_EXIST(-1) for unique variables (e.g. SURFACE_WATER)
   int         _aStateVarIndices[MAX_STATE_VAR_TYPES][MAX_SV_LAYERS]; ///< lookup table for state variable indices; the index of SOIL[1] in a state_var[] array may be returned by aStateVarIndices[(int)(SOIL)][1]
+  
+  sv_over **_pStateVarOverrides;  ///< array of pointers to state variable override time series
+  int       _nStateVarOverrides;  ///< number of state variable overrides
 
   int                _nSoilVars;  ///< number of soil layer storage units
 
@@ -460,11 +472,12 @@ class CModel: public CModelABC
   void    AddStateVariables         (const sv_type             *aSV,
                                      const int                 *aLev,
                                      const int                  nSV             );
-  void    AddLandUseClass           (       CLandUseClass      *pLU             );
+  void    AddLandUseClass           (      CLandUseClass       *pLU             );
   void    AddCustomOutput           (      CCustomOutput       *pCO             );
   void    AddCustomTable            (      CCustomTable        *pTab            );
   void    AddTransientParameter     (      CTransientParam     *pTP             );
   void    AddParameterOverride      (      param_override      *pPO             );
+  void    AddStateVarOverride       (      sv_over             *pSO             );
   void    AddPropertyClassChange    (const string             HRUgroup,
                                      const class_type         tclass,
                                      const string             new_class,
@@ -514,6 +527,8 @@ class CModel: public CModelABC
 
   //called during simulation:
   //critical simulation routines (called once during each timestep):
+  void        ApplyStateOverrrides       (const optStruct   &Options,
+                                          const time_struct &tt);
   void        UpdateTransientParams      (const optStruct   &Options,
                                           const time_struct &tt);
   void        UpdateHRUForcingFunctions  (const optStruct   &Options,

src/ModelInitialize.cpp

@@ -125,6 +125,9 @@ void CModel::Initialize(const optStruct &Options)
   for (kk=0;kk<_nHRUGroups;kk++){_pHRUGroups  [kk]->Initialize(); } //disables HRUs
   for (pp=0;pp<_nSBGroups; pp++){_pSBGroups   [pp]->Initialize(); } //disables SBs and HRUs
 
+  for (j=0;j<_nStateVarOverrides; j++){
+    _pStateVarOverrides[j]->pTS->Initialize(Options.julian_start_day, Options.julian_start_year, Options.duration, Options.timestep,true,Options.calendar);
+  }
   InitializeParameterOverrides();
 
   // Forcing grids are not "Initialized" here because the derived data have to be populated everytime a new chunk is read

src/ParseInput.cpp

@@ -330,6 +330,7 @@ bool ParseMainInputFile (CModel     *&pModel,
   Options.assimilate_stage        =false;
   Options.assim_method            =DA_RAVEN_DEFAULT;
   Options.assimilation_start      =-1.0; //start before simulation
+  Options.sv_override_endtime     =ALMOST_INF;
   Options.time_zone               =0;
   Options.rvl_read_frequency      =0.0; //do not read at all
   Options.custom_interval         =1.0; //daily
@@ -455,9 +456,9 @@ bool ParseMainInputFile (CModel     *&pModel,
     else if  (!strcmp(s[0],":SnapshotHydrograph"        )){code=66; }
     else if  (!strcmp(s[0],":SuppressWarnings"          )){code=68; }
     else if  (!strcmp(s[0],":SuppressOutput"            )){code=69; }
-    else if  (!strcmp(s[0],":DontWriteWatershedStorage" )){code=70; }//*//avoid writing WatershedStorage.csv
+    else if  (!strcmp(s[0],":DontWriteWatershedStorage" )){code=70; }//avoid writing WatershedStorage.csv
     else if  (!strcmp(s[0],":EvaluationMetrics"         )){code=71; }
-    else if  (!strcmp(s[0],":EvaluationTime"            )){code=72; }//After StartDate or JulianStartDay and JulianStartYear commands
+    else if  (!strcmp(s[0],":EvaluationTime"            )){code=72; }//After StartDate commands
     else if  (!strcmp(s[0],":EvaluationPeriod"          )){code=73; }
     else if  (!strcmp(s[0],":SuppressOutputICs"         )){code=75; }
     else if  (!strcmp(s[0],":WaterYearStartMonth"       )){code=76; }
@@ -488,6 +489,7 @@ bool ParseMainInputFile (CModel     *&pModel,
     else if  (!strcmp(s[0],":FEWSParamInfoFile"         )){code=111;}
     else if  (!strcmp(s[0],":FEWSBasinStateInfoFile"    )){code=112;}
     else if  (!strcmp(s[0],":TimeOfConcentrationMethod" )){code=113;}
+    else if  (!strcmp(s[0],":StateOverrideEndTime"      )){code=114;}//AFTER :StartDate,:Calendar commands
 
     else if  (!strcmp(s[0],":WriteGroundwaterHeads"     )){code=510;}//GWMIGRATE -TO REMOVE
     else if  (!strcmp(s[0],":WriteGroundwaterFlows"     )){code=511;}//GWMIGRATE -TO REMOVE
@@ -1959,6 +1961,18 @@ bool ParseMainInputFile (CModel     *&pModel,
       else {ExitGracefully("ParseInput :TimeOfConcentrationMethod: Unrecognized method",BAD_DATA_WARN);}
       break;
     }
+   case(114):  //--------------------------------------------
+    {/*:StateOverrideEndTime [yyyy-mm-dd] [00:00:00]*///AFTER StartDate or JulianStartDay and JulianStartYear commands
+      if(Options.noisy) { cout << "State Override End Time" << endl; }
+      if(Len<3) { ImproperFormatWarning(":StateOverrideEndTime",p,Options.noisy); break; }
+
+      time_struct tt;
+      tt = DateStringToTimeStruct(s[1],s[2],Options.calendar);
+      Options.sv_override_endtime=TimeDifference(Options.julian_start_day,Options.julian_start_year,tt.julian_day,tt.year,Options.calendar);
+      if(Options.forecast_shift!=0) { Options.sv_override_endtime+=Options.forecast_shift; }
+
+      break;
+    }
     case(160):  //--------------------------------------------
     {/*:rvh_Filename [filename.rvh]*/
       if(Options.noisy) { cout <<"rvh filename: "<<s[1]<<endl; }

src/ParseTimeSeriesFile.cpp

@@ -138,6 +138,7 @@ bool ParseTimeSeriesFile(CModel *&pModel, const optStruct &Options)
     else if  (!strcmp(s[0],":AssimilateStreamflow"        )){code=101;}
     else if  (!strcmp(s[0],":SpecifyGaugeForHRU"          )){code=102;}
     else if  (!strcmp(s[0],":AssimilateReservoirStage"    )){code=103;}
+    else if  (!strcmp(s[0],":OverrideState"               )){code=104;}
     //-----------------TRANSPORT--------------------------------
     else if  (!strcmp(s[0],":FixedConcentrationTimeSeries")){code=300; const_name=s[1];}
     else if  (!strcmp(s[0],":FixedTemperatureTimeSeries"  )){code=300; const_name="TEMPERATURE";}
@@ -1291,6 +1292,32 @@ bool ParseTimeSeriesFile(CModel *&pModel, const optStruct &Options)
       }
       break;
     }
+    case(104)://----------------------------------------------
+    {/*:OverrideState  [state_var] [HRUGroup]
+       {yyyy-mm-dd} {hh:mm:ss.0} {double timestep} {int nMeasurements}
+       {double val} x nMeasurements [mm]
+       :EndOverrideState
+     */
+      if (Options.noisy) { cout <<"Override simulated state"<<endl; }
+
+      int kk=DOESNT_EXIST;
+      string grpname=to_string(s[2]);
+      CHRUGroup *pHRUgrp=pModel->GetHRUGroup(grpname);
+      if (pHRUgrp==NULL){
+        ExitGracefully("ParseTimeSeriesFile: bad HRU group name provided in :OverrideState command",BAD_DATA_WARN);
+        break;
+      }
+      kk=pHRUgrp->GetGlobalIndex();
+      pTimeSer=CTimeSeries::Parse(p,false,"OverrideState_"+to_string(s[0])+"_"+grpname, kk, "none", Options);
+      int layer_ind;
+      sv_type typ=pModel->GetStateVarInfo()->StringToSVType(s[1],layer_ind,false);
+
+      int i=pModel->GetStateVarIndex(typ,layer_ind);
+      sv_over *pSO=new sv_over(pTimeSer,kk,i);
+      pModel->AddStateVarOverride(pSO);
+      
+      break;
+    }
     case (300)://----------------------------------------------
     {/*:FixedConcentrationTimeSeries
        :FixedConcentrationTimeSeries [string constit_name] [string state_var (storage compartment)] {optional HRU Group name}

src/RavenInclude.h

@@ -1193,6 +1193,7 @@ struct optStruct
   int              nNetCDFattribs;            ///< size of array of NetCDF attributes
   int              NetCDF_chunk_mem;          ///< [MB] size of memory chunk for each forcing grid
   bool             in_bmi_mode;               ///< true if in BMI mode (no rvt files, no end time)
+  double           sv_override_endtime;       ///< model time [d] after which state variable overrides are disabled (default: 1e99) 
 };
 
 ///////////////////////////////////////////////////////////////////

src/RavenMain.cpp

@@ -136,6 +136,7 @@ int main(int argc, char* argv[])
     for(t=t_start; t<Options.duration-TIME_CORRECTION; t+=Options.timestep)  // in [d]
     {
       pModel->UpdateTransientParams      (Options,tt);
+      pModel->ApplyStateOverrrides       (Options,tt);
       pModel->RecalculateHRUDerivedParams(Options,tt);
       pModel->GetEnsemble()->StartTimeStepOps(pModel,Options,tt,e);
       pModel->UpdateHRUForcingFunctions  (Options,tt);

src/Raven_BMI.cpp

@@ -215,6 +215,7 @@ void CRavenBMI::Initialize(std::string config_file)
 void CRavenBMI::Update()
 {
   pModel->UpdateTransientParams      (Options,tt);
+  pModel->ApplyStateOverrrides       (Options,tt);
   pModel->RecalculateHRUDerivedParams(Options,tt);
   pModel->UpdateHRUForcingFunctions  (Options,tt);
   pModel->UpdateDiagnostics          (Options,tt);