v494 - fixed timelag in :DemandExpression; removed input ordering issue; POTMELT_DD_FREEZE;

added _pUserTimeSeries to remove ordering issue where demands have to be specified before user DVs (DemandOptimization.h/.cpp)
 - removed DemandsAreInitialized() routine and member - much cleaner (DemandOptimization.cpp/.h/ParseManagementFile.cpp)
 - support :DeclareWorkflowVar before :DefineWorkflowVar (ParseManagementFile.cpp)
 - added PrepDemandProblem() to Solvers.cpp

support buffer height of zero in Lake Freezing (FrozenLake.cpp)

new potential melt routing POTMELT_DD_FREEZE (ModelParamCheck.cpp;PotentialMelt.cpp;ParseInput.cpp;RavenInclude.h)

fixed :BedTemperature reading bug in ParseInitialConditionFile.cpp (found by M. Krassovski) (ParseInitialConditionFile.cpp)

added potential melt params to RavenAlgParams.dat

Author: James Craig

src/DemandOptimization.cpp

@@ -39,6 +39,7 @@ CDemandOptimizer::CDemandOptimizer(CModel *pMod)
   _nDecisionVars=0;
   _pDecisionVars=NULL;
   _nUserDecisionVars=0;
+  _pUserDecisionVars=NULL;
 
   _nSlackVars=0;
   _aSlackValues=NULL;
@@ -68,8 +69,6 @@ CDemandOptimizer::CDemandOptimizer(CModel *pMod)
   _nGoals=0;
   _pGoals=NULL;
 
-  _demands_initialized=false;
-
   _do_debug_level=0;//no debugging
 
   _nSolverResiduals=0;
@@ -86,6 +85,7 @@ CDemandOptimizer::~CDemandOptimizer()
     cout<<"DESTROYING DEMAND OPTIMIZER"<<endl;
   }
   for (int i=0;i<_nDecisionVars;    i++){delete _pDecisionVars[i];    }delete [] _pDecisionVars;
+  for (int i=0;i<_nUserDecisionVars;i++){                             }delete [] _pUserDecisionVars;
   for (int i=0;i<_nUserTimeSeries;  i++){delete _pUserTimeSeries[i];  }delete [] _pUserTimeSeries;
   for (int i=0;i<_nUserLookupTables;i++){delete _pUserLookupTables[i];}delete [] _pUserLookupTables;
   for (int i=0;i<_nWorkflowVars;    i++){delete _pWorkflowVars[i];    }delete [] _pWorkflowVars;
@@ -175,13 +175,6 @@ int CDemandOptimizer::GetNumDemandGroups() const
   return _nDemandGroups;
 }
 //////////////////////////////////////////////////////////////////
-/// \brief returns true if demands have been initialized
-//
-bool CDemandOptimizer::DemandsAreInitialized() const
-{
-  return _demands_initialized;
-}
-//////////////////////////////////////////////////////////////////
 /// \brief returns debug level
 //
 int CDemandOptimizer::GetDebugLevel() const
@@ -301,18 +294,29 @@ void  CDemandOptimizer::AddWaterDemand(CDemand* pDem)
 //
 void CDemandOptimizer::AddDecisionVar(const decision_var* pDV)
 {
-  if (VariableNameExists(pDV->name)) {
+  if (VariableNameExists(pDV->name) && (pDV->dvar_type!=DV_USER)) { //otherwise, captures self when user vars added to master list
     string warn="CDemandOptimizer::AddDecisionVar: decision variable name "+pDV->name+" is already in use.";
     ExitGracefully(warn.c_str(),BAD_DATA_WARN);
   }
 
   if (!DynArrayAppend((void**&)(_pDecisionVars),(void*)(pDV),_nDecisionVars)){
-   ExitGracefully("CDemandOptimizer::AddDecisionVar: adding NULL DV",BAD_DATA);}
+    ExitGracefully("CDemandOptimizer::AddDecisionVar: adding NULL DV",BAD_DATA);}
 
-  if (pDV->dvar_type==DV_USER){
-    _nUserDecisionVars++;
+}
+//////////////////////////////////////////////////////////////////
+/// \brief adds USER decision variable structure to _pUserDecisionVars member aarray (later added to _pDecisionVars array in RVM initialize)
+/// \params pDV [in] - user-specified decision variable to be added
+//
+void CDemandOptimizer::AddUserDecisionVar(const decision_var* pDV)
+{
+  if (VariableNameExists(pDV->name)) {
+    string warn="CDemandOptimizer::AddUserDecisionVar: decision variable name "+pDV->name+" is already in use.";
+    ExitGracefully(warn.c_str(),BAD_DATA_WARN);
   }
+  if (!DynArrayAppend((void**&)(_pUserDecisionVars),(void*)(pDV),_nUserDecisionVars)){
+    ExitGracefully("CDemandOptimizer::AddUserDecisionVar: adding NULL DV",BAD_DATA);}
 }
+
 //////////////////////////////////////////////////////////////////
 /// \brief disables stage discharge curve handling for reservoir in subbasin p
 //
@@ -323,16 +327,16 @@ void CDemandOptimizer::OverrideSDCurve(const int p)
 //////////////////////////////////////////////////////////////////
 /// \brief sets bounds for user specified decision variable
 //
-void CDemandOptimizer::SetDecisionVarBounds(const string name, const double& min, const double& max)
+void CDemandOptimizer::SetUserDecisionVarBounds(const string name, const double& min, const double& max)
 {
-  for (int i = 0; i < _nDecisionVars; i++) {
-    if (_pDecisionVars[i]->name==name){
-      _pDecisionVars[i]->min=min;
-      _pDecisionVars[i]->max=max;
+  for (int i = 0; i < _nUserDecisionVars; i++) {
+    if (_pUserDecisionVars[i]->name==name){
+      _pUserDecisionVars[i]->min=min;
+      _pUserDecisionVars[i]->max=max;
       return;
     }
   }
-  string warn = "SetDecisionVarBounds: invalid decision variable name (" + name + "): must define before use";
+  string warn = "SetUserDecisionVarBounds: invalid decision variable name (" + name + "): must declare or define before use";
   WriteWarning(warn.c_str(),_pModel->GetOptStruct()->noisy);
 }
 
@@ -414,6 +418,9 @@ bool CDemandOptimizer::VariableNameExists(const string &name) const
   for (int i=0; i<_nDecisionVars; i++){
     if (_pDecisionVars[i]->name==name){return true;}
   }
+  for (int i=0; i<_nUserDecisionVars; i++){//not yet in _pDecisionVars when called
+    if (_pUserDecisionVars[i]->name==name){return true;}
+  }
   for (int i = 0; i < _nUserConstants; i++) {
     if (_aUserConstNames[i]==name){return true;}
   }
@@ -464,6 +471,8 @@ void CDemandOptimizer::Initialize(CModel* pModel, const optStruct& Options)
   // Populate ordered decision variable array _pDecisionVars[]
   //  This order has to be maintained because of how the decision variables are indexed (consistent with GetDVColumnInd())
   //  subbasin outflows -> reservoir outflows -> reservoir stages -> delivered demand -> user-specified DVs -> slack variables
+  //  demand variables added in InitializeDemands
+  //  user and slack variables added in InitalizePostRVM
   //------------------------------------------------------------------
   // add subbasin outflow DVs
   int SB_count=0;
@@ -596,15 +605,12 @@ void  CDemandOptimizer::IdentifyUpstreamDemands()
 }
 //////////////////////////////////////////////////////////////////
 /// \brief Initializes Demand decision variables
-/// \notes to be called during .rvm file read, PRIOR to declaring any user-specified DVs, goals, or constraints
-/// WHY? because of the expected ordering of decision variables in the GetDVColumnInd() routine
+/// Called once at end of RVM read PRIOR to InitializePostRVM()
 /// \params pModel [in] - pointer to model
 /// \params Options  [in] - model options structure
 //
 void CDemandOptimizer::InitializeDemands(CModel* pModel, const optStruct& Options)
 {
-  if (_demands_initialized){return;}//This routine has already been called
-
   if (Options.noisy){cout<<"CDemandOptimizer: Demand initialization..."<<endl;}
 
   // reserve memory for delivery arrays
@@ -657,8 +663,6 @@ void CDemandOptimizer::InitializeDemands(CModel* pModel, const optStruct& Option
     }
   }
   _nReturns=r;
-
-  _demands_initialized=true;
 }
 string ComparisonToString(comparison C)
 {
@@ -708,6 +712,12 @@ void CDemandOptimizer::InitializePostRVMRead(CModel* pModel, const optStruct& Op
     }
   }
 
+  //Add User-defined decision vars to _pDecisionVars array
+  //------------------------------------------------------------------
+  for (int i = 0; i < _nUserDecisionVars; i++) {
+    AddDecisionVar(_pUserDecisionVars[i]);
+  }
+
   int p;
   _nSlackVars = 0;
   CSubBasin    *pSB;
@@ -730,6 +740,7 @@ void CDemandOptimizer::InitializePostRVMRead(CModel* pModel, const optStruct& Op
   }
 
   // Convert reservoir commands to user-specified constraints
+  //   allowed here because it doesnt introduce new DVs
   //------------------------------------------------------------------
   AddReservoirConstraints();
 
@@ -1243,14 +1254,31 @@ int CDemandOptimizer::GetDVColumnInd(const dv_type typ, const int counter) const
   return 0;
 }
 //////////////////////////////////////////////////////////////////
+/// \brief prepares demand optimization problem
+/// \notes to be called every time step prior to demand updates and routing mass balance from Solvers.cpp
+/// \params pModel [in] - pointer to model
+/// \params Options  [in] - model options structure
+/// \params tt [in] - model time structure
+//
+void CDemandOptimizer::PrepDemandProblem(CModel *pModel, const optStruct &Options, const time_struct &tt)
+{
+  // update history arrays from previous timestep
+  // ----------------------------------------------------------------
+  UpdateHistoryArrays();
+
+  // evaluates value of all workflow variables for this time step
+  // ----------------------------------------------------------------
+  UpdateWorkflowVariables(tt,Options);
+}
+//////////////////////////////////////////////////////////////////
 /// \brief Solves demand optimization problem
 /// \notes to be called every time step in lieu of routing mass balance
 /// \params pModel [in] - pointer to model
 /// \params Options  [in] - model options structure
 /// \params aSBrunoff [in] - array of current amount of runoff released to each subbasin (m3)  [size:nSubBasins]
-/// \params t [in] - local model time
+/// \params tt [in] - model time structure
 //
-void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Options, const double *aSBrunoff,const time_struct &tt )
+void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Options, const double *aSBrunoff,const time_struct &tt)
 {
 #ifdef _LPSOLVE_
 
@@ -1274,14 +1302,6 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio
   double *Q_iter =new double [_pModel->GetNumSubBasins()];
   int    *lprow  =new int    [_pModel->GetNumSubBasins()]; //index of goal equation for non-linear reservoir stage discharge curve in subbasin p
 
-  // update history arrays from previous timestep
-  // ----------------------------------------------------------------
-  UpdateHistoryArrays();
-
-  // evaluates value of all workflow variables for this time step
-  // ----------------------------------------------------------------
-  UpdateWorkflowVariables(tt,Options);
-
   // instantiate linear programming solver
   // ----------------------------------------------------------------
   lp_lib::lprec *pLinProg;

src/DemandOptimization.h

@@ -221,6 +221,7 @@ class CDemandOptimizer
   int              _nSlackVars;         //< number of slack variables
 
   int             _nUserDecisionVars;   //< number of user-specified decision variables
+  decision_var  **_pUserDecisionVars;   //< array of pointers to user-specified decision vars [size: _nUserDecisionVars] (overlaps content of _pDecisionVars)
 
   int             _nUserConstants;      //< number of user-specified named constants
   string         *_aUserConstNames;     //< array of names of user-specified constants
@@ -246,8 +247,6 @@ class CDemandOptimizer
   ofstream        _GOALSAT;            //< ofstream for GoalSatisfaction.csv
   ofstream        _RESID;              //< ofstream for OptimizationResidual.csv
 
-  bool            _demands_initialized;//< true if demands have been initialized
-
   int             _do_debug_level;      //< =1 if debug info is to be printed to screen, =2 if LP matrix also printed (full debug), 0 for nothing
 
   //Called during simualtion
@@ -269,6 +268,7 @@ class CDemandOptimizer
 
 
   //Called during initialization
+  void              AddDecisionVar(const decision_var *pDV);
   bool        UserTimeSeriesExists(string TSname) const;
   void     AddReservoirConstraints();
   void     IdentifyUpstreamDemands();
@@ -295,17 +295,16 @@ class CDemandOptimizer
   CDemand      *GetWaterDemand        (const int d);
   int           GetNumWaterDemands    () const;
 
-  bool   DemandsAreInitialized() const;
-
   void   SetHistoryLength      (const int n);
   void   SetCumulativeDate     (const int julian_date, const string demandID);
   void   SetDebugLevel         (const int lev);
   void   SetDemandAsUnrestricted(const string dname);
   void   OverrideSDCurve       (const int p);
 
   void   AddGoalOrConstraint   (const managementGoal *pGoal);
-  void   AddDecisionVar        (const decision_var *pDV);
-  void   SetDecisionVarBounds  (const string name, const double &min, const double &max);
+  
+  void   AddUserDecisionVar    (const decision_var *pDV);
+  void   SetUserDecisionVarBounds(const string name, const double &min, const double &max);
   void   AddUserConstant       (const string name, const double &val);
   void   AddWorkflowVariable   (const workflowVar *pCV);
   void   AddUserTimeSeries     (const CTimeSeries *pTS);
@@ -327,6 +326,7 @@ class CDemandOptimizer
   void   Initialize            (CModel *pModel, const optStruct &Options);
   void   InitializePostRVMRead (CModel *pModel, const optStruct &Options);
   void   InitializeDemands     (CModel *pModel, const optStruct &Options);
+  void   PrepDemandProblem     (CModel *pModel, const optStruct &Options, const time_struct &tt);
   void   SolveDemandProblem    (CModel *pModel, const optStruct &Options, const double *aSBrunoff, const time_struct &tt);
 
   void   WriteOutputFileHeaders(const optStruct &Options);

src/EnergyTransport.h

@@ -26,7 +26,8 @@ class CEnthalpyModel :public CConstituentModel
   int            _mHyporheic;  ///< model soil layer corresponding to hyporheic mixing layer (default=2)
   double       *_aMinResTime;  ///< minimum residence time [d] (<1 dt) in each stream reach [size: nSubBasins]
 
-  double          *_aBedTemp;  ///< array of riverbed temperatures (state variable) [C] [size: nSubBasins]
+  double          *_aBedTemp;  ///< array of riverbed temperatures (*state variable*) [C] [size: nSubBasins]
+
   double       *_aTave_reach;  ///< time-lagged average reach water temperature [C] [size: nSubBasins]
   double   *_aSS_temperature;  ///< array of steady state target temperatures in each basin [C] [size: nSubBasins]
 

src/FrozenLake.cpp

@@ -120,7 +120,7 @@ void   CmvFrozenLake::GetRatesOfChange( const double      *state_vars,
     double SWE=state_vars[iSWE];
 
     double corr=1.0;
-    if      (SWE>buff_ht ){corr=0.0;            }
+    if      (SWE>=buff_ht){corr=0.0;            }
     else if (SWE<=0.0    ){corr=1.0;            }
     else                  {corr=1.0-SWE/buff_ht;}
 

src/Model.cpp

@@ -1737,7 +1737,7 @@ void CModel::OverrideStreamflow   (const long long SBID)
         CTimeSeries *pTS =new CTimeSeries(name,*pObs);//copy time series
         pTS->SetLocID(SBID);
 
-        //need to shift everything by one interval, since hydrographs are stored as period-ending
+        //need to shift everything by one interval, since HYDROGRAPHS are stored as period-ending
         pTS->ShiftInTime(-(pTS->GetInterval()),*_pOptStruct);
 
         //add as inflow hydrograph to downstream

src/ModelParamCheck.cpp

@@ -107,6 +107,12 @@ void CModel::AddFromPotMeltParamList(string *aP,class_type *aPC,int &nP,const po
     aP[nP]="MELT_FACTOR"; aPC[nP]=CLASS_LANDUSE; nP++;
     aP[nP]="DD_MELT_TEMP";aPC[nP]=CLASS_LANDUSE; nP++;
   }
+  else if(pot_melt == POTMELT_DD_FREEZE)
+  {
+    aP[nP]="MELT_FACTOR"; aPC[nP]=CLASS_LANDUSE; nP++;
+    aP[nP]="DD_MELT_TEMP";aPC[nP]=CLASS_LANDUSE; nP++;
+    aP[nP]="REFREEZE_FACTOR";aPC[nP]=CLASS_LANDUSE; nP++;
+  }  
   else if(pot_melt==POTMELT_RESTRICTED)
   {
     aP[nP]="MELT_FACTOR"; aPC[nP]=CLASS_LANDUSE; nP++;

src/ParseInitialConditionFile.cpp

@@ -845,7 +845,7 @@ bool ParseInitialConditionsFile(CModel *&pModel, const optStruct &Options)
         }
         else if (!strcmp(s[0], ":BedTemperature"))
         {
-          if (Len >= 3) {
+          if (Len >= 2) {
             CEnthalpyModel *pEnth=(CEnthalpyModel*)(pConstit);
             pEnth->SetBedTemperature (p,s_to_d(s[1]));
           }

src/ParseInput.cpp

@@ -3733,6 +3733,7 @@ potmelt_method ParsePotMeltMethod(const string s)
 {
   string tmp=StringToUppercase(s);
   if      (!strcmp(tmp.c_str(),"POTMELT_DEGREE_DAY")){return POTMELT_DEGREE_DAY;}
+  else if (!strcmp(tmp.c_str(),"POTMELT_DD_FREEZE" )){return POTMELT_DD_FREEZE;}
   else if (!strcmp(tmp.c_str(),"POTMELT_EB"        )){return POTMELT_EB;}
   else if (!strcmp(tmp.c_str(),"POTMELT_RESTRICTED")){return POTMELT_RESTRICTED;}
   else if (!strcmp(tmp.c_str(),"POTMELT_DD_RAIN"   )){return POTMELT_DD_RAIN;}