V506_lake_upgrade

.gitignore

@@ -42,5 +42,8 @@ Makefile
 Raven
 Raven_errors.txt
 
+# Directories
+test/
+
 # JetBrains
 .idea

src/EnergyTransport.h

@@ -60,8 +60,8 @@ class CEnthalpyModel :public CConstituentModel
   double GetWaterTemperature     (const double *state_vars, const int iWater) const;
 
   double GetEnergyLossesInTransit(const int p,double &Q_sens,double &Q_GW) const;
+  double GetEnergyLossesFromLake (const int p,double &Q_sens,double &Q_cond,double &Q_lat,double &Q_rad_in,double &Q_lw_in, double &Q_lw_out, double &Q_rain, double &Q_adv, double &kdiff) const;
   double GetEnergyLossesFromReach(const int p,double &Q_sens,double &Q_cond,double &Q_lat,double &Q_GW,double &Q_rad_in,double &Q_lw_in, double &Q_lw_out,double &Q_lateral, double &Q_fric, double &Tave) const;
-  double GetEnergyLossesFromLake (const int p,double &Q_sens,double &Q_cond,double &Q_lat,double &Q_rad_in,double &Q_lw_in, double &Q_lw_out, double &Q_rain) const;
 
   double GetOutflowIceFraction   (const int p) const;
   double GetAvgLatentHeatFlux    () const;

src/GlobalParams.cpp

@@ -156,6 +156,31 @@ void CGlobalParams::AutoCalculateGlobalParams(const global_struct &Gtmp, const g
   {
     G.GAMMA_SCALE_multiplier=1.0;//default=nothing (no warning needed)
   }
+
+  autocalc=SetCalculableValue(G.mix_depth_coef,Gtmp.mix_depth_coef,Gtemplate.mix_depth_coef);
+  if (autocalc)
+  {
+    G.mix_depth_coef=0.25;
+    warn="The required global parameter MIX_DEPTH_COEF was autogenerated with value "+to_string(G.mix_depth_coef);
+    if (chatty){WriteAdvisory(warn,false);}
+  }
+
+  autocalc=SetCalculableValue(G.mix_depth_expn,Gtmp.mix_depth_expn,Gtemplate.mix_depth_expn);
+  if (autocalc)
+  {
+    G.mix_depth_expn=0.45;
+    warn="The required global parameter MIX_DEPTH_EXPN was autogenerated with value "+to_string(G.mix_depth_expn);
+    if (chatty){WriteAdvisory(warn,false);}
+  }
+
+  autocalc=SetCalculableValue(G.lake_lyr_conv_coeff_max,Gtmp.lake_lyr_conv_coeff_max,Gtemplate.lake_lyr_conv_coeff_max);
+  if (autocalc)
+  {
+    G.lake_lyr_conv_coeff_max=2.5;
+    warn="The required global parameter KDIFF_MAX was autogenerated with value "+to_string(G.lake_lyr_conv_coeff_max);
+    if (chatty){WriteAdvisory(warn,false);}
+  }
+
   autocalc=SetCalculableValue(G.reference_flow_mult,Gtmp.reference_flow_mult,Gtemplate.reference_flow_mult);
   if (autocalc)
   {
@@ -417,6 +442,11 @@ void CGlobalParams::InitializeGlobalParameters(global_struct &g, bool is_templat
   g.TIME_TO_PEAK_multiplier =DefaultParameterValue(is_template,true);
   g.GAMMA_SHAPE_multiplier  =DefaultParameterValue(is_template,true);
   g.GAMMA_SCALE_multiplier  =DefaultParameterValue(is_template,true);
+
+  g.mix_depth_coef          =DefaultParameterValue(is_template,true);
+  g.mix_depth_expn          =DefaultParameterValue(is_template,true);
+  g.lake_lyr_conv_coeff_max =DefaultParameterValue(is_template,true);
+
 
   //model-specific parameters
   g.avg_annual_snow     =DefaultParameterValue(is_template,false);
@@ -564,6 +594,11 @@ void  CGlobalParams::SetGlobalProperty (global_struct &G,
   else if (!name.compare("HBVEC_LAPSE_RATE"    )){G.HBVEC_lapse_rate=value; }
   else if (!name.compare("HBVEC_LAPSE_UPPER"   )){G.HBVEC_lapse_upper=value; }
   else if (!name.compare("HBVEC_LAPSE_ELEV"    )){G.HBVEC_lapse_elev=value; }
+
+  else if (!name.compare("MIX_DEPTH_COEF"      )){G.mix_depth_coef=value;}
+  else if (!name.compare("MIX_DEPTH_EXPN"      )){G.mix_depth_expn=value;}
+  else if (!name.compare("KDIFF_MAX"           )){G.lake_lyr_conv_coeff_max=value;}
+
   else{
     WriteWarning("CGlobalParams::SetGlobalProperty: Unrecognized/invalid global parameter name ("+name+") in .rvp file",false);
 
@@ -687,6 +722,11 @@ double CGlobalParams::GetGlobalProperty(const global_struct &G, string  param_na
   else if (!name.compare("HBVEC_LAPSE_RATE"     )){return G.HBVEC_lapse_rate; }
   else if (!name.compare("HBVEC_LAPSE_UPPER"    )){return G.HBVEC_lapse_upper; }
   else if (!name.compare("HBVEC_LAPSE_ELEV"     )){return G.HBVEC_lapse_elev; }
+
+  else if (!name.compare("MIX_DEPTH_COEF"       )){return G.mix_depth_coef;}
+  else if (!name.compare("MIX_DEPTH_EXPN"       )){return G.mix_depth_expn;}
+  else if (!name.compare("KDIFF_MAX"            )){return G.lake_lyr_conv_coeff_max;}
+
   else{
     if (strict){
       string msg="CGlobalParams::GetParameter: Unrecognized/invalid global parameter name in .rvp file: "+name;

src/Makefile

@@ -25,9 +25,9 @@ LDFLAGS  :=
 # OPTION 0) some compilers require the c++11 flag, some may not
 CXXFLAGS += -std=c++11 -fPIC
 
-# OPTION 1) include netcdf - uncomment following two commands (assumes netCDF path = /usr/local):
-#CXXFLAGS += -Dnetcdf
-#LDLIBS   += -L/usr/local -lnetcdf
+# OPTION 1) include netcdf - uncomment following two commands (assumes netCDF path = /usr/local): 
+CXXFLAGS += -Dnetcdf     
+LDLIBS   += -L/usr/local -lnetcdf  
 
 # OPTION 1b) include netcdf - for newer MacOS with Apple Silicon (use with option 1 also uncommented)
 #CXXFLAGS += -I/opt/homebrew/include

src/ParseHRUFile.cpp

@@ -425,7 +425,7 @@ bool ParseHRUPropsFile(CModel *&pModel, const optStruct &Options, bool terrain_r
                 //in*=CGlobalParams::GetParameter("TOC_MULTIPLIER");    // use to be this; but has its own multiplier now; maybe set default to TOC_MULTIPLIER??
                 in *= pModel->GetGlobalParams()->GetParameter("TIME_TO_PEAK_MULTIPLIER");
               }
-	            if(!aParamStrings[i].compare("GAMMA_SHAPE") && (in!=AUTO_COMPUTE) && (in!=USE_TEMPLATE_VALUE)){
+              if(!aParamStrings[i].compare("GAMMA_SHAPE") && (in!=AUTO_COMPUTE) && (in!=USE_TEMPLATE_VALUE)){
                 in *= pModel->GetGlobalParams()->GetParameter("GAMMA_SHAPE_MULTIPLIER");
               }
               if(!aParamStrings[i].compare("GAMMA_SCALE") && (in!=AUTO_COMPUTE) && (in!=USE_TEMPLATE_VALUE)){
@@ -1915,15 +1915,15 @@ CReservoir *ReservoirParse(CParser *p,string name,const CModel *pModel,long long
   //------------------------------------------------------------------------------------
   if((type==CURVE_POWERLAW) || (type==CURVE_LINEAR))
   {
-    pRes=new CReservoir(name,SBID,a_V,b_V,a_Q,b_Q,a_A,b_A,crestht,max_depth);
+    pRes=new CReservoir(name,SBID,a_V,b_V,a_Q,b_Q,a_A,b_A,crestht,max_depth,pModel);
   }
   else if(type==CURVE_DATA)
   {
-    pRes=new CReservoir(name,SBID,aQ_ht,aQ,aQund,aA,aV,NQ);//presumes aQ_ht=aV_ht=aA_ht; NA=NV=NQ
+    pRes=new CReservoir(name,SBID,aQ_ht,aQ,aQund,aA,aV,NQ,pModel);//presumes aQ_ht=aV_ht=aA_ht; NA=NV=NQ
   }
   else if(type==CURVE_VARYING)
   {
-    pRes=new CReservoir(name,SBID,nDates,aDates,aQ_ht,aQQ,aQund,aA,aV,NQ);//presumes aQ_ht=aV_ht=aA_ht; NA=NV=NQ
+    pRes=new CReservoir(name,SBID,nDates,aDates,aQ_ht,aQQ,aQund,aA,aV,NQ,pModel);//presumes aQ_ht=aV_ht=aA_ht; NA=NV=NQ
   }
   else if(type==CURVE_LAKE)
   {
@@ -1940,7 +1940,7 @@ CReservoir *ReservoirParse(CParser *p,string name,const CModel *pModel,long long
      }
     }
 
-    pRes=new CReservoir(name,SBID,weircoeff,cwidth,crestht,lakearea,max_depth);
+    pRes=new CReservoir(name,SBID,weircoeff,cwidth,crestht,lakearea,max_depth,pModel);
 
   }
   else {

src/ParsePropertyFile.cpp

@@ -167,6 +167,15 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
 
   aP [0]="GAMMA_SCALE_MULTIPLIER"; aPC[0]=CLASS_GLOBAL;
   AddToMasterParamList(aPmaster, aPCmaster,nPmaster, aP, aPC, 1);
+
+  aP [0]="MIX_DEPTH_COEF"; aPC[0]=CLASS_GLOBAL;
+  AddToMasterParamList(aPmaster, aPCmaster,nPmaster, aP, aPC, 1);
+
+  aP [0]="MIX_DEPTH_EXPN"; aPC[0]=CLASS_GLOBAL;
+  AddToMasterParamList(aPmaster, aPCmaster,nPmaster, aP, aPC, 1);
+
+  aP [0]="KDIFF_MAX"; aPC[0]=CLASS_GLOBAL;
+  AddToMasterParamList(aPmaster, aPCmaster,nPmaster, aP, aPC, 1);
 
   //Throw warning if NULL Terrain but terrain parameter is needed
   terrain_required=false;

src/Properties.h

@@ -396,6 +396,7 @@ struct global_struct
 
   double           avg_annual_snow;          ///< [mm] avg annual snow as SWE
   double           avg_annual_runoff;        ///< [mm] avg annual runoff from basin
+
   double           reference_flow_mult;      ///< [-] multiplier for reference flow relative to annual mean flow
                                              ///< (10 for flood modelling, 0.5 for low-flow conditions)
 
@@ -426,6 +427,9 @@ struct global_struct
   double           TIME_TO_PEAK_multiplier;  ///< [0..1+] time to peak multiplier
   double           GAMMA_SHAPE_multiplier;   ///< [0..1+] Gamma shape multiplier
   double           GAMMA_SCALE_multiplier;   ///< [0..1+] Gamma scale multiplier
+  double           mix_depth_coef;           ///< [m] Reservoir mixing depth coefficient
+  double           lake_lyr_conv_coeff_max;  ///< [MJ/m2/d/K] Upper limit for convection coefficient between lake layers
+  double           mix_depth_expn;           ///< [-] Reservoir mixing depth exponent
   double           assimilation_fact;        ///< [0..1] assimilation factor (0=no assimilation to 1= full replacement)
   double           assim_upstream_decay;     ///< [1/km] assimilation upstream decay factor (0= overrides everything upstream, large- observation influence decays quickly with distance from gauge) [~0.01]
   double           assim_time_decay;         ///< [1/d] assimilation temporal decay factor (0=diminishes in future, 0.1 - diminshes in  3 days) [0.2]

src/RavenInclude.h

@@ -181,7 +181,7 @@ const double  TC_CLAY                 =0.216;
 const double  TC_ORGANIC              =0.0216;                                  ///< [MJ/m/d/K] Thermal conductivity of organic matter (0.25 W/m/K)
 const double  TC_AIR                  =0.00199;                                 ///< [MJ/m/d/K] Thermal conductivity of air (0.023 W/m/K)
 
-const double  COM_WATER		            =4.58e-10;                                ///< [1/Pa] Compressiblity of Water
+const double  COM_WATER               =4.58e-10;                                ///< [1/Pa] Compressiblity of Water
 const double  COM_ICE                 =4.58e-10;                                ///< [1/Pa] Compressiblity of Ice
 const double  HCP_WATER               =4.186;                                   ///< [MJ/m3/K] Volumetric Heat Capacity of Water
 const double  HCP_ICE                 =1.938;                                   ///< [MJ/m3/K] Volumetric Heat Capacity of Ice