Deflate

src/Assimilate.cpp

@@ -1,6 +1,6 @@
 /*----------------------------------------------------------------
   Raven Library Source Code
-  Copyright (c) 2008-2026 the Raven Development Team
+  Copyright (c) 2008-2025 the Raven Development Team
   ----------------------------------------------------------------*/
 #include "RavenInclude.h"
 #include "Model.h"
@@ -25,6 +25,8 @@ void CModel::InitializeDataAssimilation(const optStruct &Options)
   // Initialize Streamflow assimilation
   if(Options.assimilate_flow)
   {
+    _aDAscale     =new double[_nSubBasins];
+    _aDAscale_last=new double[_nSubBasins];
     _aDAQadjust   =new double[_nSubBasins];
     _aDADrainSum  =new double[_nSubBasins];
     _aDADownSum   =new double[_nSubBasins];
@@ -35,6 +37,8 @@ void CModel::InitializeDataAssimilation(const optStruct &Options)
     _aDAobsQ2     =new double[_nSubBasins];
     _aDASinceLastBlank=new double [_nSubBasins];
     for(int p=0; p<_nSubBasins; p++) {
+      _aDAscale     [p]=1.0;
+      _aDAscale_last[p]=1.0;
       _aDAQadjust   [p]=0.0;
       _aDADrainSum  [p]=0.0;
       _aDADownSum   [p]=0.0;
@@ -77,15 +81,14 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim
     double Qobs2,Qmod,Qmodlast;
     double alpha = _pGlobalParams->GetParams()->assimilation_fact;
 
+    //alpha*=(1.0-exp(-0.06*_aDASinceLastBlank[p])); //shock/oscillation prevention
+
     Qobs    = _aDAobsQ[p];
     Qobs2   = _aDAobsQ2[p];
     Qmod    = _pSubBasins[p]->GetOutflowRate();
     Qmodlast= _pSubBasins[p]->GetLastOutflowRate();
-    if (_pSubBasins[p]->GetReservoir()!=NULL){
-      Qmod    = _pSubBasins[p]->GetReservoir()->GetOutflowRate   (false); //unadjusted flows
-      Qmodlast= _pSubBasins[p]->GetReservoir()->GetOldOutflowRate(false);
-    }
     if((Qmod>PRETTY_SMALL) && (Qobs!=RAV_BLANK_DATA)){
+      _aDAscale  [p]=1.0+alpha*((Qobs-Qmod)/Qmod); //if alpha = 1, Q=Qobs in observation basin
       //_aDAQadjust[p]=alpha*(Qobs-Qmod); //Option A: end of time step flow
       //_aDAQadjust[p]=0.5*alpha*(2.0*Qobs-Qmodlast-Qmod);//Option B: mean flow - rapidly oscillatory Q - Qmean is perfect
       if (Qobs2 == RAV_BLANK_DATA) { //Option C: aim for midpoint of two observation flows (this is the way)
@@ -96,14 +99,19 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim
       }
     }
     else {
+      _aDAscale  [p]=1.0;
       _aDAQadjust[p]=0.0;
     }
+    _aDAscale_last[p]=1.0;//no need to scale previous
   }
 
-  // actually adjust flows
+  // actually scale flows
   //---------------------------------------------------------------------
   double mass_added=0.0;
-  if (Options.assim_method==DA_ECCC) {
+  if (Options.assim_method==DA_RAVEN_DEFAULT){
+    mass_added=_pSubBasins[p]->ScaleAllFlows(_aDAscale[p]/_aDAscale_last[p],_aDAoverride[p],Options.timestep,tt.model_time);
+  }
+  else if (Options.assim_method==DA_ECCC) {
     if(_aDAoverride[p]){
       mass_added=_pSubBasins[p]->AdjustAllFlows(_aDAQadjust[p],true,true,Options.timestep,tt.model_time);
     }
@@ -133,6 +141,7 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
   bool   ObsExists; //observation available in THIS basin
 
   for(p=0; p<_nSubBasins; p++) {
+    _aDAscale_last[p]=_aDAscale[p];
     _aDADrainSum  [p]=0.0;
     _aDAlength    [p]=0.0;//reboot every timestep
   }
@@ -155,7 +164,7 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
       {
         if(IsContinuousFlowObs2(_pObservedTS[i],_pSubBasins[p]->GetID()))//flow observation is available and linked to this subbasin
         {
-          Qobs  = _pObservedTS[i]->GetSampledValue(nn);   //mean timestep flow
+          Qobs = _pObservedTS[i]->GetSampledValue(nn); //mean timestep flow
           Qobs2 = _pObservedTS[i]->GetSampledValue(nn+1); //mean timestep flow
 
           //override initial conditions directly
@@ -189,6 +198,7 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
     if (ObsExists) {
       if((Qobs!=RAV_BLANK_DATA) && (tt.model_time<t_observationsOFF))
       {
+        //_aDAscale[p] calculated live in AssimilationOverride when up-to-date modelled flow available
         _aDAlength   [p]=0.0;
         _aDAtimesince[p]=0.0;
         _aDAoverride [p]=true;
@@ -198,6 +208,7 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
       }
       else
       { //found a blank or zero flow value
+        _aDAscale    [p]=_aDAscale[p];//same adjustment as before - scaling persists
         if(pdown!=DOESNT_EXIST) {
           _aDAlength   [p]+=_pSubBasins  [pdown]->GetReachLength(); //length propagates from below
         }
@@ -220,11 +231,13 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
           (_pSubBasins[p]->GetReservoir()==NULL) &&
           (!_aDAoverride[p]) &&
           (_pSubBasins[p]->IsEnabled()) && (_pSubBasins[pdown]->IsEnabled())) {
+        _aDAscale      [p]= _aDAscale    [pdown];
         _aDAlength     [p]+=_pSubBasins  [pdown]->GetReachLength();
         _aDAtimesince  [p]= _aDAtimesince[pdown];
         _aDAoverride   [p]=false;
       }
       else{ //Nothing downstream or reservoir present in this basin, no assimilation
+        _aDAscale    [p]=1.0;
         _aDAlength   [p]=0.0;
         _aDAtimesince[p]=0.0;
         _aDAoverride [p]=false;
@@ -270,6 +283,15 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
       _aDADownSum[p]=_aDADownSum[pdown];
     }
   }
+
+  // Apply time and space correction factors
+  //----------------------------------------------------------------
+  double time_fact = _pGlobalParams->GetParams()->assim_time_decay;
+  double distfact  = _pGlobalParams->GetParams()->assim_upstream_decay/M_PER_KM; //[1/km]->[1/m]
+  for(p=0; p<_nSubBasins; p++)
+  {
+    _aDAscale[p] =1.0+(_aDAscale[p]-1.0)*exp(-distfact*_aDAlength[p])*exp(-time_fact*_aDAtimesince[p]);
+  }
 }
 /////////////////////////////////////////////////////////////////
 /// \brief Calculates updated DA scaling coefficients for all subbasins for forthcoming timestep

src/ChannelXSect.cpp

@@ -260,17 +260,17 @@ double      CChannelXSect::GetNPoints()          const { return _nPoints; }
 /// \brief Returns Riverbed slope
 /// \return Riverbed slope [m/m]
 //
-double      CChannelXSect::GetBedslope   ()             const { return _bedslope;}
+double      CChannelXSect::GetBedslope()         const { return _bedslope;}
 
-double      CChannelXSect::GetAQAt       (const int i)  const { return _aQ[i];}
+double      CChannelXSect::GetAQAt(const int i)  const { return _aQ[i];}
 
-double      CChannelXSect::GetAStageAt   (const int i)  const { return _aStage[i];}
+double      CChannelXSect::GetAStageAt(const int i)  const { return _aStage[i];}
 
 double      CChannelXSect::GetATopWidthAt(const int i)  const { return _aTopWidth[i];}
 
-double      CChannelXSect::GetAXAreaAt   (const int i)  const { return _aXArea[i];}
+double      CChannelXSect::GetAXAreaAt(const int i)  const { return _aXArea[i];}
 
-double      CChannelXSect::GetAPerimAt   (const int i)  const { return _aPerim[i];}
+double      CChannelXSect::GetAPerimAt(const int i)  const { return _aPerim[i];}
 
 /*****************************************************************
    Rating Curve Interpolation functions
@@ -381,8 +381,8 @@ double  CChannelXSect::GetCelerity(const double &Q, const double &SB_slope,const
 
 //////////////////////////////////////////////////////////////////
 /// \brief Returns diffusivity of channel [m/s]
-/// \param Q [in] Channel flowrate [m3/s]
-/// \param bedslope [in] Slope of riverbed [m/m]
+/// \param &Q [in] Channel flowrate [m3/s]
+/// \param &bedslope [in] Slope of riverbed [m/m]
 /// \return Diffusivity of channel at specified flowrate and slope [m3/s]
 //
 double CChannelXSect::GetDiffusivity(const double &Q, const double &SB_slope, const double &SB_n) const
@@ -398,7 +398,7 @@ double CChannelXSect::GetDiffusivity(const double &Q, const double &SB_slope, co
 
 //////////////////////////////////////////////////////////////////
 /// \brief checks to see if reference flow is in excess of channel maximum; throws warning if it is
-/// \param Qref [in] Reference flowrate [m3/s]
+/// \param &Qref [in] Reference flowrate [m3/s]
 /// \param SB_slope [in] subbasin slope (or AUTO_COMPUTE, if channel slope to be used)
 /// \param SB_n [in] subbasin mannings  (or AUTO_COMPUTE, if channel mannings to be used)
 /// \param SBID [in] subbasin identifier
@@ -454,15 +454,15 @@ void CChannelXSect::GetPropsFromProfile(const double &elev,
       Pi=sqrt(dx*dx+(zu-zl)*(zu-zl));
       if ((i>0) && (_aElev[i-1]>_aElev[i]) && (_aX[i-1]==_aX[i])) //handles straight adjacent sides (left)
       {
-        if(elev<=_aElev[i-1]){Pi+=(elev       -_aElev[i]);}
-        else                 {Pi+=(_aElev[i-1]-_aElev[i]);}
+        if(elev<=_aElev[i-1]){Pi+=(elev      -_aElev[i]);}
+        else                {Pi+=(_aElev[i-1]-_aElev[i]);}
       }
       if ((i<(_nSurveyPts-2)) && (_aElev[i+2]>_aElev[i+1]) && (_aX[i+2]==_aX[i+1])) //handles straight adjacent sides (right)
       {
-        if(elev<=_aElev[i+2]){Pi+=(elev       -_aElev[i+1]);}
-        else                 {Pi+=(_aElev[i+2]-_aElev[i+1]);}
+        if(elev<=_aElev[i+2]){Pi+=(elev      -_aElev[i+1]);}
+        else                {Pi+=(_aElev[i+2]-_aElev[i+1]);}
       }
-      if (i==0)              {Pi+=(elev-_aElev[i]  );}
+      if (i==0)             {Pi+=(elev-_aElev[i]  );}
       if (i==_nSurveyPts-2)  {Pi+=(elev-_aElev[i+1]);}
     }
     else  //partially wet part of profile (includes riverbank)

src/CustomOutput.cpp

@@ -217,7 +217,49 @@ void CCustomOutput::SetHistogramParams(const double minv,const double maxv, cons
   _hist_max=maxv;
   _nBins=numBins;
 }
+//////////////////////////////////////////////////////////////////
+/// \brief Compute the approximate number of time steps for the output array. This is not meant to be exact. Do not allocate memory based on this value.
+/// \param Options [in] Global model options information
+/// \return Number of time steps (int)
+int CCustomOutput::ApproximateNumTimeSteps(const optStruct &Options) const
+{
+  int nsteps = 0;  // return value
+  time_struct start, end;
+
+  JulianConvert(0, Options.julian_start_day, Options.julian_start_year, Options.calendar, start);
+  JulianConvert(Options.duration, Options.julian_start_day, Options.julian_start_year, Options.calendar, end);
 
+  switch(_timeAgg) {
+    case YEARLY:
+      nsteps = int(end.year - start.year);
+      break;
+    case WATER_YEARLY:
+    // wateryr_mo: Starting month of the water year
+      nsteps = int(end.year - start.year);
+      if ((start.month < end.month) && (start.month < Options.wateryr_mo) && (end.month >= Options.wateryr_mo)) nsteps++;
+      if ((start.month > end.month) && (start.month >= Options.wateryr_mo) && (end.month < Options.wateryr_mo)) nsteps--;
+      break;
+    case MONTHLY:
+      nsteps = int((end.year - start.year) * 12 + (end.month - start.month));
+      break;
+    case DAILY:
+      nsteps = int(ceil(Options.duration));
+      break;
+    case EVERY_NDAYS:
+      nsteps = int(ceil(Options.duration / Options.custom_interval));
+      break;
+    case EVERY_TSTEP:
+      nsteps = int(ceil(Options.duration / Options.timestep));
+      break;
+    case ENTIRE_SIM:
+      nsteps = 1;
+      break;
+    default:
+      nsteps = 0;
+      break;
+  }
+  return max(0, nsteps); // Ensure non-negative and account for partial time step at end of simulation
+}
 ///////////////////////////////////////////////////////////////////
 /// \brief Allocates memory and initialize data storage of a CCustomOutput object
 /// \remarks Called prior to simulation. Determines size of and allocates memory for (member) data[][] array needed in statistical calculations
@@ -554,6 +596,7 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options)
 
   int         retval;                                // error value for NetCDF routines
   size_t      start[1], count[1];                    // determines where and how much will be written to NetCDF
+  size_t      chunksize2[2];                          // Chunksize (time, data)
   string      tmp,tmp2,tmp3,tmp4;
 
   bool cant_support=(_aggstat==AGG_RANGE || _aggstat==AGG_95CI || _aggstat==AGG_QUARTILES || _aggstat==AGG_HISTOGRAM);
@@ -574,12 +617,21 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options)
   // time
   // ----------------------------------------------------------
   // (a) Define the DIMENSIONS. NetCDF will hand back an ID for each.
+
   retval = nc_def_dim(_netcdf_ID, "time", NC_UNLIMITED, &time_dimid);              HandleNetCDFErrors(retval);
 
   // (b) Define the time variable.
   dimids1[0] = time_dimid;
   retval = nc_def_var(_netcdf_ID, "time", NC_DOUBLE, ndims1,dimids1, &varid_time); HandleNetCDFErrors(retval);
 
+  // Enable deflate compression for time variable (shuffle, zlib, deflate_level)
+  retval = nc_def_var_deflate(_netcdf_ID, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval);
+
+  // Set chunksize to len(time)
+  chunksize2[0] = ApproximateNumTimeSteps(Options) + 1;
+  retval = nc_def_var_chunking(_netcdf_ID, varid_time, NC_CHUNKED, &chunksize2[0]); HandleNetCDFErrors(retval);
+
+
   // (c) Assign units attributes to the netCDF VARIABLES.
   //     --> converts start day into "hours since YYYY-MM-DD HH:MM:SS"
   char  starttime[200]; // start time string in format 'days since YYY-MM-DD HH:MM:SS'
@@ -633,6 +685,13 @@ void CCustomOutput::WriteNetCDFFileHeader(const optStruct &Options)
     dimids2[1] = ndata_dimid;
     retval = nc_def_var(_netcdf_ID, netCDFtag.c_str(), NC_DOUBLE, ndims2, dimids2, &varid_data);    HandleNetCDFErrors(retval);
 
+    // Enable deflate compression for data variable
+    retval = nc_def_var_deflate(_netcdf_ID, varid_data, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval);
+
+    // Set chunksizes for data variable (time, ndata)
+    chunksize2[1] = max((size_t)1, min(_nData, (size_t)(NETCDF_CHUNKSIZE_MB * 1024 * 1024 / sizeof(double) / chunksize2[0]))); // Ensure at least one basin per chunk
+    retval = nc_def_var_chunking(_netcdf_ID, varid_data, NC_CHUNKED, chunksize2); HandleNetCDFErrors(retval);
+
     //(f) set some attributes to variable _netCDFtag
     tmp=_timeAggStr+" "+_statStr+" "+_varName+" "+_spaceAggStr;
     tmp2=_varUnits;
@@ -1267,4 +1326,4 @@ CCustomOutput *CCustomOutput::ParseCustomOutputCommand(char *s[MAXINPUTITEMS], c
   }
 
   return pCustom;
-}
+}

src/CustomOutput.h

@@ -126,7 +126,10 @@ class CCustomOutput
 
   void        SetHistogramParams(const double min,const double max, const int numBins);
 
+  int         ApproximateNumTimeSteps(const optStruct &Options) const;
+
   spatial_agg GetSpatialAgg() const;
+
   void        InitializeCustomOutput(const optStruct &Options);
 
   void        WriteFileHeader  (const optStruct &Options);

src/Makefile

@@ -15,6 +15,7 @@
 #            Feb 2025    James Craig and Bohao Tang - support for new Mac (Applie Silicon)
 ###
 # appname := raven_rev$(shell svnversion -n . | sed -e 's/:/-/g')_$(shell uname -o).exe
+# To debug, add -g to CXXFLAGS and LDFLAGS
 appname := Raven.exe
 
 CXX      := g++

src/Model.cpp

@@ -58,6 +58,8 @@ CModel::CModel(const int        nsoillayers,
   _aOrderedSBind  =NULL;
   _aDownstreamInds=NULL;
 
+  _aDAscale       =NULL; //Initialized in InitializeDataAssimilation
+  _aDAscale_last  =NULL;
   _aDAQadjust     =NULL;
   _aDADrainSum    =NULL;
   _aDADownSum     =NULL;
@@ -214,6 +216,8 @@ CModel::~CModel()
   delete [] _aOutputTimes;   _aOutputTimes=NULL;
   delete [] _aObsIndex;      _aObsIndex=NULL;
 
+  delete [] _aDAscale;       _aDAscale=NULL;
+  delete [] _aDAscale_last;  _aDAscale_last=NULL;
   delete [] _aDAQadjust;     _aDAQadjust=NULL;
   delete [] _aDADrainSum;    _aDADrainSum=NULL;
   delete [] _aDADownSum;     _aDADownSum=NULL;

src/Model.h

@@ -1,6 +1,6 @@
 /*----------------------------------------------------------------
   Raven Library Source Code
-  Copyright (c) 2008-2026 the Raven Development Team
+  Copyright (c) 2008-2025 the Raven Development Team
   ----------------------------------------------------------------*/
 #ifndef MODEL_H
 #define MODEL_H
@@ -143,6 +143,7 @@ class CModel: public CModelABC
   int          _nAggDiagnostics;  ///< number of aggregated diagnostics
 
   //Data Assimilation
+  double             *_aDAscale; ///< array of data assimilation flow scaling parameters [size: _nSubBasins] (NULL w/o DA)
   double            *_aDAlength; ///< array of downstream distance to nearest DA observation [m] [size: _nSubBasins] (NULL w/o DA)
   double           *_aDAQadjust; ///< array of flow adjustments [m3/s] [size: _nSubBasins]
   double          *_aDADrainSum; ///< sum of assimilated drainage areas upstream of a subbasin outlet [km2] [size: _nSubBasins]
@@ -152,6 +153,7 @@ class CModel: public CModelABC
   double              *_aDAobsQ; ///< array of observed flow values in basins [size: _nSubBasins]  (NULL w/o DA)
   double             *_aDAobsQ2; ///< array of observed flow values in next time step [size: _nSubBasins]   (NULL w/o DA)
   double    *_aDASinceLastBlank; ///< array of time steps since last blank value occurred [size: _nSubBasins] (NULL w/o DA)
+  double        *_aDAscale_last; ///< array of scale factors from previous time step  [size: _nSubBasins]  (NULL w/o DA)
 
   force_perturb**_pPerturbations; ///< array of pointers to perturbation data; defines which forcing functions to perturb and how [size: _nPerturbations]
   int            _nPerturbations; ///< number of forcing functions to perturb

src/ModelInitialize.cpp

@@ -70,6 +70,7 @@ void CModel::Initialize(const optStruct &Options)
       }
     }
   }
+
   // reserve memory for mass balance arrays
   //--------------------------------------------------------------
   _aCumulativeBal = new double * [_nHydroUnits];

src/ParseInitialConditionFile.cpp

@@ -647,7 +647,7 @@ bool ParseInitialConditionsFile(CModel *&pModel, const optStruct &Options)
             pBasin->GetReservoir()->SetControlFlow(s_to_i(s[1]),s_to_d(s[2]),s_to_d(s[3]));
           }
         }
-        else if(!strcmp(s[0],":ResDAadj"))
+        else if(!strcmp(s[0],":ResDAscale"))
         {
           if(Len>=3) {
             pBasin->GetReservoir()->SetDataAssimFactors(s_to_d(s[1]),s_to_d(s[2]));