diff --git a/src/Assimilate.cpp b/src/Assimilate.cpp
index 73fedc8..a33a570 100644
--- a/src/Assimilate.cpp
+++ b/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
diff --git a/src/ChannelXSect.cpp b/src/ChannelXSect.cpp
index d950a88..5851632 100644
--- a/src/ChannelXSect.cpp
+++ b/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)
diff --git a/src/CustomOutput.cpp b/src/CustomOutput.cpp
index a2c7d60..6cad104 100644
--- a/src/CustomOutput.cpp
+++ b/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;
-}
+}
diff --git a/src/CustomOutput.h b/src/CustomOutput.h
index d9fb177..461a98f 100644
--- a/src/CustomOutput.h
+++ b/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);
diff --git a/src/Makefile b/src/Makefile
index d897009..738277a 100644
--- a/src/Makefile
+++ b/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++
diff --git a/src/Model.cpp b/src/Model.cpp
index 0c4fe51..ca92b36 100644
--- a/src/Model.cpp
+++ b/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;
diff --git a/src/Model.h b/src/Model.h
index b55659b..a6ea20a 100644
--- a/src/Model.h
+++ b/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
diff --git a/src/ModelInitialize.cpp b/src/ModelInitialize.cpp
index 6a69e75..93fff82 100644
--- a/src/ModelInitialize.cpp
+++ b/src/ModelInitialize.cpp
@@ -70,6 +70,7 @@ void CModel::Initialize(const optStruct &Options)
       }
     }
   }
+
   // reserve memory for mass balance arrays
   //--------------------------------------------------------------
   _aCumulativeBal = new double * [_nHydroUnits];
diff --git a/src/ParseInitialConditionFile.cpp b/src/ParseInitialConditionFile.cpp
index db06b12..511d418 100644
--- a/src/ParseInitialConditionFile.cpp
+++ b/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]));
diff --git a/src/ParseInput.cpp b/src/ParseInput.cpp
index 0983819..9087a5f 100644
--- a/src/ParseInput.cpp
+++ b/src/ParseInput.cpp
@@ -332,7 +332,7 @@ bool ParseMainInputFile (CModel     *&pModel,
   Options.aNetCDFattribs          =NULL;
   Options.assimilate_flow         =false;
   Options.assimilate_stage        =false;
-  Options.assim_method            =DA_ECCC;
+  Options.assim_method            =DA_RAVEN_DEFAULT;
   Options.assimilation_start      =-1.0; //start before simulation
   Options.sv_override_endtime     =ALMOST_INF;
   Options.time_zone               =0;
@@ -1812,7 +1812,8 @@ bool ParseMainInputFile (CModel     *&pModel,
     {/*:AssimilationMethod [method]*/
       if(Options.noisy) { cout << "Assimilation Method" << endl; }
       if (Len<2){ImproperFormatWarning(":AssimilationMethod",p,Options.noisy); break;}
-      if (!strcmp(s[1],"DA_ECCC"                 )){Options.assim_method=DA_ECCC;}
+      if      (!strcmp(s[1],"DA_RAVEN_DEFAULT"        )){Options.assim_method=DA_RAVEN_DEFAULT;}
+      else if (!strcmp(s[1],"DA_ECCC"                 )){Options.assim_method=DA_ECCC;}
       break;
     }
     case(97):  //--------------------------------------------
@@ -3729,6 +3730,12 @@ bool ParseMainInputFile (CModel     *&pModel,
   ExitGracefullyIf(Options.duration<0,
                    "ParseMainInputFile::Model duration less than zero. Make sure :EndDate is after :StartDate.",BAD_DATA_WARN);
 
+  // Compute the size of the output's time dimension
+  Options.n_out_time = (int)(floor(ceil((Options.duration + TIME_CORRECTION)/Options.timestep))/Options.output_interval);
+  if (Options.n_out_time==0) {
+    WriteAdvisory("ParseMainInputFile::Number of output time steps is zero. Check :Duration, :Timestep, and :OutputInterval commands.",BAD_DATA);
+  }
+
   if((Options.nNetCDFattribs>0) && (Options.output_format!=OUTPUT_NETCDF)){
     WriteAdvisory("ParseMainInputFile: NetCDF attributes were specified but output format is not NetCDF.",Options.noisy);
   }
diff --git a/src/ParseLib.h b/src/ParseLib.h
index 38d0375..5ffedce 100644
--- a/src/ParseLib.h
+++ b/src/ParseLib.h
@@ -28,7 +28,7 @@ typedef const double * const * const Ironclad2DArray;     ///< unmodifiable 2d a
 //************************************************************************
 const int    MAXINPUTITEMS  =   1000;  ///< maximum delimited input items per line
 const int    MAXCHARINLINE  =  10000;  ///< maximum characters in line
-const bool   parserdebug=false;        ///< turn to true for debugging of parser
+const bool   parserdebug=false;   ///< turn to true for debugging of parser
 
 ///////////////////////////////////////////////////////
 /// \brief Types of parsing errors
diff --git a/src/ParseManagementFile.cpp b/src/ParseManagementFile.cpp
index 1c0077a..93835b7 100644
--- a/src/ParseManagementFile.cpp
+++ b/src/ParseManagementFile.cpp
@@ -189,7 +189,6 @@ void ParseManagementFile(CModel *&pModel,const optStruct &Options)
     else if(!strcmp(s[0],":UserTimeSeries"))              { code=70; }
 
     else if(!strcmp(s[0],":NonlinearVariable"))           { code=80; }
-    else if(!strcmp(s[0],":NonLinearVariable"))           { code=80; }
     else if(!strcmp(s[0],":MaxNLSolverIterations"))       { code=81; }
     else if(!strcmp(s[0],":NLSolverTolerance"))           { code=82; }
     else if(!strcmp(s[0],":NLRelaxationCoeff"))           { code=83; }
diff --git a/src/ParsePropertyFile.cpp b/src/ParsePropertyFile.cpp
index 3911732..fa4de19 100644
--- a/src/ParsePropertyFile.cpp
+++ b/src/ParsePropertyFile.cpp
@@ -9,8 +9,6 @@
 #include "GlobalParams.h"
 #include "SoilAndLandClasses.h"
 #include <string>
-
-
 struct val_alias
 {
   string tag;
@@ -25,10 +23,9 @@ double AutoOrDoubleOrAlias(const string s, val_alias **aAl,const int nAl)
   }
   return AutoOrDouble(s);
 }
-bool  ParsePropArray        (CParser *p,int *&indices,double **&properties,
-                             int &num_read,string *tags,const int line_length,const int max_classes,
-                             val_alias **pAliases,  const int nAliases);
-void  DeletePropArray       (int *indices,double **properties,int num_read);
+bool ParsePropArray(CParser *p,int *indices,double **properties,
+                    int &num_read,string *tags,const int line_length,const int max_classes,
+                    val_alias **pAliases,  const int nAliases);
 void  RVPParameterWarning   (string *aP, class_type *aPC, int &nP, CModel* pModel);
 void  CreateRVPTemplate     (string *aP, class_type *aPC, int &nP, const optStruct &Options);
 void  ImproperFormatWarning(string command,CParser *p,bool noisy);
@@ -63,14 +60,14 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
   global_struct     parsed_globals;
 
   int               num_parsed_veg=0;
-  CVegetationClass **pVegClasses=NULL;
-  veg_struct       **parsed_veg =NULL;
-  string            *vegtags    =NULL;
+  CVegetationClass *pVegClasses[MAX_VEG_CLASSES];
+  veg_struct        parsed_veg [MAX_VEG_CLASSES];
+  string            vegtags    [MAX_VEG_CLASSES];
 
   int               num_parsed_lult = 0;  // counter for the number of LULT classes parsed
-  surface_struct  **parsed_surf=NULL;
-  CLandUseClass   **pLUClasses =NULL;
-  string           *lulttags   =NULL;
+  surface_struct    parsed_surf [MAX_LULT_CLASSES];
+  CLandUseClass    *pLUClasses  [MAX_LULT_CLASSES];
+  string            lulttags    [MAX_LULT_CLASSES];
 
   int               num_parsed_soils=0;
   CSoilClass      **pSoilClasses=NULL;
@@ -83,14 +80,14 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
   string            terraintags [MAX_TERRAIN_CLASSES];
 
   int               num_parsed_profiles=0;
-  CSoilProfile     **pProfiles=NULL;
+  CSoilProfile     *pProfiles   [MAX_SOIL_PROFILES];
 
   int               num_parsed_aqstacks=0;
 
   bool              invalid_index;
   int               num_read;
-  int              *indices=NULL;
-  double          **properties=NULL;
+  int              *indices;
+  double          **properties;
   string            aParamStrings[MAXINPUTITEMS];
   int               nParamStrings=0;
 
@@ -100,15 +97,26 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
   pModel->GetGlobalParams()->InitializeGlobalParameters(global_template,true);
   pModel->GetGlobalParams()->InitializeGlobalParameters(parsed_globals,false);
 
+  CVegetationClass::InitializeVegetationProps ("[DEFAULT]",parsed_veg[0],true);//zero-index canopy is template
+  vegtags    [0]="[DEFAULT]";
+  num_parsed_veg++;
 
-  AddNewVegClass (pVegClasses,  parsed_veg,  vegtags, num_parsed_veg,"[DEFAULT]",true,pModel);
   AddNewSoilClass(pSoilClasses, parsed_soils, soiltags, num_parsed_soils, pModel->GetTransportModel()->GetNumConstituents(), "[DEFAULT]", true, pModel);//zero-index soil is template
-  AddNewLULTClass(pLUClasses,   parsed_surf,  lulttags, num_parsed_lult,"[DEFAULT]",true,pModel);
+
+  CLandUseClass::InitializeSurfaceProperties("[DEFAULT]", parsed_surf[0], true); // zero-index LULT is template
+  lulttags    [0]="[DEFAULT]";
+  num_parsed_lult++;
 
   CTerrainClass::InitializeTerrainProperties(parsed_terrs[0],true);//zero-index terrain is template
   terraintags [0]="[DEFAULT]";
   num_parsed_terrs++;
 
+  indices   =new int     [MAX_NUM_IN_CLASS];
+  properties=new double *[MAX_NUM_IN_CLASS];
+  for (int i=0;i<MAX_NUM_IN_CLASS;i++){
+    properties[i]=new double [MAX_PROPERTIES_PER_LINE];
+  }
+
   const int MAX_PARAMS=100;
   int         nP;
   string     *aP  =new string [MAX_PARAMS];
@@ -183,17 +191,17 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
     }
     else if (aPCmaster[p]==CLASS_VEGETATION)
     {
-      val=CVegetationClass::GetVegetationProperty(*parsed_veg[0],aPmaster[p]);
+      val=CVegetationClass::GetVegetationProperty(parsed_veg[0],aPmaster[p]);
       if (val==NOT_NEEDED_AUTO){val=AUTO_COMPUTE;}
       else if (val==NOT_NEEDED){val=NOT_SPECIFIED;}
-      CVegetationClass::SetVegetationProperty    (*parsed_veg[0],aPmaster[p],val);
+      CVegetationClass::SetVegetationProperty    (parsed_veg[0],aPmaster[p],val);
     }
     else if (aPCmaster[p]==CLASS_LANDUSE   )
     {
-      val = CLandUseClass::GetSurfaceProperty    (*parsed_surf[0], aPmaster[p]);
+      val = CLandUseClass::GetSurfaceProperty(parsed_surf[0], aPmaster[p]);
       if (val == NOT_NEEDED_AUTO){val = AUTO_COMPUTE;}
       else if (val == NOT_NEEDED){val = NOT_SPECIFIED;}
-      CLandUseClass::SetSurfaceProperty          (*parsed_surf[0], aPmaster[p], val);
+      CLandUseClass::SetSurfaceProperty(parsed_surf[0], aPmaster[p], val);
     }
     else if (aPCmaster[p]==CLASS_TERRAIN   )
     {
@@ -443,10 +451,10 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
         if      (IsComment(s[0], Len)){}//comment line
         else if (Len>=2)
         {
-          //dynamically add to stack
-          int tmp=num_parsed_profiles;
-          CSoilProfile *pNewSoil= new CSoilProfile(s[0],pModel);
-          DynArrayAppend((void**&)(pProfiles),(void*)(pNewSoil),tmp);
+          if (num_parsed_profiles>=MAX_SOIL_PROFILES-1){
+            ExitGracefully("ParseClassPropertiesFile: exceeded maximum # of soil profiles",BAD_DATA);}
+
+          pProfiles[num_parsed_profiles] = new CSoilProfile(s[0], pModel);
 
           int nhoriz=s_to_i(s[1]);
           ExitGracefullyIf(nhoriz<0,
@@ -529,7 +537,6 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
           CSoilClass::SetSoilProperty(*parsed_soils[indices[i]],aParamStrings[j+1],properties[i][j]);
         }
       }
-      DeletePropArray(indices,properties,num_read);
       bool found_in_master;
       for (int j=0;j<nParamStrings-1;j++){
         found_in_master=false;
@@ -560,13 +567,17 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
         else if (!strcmp(s[0],":Units")){} //units are explicit within Raven - useful for GUIs
         else if ((Len==3) || (Len==4))
         {
+          if (num_parsed_lult>=MAX_LULT_CLASSES-1){
+            ExitGracefully("ParseClassPropertiesFile: exceeded maximum # of LU/LT classes",BAD_DATA);}
 
-          AddNewLULTClass(pLUClasses,parsed_surf,lulttags,num_parsed_lult,s[0],false,pModel);
-
-          parsed_surf[num_parsed_lult-1]->impermeable_frac = s_to_d(s[1]);
+          CLandUseClass::InitializeSurfaceProperties(s[0], parsed_surf[num_parsed_lult], false);
+          lulttags  [num_parsed_lult]   = s[0];
+          pLUClasses[num_parsed_lult-1] = new CLandUseClass(s[0], pModel);
+          parsed_surf[num_parsed_lult].impermeable_frac = s_to_d(s[1]);
           if (Len>=3){
-            parsed_surf[num_parsed_lult-1]->forest_coverage = s_to_d(s[2]);
+            parsed_surf[num_parsed_lult].forest_coverage = s_to_d(s[2]);
           }
+          num_parsed_lult++;
         }
         else{
           ImproperFormatWarning(":LandUseClasses",p,Options.noisy);
@@ -618,12 +629,11 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
       for (int i=0; i<num_read; i++)
       {
         for (int j=0;j<nParamStrings-1;j++) {
-          CLandUseClass::SetSurfaceProperty(*parsed_surf[indices[i]],
+          CLandUseClass::SetSurfaceProperty(parsed_surf[indices[i]],
                                             aParamStrings[j+1],
                                             properties[i][j]);
         }
       }
-      DeletePropArray(indices,properties,num_read);
       bool found_in_master;
       for (int j=0;j<nParamStrings-1;j++){
         found_in_master=false;
@@ -665,11 +675,17 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
         else if (!strcmp(s[0],":Units")){} //units are explicit within Raven - useful for GUIs
         else if (Len==4)
         {
-          AddNewVegClass(pVegClasses,parsed_veg,vegtags,num_parsed_veg,s[0],false,pModel);
+          if (num_parsed_veg>=MAX_VEG_CLASSES-1){
+            ExitGracefully("ParseClassPropertiesFile: exceeded maximum # of vegetation classes",BAD_DATA);}
+
+          CVegetationClass::InitializeVegetationProps(s[0], parsed_veg [num_parsed_veg], false);
+          pVegClasses[num_parsed_veg-1] = new CVegetationClass(s[0], pModel);
+          vegtags    [num_parsed_veg] = s[0];
+          parsed_veg [num_parsed_veg].max_height    = s_to_d(s[1]);
+          parsed_veg [num_parsed_veg].max_LAI       = s_to_d(s[2]);
+          parsed_veg [num_parsed_veg].max_leaf_cond = s_to_d(s[3]);
 
-          parsed_veg [num_parsed_veg-1]->max_height    = s_to_d(s[1]);
-          parsed_veg [num_parsed_veg-1]->max_LAI       = s_to_d(s[2]);
-          parsed_veg [num_parsed_veg-1]->max_leaf_cond = s_to_d(s[3]);
+          num_parsed_veg++;
         }
         else{
           ImproperFormatWarning(":VegetationClasses",p,Options.noisy);  break;
@@ -691,10 +707,9 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
                        "ParseClassPropertiesFile: Invalid vegetation code in SeasonalCanopyLAI command",BAD_DATA);
       for (int i=0;i<num_read;i++){
         for (int mon=0;mon<12;mon++){
-          parsed_veg[indices[i]]->relative_LAI[mon]=properties[i][mon];
+          parsed_veg[indices[i]].relative_LAI[mon]=properties[i][mon];
         }
       }
-      DeletePropArray(indices,properties,num_read);
       break;
     }
     case(202):  //----------------------------------------------
@@ -708,10 +723,9 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
                        "ParseClassPropertiesFile: Invalid vegetation code in SeasonalCanopyHeight command",BAD_DATA);
       for (int i=0;i<num_read;i++){
         for (int mon=0;mon<12;mon++){
-          parsed_veg[indices[i]]->relative_ht[mon]=properties[i][mon];
+          parsed_veg[indices[i]].relative_ht[mon]=properties[i][mon];
         }
       }
-      DeletePropArray(indices,properties,num_read);
       break;
     }
     case(203):  //----------------------------------------------
@@ -761,13 +775,11 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
       {
         for (int j=0;j<nParamStrings-1;j++)
         {
-          CVegetationClass::SetVegetationProperty(*parsed_veg   [indices[i]],
+          CVegetationClass::SetVegetationProperty(parsed_veg   [indices[i]],
                                                   aParamStrings[j+1],
                                                   properties   [i][j]);
         }
       }
-      DeletePropArray(indices,properties,num_read);
-
       bool found_in_master;
       for (int j=0;j<nParamStrings-1;j++){
         found_in_master=false;
@@ -910,12 +922,15 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
       break;
     }
     case(501):  //----------------------------------------------
-    {/*
+    {/*ChannelRatingCurves
        :ChannelRatingCurves [optional string name]
          :Bedslope {double slope}
-         :StageRelations
-           {double depth double area double topwidth double discharge [double perim] }x num curve points
-         :EndStageRelations
+           :SurveyPoints
+           {double x double bed_elev}x num survey points
+         :EndSurveyPoints
+         :RoughnessZones
+           {double x_z double mannings_n}xnum roughness zones
+         :EndRoughnessZones
        :EndChannelProfile*/
       string tag;
       double slope=0.0;
@@ -1067,7 +1082,6 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
           CTerrainClass::SetTerrainProperty(parsed_terrs[indices[i]],aParamStrings[j+1],properties[i][j]);
         }
       }
-      DeletePropArray(indices,properties,num_read);
       break;
     }
     case(700):  //----------------------------------------------
@@ -1514,16 +1528,16 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
   pModel->GetGlobalParams()->AutoCalculateGlobalParams(parsed_globals,global_template);
 
   for (int c=1;c<num_parsed_veg;c++){
-    pVegClasses [c]->AutoCalculateVegetationProps  (*parsed_veg[c],  *parsed_veg[0]);
+    pVegClasses [c-1]->AutoCalculateVegetationProps  (parsed_veg[c],parsed_veg[0]);
   }
   for (int c=1;c<num_parsed_soils;c++){
-    pSoilClasses[c]->AutoCalculateSoilProps        (*parsed_soils[c],*parsed_soils[0],pModel->GetTransportModel()->GetNumConstituents());
+    pSoilClasses[c]->AutoCalculateSoilProps          (*parsed_soils[c],*parsed_soils[0],pModel->GetTransportModel()->GetNumConstituents());
   }
   for (int c=1;c<num_parsed_lult;c++) {
     pLUClasses  [c]->AutoCalculateLandUseProps     (*parsed_surf[c], *parsed_surf[0]);
   }
   for (int c=1;c<num_parsed_terrs;c++){
-    pTerrClasses[c-1]->AutoCalculateTerrainProps   ( parsed_terrs[c], parsed_terrs[0]);
+    pTerrClasses[c-1]->AutoCalculateTerrainProps    (parsed_terrs[c],parsed_terrs[0]);
   }
 
   if (!Options.silent){
@@ -1593,8 +1607,8 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
 ///     CLASS_TAGM, v1, v2, v3, v4,... \n
 ///   :EndProperties
 /// \param *p [in] File parsing object
-/// \param *indices [out] Array storing index number of class (with reference to tags[] array), dynamically generated
-/// \param **properties [out] 2D array of model properties, dynamically generated
+/// \param *indices [out] Index number of class (with reference to tags[] array)
+/// \param **properties [out] array of model properties
 /// \param &num_read [out] Number of rows read
 /// \param *tags [in] Array of tag names
 /// \param line_length [in] Number of expected columns in array
@@ -1602,8 +1616,8 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
 /// \return True if operation is successful
 //
 bool ParsePropArray(CParser          *p,           //parser
-                    int             *&indices,     //output: index number of class (with reference to tags[] array)
-                    double         **&properties,  //output: array of properties
+                    int              *indices,     //output: index number of class (with reference to tags[] array)
+                    double          **properties,  //output: array of properties
                     int              &num_read,    //output: number of rows read
                     string           *tags,        //array of tag names
                     const int         line_length, //number of expected columns in array
@@ -1631,21 +1645,10 @@ bool ParsePropArray(CParser          *p,           //parser
         cout <<"bad tag:"<<s[0]<<endl;
         return true;
       }//invalid class index
+      indices[num_read]=index;
 
-      int     *tmpind       =new int     [num_read+1];
-      double **tmpproperties=new double *[num_read+1];
-
-      for (int i=0;i<num_read;i++){
-        tmpind       [i]=indices[i];
-        tmpproperties[i]=new double [line_length];
-        for (int j=0;j<line_length;j++){
-          tmpproperties[i][j]=properties[i][j];
-        }
-      }
-      tmpind       [num_read]=index;
-      tmpproperties[num_read]=new double[line_length];
-      for(int j=1;j<line_length;j++) {
-        tmpproperties[num_read][j-1]=AutoOrDoubleOrAlias(s[j],pAliases,nAliases);
+      for (int j=1;j<line_length;j++){
+        properties[num_read][j-1]=AutoOrDoubleOrAlias(s[j],pAliases,nAliases);
       }
       DeletePropArray(indices,properties,num_read);
 
@@ -1661,13 +1664,6 @@ bool ParsePropArray(CParser          *p,           //parser
   }
   return false;
 }
-void DeletePropArray(int *indices,double **properties,int num_read)
-{
-  if (num_read==0){return;}
-  delete [] indices;
-  for (int i=0;i<num_read;i++){delete [] properties[i];}
-  delete [] properties;
-}
 ///////////////////////////////////////////////////////////////////
 /// \brief Given list of nP needed parameters aP[] of type aPC[], warns user if some are not found associated with classes
 ///
@@ -1802,6 +1798,7 @@ void AddNewSoilClass(CSoilClass **&pSoilClasses,
   soiltags=tmpSoil;
   soiltags[num_parsed_soils] = name;
 
+  //cout<<" SOIL CLASS ADDED: "<<soiltags[num_parsed_soils]<<endl;
   num_parsed_soils++;
 }
 ///////////////////////////////////////////////////////////////////
diff --git a/src/RavenInclude.h b/src/RavenInclude.h
index 1ac7e4e..6364597 100644
--- a/src/RavenInclude.h
+++ b/src/RavenInclude.h
@@ -252,6 +252,8 @@ const double  USE_TEMPLATE_VALUE      =-55555.5;
 const double  NOT_NEEDED              =-66666.6;                                ///< arbitrary value indicating that a non-auto parameter is not needed for the current model configuration
 const double  NOT_NEEDED_AUTO         =-77777.7;                                ///< arbitrary value indicating that a autogeneratable parameter is not needed for the current model configuration
 const double  NETCDF_BLANK_VALUE      =-9999.0;                                 ///< NetCDF flag for blank value
+const int     NETCDF_DEFLATE_LEVEL    = 6;                                      ///< NetCDF deflate level (compression 1-9)
+const int     NETCDF_CHUNKSIZE_MB     = 10;                                     ///< NetCDF chunk size in MB (must be >1 for compression to work)
 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
@@ -335,6 +337,9 @@ const int     MAX_SV_LAYERS       =160;         ///< Max number of layers per st
 const int     MAX_SOILLAYERS      =50;          ///< Max number of soil layers in profile
 const int     MAX_STATE_VAR_TYPES =100;         ///< Max number of *types* of state variables in model
 const int     MAX_STATE_VARS      =500;         ///< Max number of simulated state variables manipulable by one process (CAdvection worst offender)
+const int     MAX_SOIL_PROFILES   =200;         ///< Max number of soil profiles
+const int     MAX_VEG_CLASSES     =200;         ///< Max number of vegetation classes
+const int     MAX_LULT_CLASSES    =200;         ///< Max number of lult classes
 const int     MAX_TERRAIN_CLASSES =50;          ///< Max number of terrain classes
 const int     MAX_SURVEY_PTS      =500;         ///< Max number of survey points
 const int     MAX_CONSTITUENTS    =10;          ///< Max number of transported constituents
@@ -847,6 +852,7 @@ enum toc_method
 
 enum assimtype
 {
+  DA_RAVEN_DEFAULT, ///< multiplicative scaling assimilation upstream propagation
   DA_ECCC           ///< additive assimilation upstream propagation
 };
 ////////////////////////////////////////////////////////////////////
@@ -1064,6 +1070,7 @@ struct optStruct
   double           convergence_crit;          ///< convergence criteria
   double           max_iterations;            ///< maximum number of iterations for iterative solver method
   double           timestep;                  ///< numerical method timestep (in days)
+  size_t           n_out_time;                ///< size of output time dimension
   double           output_interval;           ///< write to output file every x number of timesteps
   ensemble_type    ensemble;                  ///< ensemble type (or ENSEMBLE_NONE if single model)
   string           external_script;           ///< call to external script/.exe once per timestep (or "" if none)
diff --git a/src/Reservoir.cpp b/src/Reservoir.cpp
index aa36466..f6838d3 100644
--- a/src/Reservoir.cpp
+++ b/src/Reservoir.cpp
@@ -93,8 +93,8 @@ void CReservoir::BaseConstructor(const string Name,const long long SBID)
   _assimilate_stage=false;
   _assim_blank=true;
   _pObsStage=NULL;
-  _DAadjust=0.0;
-  _DAadjust_last=0.0;
+  _DAscale=1.0;
+  _DAscale_last=1.0;
 
   _dry_timesteps=0;
 }
@@ -534,7 +534,7 @@ double  CReservoir::GetControlOutflow    (const int i) const { return _aQstruct[
 //
 double  CReservoir::GetIntegratedOutflow(const double& tstep) const
 {
-  return 0.5*(max(_Qout+_DAadjust,0.0)+max(_Qout_last+_DAadjust_last,0.0))*(tstep*SEC_PER_DAY); //integrated
+  return 0.5*(_DAscale*_Qout+_DAscale_last*_Qout_last)*(tstep*SEC_PER_DAY); //integrated
 }
 
 //////////////////////////////////////////////////////////////////
@@ -1149,10 +1149,10 @@ void CReservoir::SetDemandMultiplier(const double &value)
 //////////////////////////////////////////////////////////////////
 /// \brief sets data assimilation scale factors (read from .rvc file)
 //
-void CReservoir::SetDataAssimFactors(const double& da_adj,const double& da_adj_last)
+void CReservoir::SetDataAssimFactors(const double& da_scale,const double& da_scale_last)
 {
-  _DAadjust     =da_adj;
-  _DAadjust_last=da_adj_last;
+  _DAscale     =da_scale;
+  _DAscale_last=da_scale_last;
 }
 
 //////////////////////////////////////////////////////////////////
@@ -1229,19 +1229,39 @@ void  CReservoir::DisableOutflow()
   for (int i=0;i<_Np;i++){_aQ[i]=0.0;_aQunder[i]=0.0;}
 }
 //////////////////////////////////////////////////////////////////
-/// \brief adjust all internal flows by adjustment amount (for assimilation/nudging)
+/// \brief scales all internal flows by scale factor (for assimilation/nudging)
 /// \remark Messes with mass balance something fierce!
 ///
 /// \return mass added to system [m3]
 //
-double CReservoir::AdjustFlow(const double& Qadjust, const bool overriding, const double& tstep, const double& t)
+double CReservoir::ScaleFlow(const double& scale,const bool overriding, const double& tstep, const double &t)
 {
-  _DAadjust_last=_DAadjust;
-  _DAadjust=Qadjust;
+  double va=0.0; //volume added
+  double sf=(scale-1.0)/scale;
+
+  _DAscale=scale;
 
   //Estimate volume added through scaling
+  va+=0.5*(_Qout_last+_Qout)*sf*tstep*SEC_PER_DAY;
+
+  return va;
+}
+double CReservoir::AdjustFlow(const double& Qadjust, const bool overriding, const double& tstep, const double& t)
+{
+
+  _DAscale=1.0;
+  _DAscale_last=1.0;
+
+  double scale=(_Qout+Qadjust)/_Qout;
+  if (_Qout==0.0){scale=1.0;}
+
   double va=0.0; //volume added
-  va+=0.5*(_DAadjust_last+_DAadjust)*tstep*SEC_PER_DAY;
+  double sf=(scale-1.0)/scale;
+
+  _DAscale=scale;
+
+  //Estimate volume added through scaling
+  va+=0.5*(_Qout_last+_Qout)*sf*tstep*SEC_PER_DAY;
 
   return va;
 }
@@ -1276,6 +1296,9 @@ void  CReservoir::UpdateStage(const double &new_stage,const double &res_outflow,
     _aQstruct_last[i]=_aQstruct[i];
     _aQstruct     [i]=aQstruct_new[i];
   }
+
+  _DAscale_last=_DAscale;
+  _DAscale   =1.0;
 }
 //////////////////////////////////////////////////////////////////
 /// \brief returns AET [mm/d], meant to be called at end of time step
@@ -1856,8 +1879,8 @@ void CReservoir::WriteToSolutionFile (ofstream &RVC) const
 {
   RVC<<"    :ResFlow, "<<_Qout<<","<<_Qout_last<<endl;
   RVC<<"    :ResStage, "<<_stage<<","<<_stage_last<<endl;
-  if((_DAadjust_last!=0.0) || (_DAadjust!=0.0)){
-    RVC<<"    :ResDAadj, "<<_DAadjust<<","<<_DAadjust_last<<endl;
+  if(_DAscale_last!=1.0) {
+    RVC<<"    :ResDAscale, "<<_DAscale<<","<<_DAscale_last<<endl;
   }
   for (int i = 0; i < _nControlStructures; i++) {
     RVC<<"    :ControlFlow, "<<i<<","<<_aQstruct[i]<<","<<_aQstruct_last[i]<<endl;
diff --git a/src/Reservoir.h b/src/Reservoir.h
index 26bc365..3576f8c 100644
--- a/src/Reservoir.h
+++ b/src/Reservoir.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
   ----------------------------------------------------------------
   Reservoir.h
   ------------------------------------------------------------------
@@ -107,16 +107,14 @@ class CReservoir
   double       _stage_last;          ///< stage at beginning of current time step [m]
   double       _Qout;                ///< outflow corresponding to current stage [m3/s]
   double       _Qout_last;           ///< outflow at beginning of current time step [m3/s]
-  double      *_aQstruct;            ///< array of flows from control structures at end of time step [m3/s]
-  double      *_aQstruct_last;       ///< array of flows from control structures at start of time step [m3/s]
-  double      *_aQdelivered;         ///< amount of water demand delivered for each demand [m3/s] (for management optimization)
-  double      *_aQreturned;          ///< amount of water returned for each demand [m3/s]
-
-  //diagnostic/flux variables
   double       _MB_losses;           ///< losses over current time step [m3]
   double       _AET;                 ///< losses through AET only [m3]
   double       _Precip;              ///< gains through precipitation [m3]
   double       _GW_seepage;          ///< losses to GW only [m3] (negative for GW gains)
+  double      *_aQstruct;            ///< array of flows from control structures at end of time step [m3/s]
+  double      *_aQstruct_last;       ///< array of flows from control structures at start of time step [m3/s]
+  double      *_aQdelivered;         ///< amount of water demand delivered for each demand [m3/s] (for management optimization)
+  double      *_aQreturned;          ///< amount of water returned for each demand [m3/s]
 
   res_constraint _constraint;        ///< current constraint type
 
@@ -142,14 +140,14 @@ class CReservoir
   int                 _nControlStructures; ///< number of outflow control structures
   CControlStructure **_pControlStructures; ///< Array of pointer to outflow control structures
 
-  //GW parameters :
+  //GW information :
   double       _seepage_const;       ///< seepage constant [m3/s/m] for groundwater losses Q=k*(h-h_loc)
   double       _local_GW_head;       ///< local head [m] (same  for groundwater losses Q=k*(h-h_loc)
 
   void       BaseConstructor(const string Name,const long long SBID); //because some versions of c++ don't like delegating constructors
 
-  double     GetVolume     (const double &ht) const;
-  double     GetArea       (const double &ht) const;
+  double     GetVolume (const double &ht) const;
+  double     GetArea   (const double &ht) const;
   double     GetWeirOutflow(const double &ht, const double &adj) const;
 
   double     GetDZTROutflow(const double &V,const double &Qin,const time_struct &tt,const optStruct &Options) const;
@@ -182,8 +180,8 @@ class CReservoir
 
   double            GetStorage               () const; //[m3]
   double            GetOldStorage            () const; //[m3]
-  double            GetOutflowRate           (const bool adjusted=true) const; //[m3/s]
-  double            GetOldOutflowRate        (const bool adjusted=true) const; //[m3/s]
+  double            GetOutflowRate           () const; //[m3/s]
+  double            GetOldOutflowRate        () const; //[m3/s]
   double            GetIntegratedOutflow     (const double &tstep) const; //[m3]
   double            GetControlOutflow        (const int i) const; //[m3]
   double            GetIntegratedControlOutflow(const int i, const double& tstep) const; //[m3]
@@ -280,6 +278,7 @@ class CReservoir
   void              InitializePostRVM        (const optStruct &Options);
 
   //Called during simulation:
+
   void              UpdateDemands            (const optStruct& Options, const time_struct& tt);
   double            RouteWater               (const double      &Qin_old,
                                               const double      &Qin_new,
@@ -298,6 +297,7 @@ class CReservoir
   void              WriteToSolutionFile      (ofstream &OUT) const;
   void              UpdateReservoir          (const time_struct &tt, const optStruct &Options);
   void              UpdateMassBalance        (const time_struct &tt, const double &tstep, const optStruct &Options);
+  double            ScaleFlow                (const double &scale, const bool overriding,const double &tstep,const double &t);
   double            AdjustFlow               (const double &Qadjust, const bool overriding,const double &tstep,const double &t);
   void              AddToDeliveredDemand     (const int ii, const double &Q);
   void              RecordReturnFlow         (const int ii, const double& Qret);
diff --git a/src/StandardOutput.cpp b/src/StandardOutput.cpp
index 75ee8a6..c0580b4 100644
--- a/src/StandardOutput.cpp
+++ b/src/StandardOutput.cpp
@@ -1321,6 +1321,7 @@ void CModel::WriteMinorOutput(const optStruct &Options,const time_struct &tt)
       for (int p = 0; p < _nSubBasins; p++) {
         if((_pSubBasins[p]->IsGauged()) && (_pSubBasins[p]->IsEnabled())){
           if (Options.assim_method==DA_ECCC){ASSIM<<","<<_aDAQadjust[p]<<" ";}
+          else                              {ASSIM<<","<<_aDAscale  [p]<<" ";}
         }
       }
       ASSIM<<endl;
@@ -2020,6 +2021,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options)
   string      tmpFilename;
   int         p;                                     // loop over all sub-basins
   string      tmp,tmp2,tmp3;
+  size_t      chunksize_time;                        // Size of output time dimension
 
   // initialize all potential file IDs with -9 == "not existing and hence not opened"
   _HYDRO_ncid    = -9;   // output file ID for Hydrographs.nc         (-9 --> not opened)
@@ -2058,7 +2060,14 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options)
 
   /// Define the time variable. Assign units attributes to the netCDF VARIABLES.
   dimids1[0] = time_dimid;
+  // Set chunksize to the number of time steps
+  chunksize_time = Options.n_out_time;
+
   retval = nc_def_var(_HYDRO_ncid, "time", NC_DOUBLE, ndims1,dimids1, &varid_time); HandleNetCDFErrors(retval);
+  // Enable compression and chunking
+  retval = nc_def_var_deflate(_HYDRO_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval);
+  retval = nc_def_var_chunking(_HYDRO_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval);
+
   retval = nc_put_att_text(_HYDRO_ncid, varid_time, "units"   ,      strlen(starttime)  , starttime);   HandleNetCDFErrors(retval);
   retval = nc_put_att_text(_HYDRO_ncid, varid_time, "calendar",      strlen("gregorian"), "gregorian"); HandleNetCDFErrors(retval);
   retval = nc_put_att_text(_HYDRO_ncid, varid_time, "standard_name", strlen("time"),      "time");      HandleNetCDFErrors(retval);
@@ -2080,7 +2089,7 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options)
     // (b) create dimension "nbasins"
     retval = nc_def_dim(_HYDRO_ncid, "nbasins", nSim, &nbasins_dimid);                             HandleNetCDFErrors(retval);
 
-    // (c) create variable  and set attributes for"basin_name"
+    // (c) create variable  and set attributes for "basin_name"
     dimids1[0] = nbasins_dimid;
     retval = nc_def_var(_HYDRO_ncid, "basin_name", NC_STRING, ndims1, dimids1, &varid_bsim);       HandleNetCDFErrors(retval);
     tmp ="ID of sub-basins with simulated outflows";
@@ -2145,6 +2154,10 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options)
     /// Define the time variable. Assign units attributes to the netCDF VARIABLES.
     dimids1[0] = time_dimid;
     retval = nc_def_var     (_RESSTAGE_ncid,"time",NC_DOUBLE,ndims1,dimids1,&varid_time);           HandleNetCDFErrors(retval);
+    // Enable compression and chunking
+    retval = nc_def_var_deflate(_RESSTAGE_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval);
+    retval = nc_def_var_chunking(_RESSTAGE_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval);
+
     retval = nc_put_att_text(_RESSTAGE_ncid,varid_time,"units",strlen(starttime),starttime);        HandleNetCDFErrors(retval);
     retval = nc_put_att_text(_RESSTAGE_ncid,varid_time,"calendar",strlen("gregorian"),"gregorian"); HandleNetCDFErrors(retval);
     retval = nc_put_att_text(_RESSTAGE_ncid,varid_time,"standard_name",strlen("time"),"time");      HandleNetCDFErrors(retval);
@@ -2224,9 +2237,13 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options)
     // Define the DIMENSIONS. NetCDF will hand back an ID
     retval = nc_def_dim(_STORAGE_ncid, "time", NC_UNLIMITED, &time_dimid);  HandleNetCDFErrors(retval);
 
-    /// Define the time variable.
+     /// Define the time variable.
     dimids1[0] = time_dimid;
     retval = nc_def_var(_STORAGE_ncid, "time", NC_DOUBLE, ndims1,dimids1, &varid_time); HandleNetCDFErrors(retval);
+    // Enable compression and chunking
+    retval = nc_def_var_deflate(_STORAGE_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval);
+    retval = nc_def_var_chunking(_STORAGE_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval);
+
     retval = nc_put_att_text(_STORAGE_ncid, varid_time, "units"   , strlen(starttime)  , starttime);   HandleNetCDFErrors(retval);
     retval = nc_put_att_text(_STORAGE_ncid, varid_time, "calendar", strlen("gregorian"), "gregorian"); HandleNetCDFErrors(retval);
 
@@ -2273,8 +2290,13 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options)
     // time vector
     // ----------------------------------------------------------
     retval = nc_def_dim(_FORCINGS_ncid,"time",NC_UNLIMITED,&time_dimid);  HandleNetCDFErrors(retval);
+
     dimids1[0] = time_dimid;
     retval = nc_def_var(_FORCINGS_ncid,"time",NC_DOUBLE,ndims1,dimids1,&varid_time); HandleNetCDFErrors(retval);
+    // Enable compression and chunking
+    retval = nc_def_var_deflate(_FORCINGS_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval);
+    retval = nc_def_var_chunking(_FORCINGS_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval);
+
     retval = nc_put_att_text(_FORCINGS_ncid,varid_time,"units",strlen(starttime),starttime);   HandleNetCDFErrors(retval);
     retval = nc_put_att_text(_FORCINGS_ncid,varid_time,"calendar",strlen("gregorian"),"gregorian"); HandleNetCDFErrors(retval);
 
@@ -2321,8 +2343,15 @@ void CModel::WriteNetcdfStandardHeaders(const optStruct &Options)
     // time vector
     // ----------------------------------------------------------
     retval = nc_def_dim(_RESMB_ncid,"time",NC_UNLIMITED,&time_dimid);  HandleNetCDFErrors(retval);
+
+
     dimids1[0] = time_dimid;
     retval = nc_def_var(_RESMB_ncid,"time",NC_DOUBLE,ndims1,dimids1,&varid_time);                HandleNetCDFErrors(retval);
+
+     // Enable compression and chunking
+    retval = nc_def_var_deflate(_RESMB_ncid, varid_time, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval);
+    retval = nc_def_var_chunking(_RESMB_ncid, varid_time, NC_CHUNKED, &chunksize_time); HandleNetCDFErrors(retval);
+
     retval = nc_put_att_text(_RESMB_ncid,varid_time,"units",strlen(starttime),starttime);        HandleNetCDFErrors(retval);
     retval = nc_put_att_text(_RESMB_ncid,varid_time,"calendar",strlen("gregorian"),"gregorian"); HandleNetCDFErrors(retval);
 
@@ -3017,6 +3046,7 @@ int NetCDFAddMetadata(const int fileid,const int time_dimid,string shortname,str
   int retval;
   int dimids[1];
   dimids[0] = time_dimid;
+  size_t chunksize[1];
 
   static double fill_val[] = {NETCDF_BLANK_VALUE};
   static double miss_val[] = {NETCDF_BLANK_VALUE};
@@ -3024,6 +3054,13 @@ int NetCDFAddMetadata(const int fileid,const int time_dimid,string shortname,str
   // (a) create variable precipitation
   retval = nc_def_var(fileid,shortname.c_str(),NC_DOUBLE,1,dimids,&varid); HandleNetCDFErrors(retval);
 
+  // Set compression
+  retval = nc_def_var_deflate(fileid, varid, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval);
+  // Get the chunksize of the time dimension
+  retval = nc_inq_var_chunking(fileid, time_dimid, NULL, &chunksize[0]); HandleNetCDFErrors(retval);
+  // Set chunksize to number of time steps, or 1 if time dimension is unlimited
+  retval = nc_def_var_chunking(fileid, varid, NC_CHUNKED, chunksize); HandleNetCDFErrors(retval);
+
   // (b) add attributes to variable
   retval = nc_put_att_text  (fileid,varid,"units",units.length(),units.c_str());              HandleNetCDFErrors(retval);
   retval = nc_put_att_text  (fileid,varid,"long_name",longname.length(),longname.c_str());    HandleNetCDFErrors(retval);
@@ -3048,6 +3085,7 @@ int NetCDFAddMetadata2D(const int fileid,const int time_dimid,int nbasins_dimid,
   int    retval;
   int    dimids2[2];
   string tmp;
+  size_t chunksize2[2];
 
   static double fill_val[] = {NETCDF_BLANK_VALUE};
   static double miss_val[] = {NETCDF_BLANK_VALUE};
@@ -3057,6 +3095,17 @@ int NetCDFAddMetadata2D(const int fileid,const int time_dimid,int nbasins_dimid,
 
   // (a) create variable
   retval = nc_def_var(fileid,shortname.c_str(),NC_DOUBLE,2,dimids2,&varid); HandleNetCDFErrors(retval);
+  // Set compression
+  retval = nc_def_var_deflate(fileid, varid, 1, 1, NETCDF_DEFLATE_LEVEL); HandleNetCDFErrors(retval);
+  // Set time chunksize to number of time steps
+  retval = nc_inq_var_chunking(fileid, time_dimid, NULL, &chunksize2[0]); HandleNetCDFErrors(retval);
+   // Set nbasins chunksize to number ensuring that chunks have approximately 10 MB of data (assuming double precision)
+  retval = nc_inq_dimlen(fileid, nbasins_dimid, &chunksize2[1]); HandleNetCDFErrors(retval);
+
+
+  chunksize2[1] = max((size_t)1, min(chunksize2[1], (size_t)(NETCDF_CHUNKSIZE_MB * 1024 * 1024 / sizeof(double) / chunksize2[0]))); // Ensure at least one basin per chunk
+  // Set chunksize
+  retval = nc_def_var_chunking(fileid, varid, NC_CHUNKED, chunksize2); HandleNetCDFErrors(retval);
 
   tmp = "basin_name";
 
diff --git a/src/SubBasin.cpp b/src/SubBasin.cpp
index 22a849d..7eb44ae 100644
--- a/src/SubBasin.cpp
+++ b/src/SubBasin.cpp
@@ -57,7 +57,6 @@ CSubBasin::CSubBasin( const long long      Identifier,
   _GW_exch_coeff     =0.0;
   _bed_conductivity  =0.0;
   _bed_thickness     =0.5; //m
-  _mean_slope        =0.0;
 
   _t_conc            =AUTO_COMPUTE;
   _t_peak            =AUTO_COMPUTE;
@@ -1500,6 +1499,55 @@ void CSubBasin::SetUnusableFlowPercentage(const double &val)
     "CSubBasin:SetUnusableFlowPercentage: invalid value for :UnusableFlowPercentage (must be between zero and one).",BAD_DATA_WARN);
   _unusable_flow_pct=val;
 }
+
+//////////////////////////////////////////////////////////////////
+/// \brief scales all internal flows by scale factor (for assimilation/nudging)
+/// \remark Messes with mass balance something fierce!
+/// \param &scale [in] Flow scaling factor
+/// \param &overriding [in] True if Qlast should be scaled (overriding); false for no-data scaling
+/// \param &tstep [in] time step [d]
+/// \param &t [in] - current model time
+/// \return volume added to system [m3]
+//
+double CSubBasin::ScaleAllFlows(const double &scale, const bool overriding, const double &tstep, const double &t)
+{
+  double va=0.0; //volume added [m3]
+  double sf=(scale-1.0)/scale;
+
+  if(!overriding)
+  {
+    for(int n=0;n<_nQlatHist;n++) {
+      _aQlatHist[n]*=scale;  upperswap(_aQlatHist[n],0.0);
+      va+=_aQlatHist[n]*sf*tstep*SEC_PER_DAY;
+    }
+    for(int n=0;n<_nQinHist; n++) {
+      _aQinHist[n]*=scale; upperswap(_aQinHist[n],0.0);
+      va+=_aQinHist[n]*sf*tstep*SEC_PER_DAY;
+    }
+    for(int i=0;i<_nSegments;i++) {
+      _aQout[i]*=scale; upperswap(_aQout[i],0.0);
+      va+=0.5*(_aQout[i]+_aQout[i+1])*sf*tstep*SEC_PER_DAY;
+    }
+  }
+
+  if((overriding) && (_pReservoir==NULL)) {
+    for (int i=0;i<_nSegments;i++)
+    {
+      _aQout[i]*=scale;  upperswap(_aQout[i],0.0);
+      va+=0.5*(_aQout[i]+_aQout[i+1])*sf*tstep*SEC_PER_DAY;
+    }
+  }
+
+  if(_pReservoir!=NULL) {
+    va+=_pReservoir->ScaleFlow(scale,overriding,tstep,t);
+  }
+
+  //Estivate volume added through scaling
+  _channel_storage*=scale; va+=_channel_storage*sf;
+  _rivulet_storage*=scale; va+=_rivulet_storage*sf;
+  return va;
+}
+
 //////////////////////////////////////////////////////////////////
 /// \brief adjusts all internal flows by corresponding magnitude (for assimilation/nudging)
 /// \remark Messes with mass balance something fierce!
@@ -1540,7 +1588,7 @@ double CSubBasin::AdjustAllFlows(const double &adjust, const bool overriding,con
       va+=adjust*tstep*SEC_PER_DAY;
     }
   }
-  else if((overriding) && (assimsite)) //just update Qout
+  else if((overriding) && (assimsite)) //just update Qin
   {
     for (int i=0;i<_nSegments;i++)
     {
diff --git a/src/SubBasin.h b/src/SubBasin.h
index d51e81b..48115fd 100644
--- a/src/SubBasin.h
+++ b/src/SubBasin.h
@@ -317,6 +317,7 @@ class CSubBasin
   void            SetGauged                (const bool isgauged);
   void            Disable                  ();
   void            Enable                   ();
+  double          ScaleAllFlows            (const double &scale_factor, const bool scale_last, const double &tstep, const double &t);
   double          AdjustAllFlows           (const double &adjustment, const bool overriding, const bool assimsite, const double &tstep, const double &t);
   void            SetUnusableFlowPercentage(const double &val);
   void            IncludeInAssimilation    ();