v515 - Lateral flush/diversion with several recipient HRUs

Updated CmvLatFlush::Initialize() to handle more than one recipient HRUS (LatFlush.cpp)

cleaned up CustomOutput.cpp/.h to conform to code specs for member vars

cosmetic changes to other files

Author: James Craig

src/CustomOutput.cpp

@@ -90,15 +90,15 @@ class CCustomOutput
   string       _statStr;    ///< statistic type string
   string       _spaceAggStr;///< spatial aggregation type string
 
-  double     **data;        ///< stores accumulated data for each HRU,Basin, or WShed (size:[num_store][num_data])
-  int          num_data;    ///< number of data points
-  int          num_store;   ///< number of data items needed for each HRU, Basin or WShed
+  double     **_aData;      ///< stores accumulated data for each HRU,Basin, or WShed (size:[_nDataItems][_nData])
+  int          _nData;      ///< number of data points
+  int          _nDataItems; ///< number of data items needed for each HRU, Basin or WShed
                             //(e.g., =2 if max and min are both tracked)
-  int         _time_index;  ///< index tracking current output line (e.g., 3=3 years/months/days passed, dependent upon _timeAgg
+  int          _time_index; ///< index tracking current output line (e.g., 3=3 years/months/days passed, dependent upon _timeAgg
 
   int          _count;      ///< counts accumulated data (# of timesteps since last output dump)
 
-  int          kk_only;     ///< index of HRUGroup for which output is generated when spaceAgg==BY_SELECT_HRUS
+  int          _kk_only;     ///< index of HRUGroup for which output is generated when spaceAgg==BY_SELECT_HRUS
 
   const CModel *pModel;     ///< Reference to model
 

src/CustomOutput.h

@@ -2254,7 +2254,7 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
       else if (retval==ACCURACYERROR){code="ACCURACYERROR "; }
 
       cout<<"=========================================="<<endl;
-      cout<<"LP SOLVE CANNOT SOLVE OPTIMZIATION PROBLEM"<<endl;
+      cout<<"LP SOLVE CANNOT SOLVE OPTIMIZATION PROBLEM"<<endl;
       cout<<"=========================================="<<endl;
       cout<<"Diagnostics:"<<endl;
       cout<<"  lp_lib::solve error code: "<<retval<<" ("<<code<<")"<<endl;

src/DemandOptimization.cpp

@@ -1454,6 +1454,7 @@ CDiagPeriod::CDiagPeriod(string name,string startdate,string enddate,comparison
 
   tt=DateStringToTimeStruct(enddate  ,"00:00:00",Options.calendar);
   _t_end   = TimeDifference(Options.julian_start_day,Options.julian_start_year,tt.julian_day,tt.year,Options.calendar);
+
   if (_t_end<=_t_start){
     string warn;
     warn="CDiagPeriod: :EvaluationPeriod "+_name+": startdate after enddate. ("+to_string(_t_start)+">="+to_string(_t_end);

src/Diagnostics.cpp

@@ -58,8 +58,7 @@ class CForcingGrid //: public CForcingGridABC
 
   double     **_aVal;                        ///< Array of magnitudes of pulses (variable units)
   ///                                        ///< [size _ChunkSize, _nNonZeroWeightedGridCells]
-  ///                                        ///< time steps are in model resolution (means original input data are
-  ///                                        ///< already aggregated to match model resolution)
+  ///                                        ///< time steps are in model resolution (means original input data are already aggregated to match model resolution)
 
   double     **_GridWeight;                  ///< Sparse array of weights for each HRU for a list of cells
   //                                         ///< Dimensions : [_nHydroUnits][_nWeights[k]] (variable)

src/ForcingGrid.h

@@ -62,7 +62,7 @@ void CmvLatFlush::Initialize()
   if(_constrain_to_SBs)
   {
     //sift through all HRUs
-    for(int p=0;p<_pModel->GetNumSubBasins();p++)
+    /*for (int p = 0; p<_pModel->GetNumSubBasins(); p++)
     {
       //find 'to' HRU (only one allowed per SB)
       int kToSB=DOESNT_EXIST;
@@ -88,12 +88,45 @@ void CmvLatFlush::Initialize()
           q++;
         }
       }
+    }*/
+    //Alternate approach: multiple recipient stores
+    bool source_sink_issue=false;
+    int k1,k2;
+    for(int p=0;p<_pModel->GetNumSubBasins();p++)
+    {
+      for(int ks=0; ks<_pModel->GetSubBasin(p)->GetNumHRUs(); ks++)
+      {
+        k1=_pModel->GetSubBasin(p)->GetHRU(ks)->GetGlobalIndex();
+        for(int ks2=0; ks2<_pModel->GetSubBasin(p)->GetNumHRUs(); ks2++)
+        {
+          k2=_pModel->GetSubBasin(p)->GetHRU(ks2)->GetGlobalIndex();
+
+          if ((_pModel->IsInHRUGroup(k1,toHRUGrp)) && (_pModel->IsInHRUGroup(k2,fromHRUGrp)))
+          {
+            if (k1!=k2){
+              kTo  [q]=k1;
+              kFrom[q]=k2;
+              q++;
+              ExitGracefullyIf(q>=MAX_LAT_CONNECTIONS,"Number of HRU connections in LateralFlush or LateralDivert exceeds MAX_LAT_CONNECTIONS limit.",BAD_DATA);
+            }
+            else {
+              //may want to throw warning if source is also sink
+              source_sink_issue=true;
+            }
+          }
+        }
+      }
     }
+
     nConn=q;
     if(nConn==0) {
       string warning="CmvLatFlush::Initialize: no connections found between from and to HRU groups in any of the basins. if INTERBASIN flag not used, only transfer within basins is supported. ";
       WriteWarning(warning,true);
     }
+    if (source_sink_issue) {
+      string warning="CmvLatFlush::Initialize: one or more HRUs belongs to both source and recipient HRU groups in LateralFlush or LateralDivert process. This may have unintentended consequences. ";
+      WriteWarning(warning,true);
+    }
   }
   else //!constrain_to_SBs
   {
@@ -102,7 +135,7 @@ void CmvLatFlush::Initialize()
     for(k=0;k<_pModel->GetNumHRUs();k++) {
       if(_pModel->IsInHRUGroup(k,toHRUGrp)) {
         ExitGracefullyIf(kToSB!=DOESNT_EXIST,
-          "LatFlush::Initialize - only one HRU in the model can recieve flush output. More than one recipient HRU found.",BAD_DATA_WARN);
+          "LatFlush::Initialize - only one HRU in the model can recieve flush output if INTERBASIN is used. More than one recipient HRU found.",BAD_DATA_WARN);
         kToSB=k;
       }
     }
@@ -128,9 +161,6 @@ void CmvLatFlush::Initialize()
     _kTo     [q]    =kTo[q];
     _iFromLat[q]    =_iFlushFrom;
     _iToLat  [q]    =_iFlushTo;
-    /*cout << "latflush connections: " << _pModel->GetHydroUnit(_kFrom[q])->GetID() << " " << _pModel->GetHydroUnit(_kTo[q])->GetID() << " ";
-    cout<<_pModel->GetStateVarInfo()->SVTypeToString(_pModel->GetStateVarType(_iFromLat[q]),-1)<<" ";
-    cout<<_pModel->GetStateVarInfo()->SVTypeToString(_pModel->GetStateVarType(_iToLat[q]),-1)<<endl;*/
   }
   delete [] kFrom;
   delete [] kTo;

src/LatFlush.cpp

@@ -38,7 +38,7 @@ namespace MFUSG
 		// Routines to work with MFUSG data
 		void   mask_iunit            (int *export_iunit); // See Extern.f90 in MFUSGLib for iunit/package table (CUNIT)
 		void   add_to_flow_eq        (int *node, double *AMAT_change, double *RHS_change);
-	    void   add_to_flow_budget    (int *nnodes, int *rvn_nodes, double *rates, char *pckg_name);
+        void   add_to_flow_budget    (int *nnodes, int *rvn_nodes, double *rates, char *pckg_name);
         void   first_active_below    (int *node);
 		void   set_irvncb            (int *unit_no);
 		void   set_nrchop            (int *new_NRCHOP);
@@ -58,4 +58,4 @@ namespace MFUSG
 
 }
 #endif
-#endif
+#endif

src/MFUSGpp.h

@@ -249,84 +249,70 @@ CModel::~CModel()
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns pointer to global parameters object
-///
 /// \return Pointer to global parameters object
 //
-CGlobalParams* CModel::GetGlobalParams() const {
-  return _pGlobalParams;
-}
+CGlobalParams* CModel::GetGlobalParams() const { return _pGlobalParams;}
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns number of sub basins in model
-///
 /// \return Integer number of sub basins
 //
 int CModel::GetNumSubBasins   () const{return _nSubBasins;}
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns number of SB groups
-///
 /// \return Integer number of SB groups
 //
 int CModel::GetNumSubBasinGroups() const { return _nSBGroups; }
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns number of HRUs in model
-///
 /// \return Integer number of HRUs
 //
 int CModel::GetNumHRUs        () const{return _nHydroUnits;}
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns number of HRU groups
-///
 /// \return Integer number of HRU groups
 //
 int CModel::GetNumHRUGroups   () const{return _nHRUGroups;}
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns number of gauges in model
-///
 /// \return Integer number of gauges in model
 //
 int CModel::GetNumGauges      () const{return _nGauges;}
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns number of gridded forcings in model
-///
 /// \return Integer number of gridded forcings in model
 //
 int CModel::GetNumForcingGrids () const{return _nForcingGrids;}
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns number of state variables per HRU in model
-///
 /// \return Integer number of state variables per HRU in model
 //
 int CModel::GetNumStateVars   () const{return _nStateVars;}
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns number of soil layers
-///
 /// \return Integer number of soil layers
 //
 int CModel::GetNumSoilLayers  () const{return _nSoilVars;}
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns number of hydrologic processes simulated by model
-///
 /// \return Integer number of hydrologic processes simulated by model
 //
 int CModel::GetNumProcesses   () const{return _nProcesses;}
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns total modeled watershed area
-///
 /// \return total modeled watershed area [km2]
 //
 double CModel::GetWatershedArea () const{return _WatershedArea;}
 
-
 //////////////////////////////////////////////////////////////////
 /// \brief Returns number of observation time series
 /// \return number of observation time series
@@ -371,13 +357,11 @@ const CTimeSeriesABC* CModel::GetSimulatedTS(const int i) const {
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns specific hydrologic process denoted by parameter
-///
 /// \param j [in] Process index
 /// \return pointer to hydrologic process corresponding to passed index j
 //
 CHydroProcessABC *CModel::GetProcess(const int j) const
 {
-
 #ifdef _STRICTCHECK_
   ExitGracefullyIf((j<0) || (j>=_nProcesses),"CModel GetProcess::improper index",BAD_DATA);
 #endif
@@ -386,7 +370,6 @@ CHydroProcessABC *CModel::GetProcess(const int j) const
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns specific gauge denoted by index
-///
 /// \param g [in] Gauge index
 /// \return pointer to gauge corresponding to passed index g
 //
@@ -400,7 +383,6 @@ CGauge *CModel::GetGauge(const int g) const
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns specific forcing grid denoted by index
-///
 /// \param f [in] Forcing Grid index
 /// \return pointer to gauge corresponding to passed index g
 //
@@ -417,7 +399,6 @@ CForcingGrid *CModel::GetForcingGrid(const forcing_type &ftyp) const
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns specific HRU denoted by index k
-///
 /// \param k [in] HRU index
 /// \return pointer to HRU corresponding to passed index k
 //
@@ -431,7 +412,6 @@ CHydroUnit *CModel::GetHydroUnit(const int k) const
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns specific HRU with HRU identifier HRUID
-///
 /// \param HRUID [in] HRU identifier
 /// \return pointer to HRU corresponding to passed ID HRUID, NULL if no such HRU exists
 //
@@ -449,7 +429,6 @@ CHydroUnit *CModel::GetHRUByID(const long long int HRUID) const
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns specific HRU group denoted by parameter kk
-///
 /// \param kk [in] HRU group index
 /// \return pointer to HRU group corresponding to passed index kk
 //
@@ -463,7 +442,6 @@ CHRUGroup  *CModel::GetHRUGroup(const int kk) const
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns specific HRU group denoted by string parameter
-///
 /// \param name [in] String name of HRU group
 /// \return pointer to HRU group corresponding to passed name, or NULL if this group doesn't exist
 //
@@ -513,7 +491,6 @@ CSubBasin  *CModel::GetSubBasin(const int p) const
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns index of subbasin downstream from subbasin referred to by index
-///
 /// \param p [in] List index for accessing subbasin
 /// \return downstream subbasin index, if input index is valid; -1 if there is no downstream basin
 //
@@ -527,7 +504,6 @@ int         CModel::GetDownstreamBasin(const int p) const
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns subbasin object corresponding to passed subbasin ID
-///
 /// \param SBID [in] long long integer sub basin ID
 /// \return pointer to Sub basin object corresponding to passed ID, if ID is valid
 //
@@ -546,7 +522,6 @@ CSubBasin  *CModel::GetSubBasinByID(const long long SBID) const
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns sub basin index corresponding to passed subbasin ID
-///
 /// \param SBID [in] Integer subbasin ID
 /// \return SubBasin index corresponding to passed ID, if ID is valid
 //
@@ -564,7 +539,6 @@ int         CModel::GetSubBasinIndex(const long long SBID) const
 }
 //////////////////////////////////////////////////////////////////
 /// \brief Returns array of pointers to subbasins upstream of subbasin SBID, including that subbasin
-///
 /// \param SBID [in] long long int subbasin ID
 /// \param nUpstream [out] size of array of pointers of subbasins
 /// \return array of pointers to subbasins upstream of subbasin SBID, including that subbasin
@@ -774,7 +748,6 @@ int         CModel::GetStateVarLayer(const int ii) const
 
 //////////////////////////////////////////////////////////////////
 /// \brief Checks if state variable passed exists in model
-///
 /// \param typ [in] State variable type
 /// \return Boolean indicating whether state variable exists in model
 //
@@ -785,14 +758,12 @@ bool        CModel::StateVarExists(sv_type typ) const
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns lake storage variable index
-///
 /// \return Integer index of lake storage variable
 //
 int         CModel::GetLakeStorageIndex() const{return _lake_sv;}
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns gauge index of gauge with specified name
-///
 /// \return Integer index of gauge
 /// \param name [in] specified name
 //
@@ -806,7 +777,6 @@ int  CModel::GetGaugeIndexFromName (const string name) const
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns forcing grid index of forcing grid with specified type
-///
 /// \return Integer index of forcing grid
 /// \param name [in] specified type
 //
@@ -821,7 +791,6 @@ int  CModel::GetForcingGridIndexFromType (const forcing_type &typ) const
 //////////////////////////////////////////////////////////////////
 /// \brief Returns current mass/energy flux (mm/d, MJ/m2/d, mg/m2/d) between two storage compartments iFrom and iTo
 /// \details required for advective transport processes
-///
 /// \param k [in] HRU index
 /// \param js [in] index of process connection (i.e., j*)
 /// \param &Options [in] Global model options information
@@ -836,7 +805,6 @@ double CModel::GetFlux(const int k, const int js, const optStruct &Options) cons
 }
 //////////////////////////////////////////////////////////////////
 /// \brief returns concentration or temperature within hru k with storage index i
-///
 /// \param k [in] HRU index
 /// \param i [in] mass/energy state variable index
 //
@@ -848,7 +816,6 @@ double CModel::GetConcentration(const int k, const int i) const
 //////////////////////////////////////////////////////////////////
 /// \brief Returns current mass/energy flow (mm-m2/d, MJ/d, mg/d) between two storage compartments iFrom and iTo
 /// \details required for advective transport processes
-///
 /// \param k [in] HRU index
 /// \param qs [in] global index of process connection (i.e., q*)
 /// \param &Options [in] Global model options information
@@ -862,7 +829,6 @@ double CModel::GetLatFlow(const int qs,const optStruct &Options) const
 }
 //////////////////////////////////////////////////////////////////
 /// \brief Returns cumulative flux to/from storage unit i
-///
 /// \param k [in] index of HRU
 /// \param i [in] index of storage compartment
 /// \param to [in] true if evaluating cumulative flux to storage compartment, false for 'from'
@@ -903,7 +869,6 @@ double CModel::GetCumulativeFlux(const int k, const int i, const bool to) const
 }
 //////////////////////////////////////////////////////////////////
 /// \brief Returns cumulative gross flux between unit iFrom and iTo in HRU k
-///
 /// \param k [in] index of HRU
 /// \param iFrom [in] index of storage compartment
 /// \param iTo [in] index of storage compartment