Borland Updates March 21, 2022

borland_source/CRHMmain.cpp

@@ -1,10 +1,10 @@
-// 04/12/21 change float to long for use_rho parameter in NewModules.h
-// in SWEslope module 03/15/21
+// 04/06/22 update description for parameters, variables and units in
+// various modules in 03/21/22
 //---------------------------------------------------------------------------
 #include <vcl.h>
 #pragma hdrstop
 
-#define CurrentVersion "04/12/21"
+#define CurrentVersion "04/06/22"
 
 #include <stdexcept>
 #include "CRHMmain.h"

borland_source/NewModules.h

@@ -1,4 +1,4 @@
-// 04/12/21 with changes to 07/02/20
+// 08/11/21 with changes to 04/12/21
 //---------------------------------------------------------------------------
 
 #ifndef OurModulesH
@@ -2657,6 +2657,7 @@ ClassVolumetric(string Name, String Version = "undefined", CRHM::LMODULE Lvl = C
 // declared parameters
 const float *soil_Depth;
 const float *Si_correction;
+const float *fallstat_correction; // 08/11/2021
 const float *soil_moist_max;
 const float *soil_rechr_max;
 const long *soil_type;

crhmcode/src/modules/ClassAyers.cpp

@@ -23,21 +23,21 @@ void ClassAyers::decl(void) {
 
   Description = "'Uses Ayers, 1959 for unfrozen soil. Snow is assumed to melt immediately on contact with the ground.'";
 
-  declvar("infil", TDim::NHRU,"Potential amount of water infiltrating the soil on each HRU", "(mm/int)", &infil);
+  declvar("infil", TDim::NHRU,"interval rainfall infiltration", "(mm/int)", &infil);
 
-  declstatdiag("cuminfil", TDim::NHRU, "cumulative potential infiltration on each HRU", "(mm)", &cuminfil);
+  declstatdiag("cuminfil", TDim::NHRU, "cumulative rainfall infiltration", "(mm)", &cuminfil);
 
-  declvar("runoff", TDim::NHRU, "rainfall runoff", "(mm/int)", &runoff);
+  declvar("runoff", TDim::NHRU, "interval rainfall runoff", "(mm/int)", &runoff);
 
   declstatdiag("cumrunoff", TDim::NHRU, "cumulative rainfall runoff", "(mm)", &cumrunoff);
 
-  declvar("snowinfil", TDim::NHRU, "infiltration", "(mm/int)", &snowinfil);
+  declvar("snowinfil", TDim::NHRU, "interval snowmelt infiltration", "(mm/int)", &snowinfil);
 
-  declstatdiag("cumsnowinfil", TDim::NHRU, "cumulative infiltration", "(mm)", &cumsnowinfil);
+  declstatdiag("cumsnowinfil", TDim::NHRU, "cumulative snowmelt infiltration", "(mm)", &cumsnowinfil);
 
-  declvar("meltrunoff", TDim::NHRU, "melt runoff", "(mm/int)", &meltrunoff);
+  declvar("meltrunoff", TDim::NHRU, "interval snowmelt runoff", "(mm/int)", &meltrunoff);
 
-  declstatdiag("cummeltrunoff", TDim::NHRU, "melt runoff", "(mm)", &cummeltrunoff);
+  declstatdiag("cummeltrunoff", TDim::NHRU, "cumulative snowmelt runoff", "(mm)", &cummeltrunoff);
 
   decllocal("melt_int", TDim::NHRU, "interval melt from snowmelD", "(mm/int)", &melt_int);
 
@@ -46,10 +46,10 @@ void ClassAyers::decl(void) {
   declparam("hru_area", TDim::NHRU, "[1]", "1e-6", "1e+09", "hru area", "(km^2)", &hru_area);
 
   declparam("texture", TDim::NHRU, "[1]", "1","4",
-     "texture: 1 - coarse/medium over coarse, 2 - medium over medium, 3 - medium/fine over fine, 4 - soil over shallow bedrock.", "(%)", &texture);
+     "texture: 1 - coarse/medium over coarse, 2 - medium over medium, 3 - medium/fine over fine, 4 - soil over shallow bedrock.", "()", &texture);
 
   declparam("groundcover", TDim::NHRU, "[1]", "1","6",
-     "groundcover: 1 - bare soil, 2 - row crop, 3 - poor pasture, 4 - small grains, 5 - good pasture, 6 - forested.", "(%)", &groundcover);
+     "groundcover: 1 - bare soil, 2 - row crop, 3 - poor pasture, 4 - small grains, 5 - good pasture, 6 - forested.", "()", &groundcover);
 
   declgetvar("*",  "net_rain", "(mm/int)", &net_rain);
 

crhmcode/src/modules/ClassCRHMCanopy.cpp

@@ -85,7 +85,7 @@ void ClassCRHMCanopy::decl(void) {
 
   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*int)", &Qnsn_Var);
+  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);
 
@@ -97,43 +97,43 @@ void ClassCRHMCanopy::decl(void) {
 
   decldiag("k", TDim::NHRU, "extinction coefficient", "()", &k);
 
-  decldiag("Tauc", TDim::NHRU, "short-wave transmissivity", "(W/m^2)", &Tauc);
+  decldiag("Tauc", TDim::NHRU, "short-wave transmissivity", "()", &Tauc);
 
   decllocal("Pa", TDim::NHRU, "Average surface pressure", "(kPa)", &Pa);
 
-  declvar("ra", TDim::NHRU, "", "(s/m)", &ra);
+  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, " direct_rain + drip", "(mm)", &cum_net_rain);
+  declstatdiag("cum_net_rain", TDim::NHRU, "cumulative direct_rain + drip", "(mm)", &cum_net_rain);
 
   declvar("Subl_Cpy", TDim::NHRU, "canopy snow sublimation", "(mm/int)", &Subl_Cpy);
 
-  declstatdiag("cum_Subl_Cpy", TDim::NHRU, "canopy snow sublimation", "(mm)", &cum_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' Thru", "(mm/int)", &direct_snow);
+  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);
 
   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 Canopy unload ", "(mm)", &cum_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);
 
@@ -143,11 +143,11 @@ void ClassCRHMCanopy::decl(void) {
 
   declvar("intcp_evap", TDim::NHRU, "HRU Evaporation from interception", "(mm/int)", &intcp_evap);
 
-  declstatdiag("cum_intcp_evap", TDim::NHRU, "HRU Evaporation from interception", "(mm)", &cum_intcp_evap);
+  declstatdiag("cum_intcp_evap", TDim::NHRU, "Cumulative HRU Evaporation from interception", "(mm)", &cum_intcp_evap);
 
-  declvar("Qsisn_Var", TDim::NHRU, "incident short-wave at surface", "(W/m^2*int)", &Qsisn_Var);
+  declvar("Qsisn_Var", TDim::NHRU, "incident short-wave at surface", "(W/m^2)", &Qsisn_Var);
 
-  declvar("Qlisn_Var", TDim::NHRU, "incident long-wave at surface", "(W/m^2*int)", &Qlisn_Var);
+  declvar("Qlisn_Var", TDim::NHRU, "incident long-wave at surface", "(W/m^2)", &Qlisn_Var);
 
 
 // parameters:
@@ -179,9 +179,9 @@ void ClassCRHMCanopy::decl(void) {
 
   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);
 
-  decldiagparam("Alpha_c", TDim::NHRU, "[0.1]", "0.05", "0.2", "canopy albedo", "()", &Alpha_c);
+  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. Suggestions are Colorado - 0.023 and Alberta - 0.038", "()", &B_canopy);
+  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 ClassCRHMCanopy::init(void) {
@@ -387,7 +387,10 @@ void ClassCRHMCanopy::run(void) {
       C1 = 1.0/(D*SvDens*Nu);
 
       Alpha = 5.0;
-      Mpm = 4.0/3.0 * M_PI * PBSM_constants::DICE * Radius*Radius*Radius *(1.0 + 3.0/Alpha + 2.0/sqr(Alpha));
+
+      // 18Mar2022: remove Gamma Distribution Correction term, *(1.0 + 3.0/Alpha + 2.0/sqr(Alpha));
+      Mpm = 4.0 / 3.0 * M_PI * PBSM_constants::DICE * Radius * Radius * Radius;
+
 
 // sublimation rate of single 'ideal' ice sphere:
 
@@ -422,8 +425,11 @@ void ClassCRHMCanopy::run(void) {
 
           double IceBulbT = hru_t[hh] - (Vi* Hs/1e6/ci);
           double Six_Hour_Divisor = Global::Freq/4.0; // used to unload over 6 hours
+
+          // weekly dimensionless unloading coefficient -> to CRHM time interval
+          const double U = -1 * log(0.678) / (24 * 7 * Global::Freq / 24); 
+          // 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)
 
-          const double c = 0.678/(24*7*24/Global::Freq); // weekly dimensionless unloading coefficient -> to CRHM time interval
 
   // determine whether canopy snow is unloaded:
 
@@ -438,14 +444,19 @@ void ClassCRHMCanopy::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]*c; // the dimensionless unloading coefficient already /interval
-            if(SUnload[hh] > Snow_load[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];
+            {
+                Snow_load[hh] -= SUnload[hh];
+            }
 
             cum_SUnload[hh] += SUnload[hh];
           }

crhmcode/src/modules/ClassCRHMCanopyClearing.cpp

@@ -85,47 +85,47 @@ void ClassCRHMCanopyClearing::decl(void) {
 
   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*int)", &Qnsn_Var);
+  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*int)", &Qsisn_Var);
+  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 short-wave at surface", "(W/m^2*int)", &Qlisn_Var);
+  declvar("Qlisn_Var", TDim::NHRU, "incident short-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", "(W/m^2)", &Tauc);
+  decldiag("Tauc", TDim::NHRU, "short-wave transmissivity", "()", &Tauc);
 
   decllocal("Pa", TDim::NHRU, "Average surface pressure", "(kPa)", &Pa);
 
-  declvar("ra", TDim::NHRU, "", "(s/m)", &ra);
+  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, " direct_rain + drip", "(mm)", &cum_net_rain);
+  declstatdiag("cum_net_rain", TDim::NHRU, "cumulative direct_rain + drip", "(mm)", &cum_net_rain);
 
   declvar("Subl_Cpy", TDim::NHRU, "canopy snow sublimation", "(mm/int)", &Subl_Cpy);
 
-  declstatdiag("cum_Subl_Cpy", TDim::NHRU, "canopy snow sublimation", "(mm)", &cum_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' Thru", "(mm/int)", &direct_snow);
+  declvar("direct_snow", TDim::NHRU, "snow 'direct' through canopy", "(mm/int)", &direct_snow);
 
   declvar("SUnload", TDim::NHRU, "unloaded canopy snow", "(mm)", &SUnload);
 
@@ -137,7 +137,7 @@ void ClassCRHMCanopyClearing::decl(void) {
 
   declvar("net_snow", TDim::NHRU, "hru_snow minus interception", "(mm/int)", &net_snow);
 
-  declstatdiag("cum_net_snow", TDim::NHRU, "Cumulative Canopy unload ", "(mm)", &cum_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);
 
@@ -147,7 +147,7 @@ void ClassCRHMCanopyClearing::decl(void) {
 
   declvar("intcp_evap", TDim::NHRU, "HRU Evaporation from interception", "(mm/int)", &intcp_evap);
 
-  declstatdiag("cum_intcp_evap", TDim::NHRU, "HRU Evaporation from interception", "(mm)", &cum_intcp_evap);
+  declstatdiag("cum_intcp_evap", TDim::NHRU, "Cumulative HRU Evaporation from interception", "(mm)", &cum_intcp_evap);
 
 
 // parameters:
@@ -181,9 +181,9 @@ void ClassCRHMCanopyClearing::decl(void) {
 
   declparam("CanopyClearing", TDim::NHRU, "[0]", "0", "1", "canopy - 0/clearing - 1", "()", &CanopyClearing);
 
-  decldiagparam("Alpha_c", TDim::NHRU, "[0.1]", "0.05", "0.2", "canopy albedo", "()", &Alpha_c);
+  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. Suggestions are Colorado - 0.023 and Alberta - 0.038", "()", &B_canopy);
+  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 ClassCRHMCanopyClearing::init(void) {
@@ -405,7 +405,7 @@ void ClassCRHMCanopyClearing::run(void) {
           C1 = 1.0/(D*SvDens*Nu);
 
           Alpha = 5.0;
-          Mpm = 4.0/3.0 * M_PI * PBSM_constants::DICE * Radius*Radius*Radius *(1.0 + 3.0/Alpha + 2.0/sqr(Alpha));
+          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:
 
@@ -441,7 +441,8 @@ void ClassCRHMCanopyClearing::run(void) {
           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 c = 0.678/(24*7*24/Global::Freq); // weekly dimensionless unloading coefficient -> to CRHM time interval
+          const float 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:
 
@@ -457,7 +458,7 @@ void ClassCRHMCanopyClearing::run(void) {
             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]*c; // the dimensionless unloading coefficient already /interval
+            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;