v531 - Assimilation + management fix

support reservoir stage assimilation with management-derived outflow in same basin (Reservoir.cpp)

bug fix:
time shift issue in CModel::GenerateAveSubdailyTempFromMinMax()

cleanup:
converted CForcingGrid::GetValue_*** functions to take start index rather than raw time (as was actually happening anyhow) (ForcingGrid.cpp/.h)

Author: James Craig

src/ForcingGrid.cpp

@@ -1843,8 +1843,8 @@ int  CForcingGrid::GetChunkSize() const{return _ChunkSize;}
 int CForcingGrid::GetnHydroUnits() const{return _nHydroUnits;}
 
 ///////////////////////////////////////////////////////////////////
-/// \brief returns time index idx corresponding to t+interval/2
-/// \return time index idx corresponding to t+interval/2
+/// \brief returns time index idx corresponding to t
+/// \return time index idx corresponding to t
 /// param t [in] model time [days]
 //
 int CForcingGrid::GetTimeIndex(const double &t) const
@@ -1977,9 +1977,9 @@ double CForcingGrid::GetValue(const int ic, const int it) const
 /// \param n      [in] Number of time steps
 /// \return Magnitude of time series data point for which t is an index
 //
-double CForcingGrid::GetValue_avg(const int ic, const double &t_idx, const int nsteps) const
+double CForcingGrid::GetValue_avg(const int ic, const int t_idx, const int nsteps) const
 {
-  int it_start=max((int)(t_idx),0);
+  int it_start=max(t_idx,0);
   int lim=min(nsteps,_ChunkSize-it_start);
   double sum = 0.0;
   for (int it=it_start; it<it_start+lim;it++){
@@ -1993,14 +1993,14 @@ double CForcingGrid::GetValue_avg(const int ic, const double &t_idx, const int n
 ///////////////////////////////////////////////////////////////////
 /// \brief Returns minimum of n timesteps of time series data point for which t is an index
 /// \param ic     [in] Index of grid cell with non-zero weighting (value between 0 and _nNonZeroWeightedGridCells)
-/// \param t      [in] Time index
-/// \param n      [in] Number of time steps
+/// \param t_idx  [in] Time index
+/// \param n      [in] Number of time intervals
 /// \return Magnitude of time series data point for which t is an index
 //
-double CForcingGrid::GetValue_min(const int ic, const double &t, const int nsteps) const
+double CForcingGrid::GetValue_min(const int ic, const int t_idx, const int nsteps) const
 {
   double min_val = ALMOST_INF ;
-  int it_start=max((int)(t),0);
+  int it_start=max(t_idx,0);
   int lim=min(nsteps,_ChunkSize-it_start);
   for (int it=it_start; it<it_start+lim;it++){
     if(_aVal[it][ic] < min_val){min_val=_aVal[it][ic];}
@@ -2011,14 +2011,14 @@ double CForcingGrid::GetValue_min(const int ic, const double &t, const int nstep
 ///////////////////////////////////////////////////////////////////
 /// \brief Returns maximum of n timesteps of time series data point for which t is an index
 /// \param ic     [in] Index of grid cell with non-zero weighting (value between 0 and _nNonZeroWeightedGridCells)
-/// \param t      [in] Time index
+/// \param t_idx  [in] Time index
 /// \param n      [in] Number of time steps
 /// \return Magnitude of time series data point for which t is an index
 //
-double CForcingGrid::GetValue_max(const int ic, const double &t, const int nsteps) const
+double CForcingGrid::GetValue_max(const int ic, const int t_idx, const int nsteps) const
 {
   double max_val = -ALMOST_INF ;
-  int it_start=max((int)(t),0);
+  int it_start=max(t_idx,0);
   int lim=min(nsteps,_ChunkSize-it_start);
   for (int it=it_start; it<it_start+lim;it++){
     if(_aVal[it][ic] > max_val){max_val=_aVal[it][ic];}

src/ForcingGrid.h

@@ -149,9 +149,9 @@ class CForcingGrid //: public CForcingGridABC
 
   // accessors
   double GetValue                   (const int ic, const int it) const;
-  double GetValue_avg               (const int ic, const double &t, const int n) const;
-  double GetValue_min               (const int ic, const double &t, const int n) const;
-  double GetValue_max               (const int ic, const double &t, const int n) const;
+  double GetValue_avg               (const int ic, const int t_idx, const int n) const;
+  double GetValue_min               (const int ic, const int t_idx, const int n) const;
+  double GetValue_max               (const int ic, const int t_idx, const int n) const;
 
   // Weighting matrix associated routines
   void   AllocateWeightArray(              const int        nHydroUnits,

src/ModelForcingGrids.cpp

@@ -52,10 +52,6 @@ CForcingGrid *CModel::ForcingCopyCreate(const CForcingGrid *pGrid,
     pTout=GetForcingGrid(typ);
   }
 
-  //int nGridHRUs=pGrid->GetnHydroUnits();
-  //int nCells   =pGrid->GetRows()*pGrid->GetCols(); //handles 3D or 2D
-  //cout<<"**ForcingGrid Copy Create : "<<ForcingToString(pGrid->GetForcingType())<<"-->"<<ForcingToString(pTout->GetForcingType())<< " is new?:"<<is_new<<"**"<<endl;
-
   return pTout;
 }
 
@@ -126,6 +122,7 @@ void CModel::GenerateGriddedPrecipVars(const optStruct &Options)
 /// \brief Creates all missing gridded temperature data based on gridded information available,
 ///        e.g when only sub-daily temperature is provided estimate daily average, minimum and maximum temperature.
 ///        data are assumed to have the same resolution and hence can be initialized together.
+///        called with each new chunk of data 
 ///
 /// \param Options [in]  major options of the model
 //
@@ -261,10 +258,10 @@ void CModel::GenerateAveSubdailyTempFromMinMax(const optStruct &Options)
     for (int it=0; it<nVals; it++) {                   // loop over all time points (nVals)
       for (int ic=0; ic<nNonZero; ic++){               // loop over non-zero grid cell indexes
         time_idx_chunk = int(floor(t+TIME_CORRECTION));
-        Tmin   = pTmin->GetValue_avg(ic, floor(t +time_shift+TIME_CORRECTION)*nValsPerDay, nValsPerDay);
-        Tmax   = pTmax->GetValue_avg(ic, floor(t +time_shift+TIME_CORRECTION)*nValsPerDay, nValsPerDay);
-        T1corr = pTave->DailyTempCorrection(t+time_shift);
-        T2corr = pTave->DailyTempCorrection(t+time_shift+Options.timestep);
+        Tmin   = pTmin->GetValue_avg(ic, floor(t-time_shift+TIME_CORRECTION)*nValsPerDay, nValsPerDay);
+        Tmax   = pTmax->GetValue_avg(ic, floor(t-time_shift+TIME_CORRECTION)*nValsPerDay, nValsPerDay);
+        T1corr = pTave->DailyTempCorrection(t);
+        T2corr = pTave->DailyTempCorrection(t+Options.timestep);
         val=pTave_daily->GetValue(ic, time_idx_chunk)+0.25*(Tmax-Tmin)*(T1corr+T2corr);
         pTave->SetValue( ic, it, val);
       }
@@ -280,8 +277,8 @@ void CModel::GenerateAveSubdailyTempFromMinMax(const optStruct &Options)
       pTave = ForcingCopyCreate(pTave_daily,F_TEMP_AVE,1.0,nVals,Options);
 
       int nNonZero  =pTave->GetNumberNonZeroGridCells();
-      for(int it=0; it<nVals; it++) {                                 // loop over time points in buffer
-        for(int ic=0; ic<nNonZero; ic++) {                                  // loop over non-zero grid cell indexes
+      for(int it=0; it<nVals; it++) {                           // loop over time points in buffer
+        for(int ic=0; ic<nNonZero; ic++) {                      // loop over non-zero grid cell indexes
           pTave->SetValue(ic,it,pTave_daily->GetValue(ic,it));  // --> just copy daily average values
         }
       }
@@ -308,18 +305,20 @@ void CModel::GenerateMinMaxAveTempFromSubdaily(const optStruct &Options)
   pTmax_daily = ForcingCopyCreate(pTave,F_TEMP_DAILY_MAX,1.0,nVals,Options);
   pTave_daily = ForcingCopyCreate(pTave,F_TEMP_DAILY_AVE,1.0,nVals,Options);
 
+  int t_idx;
   double time;
   double time_shift=Options.julian_start_day-floor(Options.julian_start_day+TIME_CORRECTION);
-  int    nsteps_in_day=int(1.0/interval);
+  int    nsteps_in_day=(int)(rvn_round(1.0/interval));
 
   for (int it=0; it<nVals; it++) {                    // loop over time points in buffer
-    time=(double)it*1.0/interval;
+    time=(double)it*nsteps_in_day;
     time=floor(time-time_shift+TIME_CORRECTION); //model time corresponding to 00:00 on day of it
+    t_idx=(int)(time);
 
     for (int ic=0; ic<pTave->GetNumberNonZeroGridCells(); ic++){         // loop over non-zero grid cell indexes
-      pTmin_daily->SetValue(ic, it, pTave->GetValue_min(ic, time, nsteps_in_day));
-      pTmax_daily->SetValue(ic, it, pTave->GetValue_max(ic, time, nsteps_in_day));
-      pTave_daily->SetValue(ic, it, pTave->GetValue_avg(ic, time, nsteps_in_day));
+      pTmin_daily->SetValue(ic, it, pTave->GetValue_min(ic, t_idx, nsteps_in_day));
+      pTmax_daily->SetValue(ic, it, pTave->GetValue_max(ic, t_idx, nsteps_in_day));
+      pTave_daily->SetValue(ic, it, pTave->GetValue_avg(ic, t_idx, nsteps_in_day));
     }
   }
 

src/RavenInclude.h

@@ -235,10 +235,10 @@ const double  WIND_EXTINCT            =2.5;
 /// \note This is always 2 for broadleaves, and ranges from 2 for flat needles to pi for cylindrical needles.
 const double  LEAF_PROJ_RAT           =2.2;                                     ///< [m2/m2] ratio of total leaf area to projected area.
 
-const double HBV_REFERENCE_ELEV       =1000;                                    ///< [masl]  (zref in HBV)
-const double HBV_PET_ELEV_CORR        =0.0005;                                  ///< [mm/m-d] PET lapse rate (ECALT in HBV-EC)
-const double HBV_PET_TEMP_CORR        =0.5;                                     ///< [-]     (ETF in HBV-EC)
-const double GLOBAL_PET_CORR          =1.0;                                     ///< [-]     (ECORR in HBV-EC)
+const double  HBV_REFERENCE_ELEV      =1000;                                    ///< [masl]  (zref in HBV)
+const double  HBV_PET_ELEV_CORR       =0.0005;                                  ///< [mm/m-d] PET lapse rate (ECALT in HBV-EC)
+const double  HBV_PET_TEMP_CORR       =0.5;                                     ///< [-]     (ETF in HBV-EC)
+const double  GLOBAL_PET_CORR         =1.0;                                     ///< [-]     (ECORR in HBV-EC)
 
 //Flag variables
 const int     DOESNT_EXIST            =-1;                                      ///< return value for nonexistent index

src/Reservoir.cpp

@@ -1390,7 +1390,12 @@ void CReservoir::UpdateReservoir(const time_struct &tt, const optStruct &Options
       double obs_stage=_pObsStage->GetSampledValue(nn);
       if(obs_stage!=RAV_BLANK_DATA) {
         _stage=obs_stage;
-        _Qout =GetWeirOutflow(_stage,weir_adj);//[m3/s]
+        if ((Options.management_optimization) && (_Qoptimized != RAV_BLANK_DATA)) {
+          _Qout=_Qoptimized; //overriding stage AND outflow
+        }
+        else{
+          _Qout =GetWeirOutflow(_stage,weir_adj);//[m3/s] //default - outflow function of assimilated stage
+        }
         _assim_blank=false;
       }
     }
@@ -1542,7 +1547,12 @@ double  CReservoir::RouteWater(const double &Qin_old,
 {
   if ((_assimilate_stage) && (!_assim_blank))
   {
-    res_outflow=_Qout;//stage from data, outflow calculated from SD curve in UpdateReservoir()
+    if ((Options.management_optimization) && (_Qoptimized != RAV_BLANK_DATA)) {
+      res_outflow=_Qoptimized;
+    }
+    else{
+      res_outflow=_Qout;//stage from data, outflow calculated from SD curve in UpdateReservoir()
+    }
     return _stage;
   }