v514 - Isotope transport up and running!

Isotopes:
-Updated arguments needed for GetAdvectionCorrection() (Advection.cpp; ConstituentModel.cpp;Transport.cpp/h)
-blanks for zero volume concentrations (ConstituentModel.cpp)
-New routine UpdateInitialConditions(), member _initSWConc (IsotopeTransport.cpp/h;ParseInitialConditionsFile.cpp;Transport.cpp/h)
-Proper handling of isotopic composition and drying out stores (IsotopeTransport.cpp)
-Handle isotopic specified fluxes (MAssRouting.cpp)
-New rvc command :InitialSurfaceWaterConcentration (ParseInitialConditionsFile.cpp)
-new overloaded Transport routine GetConcentration() (Transport.cpp/.h)

bug fix: reading snow composition (ParseInput.cpp)

cleaned up IsDirichlet() routine (ConstituentModel.cpp)

Author: James Craig

src/Advection.cpp

@@ -187,25 +187,28 @@ void   CmvAdvection::GetRatesOfChange(const double      *state_vars,
       vol =sv[iToWater];
       iF=iToWater;iT=iFromWater;
     }
-    if (vol>1e-6)//note: otherwise Q should generally be constrained to be <vol/tstep & 0.0<rates[q]<(m/tstep*corr)
+    if ((vol>1e-6) && (Q[q]!=0))//note: otherwise Q should generally be constrained to be <vol/tstep & 0.0<rates[q]<(m/tstep*corr)
     {
-      corr=pTransModel->GetAdvectionCorrection(_constit_ind,pHRU,iF,iT,mass/vol*MM_PER_METER/LITER_PER_M3);//mg/m2/mm->mg/L
+      corr=pTransModel->GetAdvectionCorrection(_constit_ind,pHRU,iF,iT,mass,vol,Q[q]);
 
       rates[q]=Q[q]*mass/vol*corr; //[mg/m2/d] or [MJ/m2/d]
       if(fabs(rates[q])>mass/tstep) { rates[q]=(Q[q]/fabs(Q[q]))*mass/tstep; }//emptying out compartment
 
-      if ((!isEnthalpy) && (!isIsotope)){ //negative enthalpy/temperature/isotopic composition allowed
+      if ((!isEnthalpy)){ //negative enthalpy/temperature allowed
         if(mass<-1e-9) { ExitGracefully("CmvAdvection - negative mass",RUNTIME_ERR); }
       }
     }
+    else {
+      rates[q]=0.0;
+    }
 
     //special consideration - atmospheric precip can have negative storage but still specified concentration or temperature
     //----------------------------------------------------------------------
     if((pModel->GetStateVarType(iF)==ATMOS_PRECIP) &&
       (pConstit->IsDirichlet(iF,k,tt,Cs,1.0-pHRU->GetForcingFunctions()->snow_frac)))
     {
       if(!isEnthalpy) {
-        Cs=pConstit->ConvertConcentration(Cs);
+        Cs=pConstit->ConvertConcentration(Cs);//[mg/L]->[mg/mm-m2] or [o/oo]->[mg/mm-m2]
       }
       else {
         Cs=pTransModel->GetEnthalpyModel()->GetDirichletEnthalpy(pHRU,Cs);
@@ -227,7 +230,7 @@ void   CmvAdvection::GetRatesOfChange(const double      *state_vars,
     if(pConstit->IsDirichlet(iFromWater,k,tt,Cs))
     {
       if(!isEnthalpy) {
-        Cs=pConstit->ConvertConcentration(Cs);//[mg/L]->[mg/mm-m2]
+        Cs=pConstit->ConvertConcentration(Cs);//[mg/L]->[mg/mm-m2] or [o/oo]->[mg/mm-m2]
       }
       else{
         Cs=pTransModel->GetEnthalpyModel()->GetDirichletEnthalpy(pHRU,Cs);
@@ -244,7 +247,7 @@ void   CmvAdvection::GetRatesOfChange(const double      *state_vars,
     if(pConstit->IsDirichlet(iToWater,k,tt,Cs))
     {
       if(!isEnthalpy) {
-        Cs=pConstit->ConvertConcentration(Cs);
+        Cs=pConstit->ConvertConcentration(Cs);//[mg/L]->[mg/mm-m2] or [o/oo]->[mg/mm-m2]
       }
       else{
         Cs=pTransModel->GetEnthalpyModel()->GetDirichletEnthalpy(pHRU,Cs);

src/ConstituentModel.cpp

@@ -99,7 +99,7 @@ bool  CConstituentModel::IsPassive() const {
 /// \param iToWater   [in] index of "to" water storage state variable
 /// \param Cs         [in] RAW constituent concentration (mg/L or MJ/L or L/L, not mg/L or C or o/oo)
 //
-double CConstituentModel::GetAdvectionCorrection(const CHydroUnit* pHRU,const int iFromWater,const int iToWater,const double& Cs) const
+double CConstituentModel::GetAdvectionCorrection(const CHydroUnit* pHRU,const int iFromWater,const int iToWater,const double& mass, const double &vol, const double &Q) const
 {
   /*sv_type fromType=_pModel->GetStateVarType(iFromWater);
   sv_type toType  =_pModel->GetStateVarType(iToWater);
@@ -171,17 +171,18 @@ double CConstituentModel::GetCatchmentTransitLosses(const int p) const
 bool  CConstituentModel::IsDirichlet(const int i_stor,const int k,const time_struct &tt,double &Cs,const double blend) const
 {
   Cs=0.0;
+  double Cs2=0.0;
   int i_source=_aSourceIndices[i_stor][k];
   if(i_source==DOESNT_EXIST)          { return false; }
   if(!_pSources[i_source]->dirichlet) { return false; }
-  Cs = _pSources[i_source]->concentration;
-  if (blend!=1.0){Cs=Cs*blend+(1-blend)*_pSources[i_source]->concentration2;}
 
-  if(Cs != DOESNT_EXIST) { return true; }
-  else {//time series
-    Cs = _pSources[i_source]->pTS->GetValue(tt.model_time);
-    return true;
-  }
+  Cs =_pSources[i_source]->concentration;
+  Cs2=_pSources[i_source]->concentration2;
+  if (blend!=1.0){Cs=(blend)*Cs+(1-blend)*Cs2; }
+
+  if(Cs == DOESNT_EXIST) {Cs = _pSources[i_source]->pTS->GetValue(tt.model_time);}//time series
+
+  return true;
 }
 //////////////////////////////////////////////////////////////////
 /// \brief returns specified mass flux for given constitutent and water storage unit at time tt
@@ -668,6 +669,13 @@ void GetUnitNames(constit_type type,bool writemass,string &kg,string &kgd,string
   }
 }
 //////////////////////////////////////////////////////////////////
+/// \brief update initial conditions
+/// \details Called prior to simulation (but after initialization) from CModel::Initialize()
+//
+void CConstituentModel::UpdateInitialConditions(const optStruct& Options) {
+  //default: does nothing
+}
+//////////////////////////////////////////////////////////////////
 /// \brief Write transport output file headers
 /// \details Called prior to simulation (but after initialization) from CModel::Initialize()
 /// \param &Options [in] Global model options information
@@ -1325,7 +1333,12 @@ void CConstituentModel::WriteMinorOutput(const optStruct &Options,const time_str
 
     if(_pTransModel->GetStorWaterIndex(ii)!=iCumPrecip)
     {
-      _OUTPUT<<","<<concentration;     //print column entry
+      if ((V < 1e-6) && (_type == ISOTOPE)) {
+        _OUTPUT<<",";
+      }
+      else {
+        _OUTPUT<<","<<concentration;     //print column entry
+      }
       currentMass+=M*(area*M2_PER_KM2); //[mg]  or [MJ]  //increment total mass in system
     }
     else {

src/DemandOptimization.cpp

@@ -412,7 +412,7 @@ void CDemandOptimizer::AddNonLinVariable(const string name, const string targetD
   nonlin_var *pNonLinVar=new nonlin_var(name,targetDV);
 
   if(!DynArrayAppend((void**&)(_pNonLinVars),(void*)(pNonLinVar),_nNonLinVars)) {
-    ExitGracefully("CConstituentModel::AddDirichletCompartment: adding NULL source",BAD_DATA);
+    ExitGracefully("CDemandOptimizer::AddNonLinVariable: adding NULL source",BAD_DATA);
   }
 }
 //////////////////////////////////////////////////////////////////

src/EnergyTransport.cpp

@@ -68,7 +68,7 @@ double CEnthalpyModel::CalculateReportingConcentration(const double &M,const dou
 //
 double CEnthalpyModel::ConvertConcentration(const double &T) const
 {
-  //specified in C, convert to MJ/m2
+  //specified in degrees C, convert to MJ/mm/m2
   double pctfroz=0.0;
   if (T<0.0){pctfroz=1.0;} //treats all 0 degree water as unfrozen
   return ConvertTemperatureToVolumetricEnthalpy(T,pctfroz)/MM_PER_METER; //[C]->[MJ/m3]*[M/mm]=[MJ/mm-m2]
@@ -546,12 +546,14 @@ void CEnthalpyModel::Initialize(const optStruct& Options)
     for(int p=0;p<_pModel->GetNumSubBasins();p++)
     {
       pBasin=_pModel->GetSubBasin(p);
-      const double *aQin=pBasin->GetInflowHistory();
+      
       for(int i=0; i<pBasin->GetNumSegments(); i++)
       {
         _aMout[p][i]=pBasin->GetOutflowArray()[_pModel->GetSubBasin(p)->GetNumSegments()-1] * SEC_PER_DAY * hv; //not really - requires outflow rate from all segments in general case. Don't have access to this. assumes nSegs=1
       }
       _aMout_last[p]=_aMout[p][0];
+
+      const double *aQin=pBasin->GetInflowHistory();
       for(int i=0; i<_nMinHist[p]; i++)
       {
         _aMinHist[p][i]=aQin[i]*SEC_PER_DAY*hv;