From a007c6931e3a4637daa101d49b6bec926c6c5c18 Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Tue, 12 Nov 2024 11:26:14 -0800 Subject: [PATCH 01/15] Added new templates to add snow interception routine --- crhmcode/Makefile | 4 +- crhmcode/prj/badlake.prj | 2 +- crhmcode/src/modules/CMakeLists.txt | 1 + .../modules/ClassCRHMCanopyVectorBased.cpp | 700 ++++++++++++++++++ .../src/modules/ClassCRHMCanopyVectorBased.h | 139 ++++ 5 files changed, 843 insertions(+), 3 deletions(-) create mode 100644 crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp create mode 100644 crhmcode/src/modules/ClassCRHMCanopyVectorBased.h 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/prj/badlake.prj b/crhmcode/prj/badlake.prj index ead41fd2f..b849c63f3 100644 --- a/crhmcode/prj/badlake.prj +++ b/crhmcode/prj/badlake.prj @@ -11,7 +11,7 @@ Macros: ###### Observations: ###### -C:\Users\jhs507\repos\crhmcode\crhmcode\obs\Badlake73_76.obs +/home/alex/Documents/code/crhmcode/crhmcode/obs/Badlake73_76.obs ###### Dates: ###### diff --git a/crhmcode/src/modules/CMakeLists.txt b/crhmcode/src/modules/CMakeLists.txt index e9a5ed534..7c4cae6ed 100644 --- a/crhmcode/src/modules/CMakeLists.txt +++ b/crhmcode/src/modules/CMakeLists.txt @@ -22,6 +22,7 @@ add_library(modules OBJECT 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..dce12ae5b --- /dev/null +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp @@ -0,0 +1,700 @@ +/** +* 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 "ClassCRHMCanopyClearingGap.h" +#include "../core/GlobalDll.h" +#include "../core/ClassCRHM.h" +#include "newmodules/SnobalDefines.h" + + +using namespace CRHM; + +ClassCRHMCanopyClearingGap* ClassCRHMCanopyClearingGap::klone(string name) const{ + return new ClassCRHMCanopyClearingGap(name); +} + +void ClassCRHMCanopyClearingGap::decl(void) { + + Description = "'Prototype all season canopy/clearing module. Calculates short, long and all-wave radiation components at the snow surface.' \ + '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_ea", "(kPa)", &hru_ea); + + declgetvar("*", "hru_snow", "(mm/int)", &hru_snow); + + declgetvar("*", "hru_rain", "(mm/int)", &hru_rain); + + declgetvar("*", "hru_evap", "(mm/int)", &hru_evap); + + declgetvar("*", "SolAng", "(r)", &SolAng); + + declgetvar("*", "cosxs", "(r)", &cosxs); + + declgetvar("*", "cosxsflat", "(r)", &cosxsflat); + + declgetvar("*", "Qdfo", "(W/m^2)", &Qdfo); + + declgetvar("*", "Albedo", "()", &Albedo); + + declputvar("*", "SWE", "(mm)", &SWE); + +// declared observations + + declobs("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 transmissivity", "", &Tauc); + + decllocal("Pa", TDim::NHRU, "Average surface pressure", "(kPa)", &Pa); + + declvar("ra", TDim::NHRU, "resistance", "(s/m)", &ra); + + declvar("drip_cpy", TDim::NHRU, "canopy drip", "(mm/int)", &drip_Cpy); + + declvar("direct_rain", TDim::NHRU, "direct rainfall through canopy", "(mm/int)", &direct_rain); + + declvar("net_rain", TDim::NHRU, " direct_rain + drip", "(mm/int)", &net_rain); + + declstatdiag("cum_net_rain", TDim::NHRU, "cumulative direct_rain + drip", "(mm)", &cum_net_rain); + + declvar("pot_subl_cpy", TDim::NHRU, "dimensionless canopy snow sublimation rate aka potential sublimation rate to be multiplied by canopy snow load", "(s-1)", &pot_subl_cpy); + + declvar("Subl_Cpy", TDim::NHRU, "canopy snow sublimation", "(mm/int)", &Subl_Cpy); + + declstatdiag("cum_Subl_Cpy", TDim::NHRU, "cumulative canopy snow sublimation", "(mm)", &cum_Subl_Cpy); + + decldiag("Pevap", TDim::NHRU, "used when ground is snow covered to calculate canopy evaporation (Priestley-Taylor)", "(mm)", &Pevap); + + declstatvar("rain_load", TDim::NHRU, "canopy rain load", "(mm)", &rain_load); + + declstatvar("Snow_load", TDim::NHRU, "canopy snow load (timetep start)", "(mm)", &Snow_load); + + declvar("direct_snow", TDim::NHRU, "snow 'direct' through canopy", "(mm/int)", &direct_snow); + + declvar("SUnload", TDim::NHRU, "unloaded canopy snow", "(mm)", &SUnload); + + declvar("SUnload_H2O", TDim::NHRU, "unloaded canopy snow as water", "(mm)", &SUnload_H2O); + + declstatdiag("cum_SUnload_H2O", TDim::NHRU, "Cumulative unloaded canopy snow as water", "(mm)", &cum_SUnload_H2O); + + declstatdiag("cum_SUnload", TDim::NHRU, "Cumulative unloaded canopy snow as snow", "(mm)", &cum_SUnload); + + declvar("net_snow", TDim::NHRU, "hru_snow minus interception", "(mm/int)", &net_snow); + + declstatdiag("cum_net_snow", TDim::NHRU, "Cumulative hru_snow minus interception", "(mm)", &cum_net_snow); + + declvar("net_p", TDim::NHRU, "total precipitation after interception", "(mm/int)", &net_p); + + decldiag("u_FHt", TDim::NHRU, "wind speed at forest top (z = FHt)", "(m/s)", &u_FHt); + + decldiag("Cc", TDim::NHRU, "Canopy coverage", "()", &Cc); + + declvar("intcp_evap", TDim::NHRU, "HRU Evaporation from interception", "(mm/int)", &intcp_evap); + + declstatdiag("cum_intcp_evap", TDim::NHRU, "Cumulative HRU Evaporation from interception", "(mm)", &cum_intcp_evap); + + +// parameters: + + declparam("inhibit_evap", TDim::NHRU, "[0]", "0", "1", "inhibit evaporation, 1 -> inhibit", "()", &inhibit_evap); + + declparam("basin_area", TDim::BASIN, "3", "1e-6", "1e+09", "total basin area", "(km^2)", &basin_area); + + declparam("hru_area", TDim::NHRU, "[1]", "1e-6", "1e+09", "hru area", "(km^2)", &hru_area); + + 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", "()", &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("LAI", TDim::NHRU, "[2.2]", "0.1", "20.0", "leaf-area-index", "()", &LAI); + + declparam("Sbar", TDim::NHRU, "[6.6]", "0.0", "100.0", "maximum canopy snow interception load", "(kg/m^2)", &Sbar); + + declparam("Zvent", TDim::NHRU, "[0.75]", "0.0", "1.0", "ventilation wind speed height (z/Ht)", "()", &Zvent); + + declparam("unload_t", TDim::NHRU, "[1.0]", "-10.0", "20.0", "if ice-bulb temp >= t : canopy snow is unloaded as snow", "(" + string(DEGREE_CELSIUS) + ")", &unload_t); + + declparam("unload_t_water", TDim::NHRU, "[4.0]", "-10.0", "20.0", "if ice-bulb temp >= t: canopy snow is unloaded as water", "(" + string(DEGREE_CELSIUS) + ")", &unload_t_water); + + declparam("CanopyClearing", TDim::NHRU, "[0]", "0", "2", "canopy - 0/clearing - 1/gap - 2", "()", &CanopyClearing); + + 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 ClassCRHMCanopyClearingGap::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 + + rain_load[hh] = 0.0; + Snow_load[hh] = 0.0; + + cum_net_snow[hh] = 0.0; + cum_net_rain[hh] = 0.0; + cum_Subl_Cpy[hh] = 0.0; + cum_intcp_evap[hh] = 0.0; + cum_SUnload[hh] = 0.0; + cum_SUnload_H2O[hh] = 0.0; + + if(Ht[hh] > Zwind[hh]){ + CRHMException TExcept(string("'" + Name + " (CanopyClearingGap)' Vegetation height greater than wind reference height, i.e. (Ht > Zwind)!").c_str(), TExcept::WARNING); + LogError(TExcept); + } + } +} + +void ClassCRHMCanopyClearingGap::run(void){ + + double Exposure, LAI_, Vf, Vf_, Kstar_H, Kd; + //double Tau; variable is unreferenced commenting out for now - jhs507 + + for (hh = 0; chkStruct(); ++hh){ + + switch (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; + } // switch + + net_rain[hh] = 0.0; + direct_rain[hh] = 0.0; + drip_Cpy[hh] = 0.0; + intcp_evap[hh] = 0.0; + + net_snow[hh] = 0.0; + direct_snow[hh] = 0.0; + SUnload[hh] = 0.0; + SUnload_H2O[hh] = 0.0; + Subl_Cpy[hh] = 0.0; + +// Canopy temperature is approximated by the air temperature. + + double T1 = hru_t[hh] + CRHM_constants::Tm; + + double rho = Pa[hh]*1000/(CRHM_constants::Rgas*T1); + + 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 (kg/kg) + + Ts[hh] = T1 + (CRHM_constants::emiss*(Qli_ - CRHM_constants::sbc*pow(T1, 4.0f)) + CRHM_constants::Ls*(q - Common::Qs(Pa[hh], T1))*rho/ra[hh])/ + (4*CRHM_constants::emiss*CRHM_constants::sbc*pow(T1, 3.0f) + (CRHM_constants::Cp + CRHM_constants::Ls*deltaX)*rho/ra[hh]); + + Ts[hh] -= CRHM_constants::Tm; + + if(Ts[hh] > 0.0 || SWE[hh] <= 0.0) + Ts[hh] = 0.0; + + switch(CanopyClearing[hh]){ + case 0: // canopy + + Exposure = Ht[hh] - Common::DepthofSnow(SWE[hh]); /* depths(m) SWE(mm) */ + if(Exposure < 0.0) + Exposure = 0.0; + + LAI_ = LAI[hh]*Exposure/Ht[hh]; + + Vf = 0.45 - 0.29*log(LAI[hh]); + + Vf_ = Vf + (1.0 - Vf)*sin((Ht[hh] - Exposure)/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]); + double limit = cosxsflat[hh]/cosxs[hh]; + if(limit > 2.0) + limit = 2.0; + Tauc[hh] = exp(-k[hh]*LAI_*limit); + } + else{ + k[hh] = 0.0; + Tauc[hh] = 0.0; + } + + Kstar_H = Qsi_*(1.0 - Alpha_c[hh] - Tauc[hh]*(1.0 - Albedo[hh])); + + Qlisn[hh] = Qli_*Vf_ + (1.0f - Vf_)*CRHM_constants::emiss_c*CRHM_constants::sbc*pow(T1, 4.0f) + B_canopy[hh]*Kstar_H; + + Qlisn_Var[hh] = Qlisn[hh]; + + Qsisn[hh] = Qsi_*Tauc[hh]; + + Qsisn_Var[hh] = Qsisn[hh]; + + 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]; + + break; + case 1: // clearing + + Qlisn[hh] = Qli_; + + Qlisn_Var[hh] = Qlisn[hh]; + + Qsisn[hh] = Qsi_; + + Qsisn_Var[hh] = Qsisn[hh]; + + 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]; + + break; + case 2: // gap + Exposure = Surrounding_Ht[hh] - Common::DepthofSnow(SWE[hh]); /* depths(m) SWE(mm) */ + if(Exposure < 0.0) + Exposure = 0.0; + + LAI_ = LAI[hh]*Exposure/Surrounding_Ht[hh]; + + Vf = 0.45 - 0.29*log(LAI[hh]); + + double Tau_d = Vf + (1.0 - Vf)*sin((Surrounding_Ht[hh] - Exposure)/Surrounding_Ht[hh]*M_PI_2); // previously Vf_ + +// calculate forest clearing sky view factor (Vgap) via Reifsnyder and Lull?s (1965) expression: + + double Vgap = sqr(sin(atan2(Gap_diameter[hh], 2.0*Surrounding_Ht[hh]))); + +// calculate beam pathlength correction (variable Gap_beam_corr) for gap: + + double Gap_beam_corr = 0; + if(Qsi_ > 0.0 && SolAng[hh] > 0.001){ + double cosxsLim = 3; + if(cosxs[hh] > 0.33) + cosxsLim = 1.0/cosxs[hh]; + + Gap_beam_corr = cosxsLim*Surrounding_Ht[hh]*(1.0/cos(SolAng[hh]) - Gap_diameter[hh]/(2.0*Surrounding_Ht[hh])/sin(SolAng[hh])); + if(Gap_beam_corr > 10.0) + Gap_beam_corr = 10.0; + else if(Gap_beam_corr < 0.0) + Gap_beam_corr = 0.0; + } +// calculate beam shortwave transmittance of the gap: + + double product = LAI[hh]*Gap_beam_corr; + if(product > 50) + product = 50; + + double Tau_b_gap = exp(-product); + + Kd = Qsi_*(1.0 - Alpha_c[hh] - Tau_b_gap*(1.0 - Albedo[hh])); + + Qlisn[hh] = Vgap*Qli_ + (1.0 - Vgap)*((Qli_*Tau_b_gap + (1.0 - Tau_b_gap)*CRHM_constants::emiss_c*CRHM_constants::sbc*pow(T1, 4.0f)) + B_canopy[hh]*Kd); + + Qlisn_Var[hh] = Qlisn[hh]; + + Qsisn[hh] = cosxs[hh]*Qdfo[hh]*Tau_b_gap + Vgap*(Qsi_ - Qdfo[hh]) + (1.0 - Vgap)*Tau_d*(Qsi_ - Qdfo[hh]); + if(Qsisn[hh] < 0.0) + Qsisn[hh] = 0.0; + + Qsisn_Var[hh] = Qsisn[hh]; + + 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-Albedo[hh]); + + Qnsn_Var[hh] = Qnsn[hh]; + + break; + } // switch + + switch(CanopyClearing[hh]){ + case 0: { // canopy +//============================================================================== +// coupled forest snow interception and sublimation routine: +// after Hedstom & Pomeroy / Parviainen & Pomeroy: +// calculate maximum canopy snow load (L*): + + if(Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0){ // handle snow interception + + if(Sbar[hh] >= 0.0){ + double RhoS = 67.92 + 51.25* exp(hru_t[hh]/2.59); + double LStar = Sbar[hh]* (0.27 + 46.0/RhoS)* LAI[hh]; + + if(Snow_load[hh] > LStar){ // after increase in temperature + direct_snow[hh] = Snow_load[hh] - LStar; + Snow_load[hh] = LStar; + } + +// calculate intercepted snowload + + if(Ht[hh] - 2.0/3.0*Zwind[hh] > 0.0) + u_FHt[hh] = hru_u[hh]*log((Ht[hh] - 2.0/3.0*Zwind[hh] )/ 0.123*Zwind[hh])/log((Zwind[hh] - 2.0/3.0*Zwind[hh] )/ 0.123*Zwind[hh]); + else + u_FHt[hh] = 0.0; + + double I1 = 0.0; + + if(hru_snow[hh] > 0.0 && LStar > 0.0){ + if(fabs(hru_snow[hh]/LStar) < 50.0){ + if (u_FHt[hh] <= 1.0) // if wind speed at canopy top > 1 m/s + I1 = (LStar-Snow_load[hh])*(1.0-exp(-Cc[hh]*hru_snow[hh]/LStar)); + else + I1 = (LStar-Snow_load[hh])*(1.0-exp(-hru_snow[hh]/LStar)); + + if(I1 <= 0) + I1 = 0; + + Snow_load[hh] += I1; + } + + // calculate canopy snow throughfall before unloading: + + direct_snow[hh] += (hru_snow[hh] - I1); + } + +// calculate snow ventilation windspeed: + +//============================================================================= + const double gamma = 1.15; + const double AlbedoIce = 0.8; // albedo of ideal ice sphere + const double Radius = 5.0e-4; // radii of single 'ideal' ice sphere in, m) + const double KinVisc = 1.88e-5; // kinematic viscosity of air (Sask. avg. value) + const double ks = 0.0114; // snow shape coefficient for jack pine + const double Fract = 0.37; // fractal dimension of intercepted snow + const double ci = 2.102e-3; // heat capacity of ice (MJ/kg/K) + const double Hs = 2.838e6; // heat of sublimation (MJ/kg) +//============================================================================== + + double xi2 = 1-Zvent[hh]; + double windExt2 = (gamma * LAI[hh] * xi2); + double uVent = u_FHt[hh] * exp(-1 * windExt2); + +// calculate sublimation of intercepted snow from ideal intercepted ice sphere (500 microns diameter): + + double Alpha, A1, B1, C1, J, D, Lamb, Mpm, Nu, Nr, SStar, Sigma2; + + double Es = 611.15f * exp(22.452f*hru_t[hh]/(hru_t[hh] + 273.0f)); // {sat pressure} + + double SvDens = Es*PBSM_constants::M/(PBSM_constants::R*(hru_t[hh] + 273.0f)); // {sat density} + + Lamb = 6.3e-4*(hru_t[hh]+273.0) + 0.0673; // thermal conductivity of atmosphere + Nr = 2.0 * Radius * uVent / KinVisc; // Reynolds number + Nu = 1.79 + 0.606 * sqrt((double) Nr); // Nusselt number + SStar = M_PI * sqr(Radius) * (1.0f - AlbedoIce) * Qsi_; // SW to snow particle !!!! changed + A1 = Lamb * (hru_t[hh] + 273) * Nu; + B1 = Hs * PBSM_constants::M /(PBSM_constants::R * (hru_t[hh] + 273.0f))- 1.0; + J = B1/A1; + Sigma2 = hru_rh[hh]/100 -1; + D = 2.06e-5* pow((hru_t[hh]+273.0f)/273.0f, -1.75f); // diffusivity of water vapour + C1 = 1.0/(D*SvDens*Nu); + + Alpha = 5.0; + Mpm = 4.0 / 3.0 * M_PI * PBSM_constants::DICE * Radius * Radius * Radius; // 18Mar2022: remove Gamma Distribution Correction term, *(1.0 + 3.0/Alpha + 2.0/sqr(Alpha)); +// sublimation rate of single 'ideal' ice sphere: + + double Vs = (2.0* M_PI* Radius*Sigma2 - SStar* J)/(Hs* J + C1)/Mpm; + + pot_subl_cpy[hh] = Vs; // export the dimensionless sublimation rate (s-1) added by alex 2023-07-21 + +// snow exposure coefficient (Ce): + + double Ce; + if ((Snow_load[hh]/LStar) <= 0.0) + Ce = 0.07; + else + Ce = ks* pow((Snow_load[hh]/LStar), -Fract); + +// calculate 'potential' canopy sublimation: + + double Vi = Vs*Ce; + +// limit sublimation to canopy snow available and take sublimated snow away from canopy snow at timestep start + + Subl_Cpy[hh] = -Snow_load[hh]*Vi*Hs*Global::Interval*24*3600/Hs; // make W/m2 + if(Subl_Cpy[hh] > Snow_load[hh]){ + Subl_Cpy[hh] = Snow_load[hh]; + Snow_load[hh] = 0.0; + } + else{ + Snow_load[hh] -= Subl_Cpy[hh]; + if(Snow_load[hh] < 0.0) + Snow_load[hh] = 0.0; + } + +// calculate 'ice-bulb' temperature of intercepted snow: + + double IceBulbT = hru_t[hh] - (Vi* Hs/1e6/ci); + double Six_Hour_Divisor = Global::Freq/4.0; // used to unload over 6 hours + + const double U = -1 * log(0.678) / (24 * 7 * Global::Freq / 24); // weekly dimensionless unloading coefficient -> to CRHM time interval + // 21Mar2022 correction: invert the term 24/Global::Freq, use unloading rate coefficient U = -log(c)/t for snow unloading determined by inverse function of c = e^(-Ut) = 0.678 based on Eq. 14 in Hedstrom and Pomeroy (1998) + // determine whether canopy snow is unloaded: + + if(IceBulbT >= unload_t_water[hh]){ + drip_Cpy[hh] = Snow_load[hh]/Six_Hour_Divisor; + SUnload_H2O[hh] = drip_Cpy[hh]; + Snow_load[hh] -= SUnload_H2O[hh]; + cum_SUnload_H2O[hh] += SUnload_H2O[hh]; + } + else if(IceBulbT >= unload_t[hh]){ + SUnload[hh] = Snow_load[hh]/Six_Hour_Divisor; + Snow_load[hh] -= SUnload[hh]; + cum_SUnload[hh] += SUnload[hh]; + } + else if(IceBulbT < unload_t[hh]){ // has to be at least one interval. Trip on half step + SUnload[hh] = Snow_load[hh] * U; // the dimensionless unloading coefficient already /interval, 21Mar2022 correction: use unloading rate coefficient U + if(SUnload[hh] > Snow_load[hh]){ + SUnload[hh] = Snow_load[hh]; + Snow_load[hh] = 0.0; + } + else + Snow_load[hh] -= SUnload[hh]; + + cum_SUnload[hh] += SUnload[hh]; + } + +// calculate total sub-canopy snow: + + net_snow[hh] = direct_snow[hh] + SUnload[hh]; + + }// handle snow + break; + } // case canopy + case 1: // clearing + case 2: // gap + net_snow[hh] = hru_snow[hh]; + net_rain[hh] = hru_rain[hh]; + break; + } + } // switch + +// calculate horizontal canopy-coverage (Cc): + + double smax, Q; // cannot be in switch structure + + switch(CanopyClearing[hh]){ + + case 0: // canopy + + Cc[hh] = 0.29 * log(LAI[hh]) + 0.55; + if(Cc[hh] <= 0.0) + Cc[hh] = 0.0; + + smax = Cc[hh]*LAI[hh]*0.2; + +// Forest rain interception and evaporation model +// 'sparse' Rutter interception model (i.e. Valente 1997): + +// calculate direct throughfall: + + if(hru_rain[hh] > 0.0){ + + direct_rain[hh] = hru_rain[hh] * (1-Cc[hh]); + + // calculate rain accumulation on canopy before evap loss: + + if (rain_load[hh] + hru_rain[hh]*Cc[hh] > smax){ + drip_Cpy[hh] += (rain_load[hh] + hru_rain[hh]*Cc[hh] - smax); + rain_load[hh] = smax; + } + else + rain_load[hh] += hru_rain[hh]*Cc[hh]; + + } // if hru_rain[hh] > 0.0 + +// calculate 'actual evap' of water from canopy and canopy storage after evaporation:: + + if(rain_load[hh] > 0.0){ + if(inhibit_evap[hh] == 0){ // use Granger when no snowcover + if(rain_load[hh] >= hru_evap[hh]*Cc[hh]){ // (evaporation in mm) + intcp_evap[hh] = hru_evap[hh]*Cc[hh]; // + rain_load[hh] -= hru_evap[hh]*Cc[hh]; + } + else{ + intcp_evap[hh] = rain_load[hh]; + rain_load[hh] = 0.0; + } + } + else{ // use Priestley-Taylor when snowcover + Q = Qsi_*86400/Global::Freq/1e6/lambda(hru_t[hh]); // convert w/m2 to mm/m2/int + + if(Qsi_ > 0.0) + Pevap[hh] = 1.26*delta(hru_t[hh])*Q/(delta(hru_t[hh]) + gamma(Pa[hh], hru_t[hh])); + else + Pevap[hh] = 0.0; + + if(rain_load[hh] >= Pevap[hh]*Cc[hh]){ // (evaporation in mm) + intcp_evap[hh] = Pevap[hh]*Cc[hh]; // check + rain_load[hh] -= Pevap[hh]*Cc[hh]; + } + else{ + intcp_evap[hh] = rain_load[hh]; // check + rain_load[hh] = 0.0; + } + } + } // if rain_load[hh] > 0.0 + +// cumulative amounts.... + + net_rain[hh] = direct_rain[hh] + drip_Cpy[hh]; + + cum_intcp_evap[hh] += intcp_evap[hh]; + cum_Subl_Cpy[hh] += Subl_Cpy[hh]; + break; + default : // clearing and gap + break; + } // switch + + net_p[hh] = net_rain[hh] + net_snow[hh]; + cum_net_rain[hh] += net_rain[hh]; + cum_net_snow[hh] += net_snow[hh]; + } // end for +} + +void ClassCRHMCanopyClearingGap::finish(bool good) { + for(hh = 0; chkStruct(); ++hh) { + LogMessageA(hh, string("'" + Name + " (CanopyClearingGap)' cum_net_snow (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_net_snow[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (CanopyClearingGap)' cum_net_rain (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_net_rain[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (CanopyClearingGap)' cum_Subl_Cpy (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_Subl_Cpy[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (CanopyClearingGap)' cum_intcp_evap (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_intcp_evap[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (CanopyClearingGap)' cum_SUnload_H2O (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_intcp_evap[hh], hru_area[hh], basin_area[0]); + LogDebug(" "); + } +} + +double ClassCRHMCanopyClearingGap::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)); +} + +double ClassCRHMCanopyClearingGap::lambda(double t) // Latent heat of vaporization (mJ/(kg DEGREE_CELSIUS)) +{ + return( 2.501 - 0.002361 * t ); +} + +double ClassCRHMCanopyClearingGap::gamma(double Pa, double t) // Psychrometric constant (kPa/DEGREE_CELSIUS) +{ + return( 0.00163 * Pa / lambda(t)); // lambda (mJ/(kg DEGREE_CELSIUS)) +} + +double ClassCRHMCanopyClearingGap::RHOa(double t, double ea, double Pa) // atmospheric density (kg/m^3) +{ + const double R = 2870; + return (1E4*Pa /(R*( 273.15 + t))*(1.0 - 0.379*(ea/Pa)) ); // +} diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h new file mode 100644 index 000000000..cbffdc792 --- /dev/null +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h @@ -0,0 +1,139 @@ +/** +* 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 ClassCRHMCanopyClearingGap:public ClassModule { + +public: + +ClassCRHMCanopyClearingGap(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 *drip_Cpy { NULL }; +double *thrufall_Cpy { NULL }; +double *netRain_Cpy { NULL }; +double *direct_rain { NULL }; +double *rain_load { NULL }; +double *Snow_load { NULL }; +double *direct_snow { NULL }; +double *SUnload { NULL }; +double *SUnload_H2O { NULL }; +double *Qnsn_Var { NULL }; +double *Qsisn_Var { NULL }; +double *Qlisn_Var { NULL }; + +double *net_rain { NULL }; +double *cum_net_rain { NULL }; +double *net_snow { NULL }; +double *cum_net_snow { NULL }; +double *net_p { NULL }; +double *intcp_evap { NULL }; +double *cum_intcp_evap { NULL }; +double *pot_subl_cpy { NULL }; +double *Subl_Cpy { NULL }; +double *cum_Subl_Cpy { NULL }; +double *cum_SUnload { NULL }; +double *cum_SUnload_H2O { NULL }; + +double *Cc { NULL }; +double *k { NULL }; +double *Tauc { NULL }; +double *Pa { NULL }; +double *ra { NULL }; +double *u_FHt { NULL }; +double *Pevap { NULL }; + +// variable inputs + +const double *hru_t { NULL }; +const double *hru_u { NULL }; +const double *hru_rh { NULL }; +const double *hru_ea { NULL }; +const double *Albedo { NULL }; + +const double *hru_snow { NULL }; +const double *hru_rain { NULL }; +const double *hru_evap { 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 *basin_area { NULL }; // [BASIN] +const double *hru_area { NULL }; +const double *hru_elev { NULL }; +const double *Ht { NULL }; +const double *LAI { NULL }; +const double *Sbar { NULL }; +const double *Z0snow { NULL }; +const double *Zref { NULL }; +const double *Zwind { NULL }; +const double *Zvent { NULL }; +const double *unload_t { NULL }; +const double *unload_t_water { 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 *inhibit_evap { NULL }; + + +void decl(void); +void init(void); +void run(void); +void finish(bool good); + +double delta(double t); // Slope of sat vap p vs t, kPa/DEGREE_CELSIUS +double gamma(double Pa, double t); // Psychrometric constant (kPa/DEGREE_CELSIUS) +double RHOa(double t, double ea, double Pa); // atmospheric density (kg/m^3) +double lambda(double t); // Latent heat of vaporization (mJ/(kg DEGREE_CELSIUS)) + +ClassCRHMCanopyClearingGap* klone(string name) const; +}; From 7edf49d3922cb40901bd11c436f39a20e87df959 Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Tue, 12 Nov 2024 16:16:52 -0800 Subject: [PATCH 02/15] Other changes to incorporate new snow interception param template. Appears to work using new files w old algorithm now to modify. --- .../modules/ClassCRHMCanopyVectorBased.cpp | 34 +++++++++---------- .../src/modules/ClassCRHMCanopyVectorBased.h | 6 ++-- .../src/modules/newmodules/NewModules.cpp | 2 ++ 3 files changed, 22 insertions(+), 20 deletions(-) diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp index dce12ae5b..93298676b 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp @@ -25,7 +25,7 @@ #include #include -#include "ClassCRHMCanopyClearingGap.h" +#include "ClassCRHMCanopyVectorBased.h" #include "../core/GlobalDll.h" #include "../core/ClassCRHM.h" #include "newmodules/SnobalDefines.h" @@ -33,11 +33,11 @@ using namespace CRHM; -ClassCRHMCanopyClearingGap* ClassCRHMCanopyClearingGap::klone(string name) const{ - return new ClassCRHMCanopyClearingGap(name); +ClassCRHMCanopyVectorBased* ClassCRHMCanopyVectorBased::klone(string name) const{ + return new ClassCRHMCanopyVectorBased(name); } -void ClassCRHMCanopyClearingGap::decl(void) { +void ClassCRHMCanopyVectorBased::decl(void) { Description = "'Prototype all season canopy/clearing module. Calculates short, long and all-wave radiation components at the snow surface.' \ 'Inputs are observations Qsi (W/m^2) and Qli (W/m^2), ' \ @@ -215,7 +215,7 @@ void ClassCRHMCanopyClearingGap::decl(void) { 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 ClassCRHMCanopyClearingGap::init(void) { +void ClassCRHMCanopyVectorBased::init(void) { nhru = getdim(TDim::NHRU); // transfers current # of HRU's to module @@ -234,13 +234,13 @@ void ClassCRHMCanopyClearingGap::init(void) { cum_SUnload_H2O[hh] = 0.0; if(Ht[hh] > Zwind[hh]){ - CRHMException TExcept(string("'" + Name + " (CanopyClearingGap)' Vegetation height greater than wind reference height, i.e. (Ht > Zwind)!").c_str(), TExcept::WARNING); + CRHMException TExcept(string("'" + Name + " (CanopyVectorBased)' Vegetation height greater than wind reference height, i.e. (Ht > Zwind)!").c_str(), TExcept::WARNING); LogError(TExcept); } } } -void ClassCRHMCanopyClearingGap::run(void){ +void ClassCRHMCanopyVectorBased::run(void){ double Exposure, LAI_, Vf, Vf_, Kstar_H, Kd; //double Tau; variable is unreferenced commenting out for now - jhs507 @@ -664,18 +664,18 @@ void ClassCRHMCanopyClearingGap::run(void){ } // end for } -void ClassCRHMCanopyClearingGap::finish(bool good) { +void ClassCRHMCanopyVectorBased::finish(bool good) { for(hh = 0; chkStruct(); ++hh) { - LogMessageA(hh, string("'" + Name + " (CanopyClearingGap)' cum_net_snow (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_net_snow[hh], hru_area[hh], basin_area[0]); - LogMessageA(hh, string("'" + Name + " (CanopyClearingGap)' cum_net_rain (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_net_rain[hh], hru_area[hh], basin_area[0]); - LogMessageA(hh, string("'" + Name + " (CanopyClearingGap)' cum_Subl_Cpy (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_Subl_Cpy[hh], hru_area[hh], basin_area[0]); - LogMessageA(hh, string("'" + Name + " (CanopyClearingGap)' cum_intcp_evap (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_intcp_evap[hh], hru_area[hh], basin_area[0]); - LogMessageA(hh, string("'" + Name + " (CanopyClearingGap)' cum_SUnload_H2O (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_intcp_evap[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (CanopyVectorBased)' cum_net_snow (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_net_snow[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (CanopyVectorBased)' cum_net_rain (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_net_rain[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (CanopyVectorBased)' cum_Subl_Cpy (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_Subl_Cpy[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (CanopyVectorBased)' cum_intcp_evap (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_intcp_evap[hh], hru_area[hh], basin_area[0]); + LogMessageA(hh, string("'" + Name + " (CanopyVectorBased)' cum_SUnload_H2O (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_intcp_evap[hh], hru_area[hh], basin_area[0]); LogDebug(" "); } } -double ClassCRHMCanopyClearingGap::delta(double t) // Slope of sat vap p vs t, kPa/DEGREE_CELSIUS +double ClassCRHMCanopyVectorBased::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)); @@ -683,17 +683,17 @@ double ClassCRHMCanopyClearingGap::delta(double t) // Slope of sat vap p vs t, k return(3549.0*exp( 21.88 * t/(t+265.5)) / sqr(t+265.5)); } -double ClassCRHMCanopyClearingGap::lambda(double t) // Latent heat of vaporization (mJ/(kg DEGREE_CELSIUS)) +double ClassCRHMCanopyVectorBased::lambda(double t) // Latent heat of vaporization (mJ/(kg DEGREE_CELSIUS)) { return( 2.501 - 0.002361 * t ); } -double ClassCRHMCanopyClearingGap::gamma(double Pa, double t) // Psychrometric constant (kPa/DEGREE_CELSIUS) +double ClassCRHMCanopyVectorBased::gamma(double Pa, double t) // Psychrometric constant (kPa/DEGREE_CELSIUS) { return( 0.00163 * Pa / lambda(t)); // lambda (mJ/(kg DEGREE_CELSIUS)) } -double ClassCRHMCanopyClearingGap::RHOa(double t, double ea, double Pa) // atmospheric density (kg/m^3) +double ClassCRHMCanopyVectorBased::RHOa(double t, double ea, double Pa) // atmospheric density (kg/m^3) { const double R = 2870; return (1E4*Pa /(R*( 273.15 + t))*(1.0 - 0.379*(ea/Pa)) ); // diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h index cbffdc792..363ac5c6f 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h @@ -20,11 +20,11 @@ #include "../core/ClassModule.h" -class ClassCRHMCanopyClearingGap:public ClassModule { +class ClassCRHMCanopyVectorBased:public ClassModule { public: -ClassCRHMCanopyClearingGap(string Name, string Version = "undefined", LMODULE Lvl = LMODULE::PROTO) : ClassModule(Name, Version, Lvl, 1000, " , QliVt_Var, QsiS_Var, QsiS_Var, QsiA_Var") {}; +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}; @@ -135,5 +135,5 @@ double gamma(double Pa, double t); // Psychrometric constant (kPa/DEGREE_CELSIUS double RHOa(double t, double ea, double Pa); // atmospheric density (kg/m^3) double lambda(double t); // Latent heat of vaporization (mJ/(kg DEGREE_CELSIUS)) -ClassCRHMCanopyClearingGap* klone(string name) const; +ClassCRHMCanopyVectorBased* klone(string name) const; }; diff --git a/crhmcode/src/modules/newmodules/NewModules.cpp b/crhmcode/src/modules/newmodules/NewModules.cpp index b5c115a35..b9e91e25b 100644 --- a/crhmcode/src/modules/newmodules/NewModules.cpp +++ b/crhmcode/src/modules/newmodules/NewModules.cpp @@ -117,6 +117,7 @@ #include "../ClassXGAyers.h" //added by Manishankar Mondal #include "../ClassCRHMCanopyClearing.h" //added by Manishankar Mondal #include "../ClassCRHMCanopyClearingGap.h" //added by Manishankar Mondal +#include "../ClassCRHMCanopyVectorBased.h" //added by Alex Cebulski #include "../ClassREWroute_stream.h" //added by Manishankar Mondal #include "../ClassICEflow.h" //added by Manishankar Mondal #include "../Classglacier.h" //added by Manishankar Mondal @@ -211,6 +212,7 @@ 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 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)); From ebddb2b3f35ddc74478d111e01f3ea0f7eabbaa4 Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Wed, 13 Nov 2024 15:24:08 -0800 Subject: [PATCH 03/15] Added vector based parameterization for snow interception, changed max canopy snow load for sublimation to constant --- .../modules/ClassCRHMCanopyVectorBased.cpp | 91 ++++++++++--------- .../src/modules/ClassCRHMCanopyVectorBased.h | 3 +- 2 files changed, 50 insertions(+), 44 deletions(-) diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp index 93298676b..0d72f99cd 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp @@ -168,6 +168,8 @@ void ClassCRHMCanopyVectorBased::decl(void) { 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) following Cebulski & Pomeroy vector based param.", "(m/s)", &u_1_third_Ht); + decldiag("Cc", TDim::NHRU, "Canopy coverage", "()", &Cc); declvar("intcp_evap", TDim::NHRU, "HRU Evaporation from interception", "(mm/int)", &intcp_evap); @@ -200,7 +202,7 @@ void ClassCRHMCanopyVectorBased::decl(void) { declparam("LAI", TDim::NHRU, "[2.2]", "0.1", "20.0", "leaf-area-index", "()", &LAI); - declparam("Sbar", TDim::NHRU, "[6.6]", "0.0", "100.0", "maximum canopy snow interception load", "(kg/m^2)", &Sbar); + declparam("Lmax", TDim::NHRU, "[50]", "0.0", "100.0", "maximum canopy snow interception load, currently just used for sublimation exposure coef. 50 based on max observed in Storck et al. 2002 and Cebulski & Pomeroy", "(kg/m^2)", &Lmax); declparam("Zvent", TDim::NHRU, "[0.75]", "0.0", "1.0", "ventilation wind speed height (z/Ht)", "()", &Zvent); @@ -421,50 +423,47 @@ void ClassCRHMCanopyVectorBased::run(void){ } // switch switch(CanopyClearing[hh]){ - case 0: { // canopy + case 0: // canopy //============================================================================== // coupled forest snow interception and sublimation routine: -// after Hedstom & Pomeroy / Parviainen & Pomeroy: -// calculate maximum canopy snow load (L*): - - if(Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0){ // handle snow interception - - if(Sbar[hh] >= 0.0){ - double RhoS = 67.92 + 51.25* exp(hru_t[hh]/2.59); - double LStar = Sbar[hh]* (0.27 + 46.0/RhoS)* LAI[hh]; - - if(Snow_load[hh] > LStar){ // after increase in temperature - direct_snow[hh] = Snow_load[hh] - LStar; - Snow_load[hh] = LStar; - } - -// calculate intercepted snowload - - if(Ht[hh] - 2.0/3.0*Zwind[hh] > 0.0) - u_FHt[hh] = hru_u[hh]*log((Ht[hh] - 2.0/3.0*Zwind[hh] )/ 0.123*Zwind[hh])/log((Zwind[hh] - 2.0/3.0*Zwind[hh] )/ 0.123*Zwind[hh]); - else - u_FHt[hh] = 0.0; - - double I1 = 0.0; - - if(hru_snow[hh] > 0.0 && LStar > 0.0){ - if(fabs(hru_snow[hh]/LStar) < 50.0){ - if (u_FHt[hh] <= 1.0) // if wind speed at canopy top > 1 m/s - I1 = (LStar-Snow_load[hh])*(1.0-exp(-Cc[hh]*hru_snow[hh]/LStar)); - else - I1 = (LStar-Snow_load[hh])*(1.0-exp(-hru_snow[hh]/LStar)); +// after Cebulski & Pomeroy 2025: + +if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) +{ // calculate increase in Snow_load and direct_snow if we are in canopy (i.e., Cc > 0) + const double k_cp = -20; // rate of increase of the sigmoidal curve below + const double v_snow = 0.8; // terminal fall velocity of snowfall taken from Isyumov, 1971 + const double alpha = 0.836; // based on Cebulski & Pomeroy plot scale observations at Fortress Mountain PWL and Forest Tower plots. + double Cp = 0; // leaf contact area (Cp) based on trajectory angle + double IP = 0; // interception efficiency (IP) + double dL = 0; // change in canopy snow load + + if (hru_u[hh] > 0 && Cc[hh] < 1) + { // increase leaf contact area (Cp) based on wind speed and canopy coverage (Cc) + double Ht_1_third = Ht[hh] * (1 / 3); + double Cp_inc = 0; + if (Ht_1_third - (2.0 / 3.0) * Zwind[hh] > 0.0) + { + u_1_third_Ht[hh] = hru_u[hh] * log((Ht_1_third - (2.0 / 3.0) * Zwind[hh]) / 0.123 * Zwind[hh]) / log((Zwind[hh] - 2.0 / 3.0 * Zwind[hh]) / 0.123 * Zwind[hh]); + double snow_traj_angle = atan(u_1_third_Ht[hh] / v_snow); // in radians + Cp_inc = (1 - Cc[hh]) / (1 + exp(-k_cp * (sin(snow_traj_angle) - (1 - Cc[hh])))); // fractional increas in leaf contact area (Cp) based on horizontal trajectory. This is modified from Cebulski & Pomeroy snow interception paper. + } + Cp = Cc[hh] + Cp_inc; // calculated leaf contact area (Cp) based on trajectory angle + } + else + { + Cp = Cc[hh]; // use leaf contact area from nadir i.e., canopy coverage + } - if(I1 <= 0) - I1 = 0; + IP = Cp * alpha; // interception efficiency (IP) + dL = IP * hru_snow[hh]; // change in canopy snow load - Snow_load[hh] += I1; - } + Snow_load[hh] += dL; // calculate canopy snow throughfall before unloading: - direct_snow[hh] += (hru_snow[hh] - I1); - } + direct_snow[hh] += (1-IP) * hru_snow[hh]; +// Pomeroy et al. (1998) sublimation routine modified by Alex Cebulski to change maximum intercepted load from Hedstrom 1998 calculation to constant based on observations from several studies showing much higher intercepted loads. // calculate snow ventilation windspeed: //============================================================================= @@ -482,6 +481,15 @@ void ClassCRHMCanopyVectorBased::run(void){ double windExt2 = (gamma * LAI[hh] * xi2); double uVent = u_FHt[hh] * exp(-1 * windExt2); +// Calculate wind speed at canopy top for sublimation calculation, moved from initial interception chunk by Alex Cebulski November 2024 + + if(Ht[hh] - 2.0/3.0*Zwind[hh] > 0.0) + u_FHt[hh] = hru_u[hh]*log((Ht[hh] - 2.0/3.0*Zwind[hh] )/ 0.123*Zwind[hh])/log((Zwind[hh] - 2.0/3.0*Zwind[hh] )/ 0.123*Zwind[hh]); + else + u_FHt[hh] = 0.0; + + double I1 = 0.0; + // calculate sublimation of intercepted snow from ideal intercepted ice sphere (500 microns diameter): double Alpha, A1, B1, C1, J, D, Lamb, Mpm, Nu, Nr, SStar, Sigma2; @@ -512,10 +520,10 @@ void ClassCRHMCanopyVectorBased::run(void){ // snow exposure coefficient (Ce): double Ce; - if ((Snow_load[hh]/LStar) <= 0.0) + if ((Snow_load[hh]/Lmax[hh]) <= 0.0) Ce = 0.07; else - Ce = ks* pow((Snow_load[hh]/LStar), -Fract); + Ce = ks* pow((Snow_load[hh]/Lmax[hh]), -Fract); // calculate 'potential' canopy sublimation: @@ -569,16 +577,13 @@ void ClassCRHMCanopyVectorBased::run(void){ // calculate total sub-canopy snow: net_snow[hh] = direct_snow[hh] + SUnload[hh]; - - }// handle snow - break; + break; } // case canopy case 1: // clearing case 2: // gap net_snow[hh] = hru_snow[hh]; net_rain[hh] = hru_rain[hh]; break; - } } // switch // calculate horizontal canopy-coverage (Cc): diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h index 363ac5c6f..bf787a736 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h @@ -79,6 +79,7 @@ double *Tauc { NULL }; double *Pa { NULL }; double *ra { NULL }; double *u_FHt { NULL }; +double *u_1_third_Ht { NULL }; double *Pevap { NULL }; // variable inputs @@ -109,7 +110,7 @@ const double *hru_area { NULL }; const double *hru_elev { NULL }; const double *Ht { NULL }; const double *LAI { NULL }; -const double *Sbar { NULL }; +const double *Lmax { NULL }; const double *Z0snow { NULL }; const double *Zref { NULL }; const double *Zwind { NULL }; From af8a27f8d24ad21bd406747ea551e9d80b92b65f Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Wed, 13 Nov 2024 16:05:03 -0800 Subject: [PATCH 04/15] Change Cc calculation to limit to 1 --- crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp index 0d72f99cd..1bb80d8b0 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp @@ -594,7 +594,7 @@ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) case 0: // canopy - Cc[hh] = 0.29 * log(LAI[hh]) + 0.55; + Cc[hh] = std::min(1.0, 0.29 * std::log(LAI[hh]) + 0.55); // update by Alex so Cc is not greater than 1 if(Cc[hh] <= 0.0) Cc[hh] = 0.0; From e662bee902535adb2dea6c115eaa7892bbecb7d1 Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Wed, 13 Nov 2024 16:53:06 -0800 Subject: [PATCH 05/15] Modified IP alpha value so that it is a parameter --- crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp | 5 +++-- crhmcode/src/modules/ClassCRHMCanopyVectorBased.h | 1 + 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp index 1bb80d8b0..43e052554 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp @@ -204,6 +204,8 @@ void ClassCRHMCanopyVectorBased::decl(void) { declparam("Lmax", TDim::NHRU, "[50]", "0.0", "100.0", "maximum canopy snow interception load, currently just used for sublimation exposure coef. 50 based on max observed in Storck et al. 2002 and Cebulski & Pomeroy", "(kg/m^2)", &Lmax); + 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.", "()", &Zvent); + declparam("Zvent", TDim::NHRU, "[0.75]", "0.0", "1.0", "ventilation wind speed height (z/Ht)", "()", &Zvent); declparam("unload_t", TDim::NHRU, "[1.0]", "-10.0", "20.0", "if ice-bulb temp >= t : canopy snow is unloaded as snow", "(" + string(DEGREE_CELSIUS) + ")", &unload_t); @@ -432,7 +434,6 @@ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) { // calculate increase in Snow_load and direct_snow if we are in canopy (i.e., Cc > 0) const double k_cp = -20; // rate of increase of the sigmoidal curve below const double v_snow = 0.8; // terminal fall velocity of snowfall taken from Isyumov, 1971 - const double alpha = 0.836; // based on Cebulski & Pomeroy plot scale observations at Fortress Mountain PWL and Forest Tower plots. double Cp = 0; // leaf contact area (Cp) based on trajectory angle double IP = 0; // interception efficiency (IP) double dL = 0; // change in canopy snow load @@ -454,7 +455,7 @@ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) Cp = Cc[hh]; // use leaf contact area from nadir i.e., canopy coverage } - IP = Cp * alpha; // interception efficiency (IP) + IP = Cp * alpha[hh]; // interception efficiency (IP) dL = IP * hru_snow[hh]; // change in canopy snow load Snow_load[hh] += dL; diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h index bf787a736..56a35b239 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h @@ -111,6 +111,7 @@ const double *hru_elev { NULL }; const double *Ht { NULL }; const double *LAI { NULL }; const double *Lmax { NULL }; +const double *alpha { NULL }; const double *Z0snow { NULL }; const double *Zref { NULL }; const double *Zwind { NULL }; From 02d6e3721544397805b17f9f11466321c9d3be76 Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Wed, 13 Nov 2024 17:38:59 -0800 Subject: [PATCH 06/15] Changed Clca to a declvar so it can be displayed in the model output --- .../src/modules/ClassCRHMCanopyVectorBased.cpp | 14 ++++++++------ crhmcode/src/modules/ClassCRHMCanopyVectorBased.h | 1 + 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp index 43e052554..224ebfd23 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp @@ -172,6 +172,8 @@ void ClassCRHMCanopyVectorBased::decl(void) { decldiag("Cc", TDim::NHRU, "Canopy coverage", "()", &Cc); + declvar("Clca", TDim::NHRU, "Leaf contact area adjusted for hydrometeor trajectory angle.", "(s-1)", &Clca); + declvar("intcp_evap", TDim::NHRU, "HRU Evaporation from interception", "(mm/int)", &intcp_evap); declstatdiag("cum_intcp_evap", TDim::NHRU, "Cumulative HRU Evaporation from interception", "(mm)", &cum_intcp_evap); @@ -434,28 +436,28 @@ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) { // calculate increase in Snow_load and direct_snow if we are in canopy (i.e., Cc > 0) const double k_cp = -20; // rate of increase of the sigmoidal curve below const double v_snow = 0.8; // terminal fall velocity of snowfall taken from Isyumov, 1971 - double Cp = 0; // leaf contact area (Cp) based on trajectory angle + Clca[hh] = 0; // leaf contact area (Clca) based on trajectory angle double IP = 0; // interception efficiency (IP) double dL = 0; // change in canopy snow load if (hru_u[hh] > 0 && Cc[hh] < 1) - { // increase leaf contact area (Cp) based on wind speed and canopy coverage (Cc) + { // increase leaf contact area (Clca) based on wind speed and canopy coverage (Cc) double Ht_1_third = Ht[hh] * (1 / 3); double Cp_inc = 0; if (Ht_1_third - (2.0 / 3.0) * Zwind[hh] > 0.0) { u_1_third_Ht[hh] = hru_u[hh] * log((Ht_1_third - (2.0 / 3.0) * Zwind[hh]) / 0.123 * Zwind[hh]) / log((Zwind[hh] - 2.0 / 3.0 * Zwind[hh]) / 0.123 * Zwind[hh]); double snow_traj_angle = atan(u_1_third_Ht[hh] / v_snow); // in radians - Cp_inc = (1 - Cc[hh]) / (1 + exp(-k_cp * (sin(snow_traj_angle) - (1 - Cc[hh])))); // fractional increas in leaf contact area (Cp) based on horizontal trajectory. This is modified from Cebulski & Pomeroy snow interception paper. + Cp_inc = (1 - Cc[hh]) / (1 + exp(-k_cp * (sin(snow_traj_angle) - (1 - Cc[hh])))); // fractional increas in leaf contact area (Clca) based on horizontal trajectory. This is modified from Cebulski & Pomeroy snow interception paper. } - Cp = Cc[hh] + Cp_inc; // calculated leaf contact area (Cp) based on trajectory angle + Clca[hh] = Cc[hh] + Cp_inc; // calculated leaf contact area (Clca) based on trajectory angle } else { - Cp = Cc[hh]; // use leaf contact area from nadir i.e., canopy coverage + Clca[hh] = Cc[hh]; // use leaf contact area from nadir i.e., canopy coverage } - IP = Cp * alpha[hh]; // interception efficiency (IP) + IP = Clca[hh] * alpha[hh]; // interception efficiency (IP) dL = IP * hru_snow[hh]; // change in canopy snow load Snow_load[hh] += dL; diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h index 56a35b239..0599bfdc4 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h @@ -74,6 +74,7 @@ double *cum_SUnload { NULL }; double *cum_SUnload_H2O { NULL }; double *Cc { NULL }; +double *Clca { NULL }; double *k { NULL }; double *Tauc { NULL }; double *Pa { NULL }; From cad88ca2fe073d91dc0fe8d2a09ab46c6597b495 Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Fri, 15 Nov 2024 10:39:27 -0800 Subject: [PATCH 07/15] Fixed bug with alpha parameter --- crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp index 224ebfd23..118f6a695 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp @@ -206,7 +206,7 @@ void ClassCRHMCanopyVectorBased::decl(void) { declparam("Lmax", TDim::NHRU, "[50]", "0.0", "100.0", "maximum canopy snow interception load, currently just used for sublimation exposure coef. 50 based on max observed in Storck et al. 2002 and Cebulski & Pomeroy", "(kg/m^2)", &Lmax); - 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.", "()", &Zvent); + 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); From f37e75927d9bf533e905974197582a82429554c2 Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Fri, 15 Nov 2024 11:12:16 -0800 Subject: [PATCH 08/15] Change integers to floats so Ht_1_third is computed correctly --- crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp index 118f6a695..b7728d4b6 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp @@ -442,7 +442,7 @@ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) if (hru_u[hh] > 0 && Cc[hh] < 1) { // increase leaf contact area (Clca) based on wind speed and canopy coverage (Cc) - double Ht_1_third = Ht[hh] * (1 / 3); + double Ht_1_third = Ht[hh] * (1.0 / 3.0); double Cp_inc = 0; if (Ht_1_third - (2.0 / 3.0) * Zwind[hh] > 0.0) { From 403f1d06f79bd7e6fc0f6b2f9b2083359d85f7f3 Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Fri, 15 Nov 2024 11:53:45 -0800 Subject: [PATCH 09/15] fixed bug with Cp_inc function --- crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp index b7728d4b6..892831f23 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp @@ -434,7 +434,7 @@ void ClassCRHMCanopyVectorBased::run(void){ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) { // calculate increase in Snow_load and direct_snow if we are in canopy (i.e., Cc > 0) - const double k_cp = -20; // rate of increase of the sigmoidal curve below + const double k_cp = 20; // rate of increase of the sigmoidal curve below const double v_snow = 0.8; // terminal fall velocity of snowfall taken from Isyumov, 1971 Clca[hh] = 0; // leaf contact area (Clca) based on trajectory angle double IP = 0; // interception efficiency (IP) @@ -444,7 +444,7 @@ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) { // increase leaf contact area (Clca) based on wind speed and canopy coverage (Cc) double Ht_1_third = Ht[hh] * (1.0 / 3.0); double Cp_inc = 0; - if (Ht_1_third - (2.0 / 3.0) * Zwind[hh] > 0.0) + if ((Ht_1_third - (2.0 / 3.0) * Zwind[hh]) > 0.0) { u_1_third_Ht[hh] = hru_u[hh] * log((Ht_1_third - (2.0 / 3.0) * Zwind[hh]) / 0.123 * Zwind[hh]) / log((Zwind[hh] - 2.0 / 3.0 * Zwind[hh]) / 0.123 * Zwind[hh]); double snow_traj_angle = atan(u_1_third_Ht[hh] / v_snow); // in radians From fb29dc5c68b37920c3a766532ae9618abf7ec328 Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Mon, 18 Nov 2024 16:11:33 -0800 Subject: [PATCH 10/15] finished first try at mass unloading switch and started working on meltwater unloading switch --- crhmcode/src/core/CRHM_constants.h | 2 + .../modules/ClassCRHMCanopyVectorBased.cpp | 143 +++++++++++++----- .../src/modules/ClassCRHMCanopyVectorBased.h | 3 + 3 files changed, 114 insertions(+), 34 deletions(-) diff --git a/crhmcode/src/core/CRHM_constants.h b/crhmcode/src/core/CRHM_constants.h index a1c1f15a2..d366c31fb 100644 --- a/crhmcode/src/core/CRHM_constants.h +++ b/crhmcode/src/core/CRHM_constants.h @@ -39,6 +39,8 @@ namespace CRHM_constants { const double emiss = 0.985; // emissivity of the atmosphere and snowpack const double emiss_c = 0.96; // emissivity of the canopy const double em = 0.622; // + const double ci = 2.102e-3; // heat capacity of ice (MJ/kg/K) + } #endif // !CRHM_CONSTANTS diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp index 892831f23..d2039febe 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp @@ -138,7 +138,7 @@ void ClassCRHMCanopyVectorBased::decl(void) { declstatdiag("cum_net_rain", TDim::NHRU, "cumulative direct_rain + drip", "(mm)", &cum_net_rain); - declvar("pot_subl_cpy", TDim::NHRU, "dimensionless canopy snow sublimation rate aka potential sublimation rate to be multiplied by canopy snow load", "(s-1)", &pot_subl_cpy); + declvar("pot_subl_cpy", TDim::NHRU, "sublimation rate coefficient for single ice spheres aka potential sublimation rate to be multiplied by canopy snow load and the exposure coeficient", "(s-1)", &pot_subl_cpy); declvar("Subl_Cpy", TDim::NHRU, "canopy snow sublimation", "(mm/int)", &Subl_Cpy); @@ -168,7 +168,7 @@ void ClassCRHMCanopyVectorBased::decl(void) { 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) following Cebulski & Pomeroy vector based param.", "(m/s)", &u_1_third_Ht); + 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); decldiag("Cc", TDim::NHRU, "Canopy coverage", "()", &Cc); @@ -216,6 +216,12 @@ void ClassCRHMCanopyVectorBased::decl(void) { declparam("CanopyClearing", TDim::NHRU, "[0]", "0", "2", "canopy - 0/clearing - 1/gap - 2", "()", &CanopyClearing); + declparam("SublimationSwitch", TDim::NHRU, "[1]", "0", "1", "Pomeroy 1998 sublimation parameterisation, off - 0, on - 1", "()", &SublimationSwitch); + + declparam("MassUnloadingSwitch", TDim::NHRU, "[1]", "0", "1", "Canopy snow mass unloading parameterisation options: Hedstrom & Pomeroy 1998 constant time based unloading - 0, Cebulski & Pomeroy exponential curves for wind induced unloading, temperature based unloading, and time based unloading - 1", "()", &MassUnloadingSwitch); + + declparam("MeltwaterSwitch", TDim::NHRU, "[0]", "0", "1", "Canopy snow meltwater drip parameterisation options: Ellis et al. (2010) and Floyd (2012) - 0, CLASS - 1", "()", &MeltwaterSwitch); + 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); @@ -432,7 +438,7 @@ void ClassCRHMCanopyVectorBased::run(void){ // coupled forest snow interception and sublimation routine: // after Cebulski & Pomeroy 2025: -if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) +if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0)) { // calculate increase in Snow_load and direct_snow if we are in canopy (i.e., Cc > 0) const double k_cp = 20; // rate of increase of the sigmoidal curve below const double v_snow = 0.8; // terminal fall velocity of snowfall taken from Isyumov, 1971 @@ -440,7 +446,7 @@ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) double IP = 0; // interception efficiency (IP) double dL = 0; // change in canopy snow load - if (hru_u[hh] > 0 && Cc[hh] < 1) + if (hru_u[hh] > 0 && Cc[hh] < 1 && Cc[hh] > 0) { // increase leaf contact area (Clca) based on wind speed and canopy coverage (Cc) double Ht_1_third = Ht[hh] * (1.0 / 3.0); double Cp_inc = 0; @@ -448,13 +454,13 @@ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) { u_1_third_Ht[hh] = hru_u[hh] * log((Ht_1_third - (2.0 / 3.0) * Zwind[hh]) / 0.123 * Zwind[hh]) / log((Zwind[hh] - 2.0 / 3.0 * Zwind[hh]) / 0.123 * Zwind[hh]); double snow_traj_angle = atan(u_1_third_Ht[hh] / v_snow); // in radians - Cp_inc = (1 - Cc[hh]) / (1 + exp(-k_cp * (sin(snow_traj_angle) - (1 - Cc[hh])))); // fractional increas in leaf contact area (Clca) based on horizontal trajectory. This is modified from Cebulski & Pomeroy snow interception paper. + Cp_inc = (1 - Cc[hh]) / (1 + exp(-k_cp * (sin(snow_traj_angle) - (1 - Cc[hh])))); // fractional increas in leaf contact area (Clca) based on horizontal trajectory. This is modified from Cebulski & Pomeroy snow interception paper. Has only been tested on forest plots with Cc of .3 and .5. } Clca[hh] = Cc[hh] + Cp_inc; // calculated leaf contact area (Clca) based on trajectory angle } else { - Clca[hh] = Cc[hh]; // use leaf contact area from nadir i.e., canopy coverage + 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) @@ -466,6 +472,20 @@ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) direct_snow[hh] += (1-IP) * hru_snow[hh]; + // Finished initial loading now start the ablation paramaterisations + + double Vi = 0; // submilation rate coefficient for exposed snow (s-1) + double IceBulbT = 0; // ice bulb temperature of canopy snow + double U = 0; // unloading rate coefficient used in Hedstrom & Pomeroy 1998 param. + +switch(SublimationSwitch[hh]){ +case 0: // do not sublimate, used for debugging or experiments, recommend using case 1 otherwise. + + std::cout << "SublimationSwitch Case 0: No canopy snow sublimation applied.\n"; + break; + +case 1: // canopy snow sublimation with Pomeroy et al. (1998) paramaterisation + // Pomeroy et al. (1998) sublimation routine modified by Alex Cebulski to change maximum intercepted load from Hedstrom 1998 calculation to constant based on observations from several studies showing much higher intercepted loads. // calculate snow ventilation windspeed: @@ -476,8 +496,7 @@ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) const double KinVisc = 1.88e-5; // kinematic viscosity of air (Sask. avg. value) const double ks = 0.0114; // snow shape coefficient for jack pine const double Fract = 0.37; // fractal dimension of intercepted snow - const double ci = 2.102e-3; // heat capacity of ice (MJ/kg/K) - const double Hs = 2.838e6; // heat of sublimation (MJ/kg) + // const double Hs = 2.838e6; // heat of sublimation (MJ/kg) changed to PBSM_constants::LATH which has same value. //============================================================================== double xi2 = 1-Zvent[hh]; @@ -506,7 +525,7 @@ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) Nu = 1.79 + 0.606 * sqrt((double) Nr); // Nusselt number SStar = M_PI * sqr(Radius) * (1.0f - AlbedoIce) * Qsi_; // SW to snow particle !!!! changed A1 = Lamb * (hru_t[hh] + 273) * Nu; - B1 = Hs * PBSM_constants::M /(PBSM_constants::R * (hru_t[hh] + 273.0f))- 1.0; + B1 = PBSM_constants::LATH * PBSM_constants::M /(PBSM_constants::R * (hru_t[hh] + 273.0f))- 1.0; J = B1/A1; Sigma2 = hru_rh[hh]/100 -1; D = 2.06e-5* pow((hru_t[hh]+273.0f)/273.0f, -1.75f); // diffusivity of water vapour @@ -516,7 +535,7 @@ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) Mpm = 4.0 / 3.0 * M_PI * PBSM_constants::DICE * Radius * Radius * Radius; // 18Mar2022: remove Gamma Distribution Correction term, *(1.0 + 3.0/Alpha + 2.0/sqr(Alpha)); // sublimation rate of single 'ideal' ice sphere: - double Vs = (2.0* M_PI* Radius*Sigma2 - SStar* J)/(Hs* J + C1)/Mpm; + double Vs = (2.0* M_PI* Radius*Sigma2 - SStar* J)/(PBSM_constants::LATH* J + C1)/Mpm; pot_subl_cpy[hh] = Vs; // export the dimensionless sublimation rate (s-1) added by alex 2023-07-21 @@ -530,11 +549,11 @@ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) // calculate 'potential' canopy sublimation: - double Vi = Vs*Ce; + Vi = Vs*Ce; // limit sublimation to canopy snow available and take sublimated snow away from canopy snow at timestep start - Subl_Cpy[hh] = -Snow_load[hh]*Vi*Hs*Global::Interval*24*3600/Hs; // make W/m2 + Subl_Cpy[hh] = -Snow_load[hh]*Vi*PBSM_constants::LATH*Global::Interval*24*3600/PBSM_constants::LATH; // make W/m2 if(Subl_Cpy[hh] > Snow_load[hh]){ Subl_Cpy[hh] = Snow_load[hh]; Snow_load[hh] = 0.0; @@ -544,38 +563,94 @@ if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0) && Cc[hh] > 0) if(Snow_load[hh] < 0.0) Snow_load[hh] = 0.0; } +} -// calculate 'ice-bulb' temperature of intercepted snow: +if(MassUnloadingSwitch[hh] == 0 || MeltwaterSwitch[hh] == 0){ + // calculate 'ice-bulb' temperature of intercepted snow: + double IceBulbT = hru_t[hh] - (Vi* PBSM_constants::LATH/1e6/CRHM_constants::ci); + const double U = -1 * log(0.678) / (24 * 7 * Global::Freq / 24); // weekly dimensionless unloading coefficient -> to CRHM time interval // 21Mar2022 correction: invert the term 24/Global::Freq, use unloading rate coefficient U = -log(c)/t for snow unloading determined by inverse function of c = e^(-Ut) = 0.678 based on Eq. 14 in Hedstrom and Pomeroy (1998) +} - double IceBulbT = hru_t[hh] - (Vi* Hs/1e6/ci); - double Six_Hour_Divisor = Global::Freq/4.0; // used to unload over 6 hours +switch(MassUnloadingSwitch[hh]){ + case 0: // This is the mass unloading portion of the latest iteration of the Hedstrom & Pomeroy 1998 unloading with modifications by Ellis et al. (2010) and Floyd (2012). Generally used with MeltwaterSwitch == 0. + + // determine whether canopy snow is unloaded as mass clumps: - const double U = -1 * log(0.678) / (24 * 7 * Global::Freq / 24); // weekly dimensionless unloading coefficient -> to CRHM time interval - // 21Mar2022 correction: invert the term 24/Global::Freq, use unloading rate coefficient U = -log(c)/t for snow unloading determined by inverse function of c = e^(-Ut) = 0.678 based on Eq. 14 in Hedstrom and Pomeroy (1998) - // determine whether canopy snow is unloaded: + if(IceBulbT < unload_t[hh]){ // has to be at least one interval. Trip on half step + SUnload[hh] = Snow_load[hh] * U; // the dimensionless unloading coefficient already /interval, 21Mar2022 correction: use unloading rate coefficient U + } + + case 1: // This is the updated mass snow unloading parameterisations from Cebulski & Pomeroy to unload based on time, wind, and air temperature. + // temperature induced unloading + const double a_T = 2.584003e-05; // Cebulski & Pomeroy coef from exponential function of unloading + drip and air temp measurements at Fortress mountain when wind speed <= 1 m/s. + const double b_T = 1.646875e-01; // Cebulski & Pomeroy coef from exponential function of unloading + drip and air temp measurements at Fortress mountain when wind speed <= 1 m/s. + + double fT = a_T * exp(b_T * hru_t[hh]); // unloading rate based on warming of snow in the canopy (s-1), still need to partition out the portion of this that is drip vs mass unloading + + // mechanical wind induced unloading + const double a_u = 5.204024e-06; // Cebulski & Pomeroy coef from exponential function of unloading + drip and wind speed measurements at Fortress mountain when air temp < -6 C. + const double b_u = 7.363594e-02; // Cebulski & Pomeroy coef from exponential function of unloading + drip and wind speed measurements at Fortress mountain when air temp < -6 C. + double Ht_mid = Ht[hh] * (1.0 / 2.0); // half of canopy height + double u_mid = 0; + double fu = 0; + if ((Ht_mid - (2.0 / 3.0) * Zwind[hh]) > 0.0) + { + u_mid = hru_u[hh] * log((Ht_mid - (2.0 / 3.0) * Zwind[hh]) / 0.123 * Zwind[hh]) / log((Zwind[hh] - 2.0 / 3.0 * Zwind[hh]) / 0.123 * Zwind[hh]); + double fu = u_mid * a_u * exp(b_u * u_mid); // unloading rate due to wind (s-1) + } + + // duration based unloading + const double a_t = 8.194345e-06; // Cebulski & Pomeroy coef from exponential function of unloading + drip and duration snow has been intercepted in the canopy at Fortress mountain when wind speed <= 1 m/s and air temperature < -6 C. + const double b_t = -1.540050e+02; // Cebulski & Pomeroy coef from exponential function of unloading + drip and duration snow has been intercepted in the canopy at Fortress mountain when wind speed <= 1 m/s and air temperature < -6 C. - if(IceBulbT >= unload_t_water[hh]){ - drip_Cpy[hh] = Snow_load[hh]/Six_Hour_Divisor; + double t_snow_in_canopy = 12 * 60 * 60; // duration snow intercepted in the canopy, set to constant of 12 hours. TODO need to track duration that snow has been intercepted in the canopy. Can do this similar to snowpack albedo calculation. + + double ft = a_t * exp(b_t * t_snow_in_canopy); + + // ablation via temperature, wind, and duration based unloading + double dt = Global::Interval * 24 * 60 * 60; // converts the interval which is a time period (i.e., time/cycles, 1 day/# obs) to timestep in seconds. + SUnload[hh] = Snow_load[hh] * (fT + fu + ft) * dt; // converts our +} + +// handle mass unloading regardless of what parameterisation is chosen +if (SUnload[hh] > Snow_load[hh]) +{ + SUnload[hh] = Snow_load[hh]; + Snow_load[hh] = 0.0; +} +else + Snow_load[hh] -= SUnload[hh]; + +cum_SUnload[hh] += SUnload[hh]; + +switch(MeltwaterSwitch[hh]) { + case 0: { // Block for case 0 + // This is the meltwater drip portion of the latest iteration of the Hedstrom & Pomeroy 1998 unloading with modifications by Ellis et al. (2010) and Floyd (2012). + double Six_Hour_Divisor = Global::Freq / 4.0; // Unload over 6 hours + + if (IceBulbT >= unload_t_water[hh]) { + drip_Cpy[hh] = Snow_load[hh] / Six_Hour_Divisor; SUnload_H2O[hh] = drip_Cpy[hh]; Snow_load[hh] -= SUnload_H2O[hh]; cum_SUnload_H2O[hh] += SUnload_H2O[hh]; - } - else if(IceBulbT >= unload_t[hh]){ - SUnload[hh] = Snow_load[hh]/Six_Hour_Divisor; + } + else if (IceBulbT >= unload_t[hh]) { + SUnload[hh] = Snow_load[hh] / Six_Hour_Divisor; Snow_load[hh] -= SUnload[hh]; cum_SUnload[hh] += SUnload[hh]; - } - else if(IceBulbT < unload_t[hh]){ // has to be at least one interval. Trip on half step - SUnload[hh] = Snow_load[hh] * U; // the dimensionless unloading coefficient already /interval, 21Mar2022 correction: use unloading rate coefficient U - if(SUnload[hh] > Snow_load[hh]){ - SUnload[hh] = Snow_load[hh]; - Snow_load[hh] = 0.0; - } - else - Snow_load[hh] -= SUnload[hh]; + } + break; + } - cum_SUnload[hh] += SUnload[hh]; - } + case 1: { // Block for case 1 + // This is the meltwater drip parameterisation from CLASS + drip_Cpy[hh] = 0; + SUnload_H2O[hh] = drip_Cpy[hh]; + Snow_load[hh] -= SUnload_H2O[hh]; + cum_SUnload_H2O[hh] += SUnload_H2O[hh]; + break; + } +} // calculate total sub-canopy snow: diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h index 0599bfdc4..9a2f8578d 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h @@ -125,6 +125,9 @@ const double *Alpha_c { NULL }; const double *B_canopy { NULL }; const long *CanopyClearing { NULL }; +const long *SublimationSwitch { NULL }; +const long *MassUnloadingSwitch { NULL }; +const long *MeltwaterSwitch { NULL }; const long *inhibit_evap { NULL }; From e4fbe81a95224a3bdc8f78df6320784915f9b9a5 Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Mon, 18 Nov 2024 18:07:02 -0800 Subject: [PATCH 11/15] !!! changed formatting to vscode recommended --- .../modules/ClassCRHMCanopyVectorBased.cpp | 775 ++++++++++-------- 1 file changed, 411 insertions(+), 364 deletions(-) diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp index d2039febe..7d7ef3f25 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp @@ -1,22 +1,22 @@ /** -* 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 + * 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 @@ -30,14 +30,15 @@ #include "../core/ClassCRHM.h" #include "newmodules/SnobalDefines.h" - using namespace CRHM; -ClassCRHMCanopyVectorBased* ClassCRHMCanopyVectorBased::klone(string name) const{ +ClassCRHMCanopyVectorBased *ClassCRHMCanopyVectorBased::klone(string name) const +{ return new ClassCRHMCanopyVectorBased(name); } -void ClassCRHMCanopyVectorBased::decl(void) { +void ClassCRHMCanopyVectorBased::decl(void) +{ Description = "'Prototype all season canopy/clearing module. Calculates short, long and all-wave radiation components at the snow surface.' \ 'Inputs are observations Qsi (W/m^2) and Qli (W/m^2), ' \ @@ -46,43 +47,39 @@ void ClassCRHMCanopyVectorBased::decl(void) { '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 + // 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); - + declgetvar("*", "QliVt_Var", "(W/m^2)", &QliVt_Var); variation_set = VARIATION_2 + VARIATION_3; - declgetvar("*", "QsiS_Var", "(W/m^2)", &QsiS_Var); + declgetvar("*", "QsiS_Var", "(W/m^2)", &QsiS_Var); variation_set = VARIATION_4; - declgetvar("*", "QsiA_Var", "(W/m^2)", &QsiA_Var); - + declgetvar("*", "QsiA_Var", "(W/m^2)", &QsiA_Var); variation_set = VARIATION_ORG; -// get variables: + // get variables: - declgetvar("*", "hru_t", "(" + string(DEGREE_CELSIUS) + ")", &hru_t); + declgetvar("*", "hru_t", "(" + string(DEGREE_CELSIUS) + ")", &hru_t); - declgetvar("*", "hru_u", "(m/s)", &hru_u); + declgetvar("*", "hru_u", "(m/s)", &hru_u); - declgetvar("*", "hru_rh", "()", &hru_rh); + declgetvar("*", "hru_rh", "()", &hru_rh); - declgetvar("*", "hru_ea", "(kPa)", &hru_ea); + declgetvar("*", "hru_ea", "(kPa)", &hru_ea); declgetvar("*", "hru_snow", "(mm/int)", &hru_snow); @@ -102,7 +99,7 @@ void ClassCRHMCanopyVectorBased::decl(void) { declputvar("*", "SWE", "(mm)", &SWE); -// declared observations + // declared observations declobs("Ts", TDim::NHRU, "snow surface temperature", "(" + string(DEGREE_CELSIUS) + ")", &Ts); @@ -120,7 +117,7 @@ void ClassCRHMCanopyVectorBased::decl(void) { declobs("Qlosn", TDim::NHRU, "reflected long-wave at surface", "(W/m^2)", &Qlosn); -// declared variables + // declared variables decldiag("k", TDim::NHRU, "extinction coefficient", "()", &k); @@ -178,8 +175,7 @@ void ClassCRHMCanopyVectorBased::decl(void) { declstatdiag("cum_intcp_evap", TDim::NHRU, "Cumulative HRU Evaporation from interception", "(mm)", &cum_intcp_evap); - -// parameters: + // parameters: declparam("inhibit_evap", TDim::NHRU, "[0]", "0", "1", "inhibit evaporation, 1 -> inhibit", "()", &inhibit_evap); @@ -193,7 +189,6 @@ void ClassCRHMCanopyVectorBased::decl(void) { 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); @@ -217,7 +212,7 @@ void ClassCRHMCanopyVectorBased::decl(void) { declparam("CanopyClearing", TDim::NHRU, "[0]", "0", "2", "canopy - 0/clearing - 1/gap - 2", "()", &CanopyClearing); declparam("SublimationSwitch", TDim::NHRU, "[1]", "0", "1", "Pomeroy 1998 sublimation parameterisation, off - 0, on - 1", "()", &SublimationSwitch); - + declparam("MassUnloadingSwitch", TDim::NHRU, "[1]", "0", "1", "Canopy snow mass unloading parameterisation options: Hedstrom & Pomeroy 1998 constant time based unloading - 0, Cebulski & Pomeroy exponential curves for wind induced unloading, temperature based unloading, and time based unloading - 1", "()", &MassUnloadingSwitch); declparam("MeltwaterSwitch", TDim::NHRU, "[0]", "0", "1", "Canopy snow meltwater drip parameterisation options: Ellis et al. (2010) and Floyd (2012) - 0, CLASS - 1", "()", &MeltwaterSwitch); @@ -227,13 +222,15 @@ void ClassCRHMCanopyVectorBased::decl(void) { 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) { +void ClassCRHMCanopyVectorBased::init(void) +{ nhru = getdim(TDim::NHRU); // transfers current # of HRU's to module - for (hh = 0; hh < nhru; ++hh) { + for (hh = 0; hh < nhru; ++hh) + { - Pa[hh] = 101.3f*pow((293.0f-0.0065f*hru_elev[hh])/293.0f, 5.26f); // kPa + Pa[hh] = 101.3f * pow((293.0f - 0.0065f * hru_elev[hh]) / 293.0f, 5.26f); // kPa rain_load[hh] = 0.0; Snow_load[hh] = 0.0; @@ -245,40 +242,44 @@ void ClassCRHMCanopyVectorBased::init(void) { cum_SUnload[hh] = 0.0; cum_SUnload_H2O[hh] = 0.0; - if(Ht[hh] > Zwind[hh]){ + 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){ +void ClassCRHMCanopyVectorBased::run(void) +{ double Exposure, LAI_, Vf, Vf_, Kstar_H, Kd; - //double Tau; variable is unreferenced commenting out for now - jhs507 + // double Tau; variable is unreferenced commenting out for now - jhs507 - for (hh = 0; chkStruct(); ++hh){ + for (hh = 0; chkStruct(); ++hh) + { - switch (variation){ - case VARIATION_ORG: - Qsi_ = Qsi[hh]; - Qli_ = Qli[hh]; + switch (variation) + { + case VARIATION_ORG: + Qsi_ = Qsi[hh]; + Qli_ = Qli[hh]; break; - case VARIATION_1: - Qsi_ = Qsi[hh]; - Qli_ = QliVt_Var[hh]; + case VARIATION_1: + Qsi_ = Qsi[hh]; + Qli_ = QliVt_Var[hh]; break; - case VARIATION_2: - Qsi_ = QsiS_Var[hh]; - Qli_ = Qli[hh]; + case VARIATION_2: + Qsi_ = QsiS_Var[hh]; + Qli_ = Qli[hh]; break; - case VARIATION_3: - Qsi_ = QsiS_Var[hh]; - Qli_ = QliVt_Var[hh]; + case VARIATION_3: + Qsi_ = QsiS_Var[hh]; + Qli_ = QliVt_Var[hh]; break; - case VARIATION_4: - Qsi_ = QsiA_Var[hh]; - Qli_ = QliVt_Var[hh]; + case VARIATION_4: + Qsi_ = QsiA_Var[hh]; + Qli_ = QliVt_Var[hh]; break; } // switch @@ -293,453 +294,497 @@ void ClassCRHMCanopyVectorBased::run(void){ SUnload_H2O[hh] = 0.0; Subl_Cpy[hh] = 0.0; -// Canopy temperature is approximated by the air temperature. + // Canopy temperature is approximated by the air temperature. double T1 = hru_t[hh] + CRHM_constants::Tm; - double rho = Pa[hh]*1000/(CRHM_constants::Rgas*T1); + double rho = Pa[hh] * 1000 / (CRHM_constants::Rgas * T1); 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; + 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 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 (kg/kg) + double q = (hru_rh[hh] / 100) * Common::Qs(Pa[hh], T1); // specific humidity (kg/kg) - Ts[hh] = T1 + (CRHM_constants::emiss*(Qli_ - CRHM_constants::sbc*pow(T1, 4.0f)) + CRHM_constants::Ls*(q - Common::Qs(Pa[hh], T1))*rho/ra[hh])/ - (4*CRHM_constants::emiss*CRHM_constants::sbc*pow(T1, 3.0f) + (CRHM_constants::Cp + CRHM_constants::Ls*deltaX)*rho/ra[hh]); + Ts[hh] = T1 + (CRHM_constants::emiss * (Qli_ - CRHM_constants::sbc * pow(T1, 4.0f)) + CRHM_constants::Ls * (q - Common::Qs(Pa[hh], T1)) * rho / ra[hh]) / + (4 * CRHM_constants::emiss * CRHM_constants::sbc * pow(T1, 3.0f) + (CRHM_constants::Cp + CRHM_constants::Ls * deltaX) * rho / ra[hh]); Ts[hh] -= CRHM_constants::Tm; - if(Ts[hh] > 0.0 || SWE[hh] <= 0.0) + if (Ts[hh] > 0.0 || SWE[hh] <= 0.0) Ts[hh] = 0.0; - switch(CanopyClearing[hh]){ - case 0: // canopy + switch (CanopyClearing[hh]) + { + case 0: // canopy - Exposure = Ht[hh] - Common::DepthofSnow(SWE[hh]); /* depths(m) SWE(mm) */ - if(Exposure < 0.0) - Exposure = 0.0; + Exposure = Ht[hh] - Common::DepthofSnow(SWE[hh]); /* depths(m) SWE(mm) */ + if (Exposure < 0.0) + Exposure = 0.0; - LAI_ = LAI[hh]*Exposure/Ht[hh]; + LAI_ = LAI[hh] * Exposure / Ht[hh]; - Vf = 0.45 - 0.29*log(LAI[hh]); + Vf = 0.45 - 0.29 * log(LAI[hh]); - Vf_ = Vf + (1.0 - Vf)*sin((Ht[hh] - Exposure)/Ht[hh]*M_PI_2); + Vf_ = Vf + (1.0 - Vf) * sin((Ht[hh] - Exposure) / 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]); - double limit = cosxsflat[hh]/cosxs[hh]; - if(limit > 2.0) - limit = 2.0; - Tauc[hh] = exp(-k[hh]*LAI_*limit); - } - else{ - k[hh] = 0.0; - Tauc[hh] = 0.0; - } + 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]); + double limit = cosxsflat[hh] / cosxs[hh]; + if (limit > 2.0) + limit = 2.0; + Tauc[hh] = exp(-k[hh] * LAI_ * limit); + } + else + { + k[hh] = 0.0; + Tauc[hh] = 0.0; + } - Kstar_H = Qsi_*(1.0 - Alpha_c[hh] - Tauc[hh]*(1.0 - Albedo[hh])); + Kstar_H = Qsi_ * (1.0 - Alpha_c[hh] - Tauc[hh] * (1.0 - Albedo[hh])); - Qlisn[hh] = Qli_*Vf_ + (1.0f - Vf_)*CRHM_constants::emiss_c*CRHM_constants::sbc*pow(T1, 4.0f) + B_canopy[hh]*Kstar_H; + Qlisn[hh] = Qli_ * Vf_ + (1.0f - Vf_) * CRHM_constants::emiss_c * CRHM_constants::sbc * pow(T1, 4.0f) + B_canopy[hh] * Kstar_H; - Qlisn_Var[hh] = Qlisn[hh]; + Qlisn_Var[hh] = Qlisn[hh]; - Qsisn[hh] = Qsi_*Tauc[hh]; + Qsisn[hh] = Qsi_ * Tauc[hh]; - Qsisn_Var[hh] = Qsisn[hh]; + Qsisn_Var[hh] = Qsisn[hh]; - Qlosn[hh] = CRHM_constants::emiss*CRHM_constants::sbc*pow(Ts[hh] + CRHM_constants::Tm, 4.0f); + 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[hh] = Qlisn[hh] - Qlosn[hh] + Qsisn[hh] * (1.0 - Albedo[hh]); - Qnsn_Var[hh] = Qnsn[hh]; + Qnsn_Var[hh] = Qnsn[hh]; - break; - case 1: // clearing + break; + case 1: // clearing - Qlisn[hh] = Qli_; + Qlisn[hh] = Qli_; - Qlisn_Var[hh] = Qlisn[hh]; + Qlisn_Var[hh] = Qlisn[hh]; - Qsisn[hh] = Qsi_; + Qsisn[hh] = Qsi_; - Qsisn_Var[hh] = Qsisn[hh]; + Qsisn_Var[hh] = Qsisn[hh]; - Qlosn[hh] = CRHM_constants::emiss*CRHM_constants::sbc*pow(Ts[hh] + CRHM_constants::Tm, 4.0f); + 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[hh] = Qlisn[hh] - Qlosn[hh] + Qsisn[hh] * (1.0 - Albedo[hh]); - Qnsn_Var[hh] = Qnsn[hh]; + Qnsn_Var[hh] = Qnsn[hh]; - break; - case 2: // gap - Exposure = Surrounding_Ht[hh] - Common::DepthofSnow(SWE[hh]); /* depths(m) SWE(mm) */ - if(Exposure < 0.0) - Exposure = 0.0; + break; + case 2: // gap + Exposure = Surrounding_Ht[hh] - Common::DepthofSnow(SWE[hh]); /* depths(m) SWE(mm) */ + if (Exposure < 0.0) + Exposure = 0.0; - LAI_ = LAI[hh]*Exposure/Surrounding_Ht[hh]; + LAI_ = LAI[hh] * Exposure / Surrounding_Ht[hh]; - Vf = 0.45 - 0.29*log(LAI[hh]); + Vf = 0.45 - 0.29 * log(LAI[hh]); - double Tau_d = Vf + (1.0 - Vf)*sin((Surrounding_Ht[hh] - Exposure)/Surrounding_Ht[hh]*M_PI_2); // previously Vf_ + double Tau_d = Vf + (1.0 - Vf) * sin((Surrounding_Ht[hh] - Exposure) / Surrounding_Ht[hh] * M_PI_2); // previously Vf_ -// calculate forest clearing sky view factor (Vgap) via Reifsnyder and Lull?s (1965) expression: + // calculate forest clearing sky view factor (Vgap) via Reifsnyder and Lull?s (1965) expression: - double Vgap = sqr(sin(atan2(Gap_diameter[hh], 2.0*Surrounding_Ht[hh]))); + double Vgap = sqr(sin(atan2(Gap_diameter[hh], 2.0 * Surrounding_Ht[hh]))); -// calculate beam pathlength correction (variable Gap_beam_corr) for gap: + // calculate beam pathlength correction (variable Gap_beam_corr) for gap: - double Gap_beam_corr = 0; - if(Qsi_ > 0.0 && SolAng[hh] > 0.001){ - double cosxsLim = 3; - if(cosxs[hh] > 0.33) - cosxsLim = 1.0/cosxs[hh]; + double Gap_beam_corr = 0; + if (Qsi_ > 0.0 && SolAng[hh] > 0.001) + { + double cosxsLim = 3; + if (cosxs[hh] > 0.33) + cosxsLim = 1.0 / cosxs[hh]; - Gap_beam_corr = cosxsLim*Surrounding_Ht[hh]*(1.0/cos(SolAng[hh]) - Gap_diameter[hh]/(2.0*Surrounding_Ht[hh])/sin(SolAng[hh])); - if(Gap_beam_corr > 10.0) - Gap_beam_corr = 10.0; - else if(Gap_beam_corr < 0.0) - Gap_beam_corr = 0.0; - } -// calculate beam shortwave transmittance of the gap: + Gap_beam_corr = cosxsLim * Surrounding_Ht[hh] * (1.0 / cos(SolAng[hh]) - Gap_diameter[hh] / (2.0 * Surrounding_Ht[hh]) / sin(SolAng[hh])); + if (Gap_beam_corr > 10.0) + Gap_beam_corr = 10.0; + else if (Gap_beam_corr < 0.0) + Gap_beam_corr = 0.0; + } + // calculate beam shortwave transmittance of the gap: - double product = LAI[hh]*Gap_beam_corr; - if(product > 50) - product = 50; + double product = LAI[hh] * Gap_beam_corr; + if (product > 50) + product = 50; - double Tau_b_gap = exp(-product); + double Tau_b_gap = exp(-product); - Kd = Qsi_*(1.0 - Alpha_c[hh] - Tau_b_gap*(1.0 - Albedo[hh])); + Kd = Qsi_ * (1.0 - Alpha_c[hh] - Tau_b_gap * (1.0 - Albedo[hh])); - Qlisn[hh] = Vgap*Qli_ + (1.0 - Vgap)*((Qli_*Tau_b_gap + (1.0 - Tau_b_gap)*CRHM_constants::emiss_c*CRHM_constants::sbc*pow(T1, 4.0f)) + B_canopy[hh]*Kd); + Qlisn[hh] = Vgap * Qli_ + (1.0 - Vgap) * ((Qli_ * Tau_b_gap + (1.0 - Tau_b_gap) * CRHM_constants::emiss_c * CRHM_constants::sbc * pow(T1, 4.0f)) + B_canopy[hh] * Kd); - Qlisn_Var[hh] = Qlisn[hh]; + Qlisn_Var[hh] = Qlisn[hh]; - Qsisn[hh] = cosxs[hh]*Qdfo[hh]*Tau_b_gap + Vgap*(Qsi_ - Qdfo[hh]) + (1.0 - Vgap)*Tau_d*(Qsi_ - Qdfo[hh]); - if(Qsisn[hh] < 0.0) - Qsisn[hh] = 0.0; + Qsisn[hh] = cosxs[hh] * Qdfo[hh] * Tau_b_gap + Vgap * (Qsi_ - Qdfo[hh]) + (1.0 - Vgap) * Tau_d * (Qsi_ - Qdfo[hh]); + if (Qsisn[hh] < 0.0) + Qsisn[hh] = 0.0; - Qsisn_Var[hh] = Qsisn[hh]; + Qsisn_Var[hh] = Qsisn[hh]; - Qlosn[hh] = CRHM_constants::emiss*CRHM_constants::sbc*pow(Ts[hh] + CRHM_constants::Tm, 4.0f); + 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-Albedo[hh]); + Qnsn[hh] = Qlisn[hh] - Qlosn[hh] + Qsisn[hh] * (1 - Albedo[hh]); - Qnsn_Var[hh] = Qnsn[hh]; + Qnsn_Var[hh] = Qnsn[hh]; - break; + break; } // switch - switch(CanopyClearing[hh]){ - case 0: // canopy -//============================================================================== -// coupled forest snow interception and sublimation routine: -// after Cebulski & Pomeroy 2025: - -if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0)) -{ // calculate increase in Snow_load and direct_snow if we are in canopy (i.e., Cc > 0) - const double k_cp = 20; // rate of increase of the sigmoidal curve below - const double v_snow = 0.8; // terminal fall velocity of snowfall taken from Isyumov, 1971 - Clca[hh] = 0; // leaf contact area (Clca) based on trajectory angle - double IP = 0; // interception efficiency (IP) - double dL = 0; // change in canopy snow load - - if (hru_u[hh] > 0 && Cc[hh] < 1 && Cc[hh] > 0) - { // increase leaf contact area (Clca) based on wind speed and canopy coverage (Cc) - double Ht_1_third = Ht[hh] * (1.0 / 3.0); - double Cp_inc = 0; - if ((Ht_1_third - (2.0 / 3.0) * Zwind[hh]) > 0.0) + switch (CanopyClearing[hh]) { - u_1_third_Ht[hh] = hru_u[hh] * log((Ht_1_third - (2.0 / 3.0) * Zwind[hh]) / 0.123 * Zwind[hh]) / log((Zwind[hh] - 2.0 / 3.0 * Zwind[hh]) / 0.123 * Zwind[hh]); - double snow_traj_angle = atan(u_1_third_Ht[hh] / v_snow); // in radians - Cp_inc = (1 - Cc[hh]) / (1 + exp(-k_cp * (sin(snow_traj_angle) - (1 - Cc[hh])))); // fractional increas in leaf contact area (Clca) based on horizontal trajectory. This is modified from Cebulski & Pomeroy snow interception paper. Has only been tested on forest plots with Cc of .3 and .5. - } - Clca[hh] = Cc[hh] + Cp_inc; // calculated leaf contact area (Clca) based on trajectory angle - } - else - { - Clca[hh] = Cc[hh]; // use leaf contact area from nadir i.e., Clca == 1 for Cc == 1 and Clca == 0 when Cc == 0 - } + case 0: // canopy + //============================================================================== + // coupled forest snow interception and sublimation routine: + // after Cebulski & Pomeroy 2025: + + if ((Snow_load[hh] > 0.0 || hru_snow[hh] > 0.0)) + { // calculate increase in Snow_load and direct_snow if we are in canopy (i.e., Cc > 0) + const double k_cp = 20; // rate of increase of the sigmoidal curve below + const double v_snow = 0.8; // terminal fall velocity of snowfall taken from Isyumov, 1971 + Clca[hh] = 0; // leaf contact area (Clca) based on trajectory angle + double IP = 0; // interception efficiency (IP) + double dL = 0; // change in canopy snow load + + if (hru_u[hh] > 0 && Cc[hh] < 1 && Cc[hh] > 0) + { // increase leaf contact area (Clca) based on wind speed and canopy coverage (Cc) + double Ht_1_third = Ht[hh] * (1.0 / 3.0); + double Cp_inc = 0; + if ((Ht_1_third - (2.0 / 3.0) * Zwind[hh]) > 0.0) + { + u_1_third_Ht[hh] = hru_u[hh] * log((Ht_1_third - (2.0 / 3.0) * Zwind[hh]) / 0.123 * Zwind[hh]) / log((Zwind[hh] - 2.0 / 3.0 * Zwind[hh]) / 0.123 * Zwind[hh]); + double snow_traj_angle = atan(u_1_third_Ht[hh] / v_snow); // in radians + Cp_inc = (1 - Cc[hh]) / (1 + exp(-k_cp * (sin(snow_traj_angle) - (1 - Cc[hh])))); // fractional increas in leaf contact area (Clca) based on horizontal trajectory. This is modified from Cebulski & Pomeroy snow interception paper. Has only been tested on forest plots with Cc of .3 and .5. + } + Clca[hh] = Cc[hh] + Cp_inc; // calculated leaf contact area (Clca) based on trajectory angle + } + else + { + 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) - dL = IP * hru_snow[hh]; // change in canopy snow load + IP = Clca[hh] * alpha[hh]; // interception efficiency (IP) + dL = IP * hru_snow[hh]; // change in canopy snow load - Snow_load[hh] += dL; + Snow_load[hh] += dL; - // calculate canopy snow throughfall before unloading: + // calculate canopy snow throughfall before unloading: - direct_snow[hh] += (1-IP) * hru_snow[hh]; + direct_snow[hh] += (1 - IP) * hru_snow[hh]; - // Finished initial loading now start the ablation paramaterisations + // Finished initial loading now start the ablation paramaterisations - double Vi = 0; // submilation rate coefficient for exposed snow (s-1) - double IceBulbT = 0; // ice bulb temperature of canopy snow - double U = 0; // unloading rate coefficient used in Hedstrom & Pomeroy 1998 param. + double Vi = 0; // submilation rate coefficient for exposed snow (s-1) + double IceBulbT = 0; // ice bulb temperature of canopy snow + double U = 0; // unloading rate coefficient used in Hedstrom & Pomeroy 1998 param. -switch(SublimationSwitch[hh]){ -case 0: // do not sublimate, used for debugging or experiments, recommend using case 1 otherwise. + switch (SublimationSwitch[hh]) + { + case 0: + { // do not sublimate, used for debugging or experiments, recommend using case 1 otherwise. - std::cout << "SublimationSwitch Case 0: No canopy snow sublimation applied.\n"; - break; + std::cout << "SublimationSwitch Case 0: No canopy snow sublimation applied.\n"; + break; + } // case 0 -case 1: // canopy snow sublimation with Pomeroy et al. (1998) paramaterisation + case 1: + { // canopy snow sublimation with Pomeroy et al. (1998) paramaterisation -// Pomeroy et al. (1998) sublimation routine modified by Alex Cebulski to change maximum intercepted load from Hedstrom 1998 calculation to constant based on observations from several studies showing much higher intercepted loads. -// calculate snow ventilation windspeed: + // Pomeroy et al. (1998) sublimation routine modified by Alex Cebulski to change maximum intercepted load from Hedstrom 1998 calculation to constant based on observations from several studies showing much higher intercepted loads. + // calculate snow ventilation windspeed: -//============================================================================= - const double gamma = 1.15; - const double AlbedoIce = 0.8; // albedo of ideal ice sphere - const double Radius = 5.0e-4; // radii of single 'ideal' ice sphere in, m) - const double KinVisc = 1.88e-5; // kinematic viscosity of air (Sask. avg. value) - const double ks = 0.0114; // snow shape coefficient for jack pine - const double Fract = 0.37; // fractal dimension of intercepted snow - // const double Hs = 2.838e6; // heat of sublimation (MJ/kg) changed to PBSM_constants::LATH which has same value. -//============================================================================== + //============================================================================= + const double gamma = 1.15; + const double AlbedoIce = 0.8; // albedo of ideal ice sphere + const double Radius = 5.0e-4; // radii of single 'ideal' ice sphere in, m) + const double KinVisc = 1.88e-5; // kinematic viscosity of air (Sask. avg. value) + const double ks = 0.0114; // snow shape coefficient for jack pine + const double Fract = 0.37; // fractal dimension of intercepted snow + // const double Hs = 2.838e6; // heat of sublimation (MJ/kg) changed to PBSM_constants::LATH which has same value. + //============================================================================== - double xi2 = 1-Zvent[hh]; + double xi2 = 1 - Zvent[hh]; double windExt2 = (gamma * LAI[hh] * xi2); double uVent = u_FHt[hh] * exp(-1 * windExt2); -// Calculate wind speed at canopy top for sublimation calculation, moved from initial interception chunk by Alex Cebulski November 2024 + // Calculate wind speed at canopy top for sublimation calculation, moved from initial interception chunk by Alex Cebulski November 2024 - if(Ht[hh] - 2.0/3.0*Zwind[hh] > 0.0) - u_FHt[hh] = hru_u[hh]*log((Ht[hh] - 2.0/3.0*Zwind[hh] )/ 0.123*Zwind[hh])/log((Zwind[hh] - 2.0/3.0*Zwind[hh] )/ 0.123*Zwind[hh]); + if (Ht[hh] - 2.0 / 3.0 * Zwind[hh] > 0.0) + u_FHt[hh] = hru_u[hh] * log((Ht[hh] - 2.0 / 3.0 * Zwind[hh]) / 0.123 * Zwind[hh]) / log((Zwind[hh] - 2.0 / 3.0 * Zwind[hh]) / 0.123 * Zwind[hh]); else u_FHt[hh] = 0.0; double I1 = 0.0; -// calculate sublimation of intercepted snow from ideal intercepted ice sphere (500 microns diameter): + // calculate sublimation of intercepted snow from ideal intercepted ice sphere (500 microns diameter): double Alpha, A1, B1, C1, J, D, Lamb, Mpm, Nu, Nr, SStar, Sigma2; - double Es = 611.15f * exp(22.452f*hru_t[hh]/(hru_t[hh] + 273.0f)); // {sat pressure} + double Es = 611.15f * exp(22.452f * hru_t[hh] / (hru_t[hh] + 273.0f)); // {sat pressure} - double SvDens = Es*PBSM_constants::M/(PBSM_constants::R*(hru_t[hh] + 273.0f)); // {sat density} + double SvDens = Es * PBSM_constants::M / (PBSM_constants::R * (hru_t[hh] + 273.0f)); // {sat density} - Lamb = 6.3e-4*(hru_t[hh]+273.0) + 0.0673; // thermal conductivity of atmosphere - Nr = 2.0 * Radius * uVent / KinVisc; // Reynolds number - Nu = 1.79 + 0.606 * sqrt((double) Nr); // Nusselt number - SStar = M_PI * sqr(Radius) * (1.0f - AlbedoIce) * Qsi_; // SW to snow particle !!!! changed + Lamb = 6.3e-4 * (hru_t[hh] + 273.0) + 0.0673; // thermal conductivity of atmosphere + Nr = 2.0 * Radius * uVent / KinVisc; // Reynolds number + Nu = 1.79 + 0.606 * sqrt((double)Nr); // Nusselt number + SStar = M_PI * sqr(Radius) * (1.0f - AlbedoIce) * Qsi_; // SW to snow particle !!!! changed A1 = Lamb * (hru_t[hh] + 273) * Nu; - B1 = PBSM_constants::LATH * PBSM_constants::M /(PBSM_constants::R * (hru_t[hh] + 273.0f))- 1.0; - J = B1/A1; - Sigma2 = hru_rh[hh]/100 -1; - D = 2.06e-5* pow((hru_t[hh]+273.0f)/273.0f, -1.75f); // diffusivity of water vapour - C1 = 1.0/(D*SvDens*Nu); + B1 = PBSM_constants::LATH * PBSM_constants::M / (PBSM_constants::R * (hru_t[hh] + 273.0f)) - 1.0; + J = B1 / A1; + Sigma2 = hru_rh[hh] / 100 - 1; + D = 2.06e-5 * pow((hru_t[hh] + 273.0f) / 273.0f, -1.75f); // diffusivity of water vapour + C1 = 1.0 / (D * SvDens * Nu); Alpha = 5.0; Mpm = 4.0 / 3.0 * M_PI * PBSM_constants::DICE * Radius * Radius * Radius; // 18Mar2022: remove Gamma Distribution Correction term, *(1.0 + 3.0/Alpha + 2.0/sqr(Alpha)); -// sublimation rate of single 'ideal' ice sphere: + // sublimation rate of single 'ideal' ice sphere: - double Vs = (2.0* M_PI* Radius*Sigma2 - SStar* J)/(PBSM_constants::LATH* J + C1)/Mpm; + double Vs = (2.0 * M_PI * Radius * Sigma2 - SStar * J) / (PBSM_constants::LATH * J + C1) / Mpm; pot_subl_cpy[hh] = Vs; // export the dimensionless sublimation rate (s-1) added by alex 2023-07-21 -// snow exposure coefficient (Ce): + // snow exposure coefficient (Ce): double Ce; - if ((Snow_load[hh]/Lmax[hh]) <= 0.0) + if ((Snow_load[hh] / Lmax[hh]) <= 0.0) Ce = 0.07; else - Ce = ks* pow((Snow_load[hh]/Lmax[hh]), -Fract); + Ce = ks * pow((Snow_load[hh] / Lmax[hh]), -Fract); -// calculate 'potential' canopy sublimation: + // calculate 'potential' canopy sublimation: - Vi = Vs*Ce; + Vi = Vs * Ce; -// limit sublimation to canopy snow available and take sublimated snow away from canopy snow at timestep start + // limit sublimation to canopy snow available and take sublimated snow away from canopy snow at timestep start - Subl_Cpy[hh] = -Snow_load[hh]*Vi*PBSM_constants::LATH*Global::Interval*24*3600/PBSM_constants::LATH; // make W/m2 - if(Subl_Cpy[hh] > Snow_load[hh]){ + Subl_Cpy[hh] = -Snow_load[hh] * Vi * PBSM_constants::LATH * Global::Interval * 24 * 3600 / PBSM_constants::LATH; // make W/m2 + if (Subl_Cpy[hh] > Snow_load[hh]) + { Subl_Cpy[hh] = Snow_load[hh]; Snow_load[hh] = 0.0; } - else{ + else + { Snow_load[hh] -= Subl_Cpy[hh]; - if(Snow_load[hh] < 0.0) + if (Snow_load[hh] < 0.0) Snow_load[hh] = 0.0; } -} - -if(MassUnloadingSwitch[hh] == 0 || MeltwaterSwitch[hh] == 0){ - // calculate 'ice-bulb' temperature of intercepted snow: - double IceBulbT = hru_t[hh] - (Vi* PBSM_constants::LATH/1e6/CRHM_constants::ci); - const double U = -1 * log(0.678) / (24 * 7 * Global::Freq / 24); // weekly dimensionless unloading coefficient -> to CRHM time interval // 21Mar2022 correction: invert the term 24/Global::Freq, use unloading rate coefficient U = -log(c)/t for snow unloading determined by inverse function of c = e^(-Ut) = 0.678 based on Eq. 14 in Hedstrom and Pomeroy (1998) -} + break; + } // case 1 + } // sublimation switch + + // Mass Unloading Section + // Enters different parameterisations for mass unloading of canopy snow based on switch. + // ============================================================================= + + if (MassUnloadingSwitch[hh] == 0 || MeltwaterSwitch[hh] == 0) + { + // calculate 'ice-bulb' temperature of intercepted snow: + IceBulbT = hru_t[hh] - (Vi * PBSM_constants::LATH / 1e6 / CRHM_constants::ci); + const double U = -1 * log(0.678) / (24 * 7 * Global::Freq / 24); // weekly dimensionless unloading coefficient -> to CRHM time interval // 21Mar2022 correction: invert the term 24/Global::Freq, use unloading rate coefficient U = -log(c)/t for snow unloading determined by inverse function of c = e^(-Ut) = 0.678 based on Eq. 14 in Hedstrom and Pomeroy (1998) + } -switch(MassUnloadingSwitch[hh]){ - case 0: // This is the mass unloading portion of the latest iteration of the Hedstrom & Pomeroy 1998 unloading with modifications by Ellis et al. (2010) and Floyd (2012). Generally used with MeltwaterSwitch == 0. - - // determine whether canopy snow is unloaded as mass clumps: + switch (MassUnloadingSwitch[hh]) + { + case 0: { // This is the mass unloading portion of the latest iteration of the Hedstrom & Pomeroy 1998 unloading with modifications by Ellis et al. (2010) and Floyd (2012). Generally used with MeltwaterSwitch == 0. - if(IceBulbT < unload_t[hh]){ // has to be at least one interval. Trip on half step + // determine whether canopy snow is unloaded as mass clumps: + if (IceBulbT < unload_t[hh]) + { // has to be at least one interval. Trip on half step SUnload[hh] = Snow_load[hh] * U; // the dimensionless unloading coefficient already /interval, 21Mar2022 correction: use unloading rate coefficient U } + break; + } // case 0 + + case 1: { // This is the updated mass snow unloading parameterisations from Cebulski & Pomeroy to unload based on time, wind, and air temperature. + // temperature induced unloading + const double a_T = 2.584003e-05; // Cebulski & Pomeroy coef from exponential function of unloading + drip and air temp measurements at Fortress mountain when wind speed <= 1 m/s. + const double b_T = 1.646875e-01; // Cebulski & Pomeroy coef from exponential function of unloading + drip and air temp measurements at Fortress mountain when wind speed <= 1 m/s. + + double fT = a_T * exp(b_T * hru_t[hh]); // unloading rate based on warming of snow in the canopy (s-1), still need to partition out the portion of this that is drip vs mass unloading + + // mechanical wind induced unloading + const double a_u = 5.204024e-06; // Cebulski & Pomeroy coef from exponential function of unloading + drip and wind speed measurements at Fortress mountain when air temp < -6 C. + const double b_u = 7.363594e-02; // Cebulski & Pomeroy coef from exponential function of unloading + drip and wind speed measurements at Fortress mountain when air temp < -6 C. + double Ht_mid = Ht[hh] * (1.0 / 2.0); // half of canopy height + double u_mid = 0; + double fu = 0; + if ((Ht_mid - (2.0 / 3.0) * Zwind[hh]) > 0.0) + { + u_mid = hru_u[hh] * log((Ht_mid - (2.0 / 3.0) * Zwind[hh]) / 0.123 * Zwind[hh]) / log((Zwind[hh] - 2.0 / 3.0 * Zwind[hh]) / 0.123 * Zwind[hh]); + double fu = u_mid * a_u * exp(b_u * u_mid); // unloading rate due to wind (s-1) + } - case 1: // This is the updated mass snow unloading parameterisations from Cebulski & Pomeroy to unload based on time, wind, and air temperature. - // temperature induced unloading - const double a_T = 2.584003e-05; // Cebulski & Pomeroy coef from exponential function of unloading + drip and air temp measurements at Fortress mountain when wind speed <= 1 m/s. - const double b_T = 1.646875e-01; // Cebulski & Pomeroy coef from exponential function of unloading + drip and air temp measurements at Fortress mountain when wind speed <= 1 m/s. - - double fT = a_T * exp(b_T * hru_t[hh]); // unloading rate based on warming of snow in the canopy (s-1), still need to partition out the portion of this that is drip vs mass unloading - - // mechanical wind induced unloading - const double a_u = 5.204024e-06; // Cebulski & Pomeroy coef from exponential function of unloading + drip and wind speed measurements at Fortress mountain when air temp < -6 C. - const double b_u = 7.363594e-02; // Cebulski & Pomeroy coef from exponential function of unloading + drip and wind speed measurements at Fortress mountain when air temp < -6 C. - double Ht_mid = Ht[hh] * (1.0 / 2.0); // half of canopy height - double u_mid = 0; - double fu = 0; - if ((Ht_mid - (2.0 / 3.0) * Zwind[hh]) > 0.0) - { - u_mid = hru_u[hh] * log((Ht_mid - (2.0 / 3.0) * Zwind[hh]) / 0.123 * Zwind[hh]) / log((Zwind[hh] - 2.0 / 3.0 * Zwind[hh]) / 0.123 * Zwind[hh]); - double fu = u_mid * a_u * exp(b_u * u_mid); // unloading rate due to wind (s-1) - } + // duration based unloading + const double a_t = 8.194345e-06; // Cebulski & Pomeroy coef from exponential function of unloading + drip and duration snow has been intercepted in the canopy at Fortress mountain when wind speed <= 1 m/s and air temperature < -6 C. + const double b_t = -1.540050e+02; // Cebulski & Pomeroy coef from exponential function of unloading + drip and duration snow has been intercepted in the canopy at Fortress mountain when wind speed <= 1 m/s and air temperature < -6 C. - // duration based unloading - const double a_t = 8.194345e-06; // Cebulski & Pomeroy coef from exponential function of unloading + drip and duration snow has been intercepted in the canopy at Fortress mountain when wind speed <= 1 m/s and air temperature < -6 C. - const double b_t = -1.540050e+02; // Cebulski & Pomeroy coef from exponential function of unloading + drip and duration snow has been intercepted in the canopy at Fortress mountain when wind speed <= 1 m/s and air temperature < -6 C. + double t_snow_in_canopy = 12 * 60 * 60; // duration snow intercepted in the canopy, set to constant of 12 hours. TODO need to track duration that snow has been intercepted in the canopy. Can do this similar to snowpack albedo calculation. - double t_snow_in_canopy = 12 * 60 * 60; // duration snow intercepted in the canopy, set to constant of 12 hours. TODO need to track duration that snow has been intercepted in the canopy. Can do this similar to snowpack albedo calculation. - - double ft = a_t * exp(b_t * t_snow_in_canopy); + double ft = a_t * exp(b_t * t_snow_in_canopy); - // ablation via temperature, wind, and duration based unloading - double dt = Global::Interval * 24 * 60 * 60; // converts the interval which is a time period (i.e., time/cycles, 1 day/# obs) to timestep in seconds. - SUnload[hh] = Snow_load[hh] * (fT + fu + ft) * dt; // converts our -} + // ablation via temperature, wind, and duration based unloading + double dt = Global::Interval * 24 * 60 * 60; // converts the interval which is a time period (i.e., time/cycles, 1 day/# obs) to timestep in seconds. + SUnload[hh] = Snow_load[hh] * (fT + fu + ft) * dt; // converts our + break; + } // case 1 + } // MassUnloadingSwitch -// handle mass unloading regardless of what parameterisation is chosen -if (SUnload[hh] > Snow_load[hh]) -{ - SUnload[hh] = Snow_load[hh]; - Snow_load[hh] = 0.0; -} -else - Snow_load[hh] -= SUnload[hh]; + // handle mass unloading regardless of what parameterisation is chosen + if (SUnload[hh] > Snow_load[hh]) + { + SUnload[hh] = Snow_load[hh]; + Snow_load[hh] = 0.0; + } + else + { + Snow_load[hh] -= SUnload[hh]; + } + + cum_SUnload[hh] += SUnload[hh]; -cum_SUnload[hh] += SUnload[hh]; + // Meltwater (drip) Section + // Enters different parameterisations for meltwater drip of canopy snow based on switch. + // ============================================================================= -switch(MeltwaterSwitch[hh]) { - case 0: { // Block for case 0 - // This is the meltwater drip portion of the latest iteration of the Hedstrom & Pomeroy 1998 unloading with modifications by Ellis et al. (2010) and Floyd (2012). - double Six_Hour_Divisor = Global::Freq / 4.0; // Unload over 6 hours + switch (MeltwaterSwitch[hh]) + { + case 0: + { // Block for case 0 + // This is the meltwater drip portion of the latest iteration of the Hedstrom & Pomeroy 1998 unloading with modifications by Ellis et al. (2010) and Floyd (2012). + double Six_Hour_Divisor = Global::Freq / 4.0; // Unload over 6 hours - if (IceBulbT >= unload_t_water[hh]) { + if (IceBulbT >= unload_t_water[hh]) + { drip_Cpy[hh] = Snow_load[hh] / Six_Hour_Divisor; SUnload_H2O[hh] = drip_Cpy[hh]; Snow_load[hh] -= SUnload_H2O[hh]; cum_SUnload_H2O[hh] += SUnload_H2O[hh]; - } - else if (IceBulbT >= unload_t[hh]) { + } + else if (IceBulbT >= unload_t[hh]) + { SUnload[hh] = Snow_load[hh] / Six_Hour_Divisor; Snow_load[hh] -= SUnload[hh]; cum_SUnload[hh] += SUnload[hh]; + } + break; } - break; - } - case 1: { // Block for case 1 - // This is the meltwater drip parameterisation from CLASS - drip_Cpy[hh] = 0; - SUnload_H2O[hh] = drip_Cpy[hh]; - Snow_load[hh] -= SUnload_H2O[hh]; - cum_SUnload_H2O[hh] += SUnload_H2O[hh]; - break; - } -} + case 1: + { // Block for case 1 + // This is the meltwater drip parameterisation from CLASS + drip_Cpy[hh] = 0; + SUnload_H2O[hh] = drip_Cpy[hh]; + Snow_load[hh] -= SUnload_H2O[hh]; + cum_SUnload_H2O[hh] += SUnload_H2O[hh]; + break; + } + } -// calculate total sub-canopy snow: + // calculate total sub-canopy snow: - net_snow[hh] = direct_snow[hh] + SUnload[hh]; - break; - } // case canopy - case 1: // clearing - case 2: // gap - net_snow[hh] = hru_snow[hh]; - net_rain[hh] = hru_rain[hh]; + net_snow[hh] = direct_snow[hh] + SUnload[hh]; break; - } // switch + } // case canopy + case 1: // clearing + case 2: // gap + net_snow[hh] = hru_snow[hh]; + net_rain[hh] = hru_rain[hh]; + break; + } // switch -// calculate horizontal canopy-coverage (Cc): + // calculate horizontal canopy-coverage (Cc): double smax, Q; // cannot be in switch structure - switch(CanopyClearing[hh]){ + switch (CanopyClearing[hh]) + { - case 0: // canopy + case 0: // canopy - Cc[hh] = std::min(1.0, 0.29 * std::log(LAI[hh]) + 0.55); // update by Alex so Cc is not greater than 1 - if(Cc[hh] <= 0.0) - Cc[hh] = 0.0; + Cc[hh] = std::min(1.0, 0.29 * std::log(LAI[hh]) + 0.55); // update by Alex so Cc is not greater than 1 + if (Cc[hh] <= 0.0) + Cc[hh] = 0.0; - smax = Cc[hh]*LAI[hh]*0.2; + smax = Cc[hh] * LAI[hh] * 0.2; -// Forest rain interception and evaporation model -// 'sparse' Rutter interception model (i.e. Valente 1997): + // Forest rain interception and evaporation model + // 'sparse' Rutter interception model (i.e. Valente 1997): -// calculate direct throughfall: + // calculate direct throughfall: - if(hru_rain[hh] > 0.0){ + if (hru_rain[hh] > 0.0) + { - direct_rain[hh] = hru_rain[hh] * (1-Cc[hh]); + direct_rain[hh] = hru_rain[hh] * (1 - Cc[hh]); - // calculate rain accumulation on canopy before evap loss: + // calculate rain accumulation on canopy before evap loss: - if (rain_load[hh] + hru_rain[hh]*Cc[hh] > smax){ - drip_Cpy[hh] += (rain_load[hh] + hru_rain[hh]*Cc[hh] - smax); - rain_load[hh] = smax; - } - else - rain_load[hh] += hru_rain[hh]*Cc[hh]; + if (rain_load[hh] + hru_rain[hh] * Cc[hh] > smax) + { + drip_Cpy[hh] += (rain_load[hh] + hru_rain[hh] * Cc[hh] - smax); + rain_load[hh] = smax; + } + else + rain_load[hh] += hru_rain[hh] * Cc[hh]; - } // if hru_rain[hh] > 0.0 + } // if hru_rain[hh] > 0.0 -// calculate 'actual evap' of water from canopy and canopy storage after evaporation:: + // calculate 'actual evap' of water from canopy and canopy storage after evaporation:: - if(rain_load[hh] > 0.0){ - if(inhibit_evap[hh] == 0){ // use Granger when no snowcover - if(rain_load[hh] >= hru_evap[hh]*Cc[hh]){ // (evaporation in mm) - intcp_evap[hh] = hru_evap[hh]*Cc[hh]; // - rain_load[hh] -= hru_evap[hh]*Cc[hh]; - } - else{ - intcp_evap[hh] = rain_load[hh]; - rain_load[hh] = 0.0; - } + if (rain_load[hh] > 0.0) + { + if (inhibit_evap[hh] == 0) + { // use Granger when no snowcover + if (rain_load[hh] >= hru_evap[hh] * Cc[hh]) + { // (evaporation in mm) + intcp_evap[hh] = hru_evap[hh] * Cc[hh]; // + rain_load[hh] -= hru_evap[hh] * Cc[hh]; } - else{ // use Priestley-Taylor when snowcover - Q = Qsi_*86400/Global::Freq/1e6/lambda(hru_t[hh]); // convert w/m2 to mm/m2/int + else + { + intcp_evap[hh] = rain_load[hh]; + rain_load[hh] = 0.0; + } + } + else + { // use Priestley-Taylor when snowcover + Q = Qsi_ * 86400 / Global::Freq / 1e6 / lambda(hru_t[hh]); // convert w/m2 to mm/m2/int - if(Qsi_ > 0.0) - Pevap[hh] = 1.26*delta(hru_t[hh])*Q/(delta(hru_t[hh]) + gamma(Pa[hh], hru_t[hh])); - else - Pevap[hh] = 0.0; + if (Qsi_ > 0.0) + Pevap[hh] = 1.26 * delta(hru_t[hh]) * Q / (delta(hru_t[hh]) + gamma(Pa[hh], hru_t[hh])); + else + Pevap[hh] = 0.0; - if(rain_load[hh] >= Pevap[hh]*Cc[hh]){ // (evaporation in mm) - intcp_evap[hh] = Pevap[hh]*Cc[hh]; // check - rain_load[hh] -= Pevap[hh]*Cc[hh]; - } - else{ - intcp_evap[hh] = rain_load[hh]; // check - rain_load[hh] = 0.0; - } + if (rain_load[hh] >= Pevap[hh] * Cc[hh]) + { // (evaporation in mm) + intcp_evap[hh] = Pevap[hh] * Cc[hh]; // check + rain_load[hh] -= Pevap[hh] * Cc[hh]; + } + else + { + intcp_evap[hh] = rain_load[hh]; // check + rain_load[hh] = 0.0; } - } // if rain_load[hh] > 0.0 + } + } // if rain_load[hh] > 0.0 -// cumulative amounts.... + // cumulative amounts.... - net_rain[hh] = direct_rain[hh] + drip_Cpy[hh]; + net_rain[hh] = direct_rain[hh] + drip_Cpy[hh]; - cum_intcp_evap[hh] += intcp_evap[hh]; - cum_Subl_Cpy[hh] += Subl_Cpy[hh]; - break; - default : // clearing and gap - break; - } // switch + cum_intcp_evap[hh] += intcp_evap[hh]; + cum_Subl_Cpy[hh] += Subl_Cpy[hh]; + break; + default: // clearing and gap + break; + } // switch net_p[hh] = net_rain[hh] + net_snow[hh]; cum_net_rain[hh] += net_rain[hh]; @@ -747,8 +792,10 @@ switch(MeltwaterSwitch[hh]) { } // end for } -void ClassCRHMCanopyVectorBased::finish(bool good) { - for(hh = 0; chkStruct(); ++hh) { +void ClassCRHMCanopyVectorBased::finish(bool good) +{ + for (hh = 0; chkStruct(); ++hh) + { LogMessageA(hh, string("'" + Name + " (CanopyVectorBased)' cum_net_snow (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_net_snow[hh], hru_area[hh], basin_area[0]); LogMessageA(hh, string("'" + Name + " (CanopyVectorBased)' cum_net_rain (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_net_rain[hh], hru_area[hh], basin_area[0]); LogMessageA(hh, string("'" + Name + " (CanopyVectorBased)' cum_Subl_Cpy (mm) (mm*hru) (mm*hru/basin): ").c_str(), cum_Subl_Cpy[hh], hru_area[hh], basin_area[0]); @@ -761,23 +808,23 @@ void ClassCRHMCanopyVectorBased::finish(bool good) { double ClassCRHMCanopyVectorBased::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)); + 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)); + return (3549.0 * exp(21.88 * t / (t + 265.5)) / sqr(t + 265.5)); } double ClassCRHMCanopyVectorBased::lambda(double t) // Latent heat of vaporization (mJ/(kg DEGREE_CELSIUS)) { - return( 2.501 - 0.002361 * t ); + return (2.501 - 0.002361 * t); } double ClassCRHMCanopyVectorBased::gamma(double Pa, double t) // Psychrometric constant (kPa/DEGREE_CELSIUS) { - return( 0.00163 * Pa / lambda(t)); // lambda (mJ/(kg DEGREE_CELSIUS)) + return (0.00163 * Pa / lambda(t)); // lambda (mJ/(kg DEGREE_CELSIUS)) } double ClassCRHMCanopyVectorBased::RHOa(double t, double ea, double Pa) // atmospheric density (kg/m^3) { const double R = 2870; - return (1E4*Pa /(R*( 273.15 + t))*(1.0 - 0.379*(ea/Pa)) ); // + return (1E4 * Pa / (R * (273.15 + t)) * (1.0 - 0.379 * (ea / Pa))); // } From 54fa98318d28d003613b73c62c0623446b6b98da Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Mon, 18 Nov 2024 19:04:54 -0800 Subject: [PATCH 12/15] put hp98/ellis unloading and drip back together within the mass unloading switch. now we are getting exactly the same as before adding the switches --- .../modules/ClassCRHMCanopyVectorBased.cpp | 188 ++++++++++-------- 1 file changed, 106 insertions(+), 82 deletions(-) diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp index 7d7ef3f25..b8bcd97b6 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp @@ -582,84 +582,17 @@ void ClassCRHMCanopyVectorBased::run(void) // Mass Unloading Section // Enters different parameterisations for mass unloading of canopy snow based on switch. // ============================================================================= - - if (MassUnloadingSwitch[hh] == 0 || MeltwaterSwitch[hh] == 0) - { - // calculate 'ice-bulb' temperature of intercepted snow: - IceBulbT = hru_t[hh] - (Vi * PBSM_constants::LATH / 1e6 / CRHM_constants::ci); - const double U = -1 * log(0.678) / (24 * 7 * Global::Freq / 24); // weekly dimensionless unloading coefficient -> to CRHM time interval // 21Mar2022 correction: invert the term 24/Global::Freq, use unloading rate coefficient U = -log(c)/t for snow unloading determined by inverse function of c = e^(-Ut) = 0.678 based on Eq. 14 in Hedstrom and Pomeroy (1998) - } switch (MassUnloadingSwitch[hh]) { - case 0: { // This is the mass unloading portion of the latest iteration of the Hedstrom & Pomeroy 1998 unloading with modifications by Ellis et al. (2010) and Floyd (2012). Generally used with MeltwaterSwitch == 0. - - // determine whether canopy snow is unloaded as mass clumps: - if (IceBulbT < unload_t[hh]) - { // has to be at least one interval. Trip on half step - SUnload[hh] = Snow_load[hh] * U; // the dimensionless unloading coefficient already /interval, 21Mar2022 correction: use unloading rate coefficient U - } - break; - } // case 0 - - case 1: { // This is the updated mass snow unloading parameterisations from Cebulski & Pomeroy to unload based on time, wind, and air temperature. - // temperature induced unloading - const double a_T = 2.584003e-05; // Cebulski & Pomeroy coef from exponential function of unloading + drip and air temp measurements at Fortress mountain when wind speed <= 1 m/s. - const double b_T = 1.646875e-01; // Cebulski & Pomeroy coef from exponential function of unloading + drip and air temp measurements at Fortress mountain when wind speed <= 1 m/s. - - double fT = a_T * exp(b_T * hru_t[hh]); // unloading rate based on warming of snow in the canopy (s-1), still need to partition out the portion of this that is drip vs mass unloading - - // mechanical wind induced unloading - const double a_u = 5.204024e-06; // Cebulski & Pomeroy coef from exponential function of unloading + drip and wind speed measurements at Fortress mountain when air temp < -6 C. - const double b_u = 7.363594e-02; // Cebulski & Pomeroy coef from exponential function of unloading + drip and wind speed measurements at Fortress mountain when air temp < -6 C. - double Ht_mid = Ht[hh] * (1.0 / 2.0); // half of canopy height - double u_mid = 0; - double fu = 0; - if ((Ht_mid - (2.0 / 3.0) * Zwind[hh]) > 0.0) - { - u_mid = hru_u[hh] * log((Ht_mid - (2.0 / 3.0) * Zwind[hh]) / 0.123 * Zwind[hh]) / log((Zwind[hh] - 2.0 / 3.0 * Zwind[hh]) / 0.123 * Zwind[hh]); - double fu = u_mid * a_u * exp(b_u * u_mid); // unloading rate due to wind (s-1) - } - - // duration based unloading - const double a_t = 8.194345e-06; // Cebulski & Pomeroy coef from exponential function of unloading + drip and duration snow has been intercepted in the canopy at Fortress mountain when wind speed <= 1 m/s and air temperature < -6 C. - const double b_t = -1.540050e+02; // Cebulski & Pomeroy coef from exponential function of unloading + drip and duration snow has been intercepted in the canopy at Fortress mountain when wind speed <= 1 m/s and air temperature < -6 C. - - double t_snow_in_canopy = 12 * 60 * 60; // duration snow intercepted in the canopy, set to constant of 12 hours. TODO need to track duration that snow has been intercepted in the canopy. Can do this similar to snowpack albedo calculation. - - double ft = a_t * exp(b_t * t_snow_in_canopy); - - // ablation via temperature, wind, and duration based unloading - double dt = Global::Interval * 24 * 60 * 60; // converts the interval which is a time period (i.e., time/cycles, 1 day/# obs) to timestep in seconds. - SUnload[hh] = Snow_load[hh] * (fT + fu + ft) * dt; // converts our - break; - } // case 1 - } // MassUnloadingSwitch - - // handle mass unloading regardless of what parameterisation is chosen - if (SUnload[hh] > Snow_load[hh]) - { - SUnload[hh] = Snow_load[hh]; - Snow_load[hh] = 0.0; - } - else - { - Snow_load[hh] -= SUnload[hh]; - } - - cum_SUnload[hh] += SUnload[hh]; - - // Meltwater (drip) Section - // Enters different parameterisations for meltwater drip of canopy snow based on switch. - // ============================================================================= - - switch (MeltwaterSwitch[hh]) - { case 0: - { // Block for case 0 - // This is the meltwater drip portion of the latest iteration of the Hedstrom & Pomeroy 1998 unloading with modifications by Ellis et al. (2010) and Floyd (2012). - double Six_Hour_Divisor = Global::Freq / 4.0; // Unload over 6 hours + { // This is the mass unloading portion of the latest iteration of the Hedstrom & Pomeroy 1998 unloading with modifications by Ellis et al. (2010) and Floyd (2012). Generally used with MeltwaterSwitch == 0. + // calculate 'ice-bulb' temperature of intercepted snow: + IceBulbT = hru_t[hh] - (Vi * PBSM_constants::LATH / 1e6 / CRHM_constants::ci); + const double U = -1 * log(0.678) / (24 * 7 * Global::Freq / 24); // weekly dimensionless unloading coefficient -> to CRHM time interval // 21Mar2022 correction: invert the term 24/Global::Freq, use unloading rate coefficient U = -log(c)/t for snow unloading determined by inverse function of c = e^(-Ut) = 0.678 based on Eq. 14 in Hedstrom and Pomeroy (1998) + double Six_Hour_Divisor = Global::Freq / 4.0; // Unload over 6 hours + // determine whether canopy snow is unloaded as meltwater or mass clumps: if (IceBulbT >= unload_t_water[hh]) { drip_Cpy[hh] = Snow_load[hh] / Six_Hour_Divisor; @@ -673,19 +606,110 @@ void ClassCRHMCanopyVectorBased::run(void) Snow_load[hh] -= SUnload[hh]; cum_SUnload[hh] += SUnload[hh]; } + else if (IceBulbT < unload_t[hh]) + { // has to be at least one interval. Trip on half step + SUnload[hh] = Snow_load[hh] * U; // the dimensionless unloading coefficient already /interval, 21Mar2022 correction: use unloading rate coefficient U + if (SUnload[hh] > Snow_load[hh]) + { + SUnload[hh] = Snow_load[hh]; + Snow_load[hh] = 0.0; + } + else + Snow_load[hh] -= SUnload[hh]; + + cum_SUnload[hh] += SUnload[hh]; + } break; - } + } // case 0 case 1: - { // Block for case 1 - // This is the meltwater drip parameterisation from CLASS - drip_Cpy[hh] = 0; - SUnload_H2O[hh] = drip_Cpy[hh]; - Snow_load[hh] -= SUnload_H2O[hh]; - cum_SUnload_H2O[hh] += SUnload_H2O[hh]; + { // This is the updated mass snow unloading parameterisations from Cebulski & Pomeroy to unload based on time, wind, and air temperature. + // temperature induced unloading + const double a_T = 2.584003e-05; // Cebulski & Pomeroy coef from exponential function of unloading + drip and air temp measurements at Fortress mountain when wind speed <= 1 m/s. + const double b_T = 1.646875e-01; // Cebulski & Pomeroy coef from exponential function of unloading + drip and air temp measurements at Fortress mountain when wind speed <= 1 m/s. + + double fT = a_T * exp(b_T * hru_t[hh]); // unloading rate based on warming of snow in the canopy (s-1), still need to partition out the portion of this that is drip vs mass unloading + + // mechanical wind induced unloading + const double a_u = 5.204024e-06; // Cebulski & Pomeroy coef from exponential function of unloading + drip and wind speed measurements at Fortress mountain when air temp < -6 C. + const double b_u = 7.363594e-02; // Cebulski & Pomeroy coef from exponential function of unloading + drip and wind speed measurements at Fortress mountain when air temp < -6 C. + double Ht_mid = Ht[hh] * (1.0 / 2.0); // half of canopy height + double u_mid = 0; + double fu = 0; + if ((Ht_mid - (2.0 / 3.0) * Zwind[hh]) > 0.0) + { + u_mid = hru_u[hh] * log((Ht_mid - (2.0 / 3.0) * Zwind[hh]) / 0.123 * Zwind[hh]) / log((Zwind[hh] - 2.0 / 3.0 * Zwind[hh]) / 0.123 * Zwind[hh]); + double fu = u_mid * a_u * exp(b_u * u_mid); // unloading rate due to wind (s-1) + } + + // duration based unloading + const double a_t = 8.194345e-06; // Cebulski & Pomeroy coef from exponential function of unloading + drip and duration snow has been intercepted in the canopy at Fortress mountain when wind speed <= 1 m/s and air temperature < -6 C. + const double b_t = -1.540050e+02; // Cebulski & Pomeroy coef from exponential function of unloading + drip and duration snow has been intercepted in the canopy at Fortress mountain when wind speed <= 1 m/s and air temperature < -6 C. + + double t_snow_in_canopy = 12 * 60 * 60; // duration snow intercepted in the canopy, set to constant of 12 hours. TODO need to track duration that snow has been intercepted in the canopy. Can do this similar to snowpack albedo calculation. + + double ft = a_t * exp(b_t * t_snow_in_canopy); + + // ablation via temperature, wind, and duration based unloading + double dt = Global::Interval * 24 * 60 * 60; // converts the interval which is a time period (i.e., time/cycles, 1 day/# obs) to timestep in seconds. + SUnload[hh] = Snow_load[hh] * (fT + fu + ft) * dt; // calculate solid snow unloading over the time interval + + if (SUnload[hh] > Snow_load[hh]) + { + SUnload[hh] = Snow_load[hh]; + Snow_load[hh] = 0.0; + } + else + { + Snow_load[hh] -= SUnload[hh]; + } + + cum_SUnload[hh] += SUnload[hh]; break; - } - } + } // case 1 + } // MassUnloadingSwitch + + // handle mass unloading regardless of what parameterisation is chosen + + // Meltwater (drip) Section + // Enters different parameterisations for meltwater drip of canopy snow based on switch. + // ============================================================================= + + // switch (MeltwaterSwitch[hh]) + // { + // case 0: + // { // Block for case 0 + // // This is the meltwater drip portion of the latest iteration of the Hedstrom & Pomeroy 1998 unloading with modifications by Ellis et al. (2010) and Floyd (2012). + // double Six_Hour_Divisor = Global::Freq / 4.0; // Unload over 6 hours + // std::cout << "IceBulbT: " << IceBulbT << "\n"; // Output the value of IceBulbT to the console + // std::cout << "Six_Hour_Divisor: " << Six_Hour_Divisor << "\n"; // Output the value of IceBulbT to the console + + // if (IceBulbT >= unload_t_water[hh]) + // { + // drip_Cpy[hh] = Snow_load[hh] / Six_Hour_Divisor; + // SUnload_H2O[hh] = drip_Cpy[hh]; + // Snow_load[hh] -= SUnload_H2O[hh]; + // cum_SUnload_H2O[hh] += SUnload_H2O[hh]; + // } + // else if (IceBulbT >= unload_t[hh]) + // { + // SUnload[hh] = Snow_load[hh] / Six_Hour_Divisor; + // Snow_load[hh] -= SUnload[hh]; + // cum_SUnload[hh] += SUnload[hh]; + // } + // break; + // } + + // case 1: + // { // Block for case 1 + // // This is the meltwater drip parameterisation from CLASS + // drip_Cpy[hh] = 0; + // SUnload_H2O[hh] = drip_Cpy[hh]; + // Snow_load[hh] -= SUnload_H2O[hh]; + // cum_SUnload_H2O[hh] += SUnload_H2O[hh]; + // break; + // } + // } // calculate total sub-canopy snow: From b729d15cbc0eb7327140be10a2675c7798403245 Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Wed, 20 Nov 2024 11:57:39 -0800 Subject: [PATCH 13/15] updated vector based to take Cc as input --- crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp | 12 ++++-------- crhmcode/src/modules/ClassCRHMCanopyVectorBased.h | 2 +- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp index b8bcd97b6..e5d6b904c 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.cpp @@ -149,9 +149,9 @@ void ClassCRHMCanopyVectorBased::decl(void) declvar("direct_snow", TDim::NHRU, "snow 'direct' through canopy", "(mm/int)", &direct_snow); - declvar("SUnload", TDim::NHRU, "unloaded canopy snow", "(mm)", &SUnload); + declvar("SUnload", TDim::NHRU, "unloaded canopy snow", "(mm/int)", &SUnload); - declvar("SUnload_H2O", TDim::NHRU, "unloaded canopy snow as water", "(mm)", &SUnload_H2O); + declvar("SUnload_H2O", TDim::NHRU, "unloaded canopy snow as water", "(mm/int)", &SUnload_H2O); declstatdiag("cum_SUnload_H2O", TDim::NHRU, "Cumulative unloaded canopy snow as water", "(mm)", &cum_SUnload_H2O); @@ -167,8 +167,6 @@ void ClassCRHMCanopyVectorBased::decl(void) 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); - decldiag("Cc", TDim::NHRU, "Canopy coverage", "()", &Cc); - declvar("Clca", TDim::NHRU, "Leaf contact area adjusted for hydrometeor trajectory angle.", "(s-1)", &Clca); declvar("intcp_evap", TDim::NHRU, "HRU Evaporation from interception", "(mm/int)", &intcp_evap); @@ -197,6 +195,8 @@ void ClassCRHMCanopyVectorBased::decl(void) 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("Lmax", TDim::NHRU, "[50]", "0.0", "100.0", "maximum canopy snow interception load, currently just used for sublimation exposure coef. 50 based on max observed in Storck et al. 2002 and Cebulski & Pomeroy", "(kg/m^2)", &Lmax); @@ -732,10 +732,6 @@ void ClassCRHMCanopyVectorBased::run(void) case 0: // canopy - Cc[hh] = std::min(1.0, 0.29 * std::log(LAI[hh]) + 0.55); // update by Alex so Cc is not greater than 1 - if (Cc[hh] <= 0.0) - Cc[hh] = 0.0; - smax = Cc[hh] * LAI[hh] * 0.2; // Forest rain interception and evaporation model diff --git a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h index 9a2f8578d..3eed3d648 100644 --- a/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h +++ b/crhmcode/src/modules/ClassCRHMCanopyVectorBased.h @@ -73,7 +73,6 @@ double *cum_Subl_Cpy { NULL }; double *cum_SUnload { NULL }; double *cum_SUnload_H2O { NULL }; -double *Cc { NULL }; double *Clca { NULL }; double *k { NULL }; double *Tauc { NULL }; @@ -110,6 +109,7 @@ const double *basin_area { NULL }; // [BASIN] const double *hru_area { NULL }; const double *hru_elev { NULL }; const double *Ht { NULL }; +const double *Cc { NULL }; const double *LAI { NULL }; const double *Lmax { NULL }; const double *alpha { NULL }; From d2d772d1620773eacd64d47dce33ed426a5e5515 Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Wed, 20 Nov 2024 12:08:09 -0800 Subject: [PATCH 14/15] ignore vscode settings.json --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 49e90425b..6ae864a46 100644 --- a/.gitignore +++ b/.gitignore @@ -6,6 +6,7 @@ *.user *.userosscache *.sln.docstates +settings.json # User-specific files (MonoDevelop/Xamarin Studio) *.userprefs From bfce2980aba7ca71b6261515c9392fdad6b81f0b Mon Sep 17 00:00:00 2001 From: Alex Cebulski Date: Sat, 23 Nov 2024 09:24:48 -0800 Subject: [PATCH 15/15] fixed floating point precision problem which was causing segfault when startdate equaled the first obs date --- crhmcode/src/core/CRHMmain.cpp | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/crhmcode/src/core/CRHMmain.cpp b/crhmcode/src/core/CRHMmain.cpp index c5926d666..89465b124 100644 --- a/crhmcode/src/core/CRHMmain.cpp +++ b/crhmcode/src/core/CRHMmain.cpp @@ -2308,7 +2308,7 @@ MMSData * CRHMmain::RunClick2Start() } } - + Global::DTmin = (int)((DTstartR - Global::DTstart)* Global::Freq); Global::DTindx = Global::DTmin; Global::DTnow = Global::DTstart + Global::Interval*((long long)Global::DTindx + 1ll); @@ -2351,13 +2351,19 @@ MMSData * CRHMmain::RunClick2Start() GoodRun = false; } - + ClassData * FileData = NULL; if (ObsFilesList->size() > 0) { FileData = ObsFilesList->begin()->second; + const double tolerance = 1e-9; // to handle floating point precision when DTstartR and FileData->Dt1 are effectively equal but have slightly different precision and then entered the first if statement below. + + if (fabs(DTstartR - FileData->Dt1) < tolerance) { // check if effectively equal + + LogMessageX("DTstartR and FileData->Dt1 are effectively equal but there still may be a small difference due to floating point precision"); + + } else if (DTstartR < FileData->Dt1 - tolerance) { - if (DTstartR < FileData->Dt1) { LogMessageX("Start Time before first Observation"); GoodRun = false; }