v508 [BENCHMARKED] - (assim + management); :RevertToSDCurve

support for flow and stage assimilation working with water management (DemandOptimization.cpp/.h)
-replaces mass balance constraints with assimilation constraints
-moved to +/- delta h representation as primary DVs in LP problem
-new routine CModel::GetObservedFlow
-CREservoir::GetObsStage()

support for RevertToSDCurve (DemandExpressionHandling.cpp; ParseManagementFile.cpp)

reservoir stage asssimilation now REQUIRES :AssimilateReservoirStage command in .rvt file [BACKWARD COMPATIBILTY]
-changes to Assimilate.cpp/ParseTimeSeriesFile.cpp/Reservoir.cpp

adjustments to CModel::UpdateDiagnostics() so that invalid data streams just get ignored (Model.cpp;ModelInitialize.cpp)
(allows you to comment out reservoirs with data without breaking model)

minor bug fixes -
 - use of dry stage rather than min stage for setting stage bounds in man model
- forece positive storage in SOILEVAP_HBV (SoilEvaporation.cpp)
- mistake in writing SubbasinProperties_Auto (StandardOutput.cpp)

Author: James Craig

src/Assimilate.cpp

@@ -57,24 +57,6 @@ void CModel::InitializeDataAssimilation(const optStruct &Options)
     }
   }
 
-  //Connect Lake assimilation observations
-  if(Options.assimilate_stage)
-  {
-    for(int i=0;i<_nObservedTS;i++) {
-      if((!strcmp(_pObservedTS[i]->GetName().c_str(),"RESERVOIR_STAGE")) && (_pObservedTS[i]->GetType() == CTimeSeriesABC::TS_REGULAR)) {
-        CSubBasin *pBasin=GetSubBasinByID(_pObservedTS[i]->GetLocID());
-        CReservoir *pRes=pBasin->GetReservoir();
-        if(pRes!=NULL) {
-          pRes->TurnOnAssimilation(_pObservedTS[i]);
-        }
-        else {
-          ExitGracefully("Reservoir stage observations assigned to basin without lake or reservoir. Data cannot be assimilated",BAD_DATA_WARN);
-        }
-      }
-    }
-    WriteAdvisory("Lake stage data assimilation will lead to mass balance error estimates in the WatershedStorage output file. This is a natural side effect of assimilation.",Options.noisy);
-  }
-
 }
 /////////////////////////////////////////////////////////////////
 /// \brief Overrides flows with observations at gauges and propagates adjustments upstream of gauges 
@@ -99,9 +81,7 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim
     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 (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;

src/DemandExpressionHandling.cpp

@@ -288,6 +288,7 @@ string DVTypeToString(dv_type t)
   if (t==DV_QOUTRES ){return "DV_QOUTRES";  }
   if (t==DV_STAGE   ){return "DV_STAGE";    }
   if (t==DV_DSTAGE  ){return "DV_DSTAGE";   }
+  if (t==DV_DSTAGE2 ){return "DV_DSTAGE2";   }
   if (t==DV_BINRES  ){return "DV_BINRES";   }
   if (t==DV_DELIVERY){return "DV_DELIVERY"; }
   if (t==DV_RETURN  ){return "DV_RETURN";   }
@@ -1319,6 +1320,17 @@ void CDemandOptimizer::AddConstraintToLP(const int ii, const int kk, lp_lib::lpr
   {
     pE= pC->pOperRegimes[kk]->pExpression; //active expression
 
+    if (pE->nGroups == 0) {
+      // special case - revert to SD curve - blank row
+      row_val[0]=0.0;
+      col_ind[0]=0;
+      RHS=0.0;
+      retval = lp_lib::set_rowex(pLinProg,lpgoalrow,1,row_val,col_ind);
+      ExitGracefullyIf(retval==0,"AddConstraintToLP::Error updating user-specified constraint/goal",RUNTIME_ERR);
+      retval = lp_lib::set_rh(pLinProg,lpgoalrow,RHS);
+      return; 
+    }
+
     RHS=0.0;
     for (int j = 0; j < pE->nGroups; j++)
     {