v521 [BENCHMARKED] - added :TimeOfConcentrationMethod support

new :TimeOfConcentrationMethod
 - parse support (ParseInput.cpp)
 - new enum added to Options (RavenInclude.h)
 - new subbasin _basin_length and _mean_slope members (SubBasin.h/.cpp)
 - new routine CSubBasin::CalculateTOC (SubBasin.h/.cpp)

bug fixes:
  added QA/QC for bad HRU index in reservoir parsing (ParseHRUFile.cpp)
  fixed destructor bug in CustomOutput.cpp

Author: James Craig

src/CustomOutput.cpp

@@ -62,7 +62,7 @@ CCustomOutput::CCustomOutput( const diagnostic    variable,
   _spaceAggStr    ="UNKNOWN";
   _filename       =filename_spec;
 
-  _nData        =1;
+  _nData        =0;
   _aData        =NULL;  // pointer to data storage for aggregation
   _nDataItems   =1;
 

src/ParseHRUFile.cpp

@@ -361,7 +361,7 @@ bool ParseHRUPropsFile(CModel *&pModel, const optStruct &Options, bool terrain_r
       }
       pRes=ReservoirParse(pp,s[1],pModel,HRUID,Options);
       pSB=pModel->GetSubBasinByID(pRes->GetSubbasinID());
-      if (HRUID!=DOESNT_EXIST){
+      if ((HRUID!=DOESNT_EXIST) && (pModel->GetHRUByID(HRUID)!=NULL)){
         pRes->SetHRU(pModel->GetHRUByID(HRUID));
       }
       if (pSB!=NULL){pSB->AddReservoir(pRes);}
@@ -1937,10 +1937,15 @@ CReservoir *ReservoirParse(CParser *p,string name,const CModel *pModel,long long
     ExitGracefullyIf(max_depth==DOESNT_EXIST,"CReservoir::Parse: :LakeDepth must be specified for lake-type reservoirs",BAD_DATA_WARN);
 
    if (HRUID != DOESNT_EXIST) {
-     double HRUarea =pModel->GetHRUByID(HRUID)->GetArea() * M2_PER_KM2;
-     if (fabs((HRUarea - lakearea) / lakearea) > 0.2) {
-       string warn="CReservoirParse: specified :LakeArea and corresponding HRU area (in .rvh file) do not seem to agree for reservoir in subbasin "+to_string(SBID);
-       WriteWarning(warn.c_str(), Options.noisy);
+     if (pModel->GetHRUByID(HRUID)!=NULL){
+       double HRUarea =pModel->GetHRUByID(HRUID)->GetArea() * M2_PER_KM2;
+       if (fabs((HRUarea - lakearea) / lakearea) > 0.2) {
+         string warn="CReservoirParse: specified :LakeArea and corresponding HRU area (in .rvh file) do not seem to agree for reservoir in subbasin "+to_string(SBID);
+         WriteWarning(warn.c_str(), Options.noisy);
+       }
+     }
+     else {
+       ExitGracefullyIf(max_depth==DOESNT_EXIST,"CReservoir::Parse: invalid HRU ID in :Reservoir command",BAD_DATA_WARN);
      }
     }
 

src/ParseInput.cpp

@@ -276,6 +276,7 @@ bool ParseMainInputFile (CModel     *&pModel,
   Options.interception_factor     =PRECIP_ICEPT_USER;
   Options.recharge                =RECHARGE_NONE;
   Options.snow_depletion          =SNOWCOV_NONE;
+  Options.TOC_method              =TOC_MCDERMOTT_PILGRIM;
   Options.direct_evap             =false;
   Options.keepUBCWMbugs           =false;
   Options.suppressCompetitiveET   =false;
@@ -486,6 +487,7 @@ bool ParseMainInputFile (CModel     *&pModel,
     else if  (!strcmp(s[0],":FEWSStateInfoFile"         )){code=110;}
     else if  (!strcmp(s[0],":FEWSParamInfoFile"         )){code=111;}
     else if  (!strcmp(s[0],":FEWSBasinStateInfoFile"    )){code=112;}
+    else if  (!strcmp(s[0],":TimeOfConcentrationMethod" )){code=113;}
 
     else if  (!strcmp(s[0],":WriteGroundwaterHeads"     )){code=510;}//GWMIGRATE -TO REMOVE
     else if  (!strcmp(s[0],":WriteGroundwaterFlows"     )){code=511;}//GWMIGRATE -TO REMOVE
@@ -1946,6 +1948,17 @@ bool ParseMainInputFile (CModel     *&pModel,
       Options.flowinfo_filename = CorrectForRelativePath(s[1], Options.rvi_filename);//with .nc extension!
       break;
     }
+    case(113):  //
+    {/*:TimeOfConcentrationMethod [method]*/
+      if (Options.noisy) { cout <<  " Time of concentration method" << endl; }
+      if (Len<2){ImproperFormatWarning(":TimeOfConcentrationMethod",p,Options.noisy); break;}
+      if      (!strcmp(s[1],"TOC_MCDERMOTT_PILGRIM")){Options.TOC_method=TOC_MCDERMOTT_PILGRIM;}
+      else if (!strcmp(s[1],"TOC_AIRPORT"          )){Options.TOC_method=TOC_AIRPORT;}
+      else if (!strcmp(s[1],"TOC_BRANSBY_WILLIAMS" )){Options.TOC_method=TOC_BRANSBY_WILLIAMS;}
+      else if (!strcmp(s[1],"TOC_WILLIAMS_1922"    )){Options.TOC_method=TOC_WILLIAMS_1922;}
+      else {ExitGracefully("ParseInput :TimeOfConcentrationMethod: Unrecognized method",BAD_DATA_WARN);}
+      break;
+    }
     case(160):  //--------------------------------------------
     {/*:rvh_Filename [filename.rvh]*/
       if(Options.noisy) { cout <<"rvh filename: "<<s[1]<<endl; }
@@ -2501,7 +2514,7 @@ bool ParseMainInputFile (CModel     *&pModel,
     }
     case(215):  //----------------------------------------------
     {/*Flush
-       :Flush RAVEN_DEFAULT [state_var from] [state_var to] {Optional percentage}*/
+      :Flush RAVEN_DEFAULT [state_var from] [state_var to] {Optional percentage}*/
       if (Options.noisy){cout <<"Flushing Process"<<endl;}
       double pct=1.0;
       if (Len<4){ImproperFormatWarning(":Flush",p,Options.noisy); break;}

src/RavenInclude.h

@@ -126,6 +126,7 @@ const double  CM3_PER_METER3          =1e6;
 const double  METER_PER_CM            =0.01;                                    ///< [cm] to [m]
 const double  M2_PER_KM2              =1e6;                                     ///< [km2] to [m2]
 const double  MM2_PER_M2              =1e6;                                     ///< [m2] to [mm2]
+const double  HECTARE_PER_KM2         =100;                                     ///< [km2] to [ha]
 const double  M_PER_KM                =1000;                                    ///< [km] to [m]
 const double  GRAMS_PER_KG            =1000;                                    ///< [kg] to [g]
 const double  MG_PER_KG               =1000000;                                 ///< [kg] to [mg]
@@ -836,6 +837,15 @@ enum overflowmode
   OVERFLOW_NATURAL  ///< uses stage discharge curve to calculate Q
 };
 
+enum toc_method
+{
+  TOC_MCDERMOTT_PILGRIM,///< uses Austrailian Rainfall and runoff guidelines [McDermott and Pilgrim (1982)] [DEFAULT]
+  TOC_WILLIAMS_1922,    ///< uses Williams (1922) [DEFAULT]
+  TOC_AIRPORT,          ///< uses Airport equation (MTO Drainage Management Manual, 1997)
+  TOC_BRANSBY_WILLIAMS, ///< uses Bransby-Williams equation (MTO Drainage Management Manual, 1997)
+  TOC_TEMEZ             ///< uses Temez (1978)
+};
+
 enum assimtype
 {
   DA_RAVEN_DEFAULT, ///< multiplicative scaling assimilation upstream propagation
@@ -1111,6 +1121,7 @@ struct optStruct
   recharge_method    recharge;                ///< aquifer/soil recharge method
   bool               direct_evap;             ///< true if PET is used to directly reduce precipitation
   snowcov_method     snow_depletion;          ///< method for handling snowcover depletion curve
+  toc_method         TOC_method;              ///< method for calculating time of concentration
 
   precip_icept_method interception_factor;    ///< method for calculating canopy interception factor
 

src/SubBasin.cpp

@@ -48,6 +48,7 @@ CSubBasin::CSubBasin( const long long      Identifier,
   _drainage_area     =0.0;
   _reach_length      =reach_len;
   _reach_length2     =reach_len;
+  _basin_length      =AUTO_COMPUTE;
   _is_headwater      =true;
   _reach_HRUindex    =DOESNT_EXIST; //default
   _hyporheic_flux    =0.0; //default
@@ -1596,7 +1597,37 @@ double    CSubBasin::CalculateBasinArea()
 
   return _basin_area;
 }
+//////////////////////////////////////////////////////////////////
+/// \brief Calculates time of concentration, in days
+//
+double CSubBasin::CalculateTOC(const toc_method method)
+{
+  double slope_pct=tan(_mean_slope)*100;
+  double A=max(_basin_area,0.001); //to handle zero area issues
 
+  if      (method==TOC_MCDERMOTT_PILGRIM) // Austrailian Rainfall and runoff guidelines [McDermott and Pilgrim (1982)]
+  {
+    return 0.76*pow(A,0.38)/HR_PER_DAY; //[d]
+  }
+  else if (method==TOC_AIRPORT) // MTO Drainage Manual 1997
+  {
+    double C=0.3; //TMP DEBUG - cant have a single rational coeff for entire basin
+    return 3.26*(1.1-C)*sqrt(_basin_length*M_PER_KM)*pow(slope_pct,-0.33)/MIN_PER_DAY;
+  }
+  else if (method == TOC_BRANSBY_WILLIAMS) //MTO Drainage Manual 1997 
+  {
+    return 0.057*(_basin_length*M_PER_KM)*pow(slope_pct,-0.2)*pow(A*HECTARE_PER_KM2,-0.1);
+  }
+  else if (method == TOC_WILLIAMS_1922) 
+  {
+    return 0.2426*_reach_length*pow(A,-0.1)*pow(slope_pct/100,-0.2)/HR_PER_DAY;
+  }
+  else if (method == TOC_TEMEZ) 
+  {
+    return 0.3*pow(_reach_length,0.76)*pow(slope_pct/100,-0.19)/HR_PER_DAY;
+  }
+  return 0.0;
+}
 //////////////////////////////////////////////////////////////////
 /// \brief Initializes SB attributes
 /// \details Calculates total subbasin area; initializes _basin_area, _drainage_area,
@@ -1652,6 +1683,20 @@ void CSubBasin::Initialize(const double    &Qin_avg,          //[m3/s] from upst
     if(_pChannel!=NULL) {
       _pChannel->CheckReferenceFlow(_Q_ref,_slope,_mannings_n,_ID);
     }
+    
+    //Calculate mean slope from HRUs
+    //------------------------------------------------------------------------
+    double area;
+    for (int k=0;k<_nHydroUnits;k++){
+      area=_pHydroUnits[k]->GetArea();
+      _mean_slope = (_pHydroUnits[k]->GetSlope()) * area / _basin_area; //[radians]
+    }
+    
+    //Estimate basin length from area if not provided
+    //------------------------------------------------------------------------
+    if (_basin_length == AUTO_COMPUTE) {
+      _basin_length = 0.5*sqrt(_basin_area); //[km] assume square watershed
+    }
 
     //Estimate reach length from area if not provided
     //------------------------------------------------------------------------
@@ -1665,7 +1710,8 @@ void CSubBasin::Initialize(const double    &Qin_avg,          //[m3/s] from upst
     ///< \ref from Williams (1922), as cited in Handbook of Hydrology, eqn. 9.4.3 \cite williams1922
     if (_t_conc==AUTO_COMPUTE)
     {
-      _t_conc=0.76*pow(max(_basin_area,0.001),0.38)/HR_PER_DAY;// [d] \ref Austrailian Rainfall and runoff guidelines [McDermott and Pilgrim (1982)]
+      _t_conc=CalculateTOC(Options.TOC_method);
+
       _t_conc *= _pModel->GetGlobalParams()->GetParams()->TOC_multiplier;
     }
     if(Options.catchment_routing==ROUTE_GAMMA_CONVOLUTION ){_t_peak=AUTO_COMPUTE;}

src/SubBasin.h

@@ -60,6 +60,9 @@ class CSubBasin
   bool             _assimilate;   ///< true if observed data in this basin should be assimilated.
   const CSubBasin    *_pDownSB;   ///< pointer to downstream subbasin instance
 
+  double           _mean_slope;   ///< mean slope of basin [rad] 
+  double         _basin_length;   ///< watershed length (travel distance to reach) [km] 
+
   //catchment routing properties
   double               _t_conc;   ///< basin time of concentration [d]
   double               _t_peak;   ///< basin time to peak [d] (<=_t_conc)
@@ -161,6 +164,8 @@ class CSubBasin
   double                    TVDTheta(double In_old,double In_new,double Out_old,double Out_new,double th_in,double dx,double tstep) const;
 
   void            UpdateRoutingHydro(const double &tstep);
+
+  double                CalculateTOC(const toc_method method);
 public:/*-------------------------------------------------------*/
   //Constructors:
   CSubBasin(const long long      ID,