diff --git a/crhmcode/Makefile b/crhmcode/Makefile
index ba09182e9..23014d6e9 100644
--- a/crhmcode/Makefile
+++ b/crhmcode/Makefile
@@ -53,8 +53,8 @@ ClassalbedoBaker.o Classalbedoobs.o Classalbedoobs2.o Classalbedoparam.o \
ClassalbedoRichard.o ClassalbedoWinstral.o ClassAnnan.o ClassAyers.o \
ClassBasin.o Classbrushintcp.o Classcalcsun.o Classcontribution.o \
Classcrack.o ClassCRHMCanopy.o ClassCRHMCanopyClearing.o \
-ClassCRHMCanopyClearingGap.o Classebsm.o Classevap.o Classevap_Resist.o \
-ClassevapD.o ClassevapD_Resist.o ClassevapX.o ClassFlowInSnow.o \
+ClassCRHMCanopyClearingGap.o ClassCRHMCanopyVectorBased.o Classebsm.o Classevap.o \
+Classevap_Resist.o ClassevapD.o ClassevapD_Resist.o ClassevapX.o ClassFlowInSnow.o \
Classfrostdepth.o Classfrozen.o ClassfrozenAyers.o Classglacier.o \
Classglacier_debris.o ClassGlobal.o ClassGreenAmpt.o ClassGreencrack.o \
ClassGrow_Crop.o ClassHMSA.o ClassHtobs.o ClassIceBulb.o ClassICEflow.o \
diff --git a/crhmcode/src/modules/CMakeLists.txt b/crhmcode/src/modules/CMakeLists.txt
index e9a5ed534..c28860b3f 100644
--- a/crhmcode/src/modules/CMakeLists.txt
+++ b/crhmcode/src/modules/CMakeLists.txt
@@ -17,11 +17,14 @@ add_library(modules OBJECT
ClassBasin.cpp
Classbrushintcp.cpp
Classcalcsun.cpp
+ ClassCanopySnowBalanceBase.cpp
+ ClassCanopySnowBalanceCRHM.cpp
Classcontribution.cpp
Classcrack.cpp
ClassCRHMCanopy.cpp
ClassCRHMCanopyClearing.cpp
ClassCRHMCanopyClearingGap.cpp
+ ClassCRHMCanopyVectorBased.cpp
Classebsm.cpp
Classevap_Resist.cpp
Classevap.cpp
diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp
new file mode 100644
index 000000000..40134b363
--- /dev/null
+++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp
@@ -0,0 +1,434 @@
+/**
+ * Copyright 2022, CRHMcode's Authors or Contributors
+ * This file is part of CRHMcode.
+ *
+ * CRHMcode is free software: you can redistribute it and/or modify it under
+ * the terms of the GNU General Public License as published by the Free Software
+ * Foundation, either version 3 of the License, or (at your option) any later
+ * version.
+ *
+ * CRHMcode is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ * See the GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License along with
+ * CRHMcode. If not, see .
+ *
+ **/
+// created by Manishankar Mondal
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include "ClassCRHMCanopyVectorBased.h"
+#include "../core/GlobalDll.h"
+#include "../core/ClassCRHM.h"
+#include "newmodules/SnobalDefines.h"
+
+using namespace CRHM;
+
+ClassCRHMCanopyVectorBased* ClassCRHMCanopyVectorBased::klone(string name) const{
+ return new ClassCRHMCanopyVectorBased(name);
+}
+
+void ClassCRHMCanopyVectorBased::decl(void)
+{
+
+ Description = "'Calculates initial snow interception after Cebulski & Pomeroy (2025, Hyrological Proc.) canopy snow ablation is handled by the CanopySnowBalance module. Assumptions are for this to be applied on a forested HRU without significant clear cuts/large gaps in the forest with width many times the canopy height. This module also does some radiation calcs for consistency with the CanopyClearingGap module which this one replaces.' \
+ 'Inputs are observations Qsi (W/m^2) and Qli (W/m^2), ' \
+ 'Inputs are observation Qsi (W/m^2) and variable QliVt_Var (W/m^2), ' \
+ 'Inputs are variable QsiS_Var (W/m^2)(slope) from Annandale and observation Qli (W/m^2), ' \
+ 'Inputs are variable QsiS_Var (W/m^2)(slope) from Annandale and variable QliVt_Var (W/m^2), ' \
+ 'Inputs are variable QsiA_Var (W/m^2)(horizontal) from Annandale and variable QliVt_Var (W/m^2).'";
+
+ // Observations
+
+ variation_set = VARIATION_0 + VARIATION_1;
+
+ declreadobs("Qsi", TDim::NHRU, "incident short-wave", "(W/m^2)", &Qsi, HRU_OBS_Q);
+
+ variation_set = VARIATION_0 + VARIATION_2;
+
+ declreadobs("Qli", TDim::NHRU, "incident long-wave", "(W/m^2)", &Qli, HRU_OBS_Q);
+
+ variation_set = VARIATION_1 + VARIATION_3 + VARIATION_4;
+
+ declgetvar("*", "QliVt_Var", "(W/m^2)", &QliVt_Var);
+
+ variation_set = VARIATION_2 + VARIATION_3;
+
+ declgetvar("*", "QsiS_Var", "(W/m^2)", &QsiS_Var);
+
+ variation_set = VARIATION_4;
+
+ declgetvar("*", "QsiA_Var", "(W/m^2)", &QsiA_Var);
+
+ variation_set = VARIATION_ORG;
+
+ // get variables:
+
+ declgetvar("*", "hru_t", "(" + string(DEGREE_CELSIUS) + ")", &hru_t);
+
+ declgetvar("*", "hru_u", "(m/s)", &hru_u);
+
+ declgetvar("*", "hru_rh", "()", &hru_rh);
+
+ declgetvar("*", "hru_snow", "(mm/int)", &hru_snow);
+
+ declgetvar("*", "hru_rain", "(mm/int)", &hru_rain);
+
+ declgetvar("*", "SolAng", "(r)", &SolAng);
+
+ declgetvar("*", "cosxs", "(r)", &cosxs);
+
+ declgetvar("*", "cosxsflat", "(r)", &cosxsflat);
+
+ declgetvar("*", "Qdfo", "(W/m^2)", &Qdfo);
+
+ declgetvar("*", "Albedo", "()", &Albedo);
+
+ declgetvar("*", "QdflatE", "(W/m^2)", &QdflatE);
+
+ declputvar("*", "SWE", "(mm)", &SWE);
+
+ // declputvar("*", "TCanSnow", "(" + string(DEGREE_CELSIUS) + ")", &TCanSnow); found better performance when just using the ice bulb temp so no need for pspNew module
+
+ // declared observations
+
+ declvar("Ts", TDim::NHRU, "snow surface temperature", "(" + string(DEGREE_CELSIUS) + ")", &Ts);
+
+ declobs("Qnsn", TDim::NHRU, "net all-wave at snow surface", "(W/m^2)", &Qnsn);
+
+ declvar("Qnsn_Var", TDim::NHRU, "net all-wave at snow surface", "(W/m^2)", &Qnsn_Var);
+
+ declobs("Qsisn", TDim::NHRU, "incident short-wave at surface", "(W/m^2)", &Qsisn);
+
+ declvar("Qsisn_Var", TDim::NHRU, "incident short-wave at surface", "(W/m^2)", &Qsisn_Var);
+
+ declobs("Qlisn", TDim::NHRU, "incident long-wave at surface", "(W/m^2)", &Qlisn);
+
+ declvar("Qlisn_Var", TDim::NHRU, "incident long-wave at surface", "(W/m^2)", &Qlisn_Var);
+
+ declobs("Qlosn", TDim::NHRU, "reflected long-wave at surface", "(W/m^2)", &Qlosn);
+
+ // declared variables
+
+ decldiag("k", TDim::NHRU, "extinction coefficient", "()", &k);
+
+ decldiag("Tauc", TDim::NHRU, "short-wave transmissity", "", &Tauc);
+
+ decllocal("Pa", TDim::NHRU, "Average surface pressure", "(kPa)", &Pa);
+
+ declvar("ra", TDim::NHRU, "resistance", "(s/m)", &ra);
+
+ declvar("throughfall_rain", TDim::NHRU, "throughfall of rain through canopy (not in contact with canopy).", "(mm/int)", &throughfall_rain);
+
+ declvar("throughfall_snow", TDim::NHRU, "throughfall of snow not in contact with the canopy", "(mm/int)", &throughfall_snow);
+
+ declvar("intercepted_snow", TDim::NHRU, "initial interception of snow before ablation kicks in.", "(mm/int)", &intercepted_snow);
+
+ declvar("intercepted_rain", TDim::NHRU, "initial interception of rain before ablation kicks in.", "(mm/int)", &intercepted_rain);
+
+ decldiag("u_FHt", TDim::NHRU, "wind speed at forest top (z = FHt)", "(m/s)", &u_FHt);
+
+ decldiag("u_1_third_Ht", TDim::NHRU, "wind speed at one-third forest height (z = 1/3*Ht) for the Cebulski & Pomeroy vector based param.", "(m/s)", &u_1_third_Ht);
+
+ declvar("Clca", TDim::NHRU, "Leaf contact area adjusted for hydrometeor trajectory angle.", "()", &Clca);
+
+ declvar("LAI_", TDim::NHRU, "Leaf area index adjusted for snow cover", "()", &LAI_);
+
+ // parameters:
+
+ declparam("hru_elev", TDim::NHRU, "[637]", "0.0", "100000.0", "altitude", "(m)", &hru_elev);
+
+ declparam("Surrounding_Ht", TDim::NHRU, "[0.1, 0.25, 1.0]", "0.001", "100.0", "surrounding canopy height", "(m)", &Surrounding_Ht);
+
+ declparam("Gap_diameter", TDim::NHRU, "[100]", "10", "1000", "representative gap diameter", "(m)", &Gap_diameter);
+
+ declparam("Ht", TDim::NHRU, "[20.0]", "0.001", "100.0", "forest/vegetation height", "(m)", &Ht);
+
+ declparam("Zref", TDim::NHRU, "[1.5]", "0.01", "100.0", "temperature measurement height", "(m)", &Zref);
+
+ declparam("Zwind", TDim::NHRU, "[10]", "0.01", "100.0", "wind measurement height", "(m)", &Zwind);
+
+ declparam("Z0snow", TDim::NHRU, "[0.01]", "0.0001", "0.01", "snow roughness length", "(m)", &Z0snow);
+
+ declparam("Cc", TDim::NHRU, "[1]", "0", "1", "canopy coverage, (1-sky view fraction)", "()", &Cc);
+
+ declparam("LAI", TDim::NHRU, "[2.2]", "0.1", "20.0", "leaf-area-index", "()", &LAI);
+
+ declparam("alpha", TDim::NHRU, "[0.836]", "0.0", "1.0", "$alpha$ is an efficiency constant which determines the fraction of snowflakes that contact the $C_p$ elements and are stored in the canopy (i.e., intercepted) before canopy snow unloading or ablation processes begin. Default is based on Cebulski & Pomeroy plot scale observations at Fortress Mountain PWL and Forest Tower plots.", "()", &alpha);
+
+ declparam("Zvent", TDim::NHRU, "[0.75]", "0.0", "1.0", "ventilation wind speed height (z/Ht)", "()", &Zvent);
+
+ declparam("CanopyClearing", TDim::NHRU, "[0]", "0", "2", "canopy - 0/clearing - 1/gap - 2", "()", &CanopyClearing);
+
+ declparam("CanopyWindSwitchIP", TDim::NHRU, "[0]", "0", "1", "Canopy wind model to use at height Zcan, 0 - for Cionco, 1 - for Prandtl-von Kármán log-linear relationship", "()", &CanopyWindSwitchIP);
+
+ decldiagparam("Alpha_c", TDim::NHRU, "[0.1]", "0.05", "0.2", "canopy albedo, used for longwave-radiation enhancement estimation", "()", &Alpha_c);
+
+ decldiagparam("B_canopy", TDim::NHRU, "[0.038]", "0.0", "0.2", "canopy enhancement parameter for longwave-radiation. Suggestions are Colorado - 0.023 and Alberta - 0.038", "()", &B_canopy);
+}
+
+void ClassCRHMCanopyVectorBased::init(void)
+{
+
+ nhru = getdim(TDim::NHRU); // transfers current # of HRU's to module
+
+ for (hh = 0; hh < nhru; ++hh)
+ {
+
+ Pa[hh] = 101.3f * pow((293.0f - 0.0065f * hru_elev[hh]) / 293.0f, 5.26f); // kPa
+
+ if (Ht[hh] > Zwind[hh])
+ {
+ CRHMException TExcept(string("'" + Name + " (CanopyVectorBased)' Vegetation height greater than wind reference height, i.e. (Ht > Zwind)!").c_str(), TExcept::WARNING);
+ LogError(TExcept);
+ }
+ }
+}
+
+void ClassCRHMCanopyVectorBased::run(void)
+{
+ for (hh = 0; chkStruct(); ++hh)
+ {
+ throughfall_rain[hh] = 0.0;
+ throughfall_snow[hh] = 0.0;
+
+ intercepted_rain[hh] = 0.0;
+ intercepted_snow[hh] = 0.0;
+
+ choose_solar_data(); // decide what solar data to use based on module variation defined in prj file
+ compute_canopy_temp(); // compute canopy temperature as input for surface snow energy balance
+
+ adjust_lai_for_snow_height(Ht[hh]);
+ compute_canopy_long_short_rads(); // compute canopy long and short wave radiation components sent to canopy snow ebal and surface snowbal modules
+
+ //==============================================================================
+ // coupled forest snow interception routine:
+ // after Cebulski & Pomeroy 2025:
+
+ if (hru_snow[hh] > 0.0)
+ { // calculate increase in Snow_load and throughfall_snow if we are in canopy (i.e., Cc > 0)
+ const double v_snow = 0.8; // terminal fall velocity of snowfall taken from Isyumov, 1971
+ Clca[hh] = 0.0; // leaf contact area (Clca) based on trajectory angle
+ double IP = 0.0; // interception efficiency (IP)
+
+ if (hru_u[hh] > 0.0 && Cc[hh] < 1.0 && Cc[hh] > 0.0)
+ { // increase leaf contact area (Clca) based on wind speed and canopy coverage (Cc)
+
+ get_wind_for_leaf_contact_calc();
+ compute_leaf_contact_area(v_snow);
+
+ } else { // do not adjust snow-leaf contact area for wind effects
+ Clca[hh] = Cc[hh]; // use leaf contact area from nadir i.e., Clca == 1 for Cc == 1 and Clca == 0 when Cc == 0
+ }
+
+ IP = Clca[hh] * alpha[hh]; // interception efficiency (IP)
+ intercepted_snow[hh] = IP * hru_snow[hh]; // change in canopy snow load
+
+ // calculate canopy snow throughfall before unloading:
+
+ throughfall_snow[hh] = (1.0 - IP) * hru_snow[hh];
+
+ // net snow / rain is computed in can snobal after ablation of intercepted snow.
+ } // end snow routine
+
+ if(hru_rain[hh] > 0){
+ // Canopy Liquid Water Routine
+ throughfall_rain[hh] = hru_rain[hh] * (1.0 - Cc[hh]);
+ intercepted_rain[hh] = hru_rain[hh] * Cc[hh];
+ }
+
+ } // end for
+}
+
+// helper functions
+
+// this function chooses which solar data to use based on the module variation defined in the prj file
+void ClassCRHMCanopyVectorBased::choose_solar_data(void)
+{
+ switch ((int)variation)
+ {
+ case VARIATION_ORG:
+ Qsi_ = Qsi[hh];
+ Qli_ = Qli[hh];
+ break;
+ case VARIATION_1:
+ Qsi_ = Qsi[hh];
+ Qli_ = QliVt_Var[hh];
+ break;
+ case VARIATION_2:
+ Qsi_ = QsiS_Var[hh];
+ Qli_ = Qli[hh];
+ break;
+ case VARIATION_3:
+ Qsi_ = QsiS_Var[hh];
+ Qli_ = QliVt_Var[hh];
+ break;
+ case VARIATION_4:
+ Qsi_ = QsiA_Var[hh];
+ Qli_ = QliVt_Var[hh];
+ break;
+ default:
+ Qsi_ = 0.0;
+ Qli_ = 0.0;
+ break;
+ }
+}
+
+// this function computes canopy temperature based on air temperature and radiation balance used in the Snobal module
+void ClassCRHMCanopyVectorBased::compute_canopy_temp(void)
+{
+ // Canopy temperature is approximated by the air temperature
+ double T1 = hru_t[hh] + CRHM_constants::Tm; // deg C -> K
+
+ double rho = Pa[hh] * 1000 / (CRHM_constants::Rgas * T1); // air density
+
+ double U1 = hru_u[hh]; // wind speed (m/s)
+ ra[hh] = (log(Zref[hh] / Z0snow[hh]) * log(Zwind[hh] / Z0snow[hh])) / sqr(CRHM_constants::kappa) / U1;
+
+ double deltaX = 0.622 * CRHM_constants::Ls * Common::Qs(Pa[hh], T1) / (CRHM_constants::Rgas * sqr(T1));
+
+ double q = (hru_rh[hh] / 100) * Common::Qs(Pa[hh], T1); // specific humidity
+
+ Ts[hh] = T1 + (CRHM_constants::emiss * (Qli_ - CRHM_constants::sbc * pow(T1, 4.0))
+ + CRHM_constants::Ls * (q - Common::Qs(Pa[hh], T1)) * rho / ra[hh])
+ / (4 * CRHM_constants::emiss * CRHM_constants::sbc * pow(T1, 3.0)
+ + (CRHM_constants::Cp + CRHM_constants::Ls * deltaX) * rho / ra[hh]);
+
+ Ts[hh] -= CRHM_constants::Tm; // back to deg C
+
+ // Apply boundary conditions
+ if (Ts[hh] > 0.0 || SWE[hh] <= 0.0)
+ Ts[hh] = 0.0;
+}
+
+// this function adjusts the effective LAI based on snow height as the bottom of the canopy may be buried in snow
+void ClassCRHMCanopyVectorBased::adjust_lai_for_snow_height(double Ht)
+{
+ // Compute canopy exposure above snow
+ double exposure = Ht - Common::DepthofSnow(SWE[hh]); // Ht in meters, SWE in mm
+ if (exposure < 0.0)
+ exposure = 0.0;
+
+ // Update effective LAI
+ LAI_[hh] = LAI[hh] * exposure / Ht;
+}
+
+// this function computes canopy longwave and shortwave radiation components sent to canopy snow ebal and surface snowbal modules
+void ClassCRHMCanopyVectorBased::compute_canopy_long_short_rads(void)
+{
+ double T1 = hru_t[hh] + CRHM_constants::Tm; // deg C -> K
+ double Vf = 1.0 - Cc[hh]; // changed from approx of LAI to direct calculation from params
+
+ double Vf_ = Vf + (1.0 - Vf) * sin((Ht[hh] - (Ht[hh] - Common::DepthofSnow(SWE[hh]))) / Ht[hh] * M_PI_2);
+
+ if (SolAng[hh] > 0.001 && cosxs[hh] > 0.001 && cosxsflat[hh] > 0.001)
+ {
+ k[hh] = 1.081 * SolAng[hh] * cos(SolAng[hh]) / sin(SolAng[hh]); // Eq. 8 from Pomeroy 2009
+ double limit = cosxsflat[hh] / cosxs[hh];
+ if (limit > 2.0)
+ limit = 2.0;
+ Tauc[hh] = exp(-k[hh] * LAI_[hh] * limit);
+ }
+ else
+ {
+ k[hh] = 0.0;
+ }
+
+ double Kstar_H = Qsi_ * (1.0 - Alpha_c[hh] - Tauc[hh] * (1.0 - Albedo[hh])); // Eq. 6 from Pomeroy 2009
+
+ // incoming longwave to the surface
+ Qlisn[hh] = Qli_ * Vf_ + (1.0f - Vf_) * CRHM_constants::emiss_c * CRHM_constants::sbc * pow(T1, 4.0f) + B_canopy[hh] * Kstar_H; // looks like modification of Eq. 10 from Pomeroy 2009
+
+ Qlisn_Var[hh] = Qlisn[hh]; // goes to snobal module
+
+ // incoming shortwave radiation to the surface
+ Qsisn[hh] = Qsi_ * Tauc[hh];
+
+ Qsisn_Var[hh] = Qsisn[hh];
+
+ // incoming longwave to the canopy mid point
+
+ // Qlw_veg_Var[hh] = Qli_ * (1.0 - (Cc[hh]/2)); // canopy snow is partially exposed to atmosphere so reduce based on view factor
+
+ // incoming shortwave radiation to canopy midpoint
+
+ // Qsw_veg_Var[hh] = Qsi_ * Tauc_50; // for PSP like solar to canopy midpoint
+
+ Qlosn[hh] = CRHM_constants::emiss * CRHM_constants::sbc * pow(Ts[hh] + CRHM_constants::Tm, 4.0f);
+
+ Qnsn[hh] = Qlisn[hh] - Qlosn[hh] + Qsisn[hh] * (1.0 - Albedo[hh]);
+
+ Qnsn_Var[hh] = Qnsn[hh];
+}
+
+// this function computes wind speed, variations are defined in the prj file using the CanopyWindSwitchIP parameter
+void ClassCRHMCanopyVectorBased::get_wind_for_leaf_contact_calc(void)
+{
+ double Ht_1_third = Ht[hh] * (1.0 / 3.0);
+
+ // wind speed used for vector based initial snow interception
+ switch ((int)CanopyWindSwitchIP[hh])
+ {
+ case 0:
+ { // increase wind speed to canopy top then bring back down using using Cionco wind model for dense canopies, use this one if wind speed is below/within canopy
+ adst_wind_cpy_top(Ht[hh], hru_u[hh], Zwind[hh], u_FHt[hh]);
+ cionco_canopy_wind_spd(Ht[hh], u_FHt[hh], Ht_1_third, u_1_third_Ht[hh]);
+ break;
+ } // end case 0
+
+ case 1:
+ { // use if forcing wind speed is above/at canopy top
+ cionco_canopy_wind_spd(Ht[hh], hru_u[hh], Ht_1_third, u_1_third_Ht[hh]);
+ break;
+ } // end case 1
+
+ case 2:
+ { // Canopy wind profile developed at Fortress sparse canopy
+ double d0 = 0.5791121; // displacement height observed at sparse forest around Fortress Forest Tower
+ double z0m = 0.4995565; // roughness length observed at above site
+
+ // wind speed used for vector based initial snow interception
+ if ((Ht_1_third - d0) > z0m){
+ double Ustar = hru_u[hh]*PBSM_constants::KARMAN/(log((Zwind[hh]-d0)/z0m));
+ u_1_third_Ht[hh] = Ustar/PBSM_constants::KARMAN * log((Ht_1_third - d0)/z0m);
+ } else {
+ u_1_third_Ht[hh] = 0.0;
+ }
+ break;
+ } // end case 2
+ } // end of switch CanopyWind
+}
+
+// this function computes leaf contact area based on hydrometeor trajectory angle as in Cebulski & Pomeroy 2025 (hp)
+void ClassCRHMCanopyVectorBased::compute_leaf_contact_area(double v_snow)
+{
+ if (u_1_third_Ht[hh] > 0.0)
+ {
+ double snow_traj_angle = atan(u_1_third_Ht[hh] / v_snow); // trajectory angle in radians
+ double b_lca = 0.91; // empirical coefficient
+ double f_theta = b_lca * pow(sin(snow_traj_angle), 2.0);
+
+ double Cp_inc = (1.0 - Cc[hh]) * f_theta; // fractional increase in leaf contact area
+ Clca[hh] = Cc[hh] + Cp_inc;
+ }
+ else
+ {
+ Clca[hh] = Cc[hh]; // nadir contact area
+ }
+}
+
+void ClassCRHMCanopyVectorBased::finish(bool good)
+{
+ for (hh = 0; chkStruct(); ++hh)
+ {
+ LogDebug(" ");
+ }
+}
\ No newline at end of file
diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h
new file mode 100644
index 000000000..afcb99657
--- /dev/null
+++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h
@@ -0,0 +1,124 @@
+/**
+* Copyright 2022, CRHMcode's Authors or Contributors
+* This file is part of CRHMcode.
+*
+* CRHMcode is free software: you can redistribute it and/or modify it under
+* the terms of the GNU General Public License as published by the Free Software
+* Foundation, either version 3 of the License, or (at your option) any later
+* version.
+*
+* CRHMcode is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty
+* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+* See the GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License along with
+* CRHMcode. If not, see .
+*
+**/
+//created by Manishankar Mondal
+
+#include "../core/ClassModule.h"
+
+class ClassCRHMCanopyVectorBased:public ClassModule {
+
+public:
+
+ClassCRHMCanopyVectorBased(string Name, string Version = "undefined", LMODULE Lvl = LMODULE::PROTO) : ClassModule(Name, Version, Lvl, 1000, " , QliVt_Var, QsiS_Var, QsiS_Var, QsiA_Var") {};
+
+double Qsi_{ 0.0 }, Qli_ {0.0};
+
+// observation variables
+
+const double *Qsi { NULL };
+const double *Qli { NULL };
+
+// declared observation variables
+
+double *Ts { NULL };
+double *Qnsn { NULL };
+double *Qsisn { NULL };
+double *Qlisn { NULL };
+double *Qlosn { NULL };
+
+// put variables
+
+double *SWE { NULL };
+
+// declared variables
+
+double *throughfall_rain { NULL };
+double *throughfall_snow { NULL };
+double *intercepted_snow { NULL };
+double *intercepted_rain { NULL };
+double *Qnsn_Var { NULL };
+double *Qsisn_Var { NULL };
+double *Qlisn_Var { NULL };
+double *Qsw_veg_Var { NULL };
+double *Qlw_veg_Var { NULL };
+
+double *Clca { NULL };
+double *k { NULL };
+double *Tauc { NULL };
+double *Pa { NULL };
+double *ra { NULL };
+double *u_FHt { NULL };
+double *u_1_third_Ht { NULL };
+
+double *LAI_ { NULL };
+
+
+// variable inputs
+
+const double *hru_t { NULL };
+const double *hru_u { NULL };
+const double *hru_rh { NULL };
+const double *Albedo { NULL };
+const double *QdflatE { NULL };
+
+const double *hru_snow { NULL };
+const double *hru_rain { NULL };
+
+const double *SolAng { NULL };
+const double *cosxs { NULL };
+const double *cosxsflat { NULL };
+const double *Qdfo { NULL };
+
+const double *QsiS_Var { NULL };
+const double *QsiA_Var { NULL };
+const double *QliVt_Var { NULL };
+
+// declared parameters:
+
+const double *hru_elev { NULL };
+const double *Ht { NULL };
+const double *Cc { NULL };
+const double *LAI { NULL };
+const double *alpha { NULL };
+const double *Z0snow { NULL };
+const double *Zref { NULL };
+const double *Zwind { NULL };
+const double *Zvent { NULL };
+const double *Surrounding_Ht { NULL };
+const double *Gap_diameter { NULL };
+const double *Alpha_c { NULL };
+const double *B_canopy { NULL };
+
+const long *CanopyClearing { NULL };
+const long *CanopyWindSwitchIP { NULL };
+
+void decl(void);
+void init(void);
+void run(void);
+void finish(bool good);
+
+ClassCRHMCanopyVectorBased* klone(string name) const;
+
+private:
+ void choose_solar_data(void);
+ void compute_canopy_temp(void);
+ void adjust_lai_for_snow_height(double Ht);
+ void compute_canopy_long_short_rads(void);
+ void compute_leaf_contact_area(double v_snow);
+ void get_wind_for_leaf_contact_calc(void);
+};
diff --git a/crhmcode/src/modules/ClassCanopySnowBalanceBase.cpp b/crhmcode/src/modules/ClassCanopySnowBalanceBase.cpp
new file mode 100644
index 000000000..779eae396
--- /dev/null
+++ b/crhmcode/src/modules/ClassCanopySnowBalanceBase.cpp
@@ -0,0 +1,1475 @@
+/**
+* Copyright 2022, CRHMcode's Authors or Contributors
+* This file is part of CRHMcode.
+*
+* CRHMcode is free software: you can redistribute it and/or modify it under
+* the terms of the GNU General Public License as published by the Free Software
+* Foundation, either version 3 of the License, or (at your option) any later
+* version.
+*
+* CRHMcode is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty
+* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+* See the GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License along with
+* CRHMcode. If not, see .
+*
+**/
+//created by Alex Cebulski
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include "ClassCanopySnowBalanceBase.h"
+#include "../core/GlobalDll.h"
+#include "../core/ClassCRHM.h"
+#include "newmodules/SnobalDefines.h"
+
+
+using namespace CRHM;
+
+void ClassCanopySnowBalanceBase::init(void) {
+
+ input_rec1 = new INPUT_REC[nhru]; // input data for start of data timestep [nhru]
+ input_rec2 = new INPUT_REC[nhru]; // input data for start of data timestep [nhru]
+ time_step = new long[nhru]; // length current timestep (sec)
+ current_time = new long[nhru]; // start time of current time step (sec)
+
+ computed = new int* [nhru]; // array of flags for each timestep;
+ input_deltas = new INPUT_REC * [nhru]; // deltas for climate-input parameters for each timestep
+ precip_info = new PRECIP_REC * [nhru]; // array of precip info adjusted for each timestep
+ tstep_info = new TSTEP_REC * [nhru]; // array of info for each timestep:
+
+ for (long ii = 0; ii < nhru; ++ii) {
+ computed[ii] = new int[4]; // [nhru] [4]
+ input_deltas[ii] = new INPUT_REC[4]; // [nhru] [4]
+ precip_info[ii] = new PRECIP_REC[4]; // [nhru] [4]
+ tstep_info[ii] = new TSTEP_REC[4]; // [nhru] [4]
+ }
+
+ for (hh = 0; chkStruct(); ++hh) {
+ P_a[hh] = 101.3f * pow((293.0f - 0.0065f * hru_elev[hh]) / 293.0f, 5.26f) * 1000.0f; // Pa
+
+ h2o_sat_veg[hh] = 0.0;
+
+ snow_h2o_veg[hh] = 0.0;
+ liq_h2o_veg[hh] = 0.0;
+ m_s_veg[hh] = snow_h2o_veg[hh] + liq_h2o_veg[hh];
+
+ cc_s_veg[hh] = 0.0;
+
+ delmelt_veg[hh] = 0.0;
+ deldrip_veg[hh] = 0.0;
+ T_sf[hh] = 0.0;
+ T_rain_veg[hh] = 0.0;
+
+ isothermal[hh] = 0;
+ vegsnowcover[hh] = 0;
+ stop_no_snow[hh] = 0;
+ if (hh < nhru)
+ {
+ current_time[hh] = 0;
+ }
+
+
+ E_s_cum[hh] = 0.0;
+
+ m_precip_cum[hh] = 0.0;
+ m_rain_cum[hh] = 0.0;
+ m_snow_cum[hh] = 0.0;
+ cmlmelt_veg[hh] = 0.0;
+
+ Ql_veg[hh] = 0.0;
+ qsub_veg[hh] = 0.0;
+ Qh_veg[hh] = 0.0;
+ Qp[hh] = 0.0;
+ Qn_veg[hh] = 0.0;
+ delta_Q_veg[hh] = 0.0;
+
+ if (hh < nhru)
+ {
+ tstep_info[hh][DATA_TSTEP].level = DATA_TSTEP;
+ tstep_info[hh][DATA_TSTEP].time_step = 24 * 3600 / Global::Freq;
+ tstep_info[hh][DATA_TSTEP].intervals = 0;
+ tstep_info[hh][DATA_TSTEP].threshold = DEFAULT_NORMAL_THRESHOLD;
+
+ tstep_info[hh][NORMAL_TSTEP].level = NORMAL_TSTEP;
+ tstep_info[hh][NORMAL_TSTEP].time_step = 24 * 3600 / Global::Freq;
+ tstep_info[hh][NORMAL_TSTEP].intervals = 1;
+ tstep_info[hh][NORMAL_TSTEP].threshold = DEFAULT_NORMAL_THRESHOLD; // 60
+
+ tstep_info[hh][MEDIUM_TSTEP].level = MEDIUM_TSTEP;
+ tstep_info[hh][MEDIUM_TSTEP].time_step = 24 * 3600 / Global::Freq / 4;
+ tstep_info[hh][MEDIUM_TSTEP].intervals = 4;
+ tstep_info[hh][MEDIUM_TSTEP].threshold = DEFAULT_MEDIUM_THRESHOLD; // 10
+
+ tstep_info[hh][SMALL_TSTEP].level = SMALL_TSTEP;
+ tstep_info[hh][SMALL_TSTEP].time_step = 24 * 3600 / Global::Freq / 60;
+ tstep_info[hh][SMALL_TSTEP].intervals = 15;
+ tstep_info[hh][SMALL_TSTEP].threshold = DEFAULT_SMALL_THRESHOLD; // 1
+ }
+
+ }
+}
+
+void ClassCanopySnowBalanceBase::finish(bool good) { // only required for local storage and final output
+
+ for (hh = 0; chkStruct(); ++hh) {
+ LogMessageA(hh, string("'" + Name + " (CanopySnowBalance)' m_precip_cum (mm) (mm*hru) (mm*hru/basin): ").c_str(), m_precip_cum[hh], hru_area[hh], basin_area[0]);
+ LogMessageA(hh, string("'" + Name + " (CanopySnowBalance)' m_rain_cum (mm) (mm*hru) (mm*hru/basin): ").c_str(), m_rain_cum[hh], hru_area[hh], basin_area[0]);
+ LogMessageA(hh, string("'" + Name + " (CanopySnowBalance)' m_snow_cum (mm) (mm*hru) (mm*hru/basin): ").c_str(), m_snow_cum[hh], hru_area[hh], basin_area[0]);
+ LogMessageA(hh, string("'" + Name + " (CanopySnowBalance)' SWE (mm) (mm*hru) (mm*hru/basin): ").c_str(), m_s_veg[hh], hru_area[hh], basin_area[0]);
+ LogDebug(" ");
+ LogMessageA(hh, string("'" + Name + " (CanopySnowBalance)' E_s_cum (mm) (mm*hru) (mm*hru/basin): ").c_str(), E_s_cum[hh], hru_area[hh], basin_area[0]);
+ LogMessageA(hh, string("'" + Name + " (CanopySnowBalance)' cmlmelt_veg (mm) (mm*hru) (mm*hru/basin): ").c_str(), cmlmelt_veg[hh], hru_area[hh], basin_area[0]);
+ LogDebug(" ");
+ }
+
+ delete[] input_rec1; // [nhru]
+ delete[] input_rec2; // [nhru]
+ delete[] time_step; // [nhru]
+ delete[] current_time; // [nhru]
+
+ for (long ii = 0; ii < nhru; ++ii) {
+ delete[] computed[ii]; // [nhru] [4]
+ delete[] input_deltas[ii]; // [nhru] [4]
+ delete[] precip_info[ii]; // [nhru] [4]
+ delete[] tstep_info[ii]; // [nhru] [4]
+ }
+
+ delete[] computed;
+ delete[] input_deltas;
+ delete[] precip_info;
+ delete[] tstep_info;
+}
+
+// Initialize canopy snowpack on this timestep, i.e., no snow load on timestep start
+void ClassCanopySnowBalanceBase::init_snow_veg(void)
+{
+
+ liq_h2o_veg[hh] += m_rain[hh];
+ snow_h2o_veg[hh] += m_snow[hh];
+ m_s_veg[hh] = liq_h2o_veg[hh] + snow_h2o_veg[hh];
+
+ T_s_veg[hh] = T_sf[hh]; // init canopy snow temperature at precip temp
+
+ cc_s_veg[hh] = _cold_content_veg(T_s_veg[hh], m_s_veg[hh]);
+
+ stop_no_snow[hh] = 0;
+ vegsnowcover[hh] = 1;
+}
+
+/*
+** NAME
+** _cold_content_veg -- calculates cold content for canopy snow
+**
+** double
+** _cold_content_veg(
+** double temp, |* temperature of canopy snow *|
+** double mass) |* specific mass of canopy snow *|
+**
+** DESCRIPTION
+** This routine calculates the cold content for canopy snow (i.e., the
+** energy required to bring its temperature to freezing) from the
+** temperature and specific mass.
+**
+** RETURN VALUE
+** The canopy snow's cold content.
+*/
+
+double ClassCanopySnowBalanceBase::_cold_content_veg(
+ double temp, // temperature of canopy snow
+ double mass) // specific mass of canopy snow
+{
+ if (temp < FREEZE)
+ return heat_stor(CP_ICE(temp), mass, (temp - FREEZE));
+ else
+ return 0.0;
+}
+
+/*
+** NAME
+** int do_data_tstep_veg(void) -- run model for 1 data timestep between 2 input records
+**
+** DESCRIPTION
+** This routine performs the model's calculations for 1 data timestep
+** between 2 input-data records which are in 'input_rec1' and
+** 'input_rec2'.
+**
+** If there's precipitation during the data timestep, the flag
+** 'precip_now_veg' used be true. Furthermore, the routine requires
+** that the following precipitation variables have been initialized:
+**
+** m_pp
+** percent_snow
+** T_pp
+**
+** This routine divides the data timestep into the appropriate number
+** of normal run timesteps. The input values for each normal timestep
+** are computed from the two input records by linear interpolation.
+**
+** RETURN VALUE
+**
+** true The model's calculations were completed.
+**
+** false An error occured, and a message explaining the error has
+** been stored with the 'usrerr' routine.
+**
+** GLOBAL VARIABLES READ
+** e_a
+** I_lw
+** in_rec
+** m_pp_data
+** m_rain_data
+** m_snow_data
+** more_pr_recs
+** precip_data
+** S_n
+** T_a
+** tstep_info
+** u
+** z_snow_data
+**
+** GLOBAL VARIABLES MODIFIED
+** precip_now_veg
+*/
+
+void ClassCanopySnowBalanceBase::do_data_tstep_veg(void) {
+
+ int level; // loop index
+
+ // Copy values from first input record into global variables.
+
+ S_n[hh] = input_rec1[hh].S_n;
+ I_lw[hh] = input_rec1[hh].I_lw;
+ T_a[hh] = input_rec1[hh].T_a;
+ e_a[hh] = input_rec1[hh].e_a;
+ u[hh] = input_rec1[hh].u;
+
+ // Compute deltas for the climate input parameters over the CRHM data timestep.
+
+ input_deltas[hh][DATA_TSTEP].S_n = input_rec2[hh].S_n - input_rec1[hh].S_n;
+ input_deltas[hh][DATA_TSTEP].I_lw = input_rec2[hh].I_lw - input_rec1[hh].I_lw;
+ input_deltas[hh][DATA_TSTEP].T_a = input_rec2[hh].T_a - input_rec1[hh].T_a;
+ input_deltas[hh][DATA_TSTEP].e_a = input_rec2[hh].e_a - input_rec1[hh].e_a;
+ input_deltas[hh][DATA_TSTEP].u = input_rec2[hh].u - input_rec1[hh].u;
+
+ // If there is precipitation, then compute the amount of rain & snow in it.
+
+ if (precip_now_veg[hh]) {
+ precip_info[hh]->m_pp = m_precip[hh];
+ precip_info[hh]->m_snow = m_snow[hh];
+ precip_info[hh]->m_rain = m_rain[hh];
+
+ // Mixed snow and rain
+
+ if ((precip_info[hh]->m_snow > 0.0) && (precip_info[hh]->m_rain > 0.0)) {
+ T_sf[hh] = FREEZE;
+ h2o_sat_veg_snow[hh] = 1.0;
+ T_rain_veg[hh] = T_pp[hh]; // TODO try changing to harder ice bulb temp
+ }
+
+ // Snow only
+
+ else if (precip_info[hh]->m_snow > 0.0) {
+ if (T_pp[hh] < FREEZE) { /* Cold snow */
+ T_sf[hh] = T_pp[hh];
+ h2o_sat_veg_snow[hh] = 0.0;
+ }
+ else { /* Warm snow */
+ T_sf[hh] = FREEZE;
+ h2o_sat_veg_snow[hh] = 1.0;
+ }
+ }
+
+ // Rain only
+
+ else if (precip_info[hh]->m_rain > 0.0) {
+ T_rain_veg[hh] = T_pp[hh];
+ }
+ }
+ else {
+ precip_info[hh]->m_pp = 0.0;
+ precip_info[hh]->m_snow = 0.0;
+ precip_info[hh]->m_rain = 0.0;
+ precip_info[hh]->z_snow = 0.0;
+ }
+
+ // Clear the 'computed' flag at the other timestep levels.
+
+ for (level = NORMAL_TSTEP; level <= SMALL_TSTEP; ++level)
+ computed[hh][level] = 0;
+
+ // Divide the data timestep into normal run timesteps.
+
+ _divide_tstep_veg(tstep_info[hh]);
+}
+
+/*
+** NAME
+** int _divide_tstep_veg( -- divide a timestep into smaller timesteps
+**
+** int
+** _divide_tstep_veg(
+** TSTEP_REC *tstep; |* record of timestep to be divided *|
+**
+** DESCRIPTION
+** This routine divides a timestep into smaller timesteps. For
+** each of these smaller timestep, the routine either run the
+** model for that timestep, or further subdivides that timestep
+** into even smaller timesteps.
+**
+** RETURN VALUE
+**
+** true The timestep was successfully divided into smaller timesteps.
+**
+** false An error occured while running the model during one of the
+** smaller timesteps. An message explaining the error has
+** been stored with the 'usrerr' routine.
+**
+** GLOBAL VARIABLES READ
+** precip_now_veg
+** tstep_info
+**
+** GLOBAL VARIABLES MODIFIED
+*/
+
+int ClassCanopySnowBalanceBase::_divide_tstep_veg(TSTEP_REC* tstep) { // record of timestep to be divided
+
+ TSTEP_REC* next_lvl_tstep; // info of next level of timestep
+ INPUT_REC* curr_lvl_deltas; // -> input-deltas of current level
+ INPUT_REC* next_lvl_deltas; // -> input-deltas of next level
+ PRECIP_REC* curr_lvl_precip; // -> precip data of current level
+ PRECIP_REC* next_lvl_precip; // -> precip data of next level
+ int next_level;
+ // Fetch the record for the timestep at the next level.
+
+ next_level = tstep->level + 1;
+ next_lvl_tstep = &tstep_info[hh][next_level];
+
+ curr_lvl_deltas = &input_deltas[hh][tstep->level];
+ next_lvl_deltas = &input_deltas[hh][next_level];
+
+ curr_lvl_precip = &precip_info[hh][tstep->level];
+ next_lvl_precip = &precip_info[hh][next_level];
+
+ // If this is the first time this new level has been used during
+ // the current data timestep, then calculate its input deltas and precipitation values.
+
+ if (!computed[hh][next_level]) {
+ next_lvl_deltas->S_n = curr_lvl_deltas->S_n / next_lvl_tstep->intervals;
+ next_lvl_deltas->I_lw = curr_lvl_deltas->I_lw / next_lvl_tstep->intervals;
+ next_lvl_deltas->T_a = curr_lvl_deltas->T_a / next_lvl_tstep->intervals;
+ next_lvl_deltas->e_a = curr_lvl_deltas->e_a / next_lvl_tstep->intervals;
+ next_lvl_deltas->u = curr_lvl_deltas->u / next_lvl_tstep->intervals;
+
+ next_lvl_precip->m_pp = curr_lvl_precip->m_pp / next_lvl_tstep->intervals;
+ next_lvl_precip->m_rain = curr_lvl_precip->m_rain / next_lvl_tstep->intervals;
+ next_lvl_precip->m_snow = curr_lvl_precip->m_snow / next_lvl_tstep->intervals;
+ next_lvl_precip->z_snow = curr_lvl_precip->z_snow / next_lvl_tstep->intervals;
+
+ computed[hh][next_level] = 1;
+ }
+
+ // For each the new smaller timestep, either subdivide them if below their mass threshold, or run the model for them.
+
+ for (int ii = 0; (ii < next_lvl_tstep->intervals) && !stop_no_snow[hh]; ++ii) { // execute every interval
+
+ if ((next_level < SMALL_TSTEP) && _below_thold_veg(next_lvl_tstep->threshold)) {
+ if (!_divide_tstep_veg(next_lvl_tstep)) // call self (this routine)
+ return 0; // cannot get here unless _hl_e is coded to return error!
+ }
+ else {
+ if (!_do_tstep_veg(next_lvl_tstep)) // execute code
+ return 0; // cannot get here unless _hl_e is coded to return error!
+ }
+ } // for
+
+ return 1;
+} // _divide_tstep_veg
+
+/*
+** NAME
+** _below_thold_veg -- is a layer's mass below a threshold ?
+**
+** int
+** _below_thold_veg(
+** double threshold) |* current timestep's threshold for a
+** layer's mass *|
+**
+** DESCRIPTION
+** This routine determines if any individual layer's mass is below
+** a given threshold for the current timestep.
+**
+** RETURN VALUE
+** 1 A layer's mass is less than the threshold.
+**
+** 0 All layers' masses are greater than the threshold.
+**
+** GLOBAL VARIABLES READ
+** m_s_veg
+** m_s_veg
+**
+** GLOBAL VARIABLES MODIFIED
+*/
+
+int ClassCanopySnowBalanceBase::_below_thold_veg(double threshold) { // current timestep's threshold for a layer's mass
+
+ if (!(vegsnowcover[hh] || liq_h2o_veg[hh] > 0.0))
+ return 0;
+ else
+ return (m_s_veg[hh] < threshold);
+}
+
+/*
+** NAME
+** _do_tstep_veg -- do model calculations for 1 timestep
+**
+** int
+** _do_tstep_veg(
+** TSTEP_REC *tstep; |* timestep's record *|
+**
+** DESCRIPTION
+** This routine performs the model's calculations for a single timestep.
+** It requires that these climate variables have been initialized:
+**
+** S_n
+** I_lw
+** T_a
+** e_a
+** u
+**
+** The routine also requires the precipitation data have been adjusted
+** for the timestep, and have been stored in the array:
+**
+** precip_info
+**
+** RETURN VALUE
+**
+** true The model's calculations were completed.
+**
+** false An error occured, and a message explaining the error has
+** been stored with the 'usrerr' routine.
+**
+** GLOBAL VARIABLES READ
+** delta_Q_veg
+** delsub_veg
+** G
+** Qh_veg
+** Ql_veg
+** Qp
+** delmelt_veg
+** precip_now_veg
+** Qn_veg
+** deldrip_veg
+**
+** GLOBAL VARIABLES MODIFIED
+** current_time
+** curr_time_hrs
+** delta_Q_bar
+** e_a
+** G_bar
+** H_bar
+** liq_h2o_veg
+** I_lw
+** L_v_E_bar
+** M_bar
+** m_precip
+** m_rain
+
+** m_snow
+** R_n_bar
+** S_n
+** vegsnowcover
+** T_a
+** u
+** time_step
+** z_snow_veg
+*/
+
+int ClassCanopySnowBalanceBase::_do_tstep_veg(TSTEP_REC* tstep) // timestep's record
+{
+ time_step[hh] = tstep->time_step;
+
+ m_precip[hh] = precip_info[hh][tstep->level].m_pp;
+ m_rain[hh] = precip_info[hh][tstep->level].m_rain;
+ m_snow[hh] = precip_info[hh][tstep->level].m_snow;
+
+ delmelt_veg[hh] = 0.0;
+ delta_Q_veg[hh] = 0.0;
+
+ // Is there a vegsnowcover?
+
+ vegsnowcover[hh] = snow_h2o_veg[hh] > 1e-6; // skip energy
+
+ // determine relative measurement heights
+ if (!relative_hts[hh]) {
+ rel_z_T[hh] = z_T[hh];
+ rel_z_u[hh] = z_u[hh];
+ }
+ else {
+ rel_z_T[hh] = z_T[hh] - z_s[hh];
+ if (rel_z_T[hh] < 1.0)
+ rel_z_T[hh] = 1.0;
+ rel_z_u[hh] = z_u[hh] - z_s[hh];
+ if (rel_z_u[hh] < 1.0)
+ rel_z_u[hh] = 1.0;
+ }
+
+ // Calculate wind speed used in both the ebal and mass bal routines wrt snow intercepted in the canopy.
+ compute_canopy_snow_wind();
+
+ // Calculate energy transfer terms
+
+ if (!_e_bal_veg())
+ return 0;
+
+ // Adjust mass and calculate runoff
+
+ _mass_bal();
+
+ // Update the averages for the energy terms and the totals for mass changes since the last output.
+
+ // increment time
+ current_time[hh] += time_step[hh];
+
+ // update interval values
+ delmelt_veg_int[hh] += delmelt_veg[hh];
+ delsub_veg_int[hh] += delsub_veg[hh];
+ delevap_veg_int[hh] += delevap_veg[hh];
+ delunld_int[hh] += delunld[hh];
+ delunld_wind_int[hh] += delunld_wind[hh];
+ delunld_melt_int[hh] += delunld_melt[hh];
+ delunld_subl_int[hh] += delunld_subl[hh];
+ deldrip_veg_int[hh] += deldrip_veg[hh];
+
+ // Update the model's input parameters
+ S_n[hh] += input_deltas[hh][tstep->level].S_n;
+ I_lw[hh] += input_deltas[hh][tstep->level].I_lw;
+ T_a[hh] += input_deltas[hh][tstep->level].T_a;
+ e_a[hh] += input_deltas[hh][tstep->level].e_a;
+ u[hh] += input_deltas[hh][tstep->level].u;
+
+ return 1;
+}
+
+/*
+** NAME
+** void compute_canopy_snow_wind.
+**
+** DESCRIPTION
+** Calculate wind speed used in both the ebal and mass bal routines wrt snow intercepted in the canopy.
+**
+** RETURN VALUE
+** void
+**
+** GLOBAL VARIABLES READ
+**
+** CanopyWindSwitchCanSno[hh]
+** Ht[hh]
+** u[hh]
+** rel_z_u[hh]
+** z_u[hh]
+**
+** GLOBAL VARIABLES MODIFIED
+** u_2_3rds[hh]
+**
+*/
+void ClassCanopySnowBalanceBase::compute_canopy_snow_wind(void)
+{
+ switch ((int)CanopyWindSwitchCanSno[hh])
+ {
+ case 0:
+ { // used if wind speed is measured within canopy
+ // scale wind to top of canopy and then down to within canopy using cionco
+ double z_veg_u = Ht[hh]*(2.0/3.0);
+ double u_veg_ht = 0.0;
+
+ adst_wind_cpy_top(Ht[hh], u[hh], rel_z_u[hh], u_veg_ht);
+ cionco_canopy_wind_spd(Ht[hh], u_veg_ht, z_veg_u, u_2_3rds[hh]);
+ break;
+ } // end wind case 0
+
+ case 1:
+ { // used if wind speed is measured at or above canopy, then wind speed is scaled down to within canopy using cionco
+ double z_veg_u = Ht[hh]*(2.0/3.0);
+
+ cionco_canopy_wind_spd(Ht[hh], u[hh], z_veg_u, u_2_3rds[hh]);
+ break;
+ } // end wind case 1
+
+ case 2:
+ { // Canopy wind profile developed at Fortress sparse canopy
+ double d0_wind = 0.5791121; // displacement height observed at sparse forest around Fortress Forest Tower
+ double z0m_wind = 0.4995565; // roughness length observed at above site
+ double z_veg_u = Ht[hh]*(2.0/3.0);
+
+
+ // wind speed for wind induced unloading at 2/3 canopy height
+ if ((z_veg_u - d0_wind) > z0m_wind)
+ {
+ double Ustar = u[hh] * PBSM_constants::KARMAN / (log((z_u[hh] - d0_wind) / z0m_wind));
+ u_2_3rds[hh] = Ustar / PBSM_constants::KARMAN * log((z_veg_u - d0_wind) / z0m_wind);
+ }
+ else
+ {
+ u_2_3rds[hh] = 0.0;
+ }
+
+ break;
+ } // end wind case 2
+ } // end of switch CanopyWindSwitchCanSno
+}
+
+/*
+** NAME
+** int _e_bal_veg(void) -- calculates point energy budget for the vegsnowcover. uses a similar iteration scheme as in CLASS canopy snow.
+**
+** DESCRIPTION
+** Calculates point energy budget for the vegsnowcover and adjusts the canopy snowpack temperature.
+**
+** RETURN VALUE
+**
+** true The calculations were completed.
+**
+** false An error occured, and a message explaining the error has
+** been stored with the 'usrerr' routine.
+**
+** GLOBAL VARIABLES READ
+**
+** GLOBAL VARIABLES MODIFIED
+**
+*/
+
+int ClassCanopySnowBalanceBase::_e_bal_veg(void)
+{
+
+ if (vegsnowcover[hh] || liq_h2o_veg[hh] > 0.0)
+ {
+ double resid; // left over energy not accounted for in the ebal
+ double Tstep = 0.25; // increment to adjust canopy snow temperature by
+ int niter = 1;
+ int max_iter = 100;
+ double T0_s_veg = T_s_veg[hh];
+
+ do
+ { /* iterate on canopy snow e bal */
+
+ // calculate radiative terms
+ _net_rad_veg();
+
+ // calculate latent and sensible heat transfer
+ if (!init_turb_transfer())
+ return 0;
+
+ // calculate advection
+ _advec_veg();
+
+ // veg snowpack energy budget
+ delta_Q_veg[hh] = ((T_s_veg[hh] - T0_s_veg) * m_s_veg[hh] * CP_ICE(T_s_veg[hh])) / time_step[hh];
+
+ resid = Qn_veg[hh] + Qh_veg[hh] + Ql_veg[hh] + Qp[hh] - delta_Q_veg[hh];
+
+ // Check energy budget residual and if too large adjust canopy snow temperature. Temp adjust algorithm below is after CLASS.
+ if (fabs(resid) < 5.0 || fabs(Tstep) < 1e-2)
+ {
+ break;
+ }
+
+ if (niter == 1)
+ {
+ if (resid > 0.0)
+ { // we overshot the canopy snow surface temperature so try reducing it 1 k step at a time
+ T_s_veg[hh] += Tstep;
+ }
+ else
+ {
+ T_s_veg[hh] -= Tstep;
+ }
+ }
+ else
+ {
+ if ((resid > 0.0 && Tstep < 0.0) || (resid < 0.0 && Tstep > 0.0))
+ {
+ Tstep = -Tstep / 2.0;
+ }
+ T_s_veg[hh] += Tstep;
+ }
+ niter += 1;
+
+ // if (std::abs(T_s_veg[hh] - FREEZE) < 1.0E-6) // likely uneeded but is in CLASS
+ // {
+ // T_s_veg[hh] = FREEZE;
+ // }
+
+ // Check if max iterations were reached
+ if (niter >= max_iter)
+ {
+ string D = StandardConverterUtility::GetDateTimeInString(Global::DTnow);
+ string SS = D + "hh " + to_string(hh);
+
+ CRHMException TExcept2(SS.c_str(), TExcept::WARNING);
+ LogError(TExcept2);
+
+ CRHMException TExcept("Warning: Max iterations reached in canopy snow energy balance", TExcept::WARNING);
+ LogError(TExcept);
+ }
+ } while (niter < max_iter); // breaks inside based on resid or temp step
+
+ if (niter > max_iter)
+ { // solution has not been reached, set canopy temp equal to air temperature
+ T_s_veg[hh] = T_a[hh];
+
+ Qn_veg[hh] = 0.0;
+
+ Qh_veg[hh] = Ql_veg[hh] = qsub_veg[hh] = 0.0;
+
+ Qp[hh] = 0.0;
+
+ delta_Q_veg[hh] = 0.0;
+ }
+ }
+ else
+ {
+ Qn_veg[hh] = 0.0;
+
+ Qh_veg[hh] = Ql_veg[hh] = qsub_veg[hh] = 0.0;
+
+ Qp[hh] = 0.0;
+
+ delta_Q_veg[hh] = 0.0;
+ }
+ return 1;
+}
+
+/*
+** NAME
+** void _net_rad_veg(void) -- calculates net allwave radiation
+**
+** DESCRIPTION
+** Calculates net allwave radiation from the net solar radiation
+** incoming thermal/longwave radiation, and the snow surface
+** temperature. This is relative to a snow clump sitting on the canopy.
+**
+** GLOBAL VARIABLES READ
+** I_lw
+** S_n
+** T_s_veg
+**
+** GLOBAL VARIABLES MODIFIED
+** Qn_veg
+*/
+
+void ClassCanopySnowBalanceBase::_net_rad_veg(void)
+{
+ // Qn_veg[hh] = S_n[hh] + (SNOW_EMISSIVITY * (I_lw[hh] - STEF_BOLTZ * pow(T_s_veg[hh], 4.0f)));
+ O_LW_cpysnow[hh] = 2.0 * SNOW_EMISSIVITY * STEF_BOLTZ * pow(T_s_veg[hh], 4.0f); // since we include incoming from atmosphere and from branch elements we multiply by two
+ I_LW_cpy_2_cpy[hh] = CAN_EMISSIVITY*(1.0 - SNOW_EMISSIVITY)*O_LW_cpysnow[hh]; // canopy longwave reflected from ground and absorbed by the canopy (W m-2)
+
+ // TODO move below to a function
+ double SW_ext_in_canopy = S_n[hh] * (1-Albedo_veg[hh] - Tauc[hh]*(1-Albedo_surface[hh])); // Eq 6. from Pomeroy et al., (2009)
+ I_LW_cpy[hh] = (CAN_EMISSIVITY * STEF_BOLTZ * pow(T_a[hh], 4.0f)) + SW_to_LW_fn[hh] * SW_ext_in_canopy; // Longwave emission from the canopy as in Pomeroy et al., (2009) Eq 10 but modified relative to calculate just the canopy emission assuming canopy is at air temperature at night
+ Qn_veg[hh] = (S_n[hh]*(1-Albedo_vegsnow[hh])) + (I_lw[hh] + I_LW_cpy_2_cpy[hh] + I_LW_cpy[hh] - O_LW_cpysnow[hh]); // assumes canopy snow is getting full atmosphere input of LW on top and LW from canopy on branch may need to adjust contribution on top to also include some canopy emission
+}
+
+/*
+** NAME
+** int init_turb_transfer(void) -- calculates turbulent transfer at a point within the canopy after Essery et al., 2003.
+**
+** DESCRIPTION
+** Calculates point turbulent transfer (Qh_veg and Ql_veg) for a 1-layer
+** vegsnowcover. Ignores stability corrections since within the roughness elements.
+**
+** GLOBAL VARIABLES READ
+**
+** GLOBAL VARIABLES MODIFIED
+**
+*/
+
+int ClassCanopySnowBalanceBase::init_turb_transfer(void) {
+
+ double e_s, e_a_fix;
+ double sat_vp;
+
+// calculate saturation vapor pressure
+
+ e_s = sati(T_s_veg[hh]);
+
+ //*** error check for bad vapor pressures ***
+
+ sat_vp = sati(T_a[hh]);
+ if (e_a[hh] > sat_vp)
+ e_a_fix = sat_vp;
+ // const_cast (e_a)[hh] = sat_vp;
+ else
+ e_a_fix = e_a[hh];
+
+ // calculate Qh_veg & Ql_veg
+
+ if (calc_turb_transfer(P_a[hh], T_a[hh], rel_z_T[hh], T_s_veg[hh], e_a_fix, e_s, u_2_3rds[hh], rel_z_u[hh], Qh_veg[hh], Ql_veg[hh]))
+ return 1; // !!!!!TB
+
+ return 1;
+}
+
+/* ----------------------------------------------------------------------- */
+
+int ClassCanopySnowBalanceBase::calc_turb_transfer(
+ double press, // aiFr pressure (Pa)
+ double ta, // air temperature (K) at height za
+ double tz, // height of air temp measurement (m)
+ double ts, // surface temperature (K)
+ double ea, // vapor pressure (Pa) at height zq
+ double es, // vapor pressure (Pa) at surface
+ double u, // wind speed (m/s) at height zu
+ double uz, // height of wind speed measurement (m)
+
+ // output variables
+
+ double &CRHM_h, // sens heat flux (+ to surf) (W/m^2)
+ double &CRHM_le // latent heat flux (+ to surf) (W/m^2)
+)
+{
+ double h = 0; // sens heat flux (+ to surf) (W/m^2)
+ double e = 0; // mass flux (+ to surf) (kg/m^2/s)
+ double qa; // specific humidity at height zq (kg/kg)
+ double qs; // specific humidity at surface (kg/kg)
+ double dens; // air density (kg/m3)
+
+ double Nr; // Reynolds number
+ double NuSh; // Nusselt number (-)
+ double D; // diffusivity of water vapour in still air (m^2/s)
+ const double KinVisc = 1.88E-5; /* Kinematic viscosity of air (Sask. avg. value) */
+ // const double ks = 0.0114; // snow shape coefficient for jack pine from Pomeroy et al., 1998
+ const double ks = 0.02; // recalibrated from Essery et al., (2003)
+ // const double ks = 0.2; // recalibrated from AC
+ const double Fract = 0.37; // fractal dimension of intercepted snow
+ double Ce; // canopy fullness index to represent how exposed snow is for sublimation, more sublimation for less snow in the canopy
+ const double Radius = 0.0005; /* Ice sphere radius, metres */
+ double dice = 900.0; // density of ice from canopy module
+ double ra; // aerodynamic resistance (s/m)
+ double ri; // resistance of moisture transfer from intercepted snow (s/m)
+ double z_0; // roughness length
+ double d_0; // displacement height
+ int ier = 0; // return error code
+
+ // check for bad input
+
+ // temperatures are Kelvin
+ if (ta <= 0 || ts <= 0)
+ {
+ CRHMException TExcept("temps not K", TExcept::TERMINATE);
+ LogError(TExcept);
+ ier = -2;
+ return (ier);
+ }
+
+ // pressures must be positive
+ if (ea <= 0 || es <= 0 || press <= 0 || ea >= press || es >= press)
+ {
+
+ string D = StandardConverterUtility::GetDateTimeInString(Global::DTnow);
+ string SS = D + "hh " + to_string(hh) + " 'calc_turb_transfer' " + ". Qh_veg: " + FloatToStrF(h, TFloatFormat::ffFixed, 10, 4) + " le: " + FloatToStrF(e, TFloatFormat::ffFixed, 10, 4);
+ SS = SS + " ta: " + FloatToStrF(ta - FREEZE, TFloatFormat::ffFixed, 10, 4) +
+ ", ts: " + FloatToStrF(ts - FREEZE, TFloatFormat::ffFixed, 10, 4) +
+ ", ea: " + FloatToStrF(ea, TFloatFormat::ffFixed, 10, 4) +
+ ", es: " + FloatToStrF(es, TFloatFormat::ffFixed, 10, 4) +
+ ", u: " + FloatToStrF(u, TFloatFormat::ffFixed, 10, 4) +
+ ", m_s_veg: " + FloatToStrF(m_s_veg[hh], TFloatFormat::ffFixed, 10, 4);
+
+ CRHM_le = 0.0; // addition TB 09/23/14
+ CRHM_h = 0.0;
+
+ CRHMException TExcept2(SS.c_str(), TExcept::WARNING);
+ LogError(TExcept2);
+
+ // Manishankar made this change to keep the program running.
+ // CRHMException TExcept("calc_turb_transfer: pressures must be positive", TERMINATE);
+ CRHMException TExcept("calc_turb_transfer: pressures must be positive", TExcept::WARNING);
+ LogError(TExcept);
+
+ ier = -2;
+ return (ier);
+ }
+
+ // vapor pressures can't exceed saturation if way off stop
+ if ((es - 25.0) > sati(ts) || (ea - 25.0) > satw(ta))
+ {
+ CRHMException TExcept("calc_turb_transfer: vapor pressures can't exceed saturation", TExcept::TERMINATE);
+ LogError(TExcept);
+ ier = -2;
+ return (ier);
+ }
+ // else fix them up
+ if (es > sati(ts))
+ {
+ es = sati(ts);
+ }
+ if (ea > satw(ta))
+ {
+ ea = satw(ta);
+ }
+
+ D = 2.06E-5 * pow(ts / FREEZE, -1.75);
+
+ if (u < 0.1)
+ u = 0.1;
+ else if (u > 15)
+ u = 15;
+
+ Nr = 2.0 * Radius * u / KinVisc;
+ NuSh = 1.79 + 0.606 * sqrt(Nr);
+
+ // air density at press, virtual temp of geometric mean of air and surface
+ dens = GAS_DEN(press, MOL_AIR, VIR_TEMP(sqrt(ta * ts), sqrt(ea * es), press));
+
+ // convert vapor pressures to specific humidities, to absolute handled within the Qe calc
+ qa = SPEC_HUM(ea, press);
+ qs = SPEC_HUM(es, press);
+
+ // calculates how exposed canopy snow is to the atmosphere based on the fullness of the canopy, more exposed with less snow as more gaps / surface area
+ if ((m_s_veg[hh] / Lmax[hh]) <= 0.0) // Using Lmax instead of Lstar as gives more appropriate index of canopy fullness
+ Ce = 0.07;
+ else
+ Ce = ks * pow((m_s_veg[hh] / Lmax[hh]), -Fract); // Ce is higher when the canopy is less full with snow as more of it is exposed, TODO maybe limit snow canopy fraction to 1.0 also need to reconsider Lstar
+
+ d_0 = Ht[hh] * (2/3);
+ z_0 = Ht[hh] * 0.1;
+
+ // resitances
+ ra = 1.0/(VON_KARMAN2*VON_KARMAN2*u)*(log((tz - d_0)/z_0)*log(((uz - d_0)/z_0))); // Allen 1998 Eq. 4
+ ri = 2.0 * dice * Radius * Radius / (3.0 * Ce * m_s_veg[hh] * D * NuSh); // Eq. 28 from Essery et al., 2003
+
+ CRHM_le = (dens / (ra + ri)) * (qa - qs) * LH_SUB(ts); // Eq. 29 from Essery et al., 2003
+ CRHM_h = (dens / ra) * CP_AIR * (ta - ts); // Eq. 4 Essery et al., 2003 and Pomeroy et al., 2016
+
+ return (ier);
+}
+
+/*
+** NAME
+** void _advec_veg(void) -- calculates advected energy at a point
+**
+** DESCRIPTION
+** This routine calculates the advected energy for a 2-layer vegsnowcover
+** if there's precipitation for the current timestep.
+**
+** GLOBAL VARIABLES READ
+** m_rain
+** m_snow
+** precip_now_veg
+** T_rain_veg
+** T_s_veg
+** T_sf
+** time_step
+**
+** GLOBAL VARIABLES MODIFIED
+** Qp
+*/
+
+void ClassCanopySnowBalanceBase::_advec_veg(void) {
+
+ if (precip_now_veg[hh]) {
+ Qp[hh] = (heat_stor(CP_WATER(T_rain_veg[hh]), m_rain[hh], (T_rain_veg[hh] - T_s_veg[hh])) +
+ heat_stor(CP_ICE(T_sf[hh]), m_snow[hh], (T_sf[hh] - T_s_veg[hh]))) / time_step[hh];
+ }
+ else
+ Qp[hh] = 0.0;
+}
+
+/*
+** NAME
+** void _mass_bal(void) -- calculates point mass budget of snow intercepted in the canopy
+**
+** DESCRIPTION
+** Calculates the point mass budget for snow intercepted in the canopy (i.e., adds new precip, sublimation/evap mass removal, melt/drip, unloading)
+** model. It then solves for new snow temperatures.
+**
+** GLOBAL VARIABLES READ
+**
+** GLOBAL VARIABLES MODIFIED
+**
+*/
+
+void ClassCanopySnowBalanceBase::_mass_bal(void) {
+
+ // adjust mass and calc. runoff
+
+ // age snow by compacting snow due to time passing
+ // _time_compact(); rm for now, likely needs to be adjusted for canopy snow
+
+ // process precipitation event
+ _precip_veg();
+
+ // calculate sublimation/evaporation and adjust canopy snowpack
+
+ _subl_evap();
+
+ // calculate melt or freezing and adjust cold content
+
+ _snowmelt();
+
+ // calculate mass unloading of snow due to wind, sublimation, melting
+
+ _mass_unld();
+
+ // calculate runoff of liquid water intercepted in the canopy */
+
+ _runoff_veg();
+
+ // adjust layer temps if there was a vegsnowcover at start of the
+ // timestep and there's still snow on the ground
+
+ if(snow_h2o_veg[hh] > 0.0){
+ vegsnowcover[hh] = 1;
+ T_s_veg[hh] = new_tsno_veg(m_s_veg[hh], T_s_veg[hh], cc_s_veg[hh]);
+ } else {
+ vegsnowcover[hh] = 0;
+ }
+
+}
+
+/*
+** NAME
+** void _precip_veg(void) -- process a precipitation event
+**
+** DESCRIPTION
+** This routine processes a precipitation event, i.e., the current
+** precip record, if there's one for the current timestep. It
+** determines if the precip is rain or snow which increases the
+** vegsnowcover.
+**
+** GLOBAL VARIABLES READ
+** h2o_sat_veg_snow
+** m_rain
+** m_precip
+** max_h2o_vol_veg
+** precip_now_veg
+** vegsnowcover
+** T_sf
+** z_snow_veg
+**
+** GLOBAL VARIABLES MODIFIED
+** liq_h2o_veg
+** h2o_sat_veg
+** liq_h2o_veg
+** T_s_veg
+** T_s_veg
+*/
+
+void ClassCanopySnowBalanceBase::_precip_veg(void)
+{
+
+ if (precip_now_veg[hh]) {
+ if (vegsnowcover[hh] || liq_h2o_veg[hh] > 0.0) {
+
+ // Adjust vegsnowcover's load by new snow/rain.
+ liq_h2o_veg[hh] += m_rain[hh];
+ snow_h2o_veg[hh] += m_snow[hh];
+ m_s_veg[hh] = liq_h2o_veg[hh] + snow_h2o_veg[hh];
+ }
+ else { // no vegsnowcover
+
+ // Use snowfall, if any, to setup a new vegsnowcover.
+
+ if (m_precip[hh] > 0.0) {
+ init_snow_veg();
+ }
+ }
+
+ }
+
+}
+
+/*
+** NAME
+** void _snowmelt(void) -- calculates melt or re-freezing at a point
+**
+** DESCRIPTION
+** Calculates melting or re-freezing for point 2-layer energy balance
+** snowmelt model.
+**
+** GLOBAL VARIABLES READ
+**
+** GLOBAL VARIABLES MODIFIED
+**
+*/
+
+void ClassCanopySnowBalanceBase::_snowmelt(void) {
+
+ double Q_0; // energy available for surface melt
+ double Q_freeze; // energy used for re-freezing
+ double Q_left; // energy left after re_freezing
+ double h2o_refrozen; // amount of liquid H2O that was refrozen
+
+// If no snow in the canopy, then just exit.
+
+ if (!vegsnowcover[hh]) {
+ delmelt_veg[hh] = 0.0;
+ return;
+ }
+
+ // calculate melt or freezing, and adjust cold content
+
+ // calculate surface melt
+
+ // energy for surface melt
+ Q_0 = (delta_Q_veg[hh] * time_step[hh]) + cc_s_veg[hh];
+
+ if (Q_0 > 0.0) {
+ delmelt_veg[hh] = MELT(Q_0); // liquid water will only freeze on next timestep
+ if(snow_h2o_veg[hh] - delmelt_veg[hh] < 0.0){
+ delmelt_veg[hh] = snow_h2o_veg[hh]; // cannot melt more than the amount of snow in the canopy
+ }
+ cc_s_veg[hh] = 0.0; // all energy is used to melt snow
+ }
+ else if (Q_0 == 0.0) {
+ delmelt_veg[hh] = 0.0;
+ cc_s_veg[hh] = 0.0;
+ }
+ else {
+ delmelt_veg[hh] = 0.0;
+ cc_s_veg[hh] = Q_0;
+
+ // if (m_s_veg[hh] < 2.0) // addition TB 06/03/10, TODO mod for canopy snow
+ // cc_s_veg[hh] = 0.0;
+ }
+
+ liq_h2o_veg[hh] += delmelt_veg[hh];
+ snow_h2o_veg[hh] -= delmelt_veg[hh];
+
+ // adjust layers for re-freezing
+
+ // adjust surface layer
+
+ h2o_refrozen = 0.0;
+
+ if (cc_s_veg[hh] < 0.0) { // cannot enter here if snow is melting on current timestep since cc_s_veg == 0
+ // if liquid h2o present, calc refreezing and adj cc_s_veg, assumes liquid is at freezing point
+ if (liq_h2o_veg[hh] > 0.0) {
+ Q_freeze = liq_h2o_veg[hh] * LH_FUS(FREEZE);
+ Q_left = Q_0 + Q_freeze;
+
+ // cold content is negative considering updated Q_0, thus enough energy to freeze all liquid water
+ if (Q_left <= 0.0) {
+ h2o_refrozen = liq_h2o_veg[hh];
+ cc_s_veg[hh] = Q_left;
+ }
+ else {
+ h2o_refrozen = liq_h2o_veg[hh] - MELT(Q_left);
+ cc_s_veg[hh] = 0.0;
+ }
+ }
+ }
+
+ // Note: because of rounding errors, h2o_refrozen may not
+ // be exactly the same as liq_h2o_veg. Check for this case, and if so, then just zero out liq_h2o_veg.
+
+ if (fabs(liq_h2o_veg[hh] - h2o_refrozen) <= 1e-8) {
+ liq_h2o_veg[hh] = 0.0;
+ }
+ else {
+ liq_h2o_veg[hh] -= h2o_refrozen;
+ snow_h2o_veg[hh] += h2o_refrozen;
+
+ }
+
+ // determine if vegsnowcover is isothermal
+
+ if ((vegsnowcover[hh] == 1) && (cc_s_veg[hh] == 0.0))
+ isothermal[hh] = 1;
+ else
+ isothermal[hh] = 0;
+
+ m_s_veg[hh] = liq_h2o_veg[hh] + snow_h2o_veg[hh];
+}
+
+/*
+** NAME
+** void _mass_unld(void) -- adjust the canopy snowpack's snow storage due to mass unloading of snow liquid water is handled seperately under runoff
+**
+** DESCRIPTION
+** This is the updated mass snow unloading parameterisations from
+** Cebulski & Pomeroy 2025 to unload based on wind, snowmelt, and sublimation
+**
+** GLOBAL VARIABLES READ
+**
+** GLOBAL VARIABLES MODIFIED
+**
+*/
+
+void ClassCanopySnowBalanceBase::_mass_unld(void)
+{
+ delunld[hh] = 0.0; // need to clear this val otherwise accumulates due to += below
+
+ if (!vegsnowcover[hh]) { // no snow in the canopy exit early
+ delunld_wind[hh] = 0.0;
+ delunld_melt[hh] = 0.0;
+ delunld_subl[hh] = 0.0;
+ return;
+ }
+
+ switch ((int)MassUnloadingSwitch[hh]) // options are 0 for Cebulski & Pomeroy 2025 ablation paper; 1 for Andreadis 2009; 2 for Roesch 2001. HP98/Ellis is enabled using the original canopyclearinggap module due to major differences in how it handles canopy snowmelt, icebulb temp using the pom98 analytical sublimation, etc.
+ {
+ case 0: // This is the updated mass snow unloading parameterisations from Cebulski & Pomeroy 2025 ablation paper
+ {
+ // check maximum canopy snow load spill what is overflowing
+ if (snow_h2o_veg[hh] > Lmax[hh])
+ {
+ delunld[hh] = snow_h2o_veg[hh] - Lmax[hh];
+ }
+
+ // melt induced mass unloading of solid snow based on ratio relative to canopy snowmelt found to be function of canopy snow load for Fortress obs
+ double unld_to_melt_ratio_m = 0.15;
+ double unld_to_melt_ratio_b = -0.43;
+
+ double unld_to_melt_ratio = snow_h2o_veg[hh] * unld_to_melt_ratio_m + unld_to_melt_ratio_b; // WARNING this can go negative so handle below
+ unld_to_melt_ratio = std::max(0.0, unld_to_melt_ratio);
+ delunld_melt[hh] = delmelt_veg[hh] * unld_to_melt_ratio;
+
+ // mechanical wind induced unloading
+ // const double a_u = 1.282646e-06; // Cebulski & Pomeroy coef from exponential function of unloading as function of wind speed and canopy snow load measurements at Fortress mountain w no melt.
+ // const double b_u = 3.925391e-01; // TODO move to par file
+
+ // wind induced unloading as function of shear stress
+ const double a_tau = 0.371/(60.0*60.0); // derrived from unloading vs. shear stress relationship presented in Cebulski & Pomeroy 2025 (HESS paper) converting from per hour to per sec here.
+
+ double fu = 0.0;
+
+ if (u_2_3rds[hh] >= 0.0)
+ {
+ // fu = u_2_3rds[hh] * a_u * exp(b_u * u_2_3rds[hh]); // unloading rate due to wind (s-1)
+ double tau_mid = u_2_3rds[hh] * u_2_3rds[hh] * 0.02; // wind to tau conversion developed at forest tower using obs shear stress vs. wind speed
+ fu = a_tau * tau_mid; // unloading rate due to wind as predicted by shear stress (s-1) multiplied by canopy load later
+ }
+ else
+ {
+ fu = 0.0; // less than wind induced unloading threshold so set equal to 0.
+ }
+
+ // ablation via temperature, wind, and duration based unloading
+ // delunld[hh] = snow_h2o_veg[hh] * (fT + fu + ft) * dt; // ODE solution: calculate solid snow unloading over the time interval
+
+ // analytical solution which is more exact over longer time intervals, following from Cebulski & Pomeroy derivation of the HP98 unloading parameterisation
+ delunld_wind[hh] = snow_h2o_veg[hh] * (1 - exp(-fu * time_step[hh])); // similar analytical solution for ODE equation 30 in Cebulski & Pomeroy 2025 Wires WATER Review ,similar as steps from eq 27 to 28.
+ delunld[hh] += delunld_wind[hh];
+ delunld[hh] += delunld_melt[hh];
+
+ break;
+
+ } // end case 0
+
+ case 1:
+ { // This is the mass snow unloading parameterisation from Andreadis et al., 2009 based on the snowmelt rate. No wind induced unloading.
+
+ double Lres = 5.0; // threshold above which melt unloading starts to occur.
+
+ if (snow_h2o_veg[hh] > Lres)
+ {
+ double melt_drip_ratio_andreadis = 0.4;
+ delunld_melt[hh] = delmelt_veg[hh] * melt_drip_ratio_andreadis;
+ } else {
+ delunld_melt[hh] = 0.0; // snow load is below threshold and can only melt or sublimate
+ }
+
+ delunld[hh] += delunld_melt[hh];
+
+ break;
+
+ } // case 2
+
+ case 2:
+ { // This is the mass snow unloading parameterisation from Roesch et al., 2001 based on a temperature and wind speed function.
+
+ // check maximum canopy snow load
+ if (snow_h2o_veg[hh] > Lmax[hh])
+ {
+ delunld[hh] = snow_h2o_veg[hh] - Lmax[hh];
+ }
+
+ // Mass unloading of canopy snow due to melt
+
+ const double C1 = 1.87e5; // constant provided in Roesch et al., 2001
+ double fT = 0.0; // unloading rate per second.
+
+ if ((T_a[hh]) < CRHM_constants::Tm) // Follows equation 31 in Cebulski & Pomeroy 2025 Wires WATER Review.
+ {
+ fT = 0.0;
+ } else {
+ fT = (T_a[hh] - CRHM_constants::Tm)/C1;
+ }
+
+ delunld_melt[hh] = snow_h2o_veg[hh] * (1-exp(-(fT) * time_step[hh])); // analytical solution for ODE equation 30 in Cebulski & Pomeroy 2025 Wires WATER Review
+
+ // Mass unloading of canopy snow due to mechanical removal from wind at canopy top
+ const double C2 = 1.56e5; // wind unloading constant from Roesch et al., 2001
+ double fu = 0.0;
+ double u_ht = 0.0;
+
+ adst_wind_cpy_top(Ht[hh], u[hh], z_u[hh], u_ht);
+
+ if(u_ht >= 0.0){
+ fu = u_ht/C2;
+ } else {
+ fu = 0.0; // unloading rate due to wind (s-1)
+ }
+
+ delunld_wind[hh] = snow_h2o_veg[hh] * (1-exp(-(fu) * time_step[hh])); // analytical solution for ODE equation 30 in Cebulski & Pomeroy 2025 Wires WATER Review
+
+ delunld[hh] += delunld_wind[hh];
+ delunld[hh] += delunld_melt[hh];
+
+ break;
+
+ } // end case 3
+ } // MassUnloadingSwitch
+
+
+
+ if (delunld[hh] > snow_h2o_veg[hh])
+ {
+ delunld[hh] = snow_h2o_veg[hh];
+ snow_h2o_veg[hh] = 0.0;
+ }
+ else
+ {
+ snow_h2o_veg[hh] -= delunld[hh];
+ }
+
+ m_s_veg[hh] = liq_h2o_veg[hh] + snow_h2o_veg[hh];
+
+}
+
+/*
+** NAME
+** void _subl_evap(void) -- calculates evaporation/condensation at a point
+**
+** DESCRIPTION
+** Calculates mass lost or gained by sublimation/evaporation/condensation
+**
+** GLOBAL VARIABLES READ
+** qsub_veg
+** time_step
+**
+** GLOBAL VARIABLES MODIFIED
+** delsub_veg snow_h2o_veg liq_h2o_veg
+**
+*/
+
+void ClassCanopySnowBalanceBase::_subl_evap(void) {
+
+ if (!vegsnowcover[hh] && !(snow_h2o_veg[hh] > 0.0)) { // no h2o in the canopy
+ delsub_veg[hh] = 0.0;
+ } else {
+ // Compute total mass change due to sublimation/evaporation over the time step
+ qsub_veg[hh] = Ql_veg[hh] / LH_SUB(T_s_veg[hh]);
+ delsub_veg[hh] = qsub_veg[hh] * time_step[hh];
+
+ // Apply mass changes based on the presence of snow or liquid water
+ if (snow_h2o_veg[hh] > -delsub_veg[hh]) {
+ // Snow is present in the canopy; modify snow water equivalent
+ snow_h2o_veg[hh] += delsub_veg[hh];
+
+ } else {
+ delsub_veg[hh] = -snow_h2o_veg[hh];
+ snow_h2o_veg[hh] = 0.0;
+ }
+
+ }
+
+ if (!(liq_h2o_veg[hh] > 0.0)){
+ delevap_veg[hh] = 0.0;
+ } else {
+ // liquid water is present; modify liquid water content
+ // Compute total mass change due to sublimation/evaporation over the time step
+ delevap_veg[hh] = (Ql_veg[hh] / LH_VAP(T_s_veg[hh])) * time_step[hh];
+
+ // Apply mass changes based on the presence of snow or liquid water
+ if (liq_h2o_veg[hh] > -delevap_veg[hh]) {
+ // Snow is present in the canopy; modify snow water equivalent
+ liq_h2o_veg[hh] += delevap_veg[hh];
+
+ } else {
+ delevap_veg[hh] = -liq_h2o_veg[hh];
+ liq_h2o_veg[hh] = 0.0;
+ }
+
+ }
+ m_s_veg[hh] = liq_h2o_veg[hh] + snow_h2o_veg[hh];
+}
+
+#define BB 0.4 // changed from 0.05
+
+/*
+** NAME
+** void _runoff_veg(void) -- calculates runoff from vegsnowcover
+**
+** DESCRIPTION
+** Calculates runoff for point energy budget 2-layer snowmelt model
+**
+** GLOBAL VARIABLES READ
+** liq_h2o_veg
+** vegsnowcover
+** max_h2o_vol_veg
+**
+** GLOBAL VARIABLES MODIFIED
+** liq_h2o_veg
+** max_liq_veg_frac
+** h2o_sat_veg
+** h2o_vol_veg
+** deldrip_veg
+*/
+
+void ClassCanopySnowBalanceBase::_runoff_veg(void) {
+
+ // Determine runoff, and water left in the snow
+ if((int)MassUnloadingSwitch[hh] == 1){ // Andreadis2009
+ max_liq_veg[hh] = m_s_veg[hh] * max_liq_veg_frac[hh] + exp(-4.0) * LAI[hh] * 2.0; // times 2 for two sided lai as in Andreadis2009
+ } else { // From CRHM Canopy
+ max_liq_veg[hh] = m_s_veg[hh] * max_liq_veg_frac[hh] + Cc[hh]*LAI[hh]*0.2;
+ }
+
+ if (liq_h2o_veg[hh] > max_liq_veg[hh]) {
+
+ deldrip_veg[hh] = liq_h2o_veg[hh] - max_liq_veg[hh];
+ liq_h2o_veg[hh] = max_liq_veg[hh];
+ m_s_veg[hh] = liq_h2o_veg[hh] + snow_h2o_veg[hh];
+
+ } else {
+ deldrip_veg[hh] = 0.0;
+ }
+
+
+}
+
+/* ----------------------------------------------------------------------- */
+
+double ClassCanopySnowBalanceBase::new_tsno_veg(
+ double spm, /* layer's specific mass (kg/m^2) */
+ double t0, /* layer's last temperature (K) */
+ double ccon) /* layer's adjusted cold content (J/m^2) */
+{
+ double tsno;
+ double cp;
+ double tdif;
+
+ cp = CP_ICE(t0);
+
+ tdif = ccon / (spm * cp);
+ tsno = tdif + FREEZE;
+
+ if (tsno < MIN_SNOW_TEMP + FREEZE)
+ tsno = MIN_SNOW_TEMP + FREEZE;
+ else {
+ if (tsno > FREEZE)
+ tsno = FREEZE;
+ }
+
+ return (tsno);
+}
+/* ----------------------------------------------------------------------- */
diff --git a/crhmcode/src/modules/ClassCanopySnowBalanceBase.h b/crhmcode/src/modules/ClassCanopySnowBalanceBase.h
new file mode 100644
index 000000000..0992067f1
--- /dev/null
+++ b/crhmcode/src/modules/ClassCanopySnowBalanceBase.h
@@ -0,0 +1,255 @@
+/**
+* Copyright 2022, CRHMcode's Authors or Contributors
+* This file is part of CRHMcode.
+*
+* CRHMcode is free software: you can redistribute it and/or modify it under
+* the terms of the GNU General Public License as published by the Free Software
+* Foundation, either version 3 of the License, or (at your option) any later
+* version.
+*
+* CRHMcode is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty
+* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+* See the GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License along with
+* CRHMcode. If not, see .
+*
+**/
+#ifndef ClassCanopySnowBalanceBaseH
+#define ClassCanopySnowBalanceBaseH
+
+//created by Manishankar Mondal
+
+#include "../core/ClassModule.h"
+#include "INPUT_REC.h"
+#include "TSTEP_REC.h"
+#include "PRECIP_REC.h"
+
+class ClassCanopySnowBalanceBase : public ClassModule {
+
+public:
+
+ ClassCanopySnowBalanceBase(string Name, string Version = "undefined", LMODULE Lvl = LMODULE::PROTO) : ClassModule(Name, Version, Lvl) {};
+
+ // declared variables
+
+ // snowpack information
+
+ double* snow_h2o_veg{ NULL }; // snow h2o content as specific mass(kg/m^2)
+ double* liq_h2o_veg{ NULL }; // liquid h2o content as specific mass(kg/m^2)
+ double* m_s_veg{ NULL }; // snowcover's specific mass, liquid and snow (kg/m^2). Init by init_snow_veg.
+ double* T_s_veg{ NULL }; // average snowcover temp (K). Init by init_snow_veg
+ double* cc_s_veg{ NULL }; // snowcover's cold content (J/m^2). Init by init_snow_veg.
+ double* h2o_sat_veg{ NULL }; // % of liquid H2O saturation (relative water content, i.e., ratio of water in snowcover
+ double* h2o_vol_veg{ NULL }; // liquid h2o content as volume ratio: V_water/(V_snow - V_ice) (unitless).init_snow_veg
+ double* max_liq_veg{ NULL }; // max liquid h2o content as specific mass(kg/m^2)
+
+// energy balance info for current timestep
+
+ double* Qn_veg{ NULL }; // net allwave radiation wrt the canopy (W/m^2)
+ double* Qh_veg{ NULL }; // sensible heat xfr positive towards the canopy (W/m^2)
+ double* Ql_veg{ NULL }; // latent heat xfr positive towards the canopy (W/m^2)
+ double* Qp{ NULL }; // advected heat from precip wrt the canopy (W/m^2)
+ double* delta_Q_veg{ NULL }; // change in snowcover's energy wrt the canopy (W/m^2)
+ const double *Tauc{ NULL }; // Canopy transmittance through the entire canopy calculated in the canopy module
+
+
+// mass balance vars for current timestep
+
+ double* delmelt_veg_int{ NULL }; // specific melt (kg/m^2 or m)
+ double* delL{ NULL }; // interval change in SWE
+ double* delevap_veg_int{ NULL }; // mass of evap into air & soil from snowcover (kg/m^2*int) delunld_int
+ double* delsub_veg_int{ NULL }; // mass of evap into air & soil from snowcover (kg/m^2*int) delunld_int
+ double* delunld_int{ NULL }; // specific mass of canopy snow unloaded to subcanopy (kg/m^2*int)
+ double* delunld{ NULL }; // canopy snow unloading rate (kg/m^2*s)
+ double* delunld_wind{ NULL }; // solid snow unloading from the canopy induced by wind (kg/m^2*s)
+ double* delunld_melt{ NULL }; // canopy snow unloading rate due to melting (kg/m^2*s)
+ double* delunld_subl{ NULL }; // canopy snow unloading due to sublimation (kg/m^2*s)
+ double* delunld_wind_int{ NULL }; // solid snow unloading from the canopy induced by wind (kg/m^2*int)
+ double* delunld_melt_int{ NULL }; // canopy snow unloading rate due to melting (kg/m^2*int)
+ double* delunld_subl_int{ NULL }; // canopy snow unloading due to sublimation (kg/m^2*int)
+ double* deldrip_veg_int{ NULL }; // canopy snowmelt drainage (kg/m^2*int)
+
+ double* delmelt_veg_day{ NULL }; // daily predicted specific runoff (m/sec)
+ double* cmlmelt_veg_day{ NULL }; // daily predicted specific runoff accumulator (m/sec)
+
+ const double* hru_evap{ NULL }; // liquid water evaporated off the canopy. computed in the evap module aka "internal evap". (kg m^-2)
+
+ double* net_rain { NULL }; // throughfall rain + drip (runoff/melt drainage) of snow intercepted in the canopy (kg/m^2*int)
+ double* net_snow { NULL }; // throughfall snow + unloaded snow from the canopy (kg/m^2*int)
+ double* net_p { NULL }; // net rain + net snow (kg/m^2*int)
+
+ const double *throughfall_rain { NULL }; // throughfall of rain, not in contact with the canopy to be added with canopy drip (kg/m^2*int)
+ const double *throughfall_snow { NULL }; // throughfall of snow, not in contact with the canopy to be added with canopy snow unloading (kg/m^2*int)
+ const double *z_s { NULL }; // total surface snowcover thickness from snobal module for relative heights calculation (m)
+
+// mass balance vars for variable timestep
+
+ double* delmelt_veg{ NULL }; // specific melt (kg/m^2 or m)
+ double* qsub_veg{ NULL }; // mass flux by subl/evap (+ to surf) (kg/m^2/s)
+ double* delsub_veg{ NULL }; // mass flux by subl/evap (+ to surf) (kg/m^2/int)
+ double* delevap_veg{ NULL }; // mass flux by subl/evap (+ to surf) (kg/m^2/int)
+ double* deldrip_veg{ NULL }; // predicted specific runoff (m/sec)
+ double* E_l{ NULL }; // mass flux by evap/cond to soil (kg/m^2/s)
+
+// precipitation info adjusted for current run timestep
+
+ double* m_precip{ NULL }; // specific mass of total precip (kg/m^2)
+ double* m_rain{ NULL }; // specific mass of rain in precip (kg/m^2)
+ double* m_snow{ NULL }; // specific mass in snow in precip (kg/m^2)
+ double* T_pp{ NULL }; // precip temp (K)
+
+// precipitation info for the current DATA timestep
+
+ long* precip_now_veg{ NULL }; // precipitation occur for current timestep?
+ double* T_rain_veg{ NULL }; // rain's temp (K)
+ double* T_sf{ NULL }; // snowfall's temp (K)
+ double* h2o_sat_veg_snow{ NULL }; // snowfall's % of liquid H2O saturation
+
+// local climate-data values for the current run timestep
+ double* S_n{ NULL }; // net solar radiation (W/m^2)
+ double* I_lw{ NULL }; // incoming longwave (thermal) rad (W/m^2)
+ double* T_a{ NULL }; // air temp (K)
+ double* e_a{ NULL }; // vapor pressure (Pa)
+ double* u{ NULL }; // wind speed (m/sec)
+ const double *T_s_0 { NULL }; // temp of the surface snowpack active layer (C)
+
+
+ long* isothermal{ NULL }; // melting?
+ long* vegsnowcover{ NULL }; // snow on veg at start of current timestep?
+ long* stop_no_snow{ NULL }; //
+
+// local variables
+ double* P_a{ NULL }; // air pressure (Pa)
+ double* m_precip_cum{ NULL }; //
+ double* m_rain_cum{ NULL }; //
+ double* m_snow_cum{ NULL }; //
+ double* E_s_cum{ NULL }; //
+ double* cmlmelt_veg{ NULL }; //
+ double* Fault{ NULL }; //
+ double* I_LW_atm{ NULL }; // Downwelling longwave from the atmoshpere (W/m^2)
+ double* I_LW_gnd{ NULL }; // Upwelling longwave from the ground (W/m^2)
+ double* I_LW_cpy_2_cpy{ NULL }; // Longwave from the canopy reflected off the surface back to the canopy (W/m^2)
+ double* O_LW_cpysnow{ NULL }; // Outgoing longwave radiation emitted from the canopy snow (W/m^2)
+ double* I_LW_cpy{ NULL }; // Incoming longwave radiation emitted from the canopy (W/m^2)
+ double* u_2_3rds{ NULL }; // Wind speed at 2/3 canopy height for canopy snow energy balance and unloading (m s^-1)
+ double* rel_z_u{ NULL }; // Height of wind speed measurement relative to top of snowpack, only implemented if relative_hts[hh] == 1 (m)
+ double* rel_z_T{ NULL }; // Height of temp measurement relative to top of snowpack, only implemented if relative_hts[hh] == 1 (m)
+
+// debug variables
+/* double *Length;
+ double **Length_array;
+ double *Ustar;
+ double **Ustar_array;
+ double *e;
+ double **e_array;
+ double *h;
+ double **h_array;
+ long *ier;
+ long **ier_array;
+ long *ArrayCnt; */
+
+ // declared parameters
+
+ // measurement heights/depths
+
+ const double* hru_elev{ NULL }; // HRU elevation
+ const double* basin_area{ NULL }; // [BASIN]
+ const double* hru_area{ NULL };
+ const double *Albedo_surface{ NULL }; // albedo of surface ()
+ const double *Albedo_veg{ NULL }; // albedo of vegetation ()
+ const long* inhibit_evap{ NULL }; // inhibit evaporation, 1 -> inhibit
+ const long* relative_hts{ NULL }; // true if measurements heights, z_T and z_u, are relative to snow surface
+ // false if they are absolute heights above the ground
+ const double* z_u{ NULL }; // height of wind measurement (m)
+ const double* z_T{ NULL }; // height of air temp & vapor pressure measurement (m)
+ const double* Albedo_vegsnow{ NULL }; // albedo of snow on vegetation ()
+ const double* SW_to_LW_fn{ NULL }; // dimensionless shortwave to longwave transfer efficiency function. 0.038 from Pomeroy et al., (2009) for marmot forced through the origin, alternative value is 0.023 from Fraser site. (-)
+ const double* max_liq_veg_frac{ NULL }; // max liquid h2o content as fraction of snow mass (-)
+ const double *Cc{ NULL }; // canopy coverage, (1-sky view fraction)
+ const double *LAI{ NULL }; // LAI
+ const double *Lmax{ NULL }; // maximum canopy snow interception load, currently just used for sublimation exposure coef. 50 kg m-2 based on max observed in Storck et al. 2002, Floyd 2012 and Cebulski & Pomeroy (kg/m^2)
+ const double *Ht{ NULL }; // forest/vegetation height (m)
+ const long *CanopyWindSwitchCanSno{ NULL }; // Canopy wind model to use at height Zcan, 0 - for Cionco (dense canopy), 1 - for Prandtl-von Kármán log-linear relationship (sparse forest)".
+ const long *MassUnloadingSwitch { NULL }; // canopy snow ablation parameterization to use, 0 - Cebulski & Pomeroy 2025 ablation paper, 1- Andreadis 2009, 2 - Roesch2001 (enable HP98/Ellis2010 using original canopy clearing gap module)
+
+// void decl(void);
+
+ void init(void);
+
+ // void run(void);
+
+ void finish(bool good); // delete local storage used
+
+ void init_snow_veg(void);
+
+ double _cold_content_veg(double temp, double mass); // temperature of layer specific mass of layer
+
+ void do_data_tstep_veg(void);
+
+ int _divide_tstep_veg(TSTEP_REC* tstep); // record of timestep to be divided
+
+ int _below_thold_veg(double threshold); // current timestep's threshold for a layer's mass
+
+ int _do_tstep_veg(TSTEP_REC* tstep); // timestep's record
+
+ void compute_canopy_snow_wind(void);
+
+ int _e_bal_veg(void);
+
+ void _net_rad_veg(void);
+
+ int init_turb_transfer(void);
+
+ void _advec_veg(void);
+
+ void _mass_bal(void);
+
+ void _precip_veg(void);
+
+ void _snowmelt(void);
+
+ void _mass_unld(void);
+
+ void _subl_evap(void);
+
+ void _runoff_veg(void);
+
+ double new_tsno_veg(double spm, double t0, double ccon);
+
+ int calc_turb_transfer(double press, double ta, double rel_z_T, double ts, double ea, double es, double u, double rel_z_u,
+ double &h, double &le);
+
+ int subl_ice_sphere(double ea, double es, double ta, double ts, double u, double press);
+
+ int init_subl_ice_sphere(void);
+
+ // time step information
+
+ TSTEP_REC** tstep_info{ NULL }; // array of info for each timestep [nhru] [4]:
+ // 0 : data timestep
+ // 1 : normal run timestep
+ // 2 : medium " "
+ // 3 : small " "
+
+ long* time_step{ NULL }; // length current timestep (sec)
+ long* current_time{ NULL }; // start time of current time step (sec)
+
+ // climate-data input records
+
+ INPUT_REC* input_rec1{ NULL }; // input data for start of data timestep [nhru]
+ INPUT_REC* input_rec2{ NULL }; // " " " end " " " [nhru]
+
+ INPUT_REC** input_deltas{ NULL }; // deltas for climate-input parameters
+ // over each timestep [nhru] [4]
+
+ PRECIP_REC** precip_info{ NULL }; // array of precip info adjusted for
+ // each timestep [nhru] [4]
+
+ int** computed{ NULL }; // array of flags for each timestep;
+ // true if computed values for input
+ // deltas and precip arrays [nhru] [4]
+};
+
+#endif
diff --git a/crhmcode/src/modules/ClassCanopySnowBalanceCRHM.cpp b/crhmcode/src/modules/ClassCanopySnowBalanceCRHM.cpp
new file mode 100644
index 000000000..b74967d8b
--- /dev/null
+++ b/crhmcode/src/modules/ClassCanopySnowBalanceCRHM.cpp
@@ -0,0 +1,347 @@
+/**
+* Copyright 2022, CRHMcode's Authors or Contributors
+* This file is part of CRHMcode.
+*
+* CRHMcode is free software: you can redistribute it and/or modify it under
+* the terms of the GNU General Public License as published by the Free Software
+* Foundation, either version 3 of the License, or (at your option) any later
+* version.
+*
+* CRHMcode is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty
+* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+* See the GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License along with
+* CRHMcode. If not, see .
+*
+**/
+//created by Alex Cebulski
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include "ClassCanopySnowBalanceCRHM.h"
+#include "../core/GlobalDll.h"
+#include "../core/ClassCRHM.h"
+#include "newmodules/SnobalDefines.h"
+
+
+using namespace CRHM;
+
+ClassCanopySnowBalanceCRHM* ClassCanopySnowBalanceCRHM::klone(string name) const{
+ return new ClassCanopySnowBalanceCRHM(name);
+}
+
+void ClassCanopySnowBalanceCRHM::decl(void) {
+
+ Description = "Handles the mass and energy balance of snow intercepted in the canopy. Should be paired with the CanopyVectorBased module which does the initial loading of precip in the canopy and some radiation calculations. There is a timestep subset routine in this module following the surface snowpack routine. Standard CRHM module.' \
+ 'use Qsi (W/m^2) and Qli (W/m^2) observations,' \
+ 'use variables Qsisn_Var (W/m^2) and Qlisn_Var (W/m^2) from module CanopyClearing.' \
+ 'use observation Qsi (W/m^2) and QliVt_Var (W/m^2) from module longVt.' \
+ 'use variables QsiS_Var (W/m^2) from Annandale and QliVt_Var (W/m^2) from module longVt.'";
+
+ declstatvar("isothermal", TDim::NHRU, "melting: 0/1", "()", &isothermal);
+ declstatvar("vegsnowcover", TDim::NHRU, "snow on veg at start of current timestep: 0/1", "()", &vegsnowcover);
+
+ declvar("Qn_veg", TDim::NHRU, "net allwave radiation wrt the canopy", "(W/m^2)", &Qn_veg);
+ declvar("Qh_veg", TDim::NHRU, "sensible heat xfr wrt the canopy", "(W/m^2)", &Qh_veg);
+ declvar("Ql_veg", TDim::NHRU, "latent heat xfr wrt the canopy", "(W/m^2)", &Ql_veg);
+ declvar("Qp", TDim::NHRU, "advected heat from precip wrt the canopy", "(W/m^2)", &Qp);
+ declvar("delta_Q_veg", TDim::NHRU, "change in canopy snow's energy wrt the canopy", "(W/m^2)", &delta_Q_veg);
+
+ declstatvar("cc_s_veg", TDim::NHRU, "canopy snow's cold content", "(J/m^2)", &cc_s_veg);
+ declstatvar("liq_h2o_veg", TDim::NHRU, "canopy liquid h2o load as specific mass", "(kg/m^2)", &liq_h2o_veg);
+ declstatvar("snow_h2o_veg", TDim::NHRU, "canopy snow load as specific mass", "(kg/m^2)", &snow_h2o_veg);
+
+ decllocal("delevap_veg", TDim::NHRU, "mass of subl/evap from canopy snow/liquid (+to surface)", "(kg/m^2*s)", &delevap_veg);
+ decllocal("delsub_veg", TDim::NHRU, "mass of subl/evap from canopy snow/liquid (+to surface)", "(kg/m^2*s)", &delsub_veg);
+ decllocal("qsub_veg", TDim::NHRU, "mass flux by evap into air from active layer", "(kg/m^2*s)", &qsub_veg);
+ decllocal("delmelt_veg", TDim::NHRU, "specific melt (kg/m^2 or m)", "(kg/m^2)", &delmelt_veg);
+ decllocal("delunld_wind", TDim::NHRU, "solid snow unloading from the canopy induced by wind", "(kg/m^2*s)", &delunld_wind);
+ decllocal("delunld_melt", TDim::NHRU, "canopy snow unloading rate due to melting", "(kg/m^2*s)", &delunld_melt);
+ decllocal("delunld_subl", TDim::NHRU, "canopy snow unloading due to sublimation", "(kg/m^2*s)", &delunld_subl);
+ decllocal("delunld", TDim::NHRU, "canopy snow unloading rate", "(kg/m^2*s)", &delunld);
+ declvar("deldrip_veg", TDim::NHRU, "drip from canopy snowmelt and intercepted rainfall", "(kg/m^2*s)", &deldrip_veg);
+
+ declvar("delevap_veg_int", TDim::NHRU, "mass of evap into air from vegsnowcover", "(kg/m^2*int)", &delevap_veg_int);
+ declvar("delsub_veg_int", TDim::NHRU, "mass of sublimation into air from vegsnowcover", "(kg/m^2*int)", &delsub_veg_int);
+ declvar("delmelt_veg_int", TDim::NHRU, "specific melt (kg/m^2 or m)", "(kg/m^2*int)", &delmelt_veg_int);
+ declvar("delunld_int", TDim::NHRU, "specific mass of canopy snow unloaded to subcanopy", "(kg/m^2*int)", &delunld_int);
+ declvar("net_rain", TDim::NHRU, "throughfall rain + drip (runoff/melt drainage) of snow intercepted in the canopy", "(kg/m^2*int)", &net_rain);
+ declvar("net_snow", TDim::NHRU, "throughfall snow + unloaded snow from the canopy", "(kg/m^2*int)", &net_snow);
+ declvar("net_p", TDim::NHRU, "net rain + net snow", "(kg/m^2*int)", &net_p);
+ decllocal("delunld_wind_int", TDim::NHRU, "solid snow unloading from the canopy induced by wind", "(kg/m^2*int)", &delunld_wind_int);
+ decllocal("delunld_melt_int", TDim::NHRU, "canopy snow unloading rate due to melting", "(kg/m^2*int)", &delunld_melt_int);
+ decllocal("delunld_subl_int", TDim::NHRU, "canopy snow unloading due to sublimation", "(kg/m^2*int)", &delunld_subl_int);
+ decllocal("deldrip_veg_int", TDim::NHRU, "drip from canopy snowmelt and intercepted rainfall", "(kg/m^2*int)", &deldrip_veg_int);
+
+ declvar("delL", TDim::NHRU, "interval change in SWE", "(kg/m^2*int)", &delL);
+ declvar("delmelt_veg_day", TDim::NHRU, "daily snow melt", "(mm/d)", &delmelt_veg_day);
+ decllocal("cmlmelt_veg_day", TDim::NHRU, "daily snow melt accumulator", "(mm/d)", &cmlmelt_veg_day);
+ declvar("cmlmelt_veg", TDim::NHRU, "cumulative melt", "(mm)", &cmlmelt_veg);
+
+ declstatvar("m_s_veg", TDim::NHRU, "canopy snow's specific mass, liquid and snow", "(kg/m^2)", &m_s_veg);
+
+ declstatvar("max_liq_veg", TDim::NHRU, "maximum veg canopy snow's liquid mass", "(kg/m^2)", &max_liq_veg);
+
+ declstatvar("T_s_veg", TDim::NHRU, "average canopy snow temp", "(" + string(DEGREE_CELSIUS) + ")", &T_s_veg);
+
+ declstatvar("h2o_sat_veg", TDim::NHRU, "fraction of liquid H2O saturation (0 to 1.0)", "()", &h2o_sat_veg);
+ declvar("h2o_vol_veg", TDim::NHRU, "liquid h2o content as volume ratio: V_water/(V_snow - V_ice)", "()", &h2o_vol_veg);
+
+ declvar("h2o_sat_veg_snow", TDim::NHRU, "snowfall's % of liquid H2O saturation", "()", &h2o_sat_veg_snow);
+
+ declvar("precip_now_veg", TDim::NHRU, "precipitation occur for current timestep - 0/1", "()", &precip_now_veg);
+ declvar("T_rain_veg", TDim::NHRU, "rain's temp", "(" + string(DEGREE_CELSIUS) + ")", &T_rain_veg);
+ declvar("T_sf", TDim::NHRU, "snowfall's temp", "(" + string(DEGREE_CELSIUS) + ")", &T_sf);
+
+ decllocal("S_n_L", TDim::NHRU, "net solar radiation", "(W/m^2)", &S_n);
+ decllocal("I_lw_L", TDim::NHRU, "incoming longwave (thermal) rad ", "(W/m^2)", &I_lw);
+ decllocal("T_a_L", TDim::NHRU, "air temp", "(" + string(DEGREE_CELSIUS) + ")", &T_a);
+ decllocal("e_a_L", TDim::NHRU, "vapor pressure", "(Pa)", &e_a);
+ decllocal("u_L", TDim::NHRU, "wind speed", "(m/s)", &u);
+ decllocal("I_LW_atm", TDim::NHRU, "Downwelling longwave from the atmoshpere", "(W/m^2)", &I_LW_atm);
+ decllocal("I_LW_gnd", TDim::NHRU, "Upwelling longwave from the ground", "(W/m^2)", &I_LW_gnd);
+ decllocal("I_LW_cpy_2_cpy", TDim::NHRU, "Longwave from the canopy reflected off the surface back to the canopy", "(W/m^2)", &I_LW_cpy_2_cpy);
+ decllocal("I_LW_cpy", TDim::NHRU, "Incoming longwave radiation emitted from the canopy", "(W/m^2)", &I_LW_cpy);
+ decllocal("O_LW_cpysnow", TDim::NHRU, "Outgoing longwave radiation emitted from the canopy snow", "(W/m^2)", &O_LW_cpysnow);
+ decllocal("u_2_3rds", TDim::NHRU, "Wind speed at 2/3 canopy height for canopy snow energy balance and unloading", "(m s^-1)", &u_2_3rds);
+ decllocal("rel_z_u", TDim::NHRU, "Height of wind speed measurement relative to top of snowpack, only implemented if relative_hts[hh] == 1", "(m)", &rel_z_u);
+ decllocal("rel_z_T", TDim::NHRU, "Height of temp measurement relative to top of snowpack, only implemented if relative_hts[hh] == 1", "(m)", &rel_z_T);
+
+
+ decllocal("m_precip_L", TDim::NHRU, "specific mass of total precip", "(kg/m^2)", &m_precip);
+ declvar("rain_on_snow_veg", TDim::NHRU, "specific mass of rain in precip", "(kg/m^2)", &m_rain);
+ decllocal("m_snow_L", TDim::NHRU, "specific mass in snow in precip", "(kg/m^2)", &m_snow);
+ decllocal("T_pp_L", TDim::NHRU, "precip temp", "(" + string(DEGREE_CELSIUS) + ")", &T_pp);
+
+ decllocal("P_a", TDim::NHRU, "air pressure", "(Pa)", &P_a);
+
+ decllocal("m_precip_cum", TDim::NHRU, "cumulative specific mass of total precip", "(kg/m^2)", &m_precip_cum);
+ decllocal("m_rain_cum", TDim::NHRU, "cumulative specific mass of total rain", "(kg/m^2)", &m_rain_cum);
+ decllocal("m_snow_cum", TDim::NHRU, "cumulative specific mass of total snow", "(kg/m^2)", &m_snow_cum);
+ decllocal("E_s_cum", TDim::NHRU, "cumulative mass flux by evap into air from active layer", "(kg/m^2)", &E_s_cum);
+
+ decllocal("stop_no_snow", TDim::NHRU, "snow flag", "()", &stop_no_snow);
+ declparam("max_liq_veg_frac", TDim::NHRU, "[0.01]", "0.0001", "0.2", "max liquid h2o content as fraction of specific snow mass", "()", &max_liq_veg_frac);
+ declparam("Albedo_vegsnow", TDim::NHRU, "[0.6]", "0.6", "0.9", "Albedo_vegsnow", "()", &Albedo_vegsnow);
+ declparam("SW_to_LW_fn", TDim::NHRU, "[0.01]", "0.0001", "0.5", "dimensionless shortwave to longwave transfer efficiency function. 0.038 from Pomeroy et al., (2009) for marmot forced through the origin, alternative value is 0.023 from Fraser site.", "()", &SW_to_LW_fn);
+
+ declgetparam("*", "z_u", "()", &z_u); // height of wind measurement (m)
+ declgetparam("*", "z_T", "()", &z_T); // height of air temp & vapor pressure measurement(m)
+ declgetparam("*", "hru_elev", "()", &hru_elev); // altitude (m)
+ declgetparam("*", "basin_area", "()", &basin_area); // total basin area (km^2)
+ declgetparam("*", "hru_area", "()", &hru_area); // hru area (km^2)
+ declgetparam("*", "hru_rho_snow", "()", &rho_snow_X); // rho of falling snow (kg/m^3)
+ declgetparam("*", "Cc", "()", &Cc); // canopy coverage, (1-sky view fraction)
+ declgetparam("*", "LAI", "()", &LAI);
+ declgetparam("*", "Ht", "()", &Ht);
+ declparam("Lmax", TDim::NHRU, "[50]", "0", "100", "maximum canopy snow load", "(kg/m^2)", &Lmax);
+ declparam("CanopyWindSwitchCanSno", TDim::NHRU, "[0]", "0", "1", "Canopy wind model to use for wind induced unloading at 1/2 canopy height, 0 - for Cionco, 1 - for Prandtl-von Kármán log-linear relationship", "()", &CanopyWindSwitchCanSno);
+ declparam("MassUnloadingSwitch", TDim::NHRU, "[0]", "0", "1", "canopy snow ablation parameterization to use, 0 - Cebulski & Pomeroy 2025 ablation paper, 1- Andreadis 2009, 2 - Roesch2001 (enable HP98/Ellis2010 using original canopy clearing gap module)", "()", &MassUnloadingSwitch);
+ declgetparam("*", "relative_hts", "()", &relative_hts);
+ declgetparam("*", "inhibit_evap", "()", &inhibit_evap);
+
+ declgetparam("*", "Alpha_c", "()", &Albedo_veg); // canopy albedo
+ declgetvar("*", "Albedo", "()", &Albedo_surface);
+ declgetvar("*", "Tauc", "(r)", &Tauc);
+
+ declgetvar("*", "hru_t", "(" + string(DEGREE_CELSIUS) + ")", &T_a_X);
+ declgetvar("*", "hru_t", "(" + string(DEGREE_CELSIUS) + ")", &T_pp_X); // default precipitation temperature to air
+ declgetvar("*", "hru_ea", "(kPa)", &e_a_X);
+ declgetvar("*", "hru_u", "(m/s)", &u_X);
+ declgetvar("*", "T_s_0", "(" + string(DEGREE_CELSIUS) + ")", &T_s_0);
+ // declreadobs("obs_snow_load", TDim::NHRU, "Weighed tree canopy snow load", "(kg/m^2)", &obs_snow_load, HRU_OBS_misc);
+
+ declgetvar("*", "z_s", "(m)", &z_s); // height of surface snowpack used for relative height calculation
+ declgetvar("*", "intercepted_snow", "(kg/m^2)", &new_snow); // new snow intercepted in canopy before ablation processes have kicked in, from vector based module
+ declgetvar("*", "intercepted_rain", "(kg/m^2)", &new_rain); // new snow intercepted in canopy before ablation processes have kicked in, from vector based module
+ declgetvar("*", "throughfall_rain", "(kg/m^2)", &throughfall_rain); // throughfall of rain to be added with canopy drip
+ declgetvar("*", "throughfall_snow", "(kg/m^2)", &throughfall_snow); // throughfall of snow to be added with canopy snow unloading
+ declgetvar("*", "hru_evap", "(kg/m^2)", &hru_evap);
+
+ variation_set = VARIATION_0 + VARIATION_1;
+
+ declreadobs("Qsi", TDim::NHRU, "incident short-wave", "(W/m^2)", &Qsi, HRU_OBS_Q, true); // must check for NULL
+
+ variation_set = VARIATION_0;
+
+ declreadobs("Qli", TDim::NHRU, "incident long-wave", "(W/m^2)", &Qli, HRU_OBS_Q, true);
+
+ variation_set = VARIATION_2;
+
+ declgetvar("*", "QsiS_Var", "(W/m^2)", &Qsw_in_veg); // downwelling SW to slope simulated from annandale module or Slope_Qsi module depending on if obs Qsi are used, same as regular snobal but without the canopy tau applied
+
+ variation_set = VARIATION_1 + VARIATION_2;
+
+ declgetvar("*", "QliVt_Var", "(W/m^2)", &Qlw_out_atm); // downwelling longwave radiation from the atmosphere, adjusted for terrain emission as well
+
+ variation_set = VARIATION_ORG;
+
+}
+
+void ClassCanopySnowBalanceCRHM::init(void) {
+
+ ClassCanopySnowBalanceBase::init();
+ for (hh = 0; chkStruct(); ++hh) {
+ delmelt_veg_day[hh] = 0.0;
+ }
+}
+
+void ClassCanopySnowBalanceCRHM::run(void) { // executed every interval
+
+ long nstep = getstep() % Global::Freq;
+
+ if(getstep() == 1){ // beginning of model run. Handle initial state file problems
+ for (hh = 0; chkStruct(); ++hh) {
+ if(snow_h2o_veg[hh] <= 0)
+ vegsnowcover[hh] = 0;
+ else{
+ vegsnowcover[hh] = 1;
+ }
+ }
+ }
+
+ for (hh = 0; chkStruct(); ++hh) {
+
+ switch ((int)variation)
+ {
+ case VARIATION_ORG: // obs SW and obs LW
+ input_rec2[hh].S_n = Qsi[hh];
+ input_rec2[hh].I_lw = Qli[hh];
+
+ break;
+ case VARIATION_1: // obs SW and mod LW
+ input_rec2[hh].S_n = Qsi[hh];
+ I_LW_atm[hh] = (CAN_EMISSIVITY + (1.0 - CAN_EMISSIVITY) * (1.0 - SNOW_EMISSIVITY) * CAN_EMISSIVITY) * Qlw_out_atm[hh]; // ! downward atmospheric longwave radiation absorbed by the canopy (W m-2) from SUMM
+ input_rec2[hh].I_lw = I_LW_atm[hh];
+
+ break;
+ case VARIATION_2: // obs or mod SW (depends on annandale vs slope_qsi module selection) and mod LW
+ input_rec2[hh].S_n = Qsw_in_veg[hh]; // after CLASSIC just take the incoming solar to slope which is multiplied by 1 - canopy albedo later on. Differs slightly from class which uses incoming SW to horizontal surface where this is the SW to slope
+ I_LW_atm[hh] = (CAN_EMISSIVITY + (1.0 - CAN_EMISSIVITY) * (1.0 - SNOW_EMISSIVITY) * CAN_EMISSIVITY) * Qlw_out_atm[hh]; // ! downward atmospheric longwave radiation absorbed by the canopy (W m-2) from SUMM
+ input_rec2[hh].I_lw = I_LW_atm[hh];
+
+ break;
+ }
+
+ input_rec2[hh].T_a = T_a_X[hh] + FREEZE;
+ input_rec2[hh].e_a = e_a_X[hh]*1000;
+ input_rec2[hh].u = u_X[hh];
+
+ // handles non throughfall precip
+ m_precip[hh] = new_snow[hh] + new_rain[hh];
+ m_rain[hh] = new_rain[hh];
+ m_snow[hh] = new_snow[hh];
+
+ m_precip_cum[hh] += m_precip[hh]; // change
+ m_rain_cum[hh] += m_rain[hh];
+ m_snow_cum[hh] += m_snow[hh];
+
+ T_pp[hh] = T_pp_X[hh] + FREEZE;
+
+// clear interval values
+
+ delmelt_veg_int[hh] = 0.0;
+ delsub_veg_int[hh] = 0.0;
+ delunld_int[hh] = 0.0;
+ delunld_wind_int[hh] = 0.0;
+ delunld_melt_int[hh] = 0.0;
+ delunld_subl_int[hh] = 0.0;
+ deldrip_veg_int[hh] = 0.0;
+ delevap_veg_int[hh] = 0.0;
+
+ long Step = getstep();
+
+ // uncomment below to hop to specific time in debug
+
+ // string test = StandardConverterUtility::GetDateTimeInString(Global::DTnow);
+
+ // if (test == "10/24/2022 14:0") {
+ // // if (test == "3/26/2023 15:0") { // TOP OF THE HOUR IS ONE ZERO
+ // // if (test == "10/1/2021 0:15") {
+ // std::cout << "Breakpoint here: Date matched" << std::endl;
+ // }
+
+
+ if(getstep() > 1){ // Not first step
+
+ if (m_precip[hh] > 0.0)
+ {
+ stop_no_snow[hh] = 0;
+ precip_now_veg[hh] = true;
+ }
+ else
+ {
+ if (!(vegsnowcover[hh] || liq_h2o_veg[hh] > 0.0))
+ stop_no_snow[hh] = 1; // _do_tstep_veg will not execute
+ precip_now_veg[hh] = false;
+ }
+
+ do_data_tstep_veg(); // executes Snobal code only if new snow intercepted/snow is in canopy, will not execute for hrus w.o. canopy.
+
+ if (snow_h2o_veg[hh] < 1e-6 & snow_h2o_veg[hh] > 0.0){ // if very small amount of snow on canopy then sublimate it off, handles non-convergence of energy balance for small snow values
+
+ delsub_veg_int[hh] += snow_h2o_veg[hh];
+ snow_h2o_veg[hh] = 0.0;
+ m_s_veg[hh] = snow_h2o_veg[hh] + liq_h2o_veg[hh];
+ vegsnowcover[hh] = 0.0;
+ }
+
+ }
+
+ else if(m_precip[hh] > 0.0) {
+ CRHMException TExcept("Snobal - cannot handle precipitation during first day of model run", TExcept::WARNING);
+ LogError(TExcept);
+ }
+
+ Step = Step % Global::Freq;
+ if(Step == 1)
+ cmlmelt_veg_day[hh] = delmelt_veg_int[hh];
+ else
+ cmlmelt_veg_day[hh] += delmelt_veg_int[hh];
+
+ if(Step == 0){
+ delmelt_veg_day[hh] = cmlmelt_veg_day[hh];
+ cmlmelt_veg_day[hh] = 0.0;
+ }
+
+ input_rec1[hh].S_n = input_rec2[hh].S_n;
+ input_rec1[hh].I_lw = input_rec2[hh].I_lw;
+ input_rec1[hh].T_a = input_rec2[hh].T_a;
+ input_rec1[hh].e_a = input_rec2[hh].e_a;
+ input_rec1[hh].u = input_rec2[hh].u;
+
+ E_s_cum[hh] += delsub_veg_int[hh];
+ delL[hh] += m_s_veg[hh];
+ cmlmelt_veg[hh] += delmelt_veg_int[hh];
+
+ // net rain and snow falling to the subcanopy/open snowpack
+ net_rain[hh] = throughfall_rain[hh] + deldrip_veg_int[hh];
+ net_snow[hh] = throughfall_snow[hh] + delunld_int[hh];
+
+ net_p[hh] = net_rain[hh] + net_snow[hh];
+
+ // set assimilate observed snow load for begining of select events
+ // if(obs_snow_load[hh] < 9999){
+ // m_s_veg[hh] = obs_snow_load[hh];
+ // snow_h2o_veg[hh] = obs_snow_load[hh];
+ // liq_h2o_veg[hh] = 0.0; // could be incorrect for first timestep but will be small impact as all liq water assumed to drain on each timestep
+ // delmelt_veg_int[hh] = 0.0;
+ // delsub_veg_int[hh] = 0.0;
+ // delunld_int[hh] = 0.0;
+ // deldrip_veg[hh] = 0.0;
+ // }
+
+ }
+}
+
+void ClassCanopySnowBalanceCRHM::finish(bool good) { // only required for local storage and final output
+
+ ClassCanopySnowBalanceBase::finish(good);
+}
diff --git a/crhmcode/src/modules/ClassCanopySnowBalanceCRHM.h b/crhmcode/src/modules/ClassCanopySnowBalanceCRHM.h
new file mode 100644
index 000000000..3255b568e
--- /dev/null
+++ b/crhmcode/src/modules/ClassCanopySnowBalanceCRHM.h
@@ -0,0 +1,65 @@
+/**
+* Copyright 2022, CRHMcode's Authors or Contributors
+* This file is part of CRHMcode.
+*
+* CRHMcode is free software: you can redistribute it and/or modify it under
+* the terms of the GNU General Public License as published by the Free Software
+* Foundation, either version 3 of the License, or (at your option) any later
+* version.
+*
+* CRHMcode is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty
+* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+* See the GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License along with
+* CRHMcode. If not, see .
+*
+**/
+//created by Manishankar Mondal
+
+#include "ClassCanopySnowBalanceBase.h"
+
+class ClassCanopySnowBalanceCRHM : public ClassCanopySnowBalanceBase {
+
+ public:
+
+ ClassCanopySnowBalanceCRHM(string Name, string Version = "undefined", LMODULE Lvl = LMODULE::PROTO) : ClassCanopySnowBalanceBase(Name, Version, Lvl) {};
+
+ ClassCanopySnowBalanceCRHM* klone(string name) const;
+
+ double Qsi_{0.0};
+ double Qli_{0.0};
+
+// Observations
+
+ const double *Qsi{ NULL }; // incoming solar radiation (W/m^2)
+ const double *Qli{ NULL }; // incoming longwave (thermal) rad (W/m^2) long
+ const double *obs_snow_load { NULL };
+
+// variable climate-data inputs
+
+ const double *T_a_X{ NULL }; // air temp (C) hru_t
+ const double *T_pp_X{ NULL }; // precip temp (C) hru_t
+ const double *e_a_X{ NULL }; // vapor pressure (Pa) hru_ea
+ const double *u_X{ NULL }; // wind speed (m/sec) hru_u
+ const double *Qsw_in_veg{ NULL }; // downwelling shortwave radiation to the canopy (W m^-2)
+ const double *Qlw_out_atm{ NULL }; // downwelling longwave radiation from the atmosphere plus terrain (W m^-2)
+
+// variable precipitation inputs
+
+ const double *new_snow{ NULL }; // snow intercepted in the canopy before ablation (kg/m^2*int)
+ const double *new_rain{ NULL }; // rain intercepted in the canopy before ablation (kg/m^2*int)
+
+// parameters
+
+ const double *rho_snow_X{ NULL }; // density of snowfall (kg/m^3)
+
+ void decl(void);
+
+ void init(void);
+
+ void run(void);
+
+ void finish(bool good); // delete local storage used
+};
\ No newline at end of file
diff --git a/crhmcode/src/modules/ClassObs.cpp b/crhmcode/src/modules/ClassObs.cpp
index 0a845a419..4105c7a56 100644
--- a/crhmcode/src/modules/ClassObs.cpp
+++ b/crhmcode/src/modules/ClassObs.cpp
@@ -520,7 +520,10 @@ void Classobs::Harder(void) {
else
- ratio = 1.0/(1.0 + 2.50286*pow(0.125006, hru_icebulb));
+ ratio = 1.0/(1.0 + 2.50286*pow(0.125006, hru_icebulb));
+ if(ratio < 0.0001){ // added by AC to avoid erroneously going into mixed precip routines even if very small frac of rain
+ ratio = 0.0;
+ }
hru_snow[hh] = 0.0;
diff --git a/crhmcode/src/modules/ClassSnobalBase.cpp b/crhmcode/src/modules/ClassSnobalBase.cpp
index 9b836fb6c..ef4cf35b8 100644
--- a/crhmcode/src/modules/ClassSnobalBase.cpp
+++ b/crhmcode/src/modules/ClassSnobalBase.cpp
@@ -995,7 +995,7 @@ double ClassSnobalBase::g_snow(
void ClassSnobalBase::_advec(void) {
- if (precip_now) {
+ if (precip_now[hh]) {
M[hh] = (heat_stor(CP_WATER(T_rain[hh]), m_rain[hh], (T_rain[hh] - T_s_0[hh])) +
heat_stor(CP_ICE(T_snow[hh]), m_snow[hh], (T_snow[hh] - T_s_0[hh]))) / time_step[hh];
}
@@ -1209,9 +1209,12 @@ void ClassSnobalBase::_precip(void)
// ro_predict[hh] += m_rain[hh];
} // current precip
- // Add water in the snowcover to total liquid water.
+ // Add water in the snowcover to total liquid water.
else // no precip
h2o_total[hh] += h2o[hh];
+
+ // Uncomment below (comment lines above) to fix obvious bug and add water in the snowcover from prev timestep to total liquid water regardless of new precip or not.
+ // h2o_total[hh] += h2o[hh];
}
void ClassSnobalBase::_drift(void)
@@ -1841,145 +1844,6 @@ void ClassSnobalBase::_runoff(void) {
}
}
-double satw(
- double tk) { /* air temperature (K) */
-
- double x;
- double l10;
-
- if (tk <= 0.) {
- CRHMException TExcept("satw temperature <= 0.0", TExcept::TERMINATE);
- LogError(TExcept);
- }
-
- errno = 0;
- l10 = log(1.e1);
-
- x = -7.90298 * (BOIL / tk - 1.0f) + 5.02808f * log(BOIL / tk) / l10 -
- 1.3816e-7 * (pow(1.0e1, 1.1344e1f * (1.0 - tk / BOIL)) - 1.0f) +
- 8.1328e-3 * (pow(1.0e1, -3.49149f * (BOIL / tk - 1.0f)) - 1.0f) +
- log(SEA_LEVEL) / l10;
-
- x = pow(1.0e1f, x);
-
- if (errno) {
- CRHMException TExcept("satw: bad return from log or pow", TExcept::TERMINATE);
- LogError(TExcept);
- }
-
- return(x);
-}
-
-
-// psi-functions
-// code = SM momentum
-// SH sensible heat flux
-// SV latent heat flux
-
-static double
-psi(double zeta, // z/lo
- int code) // which psi function? (see above)
-{
- double x; // height function variable
- double result{};
-
- if (zeta > 0) // stable
- {
- if (zeta > 1)
- {
- zeta = 1;
- }
- result = -BETA_S * zeta;
- }
- else if (zeta < 0) // unstable
- {
- x = sqrt(sqrt(1.0 - BETA_U * zeta));
-
- switch (code)
- {
- case SM:
- result = 2.0 * log((1.0 + x) / 2.0) + log((1.0 + x * x) / 2.0) -
- 2.0 * atan(x) + M_PI_2;
- break;
-
- case SH:
- case SV:
- result = 2.0 * log((1.0 + x * x) / 2.0);
- break;
-
- default: // shouldn't reach
- CRHMException TExcept("psi-function code not of these: SM, SH, SV", TExcept::TERMINATE);
- LogError(TExcept);
- }
- }
- else //Zeta == 1, neutral
- {
- result = 0.0;
- }
-
- return (result);
-}
-
-
-
-double ClassSnobalBase::sati(double tk) { //* air temperature (K)
-
- double l10;
- double x;
-
- if (tk <= 0.0) {
- CRHMException TExcept("sati temperature <= 0.0", TExcept::TERMINATE);
- LogError(TExcept);
- // tk = FREEZE + 0.01;
- }
-
- if (tk > FREEZE) {
- x = satw(tk);
- return(x);
- }
-
- errno = 0;
- l10 = log(1.e1);
-
- x = pow(1.e1, -9.09718 * ((FREEZE / tk) - 1.0) - 3.56654 * log(FREEZE / tk) / l10 +
- 8.76793e-1 * (1.0 - (tk / FREEZE)) + log(6.1071) / l10);
-
- if (errno) {
- CRHMException TExcept("sati: bad return from log or pow", TExcept::TERMINATE);
- LogError(TExcept);
- }
-
- return(x * 1.e2);
-}
-/* ----------------------------------------------------------------------- */
-
-double ClassSnobalBase::ssxfr(
- double k1, /* layer 1's thermal conductivity (J / (m K sec)) */
- double k2, /* layer 2's " " */
- double t1, /* layer 1's average layer temperature (K) */
- double t2, /* layer 2's " " " */
- double d1, /* layer 1's thickness (m) */
- double d2) /* layer 2's " " */
-{
- double xfr;
-
- xfr = 2.0 * (k1 * k2 * (t2 - t1)) / ((k2 * d1) + (k1 * d2));
-
- return (xfr);
-}
-/* ----------------------------------------------------------------------- */
-
-double ClassSnobalBase::heat_stor(
- double cp, /* specific heat of layer (J/kg K) */
- double spm, /* layer specific mass (kg/m^2) */
- double tdif) /* temperature change (K) */
-{
- double stor;
-
- stor = cp * spm * tdif;
-
- return (stor);
-}
/* ----------------------------------------------------------------------- */
double ClassSnobalBase::new_tsno(
diff --git a/crhmcode/src/modules/ClassSnobalBase.h b/crhmcode/src/modules/ClassSnobalBase.h
index a3ec62a1e..8550a9915 100644
--- a/crhmcode/src/modules/ClassSnobalBase.h
+++ b/crhmcode/src/modules/ClassSnobalBase.h
@@ -242,12 +242,6 @@ class ClassSnobalBase : public ClassModule {
double new_tsno(double spm, double t0, double ccon);
- double heat_stor(double cp, double spm, double tdif);
-
- double sati(double tk);
-
- double ssxfr(double k1, double k2, double t1, double t2, double d1, double d2);
-
double efcon(double k, double t, double p);
int hle1(double press, double ta, double ts, double za, double ea, double es, double zq, double u, double zu,
diff --git a/crhmcode/src/modules/ClassSnobalCRHM.cpp b/crhmcode/src/modules/ClassSnobalCRHM.cpp
index 7ed265914..2e228f11c 100644
--- a/crhmcode/src/modules/ClassSnobalCRHM.cpp
+++ b/crhmcode/src/modules/ClassSnobalCRHM.cpp
@@ -118,7 +118,7 @@ void ClassSnobalCRHM::decl(void) {
decllocal("F_g_L", TDim::NHRU, "soil flux at depth z_g", "(W/m^2)", &F_g);
decllocal("m_precip_L", TDim::NHRU, "specific mass of total precip", "(kg/m^2)", &m_precip);
- declvar ("rain_on_snow", TDim::NHRU, "specific mass of rain in precip", "(kg/m^2)", &m_rain);
+ declvar("rain_on_snow", TDim::NHRU, "specific mass of rain in precip", "(kg/m^2)", &m_rain);
decllocal("m_snow_L", TDim::NHRU, "specific mass in snow in precip", "(kg/m^2)", &m_snow);
decllocal("m_drift_L", TDim::NHRU, "specific mass of drifting snow", "(kg/m^2)", &m_drift);
decllocal("m_subl_L", TDim::NHRU, "specific mass of drifting snow", "(kg/m^2)", &m_subl);
@@ -252,12 +252,22 @@ void ClassSnobalCRHM::run(void) { // executed every interval
SWE_change[hh] = -m_s[hh];
+ // // uncomment below to hop to specific time in debug
+
+ // string test = StandardConverterUtility::GetDateTimeInString(Global::DTnow);
+
+ // if (test == "6/15/2022 15:0") {
+ // // if (test == "3/26/2023 15:0") { // TOP OF THE HOUR IS ONE ZERO
+ // // if (test == "10/1/2021 0:15") {
+ // std::cout << "Breakpoint here: Date matched" << std::endl;
+ // }
+
switch (variation){
case VARIATION_ORG:
input_rec2[hh].S_n = Qsi[hh]*(1.0 - Albedo[hh]);
input_rec2[hh].I_lw = Qli[hh];
break;
- case VARIATION_1:
+ case VARIATION_1: // this is the default
input_rec2[hh].S_n = Qsisn_Var[hh]*(1.0 - Albedo[hh]);
input_rec2[hh].I_lw = Qlisn_Var[hh];
break;
@@ -286,16 +296,19 @@ void ClassSnobalCRHM::run(void) { // executed every interval
else
input_rec2[hh].T_g = T_g_X[hh] + CRHM_constants::Tm;
- if(m_snow_X[hh] == 0)
+ if(m_snow_X[hh] == 0){
m_snow[hh] = 0.0;
+ }
if(snowcover[hh]){
m_snow[hh] = m_snow_X[hh];
}
else {
- if(m_snow_X[hh] > 0.0)
+ if(m_snow_X[hh] > 0.0){
snow_store[hh] += m_snow_X[hh];
- m_snow[hh] = 0.0;
+ } else {
+ m_snow[hh] = 0.0;
+ }
}
if((nstep == 1 && snow_store[hh] > 0.0) || snow_store[hh] > 1.0){
diff --git a/crhmcode/src/modules/newmodules/NewModules.cpp b/crhmcode/src/modules/newmodules/NewModules.cpp
index b5c115a35..3ecd44ec7 100644
--- a/crhmcode/src/modules/newmodules/NewModules.cpp
+++ b/crhmcode/src/modules/newmodules/NewModules.cpp
@@ -130,6 +130,9 @@
#include "../Classshared.h" //added by Manishankar Mondal
#include "../ClassNOP.h" //added by Manishankar Mondal
+#include "../ClassCRHMCanopyVectorBased.h" //added by Alex Cebulski
+#include "../ClassCanopySnowBalanceBase.h" //added by Alex Cebulski
+#include "../ClassCanopySnowBalanceCRHM.h" //added by Alex Cebulski
//---------------------------------------------------------------------------
@@ -158,7 +161,7 @@ bool RELEASE = false;
void MoveModulesToGlobal(string DLLName)
{
-
+
DLLModules.AddModule(new Classshared("Shared", "10/25/10", LMODULE::CUSTOM)); // essential for parameter screen
DLLModules.AddModule(new ClassNOP("NOP", "05/20/16", LMODULE::ADVANCE)); // essential for parameter screen
DLLModules.AddModule(new Classbasin("basin", "02/24/12", LMODULE::BASIC));
@@ -211,6 +214,9 @@ void MoveModulesToGlobal(string DLLName)
DLLModules.AddModule(new ClassCRHMCanopy("Canopy", "04/05/22", LMODULE::ADVANCE));
DLLModules.AddModule(new ClassCRHMCanopyClearing("CanopyClearing", "04/05/22", LMODULE::ADVANCE));
DLLModules.AddModule(new ClassCRHMCanopyClearingGap("CanopyClearingGap", "04/05/22", LMODULE::ADVANCE));
+ DLLModules.AddModule(new ClassCRHMCanopyVectorBased("CanopyVectorBased", "11/12/24", LMODULE::ADVANCE));
+ DLLModules.AddModule(new ClassCanopySnowBalanceBase("CanopySnowBalanceBase", "02/07/25", LMODULE::ADVANCE));
+ DLLModules.AddModule(new ClassCanopySnowBalanceCRHM("CanopySnowBalanceCRHM", "02/07/25", LMODULE::ADVANCE));
DLLModules.AddModule(new ClassNeedle("NeedleLeaf", "04/05/22", LMODULE::ADVANCE));
DLLModules.AddModule(new Classwalmsley_wind("walmsley_wind", "07/30/08", LMODULE::ADVANCE));
DLLModules.AddModule(new ClassXG("XG", "04/05/22", LMODULE::ADVANCE));
diff --git a/crhmcode/src/modules/newmodules/SnobalDefines.h b/crhmcode/src/modules/newmodules/SnobalDefines.h
index c44ab33a6..7e7780c93 100644
--- a/crhmcode/src/modules/newmodules/SnobalDefines.h
+++ b/crhmcode/src/modules/newmodules/SnobalDefines.h
@@ -38,7 +38,7 @@
#define RGAS 8.31432e3
/*
-* specific humidity
+* specific humidity (-) / (kg/kg)
*
* e = vapor pressure
* P = pressure (same units as e)
@@ -148,7 +148,8 @@
/*
* thermal emissivity of snow
*/
-#define SNOW_EMISSIVITY 0.99 // changed from 0.98 03/27/15
+#define SNOW_EMISSIVITY 0.99 // changed from 0.98 03/27/15 as in JULES
+#define CAN_EMISSIVITY 0.98 // as in JULES
/*
* Macros
@@ -273,6 +274,7 @@
* Von Karman constant
*/
#define VON_KARMAN 3.5e-1
+#define VON_KARMAN2 0.41
/*
* virtual temperature, i.e. the fictitious temperature that air must
@@ -295,3 +297,197 @@
#define MIN_SNOW_TEMP -75
#endif
+
+inline double satw(double tk) { // air temperature (K)
+ double x;
+ double l10;
+
+ if (tk <= 0.) {
+ CRHMException TExcept("satw temperature <= 0.0", TExcept::TERMINATE);
+ LogError(TExcept);
+ }
+
+ errno = 0;
+ l10 = log(1.e1);
+
+ x = -7.90298 * (BOIL / tk - 1.0) + 5.02808 * log(BOIL / tk) / l10 -
+ 1.3816e-7 * (pow(1.0e1, 1.1344e1 * (1.0 - tk / BOIL)) - 1.0) +
+ 8.1328e-3 * (pow(1.0e1, -3.49149 * (BOIL / tk - 1.0)) - 1.0) +
+ log(SEA_LEVEL) / l10;
+
+ x = pow(1.0e1, x);
+
+ if (errno) {
+ CRHMException TExcept("satw: bad return from log or pow", TExcept::TERMINATE);
+ LogError(TExcept);
+ }
+
+ return x;
+}
+
+// psi-functions
+// code = SM momentum
+// SH sensible heat flux
+// SV latent heat flux
+
+inline double
+psi(double zeta, // z/lo
+ int code) // which psi function? (see above)
+{
+ double x; // height function variable
+ double result{};
+
+ if (zeta > 0) // stable
+ {
+ if (zeta > 1)
+ {
+ zeta = 1;
+ }
+ result = -BETA_S * zeta;
+ }
+ else if (zeta < 0) // unstable
+ {
+ x = sqrt(sqrt(1.0 - BETA_U * zeta));
+
+ switch (code)
+ {
+ case SM:
+ result = 2.0 * log((1.0 + x) / 2.0) + log((1.0 + x * x) / 2.0) -
+ 2.0 * atan(x) + M_PI_2;
+ break;
+
+ case SH:
+ case SV:
+ result = 2.0 * log((1.0 + x * x) / 2.0);
+ break;
+
+ default: // shouldn't reach
+ CRHMException TExcept("psi-function code not of these: SM, SH, SV", TExcept::TERMINATE);
+ LogError(TExcept);
+ }
+ }
+ else //Zeta == 1, neutral
+ {
+ result = 0.0;
+ }
+
+ return (result);
+}
+
+
+
+inline double sati(double tk) { //* air temperature (K)
+
+ double l10;
+ double x;
+
+ if (tk <= 0.0) {
+ CRHMException TExcept("sati temperature <= 0.0", TExcept::TERMINATE);
+ LogError(TExcept);
+ // tk = FREEZE + 0.01;
+ }
+
+ if (tk > FREEZE) {
+ x = satw(tk);
+ return(x);
+ }
+
+ errno = 0;
+ l10 = log(1.e1);
+
+ x = pow(1.e1, -9.09718 * ((FREEZE / tk) - 1.0) - 3.56654 * log(FREEZE / tk) / l10 +
+ 8.76793e-1 * (1.0 - (tk / FREEZE)) + log(6.1071) / l10);
+
+ if (errno) {
+ CRHMException TExcept("sati: bad return from log or pow", TExcept::TERMINATE);
+ LogError(TExcept);
+ }
+
+ return(x * 1.e2);
+}
+/* ----------------------------------------------------------------------- */
+
+inline double ssxfr(
+ double k1, /* layer 1's thermal conductivity (J / (m K sec)) */
+ double k2, /* layer 2's " " */
+ double t1, /* layer 1's average layer temperature (K) */
+ double t2, /* layer 2's " " " */
+ double d1, /* layer 1's thickness (m) */
+ double d2) /* layer 2's " " */
+{
+ double xfr;
+
+ xfr = 2.0 * (k1 * k2 * (t2 - t1)) / ((k2 * d1) + (k1 * d2));
+
+ return (xfr);
+}
+/* ----------------------------------------------------------------------- */
+
+inline double heat_stor(
+ double cp, /* specific heat of layer (J/kg K) */
+ double spm, /* layer specific mass (kg/m^2) */
+ double tdif) /* temperature change (K) */
+{
+ double stor;
+
+ stor = cp * spm * tdif;
+
+ return (stor);
+}
+
+/* ----------------------------------------------------------------------- */
+
+inline double lambda(double t) // Latent heat of vaporization (mJ/(kg DEGREE_CELSIUS))
+{
+ return( 2.501 - 0.002361 * t );
+}
+
+inline double delta(double t) // Slope of sat vap p vs t, kPa/DEGREE_CELSIUS
+{
+ if (t > 0.0)
+ return (2504.0 * exp(17.27 * t / (t + 237.3)) / sqr(t + 237.3));
+ else
+ return (3549.0 * exp(21.88 * t / (t + 265.5)) / sqr(t + 265.5));
+}
+
+inline double gamma(double Pa, double t) // Psychrometric constant (kPa/DEGREE_CELSIUS)
+{
+ return (0.00163 * Pa / lambda(t)); // lambda (mJ/(kg DEGREE_CELSIUS))
+}
+
+inline double adst_wind_cpy_top(
+ double veg_ht, /* Height of vegetation (m) */
+ double uz, /* Wind speed at height z (m/s) */
+ double z, /* Height of wind speed measurement (m) */
+
+ // output
+ double& u_veg_ht /* Wind speed at canopy top (m/s)*/
+ )
+{
+ if(z >= veg_ht){
+ u_veg_ht = uz;
+ } else if (veg_ht - 2.0 / 3.0 * z > 0.0) {
+ u_veg_ht = uz * log((veg_ht - 2.0 / 3.0 * z) / 0.123 * z) / log((z - 2.0 / 3.0 * z) / 0.123 * z);
+ } else {
+ u_veg_ht = 0.0;
+ }
+ return(u_veg_ht);
+}
+
+inline double cionco_canopy_wind_spd(
+ double veg_ht, /* Height of vegetation (m) */
+ double u_veg_ht, /* Wind speed at height veg ht (m/s) */
+ double target_ht, /* target height of output wind speed (m) */
+
+ // output
+ double &u_target_ht /* Wind speed at target height (m/s)*/
+)
+{
+ double A = 2.4338 + 3.45753 * exp(-u_veg_ht); /* Modified Cionco wind model from Parviainen & Pomeroy 2000 for mature forest*/
+ //double A = 2.97 + 3.2 * exp(-u_veg_ht); /* Modified Cionco wind model from Parviainen & Pomeroy 2000 for regen forest*/
+
+ u_target_ht = u_veg_ht * exp(A * (target_ht / veg_ht - 1.0)); /* calculates canopy windspd */
+
+ return (u_target_ht);
+}
+