v526 - Support for :Flush BY_SUBBASIN, minor assimilation change

Support for :Flush BY_SUBBASIN
-full model pointer in CmvFlush (HydroProcessABC.h;Flush.cpp)
-parsing support (ParseInput.cpp)
-new global constant BY_SUBBASIN_FLAG (RavenInclude.h)
-new SubBasin member _flush_fract (SubBasin.cpp/.h)

adjustment to AdjustAllFlows to properly correct inflows in long basin

minor edit to StandardOutput.cpp to support reset of output times in ensemble simulation

Author: James Craig

src/Flush.cpp

@@ -1,6 +1,6 @@
 /*----------------------------------------------------------------
   Raven Library Source Code
-  Copyright (c) 2008-2022 the Raven Development Team
+  Copyright (c) 2008-2025 the Raven Development Team
   ------------------------------------------------------------------
   Flush (abstract move of all water from one compartment to another)
   Split (move water to two compartments)
@@ -24,9 +24,11 @@
 CmvFlush::CmvFlush(int In_index,                       //soil water storage
                    int Out_index,
                    const double &pct,
-                   CModelABC *pModel)
+                   CModelABC *pModel,
+                   CModel *pModel_full)
   :CHydroProcessABC(FLUSH, In_index, Out_index, pModel)
 {
+  _pModel_full=pModel_full;
   _percentage=pct;
   ExitGracefullyIf(In_index==DOESNT_EXIST,
                    "CmvFlush Constructor: invalid 'from' compartment specified",BAD_DATA);
@@ -75,8 +77,13 @@ void   CmvFlush::GetRatesOfChange( const double      *storage,
                                    const optStruct   &Options,
                                    const time_struct &tt,
                                          double     *rates) const
-{
-  rates[0]=_percentage*storage[iFrom[0]]/Options.timestep;
+{ 
+  double mult=_percentage;
+  if (_percentage==BY_SUBBASIN_FLAG){
+    int p=pHRU->GetSubBasinIndex();
+    mult=_pModel_full->GetSubBasin(p)->GetFlushFract();
+  }
+  rates[0]=mult*storage[iFrom[0]]/Options.timestep;
 }
 
 //////////////////////////////////////////////////////////////////

src/HydroProcessABC.h

@@ -1,6 +1,6 @@
 /*----------------------------------------------------------------
   Raven Library Source Code
-  Copyright (c) 2008-2019 the Raven Development Team
+  Copyright (c) 2008-2025 the Raven Development Team
   ----------------------------------------------------------------*/
 
 #ifndef HYDROPROCESS_H
@@ -117,14 +117,16 @@ class CHydroProcessABC
 class CmvFlush: public CHydroProcessABC
 {
 private:/*------------------------------------------------------*/
-  double _percentage;
+  double  _percentage;
+  CModel *_pModel_full;
 
 public:/*-------------------------------------------------------*/
   //Constructors/destructors:
   CmvFlush(int          from_index,
            int          to_index,
            const double &pct,
-           CModelABC    *pM);
+           CModelABC    *pM,
+           CModel       *pM_full);
   ~CmvFlush();
 
   //inherited functions

src/ModelInitialize.cpp

@@ -775,7 +775,7 @@ void CModel::InitializeParameterOverrides()
     {
       double *address = this->_pGlobalParams->GetAddress(name);
       if (address==NULL){
-        ExitGracefully("CModel::InitializeParameterOverrides() : Invalid or unsupported global parameter name in :LocalParameterOverride command",BAD_DATA);
+        ExitGracefully("CModel::InitializeParameterOverrides() : Invalid or unsupported global parameter name in :GlobalParameterOverride command",BAD_DATA);
       }
       _pParamOverrides[i]->pxAddress=address;
 

src/ParseInput.cpp

@@ -2528,18 +2528,25 @@ 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 or 'FRACT_BY_SUBBASIN'}*/
       if (Options.noisy){cout <<"Flushing Process"<<endl;}
       double pct=1.0;
       if (Len<4){ImproperFormatWarning(":Flush",p,Options.noisy); break;}
       tmpS[0] = pStateVar->StringToSVType(s[2],tmpLev[0],true);
       tmpS[1] = pStateVar->StringToSVType(s[3],tmpLev[1],true);
       pModel->AddStateVariables(tmpS,tmpLev,2);
-      if ((Len>=5) && (s[4][0]!='#')){pct=max(min(s_to_d(s[4]),1.0),0.0);}
+      if (Len>=5){
+        if (to_string(s[4])=="FRACT_BY_SUBBASIN"){
+          pct=BY_SUBBASIN_FLAG;
+        }
+        else{
+          pct=max(min(s_to_d(s[4]),1.0),0.0);
+        }
+      }
 
       pMover = new CmvFlush(ParseSVTypeIndex(s[2], pModel, pStateVar),
                             ParseSVTypeIndex(s[3], pModel, pStateVar),
-                            pct, pModel);
+                            pct, pModel,pModel);
       AddProcess(pModel,pMover,pProcGroup);
       break;
     }

src/RavenInclude.h

@@ -254,6 +254,7 @@ const double  NETCDF_BLANK_VALUE      =-9999.0;
 const double  RAV_BLANK_DATA          =-1.2345;                                 ///< double corresponding to blank/void data item (also used in input files)
 const double  DIRICHLET_TEMP          =-9999.0;                                 ///< dirichlet concentration flag corresponding to air temperature
 const int     FROM_STATION_VAR        =-55;                                     ///< special flag indicating that NetCDF indices should be looked up from station attribute table
+const double  BY_SUBBASIN_FLAG        =64;                                      ///< special flag indicating flush percentage should be retrieved from subbasin parameter
 
 //Decision constants
 const double  HUGE_RESIST             =1e20;                                    ///< [d/mm]   essentially infinite resistance
@@ -865,7 +866,7 @@ enum assimtype
 //
 enum sv_type
 {
-  //Water Storage
+  // Water Storage Compartments 
   SURFACE_WATER,           ///< [mm] Streams & rivers: see surface_struct (REQUIRED)
   ATMOSPHERE,              ///< [mm] atmosphere : recieves water only!! (REQUIRED)
   ATMOS_PRECIP,            ///< [mm] atmosphere : provides water only!! (REQUIRED)
@@ -877,11 +878,11 @@ enum sv_type
   TRUNK,                   ///< [mm] water stored in trunks of trees
   ROOT,                    ///< [mm] water stored in roots
   GROUNDWATER,             ///< [mm] Deep groundwater
-  DEPRESSION,              ///< [mm] depression/surface storage
+  DEPRESSION,              ///< [mm] depression/surface/wetland storage
   SNOW,                    ///< [mm] frozen snow depth (mm SWE : snow water equivalent)
   NEW_SNOW,                ///< [mm] new snowfall waiting to be handled by snow balance (as SWE)
   SNOW_LIQ,                ///< [mm] liquid water content of snowpack
-  TOTAL_SWE,               ///< [mm] equivalent to SNOW[0]+SNOW[1]+...+SNOW_LIQ[0]..
+  TOTAL_SWE,               ///< [mm] equivalent to SNOW[0]+SNOW[1]+...+SNOW_LIQ[0].. (diagnostic variable)
   WETLAND,                 ///< [mm] deep wetland depression storage
   GLACIER,                 ///< [mm] Glacier melt/reservoir storage
   GLACIER_ICE,             ///< [mm] Glacier ice - typically assumed to be infinite reservoir.
@@ -893,14 +894,14 @@ enum sv_type
   // Distribution tracking variables
   MIN_DEP_DEFICIT,         ///< [mm or -1..0] Minimum depression deficit (describes deficit distribution), =-percent full if negative
 
-  // Memory variables
+  // Memory/diagnostic  variables
   CUM_INFIL,               ///< [mm] Cumulative infiltration to topsoil
   GA_MOISTURE_INIT,        ///< [mm] Initial topsoil moisture content for Green Ampt infiltration
   CUM_SNOWMELT,            ///< [mm] Cumulative snowmelt as SWE
   AET,                     ///< [mm] PET used up in given time step (diagnostic variable)
   RUNOFF,                  ///< [mm] net release of water to surface water in given times step (diagnostic variable)
 
-  //Temperature/Energy storage [C] or [MJ/m^2]
+  // Temperature/Energy storage [C] or [MJ/m^2]
   ENERGY_LOSSES,           ///< [MJ/m2] general energy losses
 
   SURFACE_WATER_TEMP,      ///< [C] Temperature of surface water
@@ -909,7 +910,7 @@ enum sv_type
   SOIL_TEMP,               ///< [C] Temperature of soil
   CANOPY_TEMP,             ///< [C] Temperature fo canopy
 
-  //Snow/Glacier variables
+  // Snow/Glacier variables
   SNOW_DEPTH,              ///< [mm] Snow depth - surrogate for density
   PERMAFROST_DEPTH,        ///< [mm] depth of permafrost
   THAW_DEPTH,              ///< [mm] depth of thaw
@@ -924,22 +925,22 @@ enum sv_type
 
   SNOW_ALBEDO,             ///< [-] Snow Surface albedo
 
-  //Crop variables
+  // Crop variables
   CROP_HEAT_UNITS,         ///< [-] cumulative crop heat units
 
-  //Transport variables
+  // Transport variables
   CONSTITUENT,             ///< chemical species [mg/m2], enthalpy [MJ/m2], or tracer [-]
   CONSTITUENT_SRC,         ///< chemical species [mg/m2], enthalpy [MJ/m2], or tracer [-] cumulative source
   CONSTITUENT_SW,          ///< chemical species [mg/m2], enthalpy [MJ/m2], or tracer [-] dumped to surface water
   CONSTITUENT_SINK,        ///< chemical species [mg/m2], enthalpy [MJ/m2], or tracer [-] cumulative sink (e.g., decay)
 
-  //Energy variables
+  // Energy variables
   LATENT_HEAT,             ///< [MJ/m2] cumulative energy lost to evaporative phase change
 
-  //Lateral exchange
+  // Lateral exchange
   LATERAL_EXCHANGE,        ///< [mm] water storage in transit from HRU awaiting lateral transfer to other HRUs
 
-  //Special - internal flags
+  // Special - internal flags
   STREAMFLOW,              ///< only used for referencing in data assimilation
   RESERVOIR_STAGE,         ///< only used for referencing in data assimilation
   UNRECOGNIZED_SVTYPE,     ///< Unrecognized type of state variable

src/StandardOutput.cpp

@@ -1343,6 +1343,9 @@ void CModel::WriteMinorOutput(const optStruct &Options,const time_struct &tt)
 
   // Write major output, if necessary
   //--------------------------------------------------------------
+  if (tt.model_time == 0.0) {//for ensembles
+    _currOutputTimeInd=0;
+  }
   if ((_nOutputTimes>0) && (_currOutputTimeInd<_nOutputTimes) && (tt.model_time>_aOutputTimes[_currOutputTimeInd]-0.5*Options.timestep))
   {
     string thishour=DecDaysToHours(tt.julian_day);

src/SubBasin.cpp

@@ -72,6 +72,7 @@ CSubBasin::CSubBasin( const long long      Identifier,
   _temperature_corr = 0.0;
   _unusable_flow_pct =0.0;
   _divert_fract      =1.0;
+  _flush_fract       =1.0;
 
   _res_disabled      =false;
   _assimilate        =false;
@@ -1085,6 +1086,12 @@ double CSubBasin::GetUnusableFlowPercentage() const {
 double CSubBasin::GetDivertFract       () const{
   return _divert_fract;
 }
+//////////////////////////////////////////////////////////////////
+/// \brief gets flush fraction (percentage of flush applied in basin)
+//
+double CSubBasin::GetFlushFract       () const{
+  return _flush_fract;
+}
 /*****************************************************************
    Manipulators
 *****************************************************************/
@@ -1137,6 +1144,7 @@ bool CSubBasin::SetBasinProperties(const string label,
   else if (!label_n.compare("GAMMA_SHAPE"   ))  {_gamma_shape=value;}
   else if (!label_n.compare("GAMMA_SCALE"   ))  {_gamma_scale=value;}
   else if (!label_n.compare("DIVERT_FRACT"  ))  {_divert_fract=value;}
+  else if (!label_n.compare("FLUSH_FRACT"   ))  {_flush_fract=value;}
 
   else if (!label_n.compare("Q_REFERENCE"   ))  {_Q_ref=value;}
   else if (!label_n.compare("MANNINGS_N"    ))  {_mannings_n=value;}
@@ -1224,6 +1232,7 @@ double CSubBasin::GetBasinProperties(const basin_props prop)
     case BP_GAMMA_SHAPE:      return _gamma_shape;
     case BP_GAMMA_SCALE:      return _gamma_scale;
     case BP_DIVERT_FRACT   :  return _divert_fract;
+    case BP_FLUSH_FRACT    :  return _flush_fract;
 
     case BP_Q_REFERENCE    : return _Q_ref;
     case BP_MANNINGS_N     : return _mannings_n;
@@ -1271,6 +1280,7 @@ double CSubBasin::GetBasinProperties(const string label) const
   else if (!label_n.compare("GAMMA_SHAPE"   ))  { return _gamma_shape;}
   else if (!label_n.compare("GAMMA_SCALE"   ))  { return _gamma_scale;}
   else if (!label_n.compare("DIVERT_FRACT"  ))  { return _divert_fract;}
+  else if (!label_n.compare("FLUSH_FRACT"   ))  { return _divert_fract;}
 
   else if (!label_n.compare("Q_REFERENCE"   ))  { return _Q_ref;}
   else if (!label_n.compare("MANNINGS_N"    ))  { return _mannings_n;}
@@ -1543,8 +1553,11 @@ double CSubBasin::AdjustAllFlows(const double &adjust, const bool overriding, co
 
   if(!overriding)
   {
+    
+    double distfact  = _pModel->GetGlobalParams()->GetParams()->assim_upstream_decay/M_PER_KM; //[1/km]->[1/m]
+    double corr=exp(-distfact*_reach_length);
     for(int n=0;n<_nQinHist; n++) {
-      _aQinHist[n]+=adjust*(_drainage_area-_basin_area)/_drainage_area;
+      _aQinHist[n]+=adjust*(_drainage_area-_basin_area)/_drainage_area*corr;
       upperswap(_aQinHist[n],0.0);
       va+=adjust*tstep*SEC_PER_DAY;
     }

src/SubBasin.h

@@ -97,6 +97,7 @@ class CSubBasin
   double            _recharge_corr;    ///< correction factor for recharge [-]
   double            _temperature_corr; ///< correction factor (additive) for temperature [-]
   double            _divert_fract;     ///< fraction of water flushed from source HRUs in basin [0..1] (:LateralDivert)
+  double            _flush_fract;      ///< fraction of water flushed from source compartment in basin [0..1] (:Flush with BY_SUBBASIN)
 
   int               _nSegments;        ///< Number of river segments used in routing( >=1)
 
@@ -224,6 +225,7 @@ class CSubBasin
   int                  GetNumWaterDemands   () const;
   double               GetUnusableFlowPercentage() const;
   double               GetDivertFract       () const;
+  double               GetFlushFract        () const;
 
   const double   *GetUnitHydrograph        () const;
   const double   *GetRoutingHydrograph     () const;