Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions borland_source/CRHMmain.cpp
Original file line number Diff line number Diff line change
@@ -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"
Expand Down
711 changes: 358 additions & 353 deletions borland_source/NewModules.cpp

Large diffs are not rendered by default.

3 changes: 2 additions & 1 deletion borland_source/NewModules.h
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// 04/12/21 with changes to 07/02/20
// 08/11/21 with changes to 04/12/21
//---------------------------------------------------------------------------

#ifndef OurModulesH
Expand Down Expand Up @@ -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;
Expand Down
614 changes: 613 additions & 1 deletion borland_source/Notes.txt

Large diffs are not rendered by default.

18 changes: 9 additions & 9 deletions crhmcode/src/modules/ClassAyers.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand All @@ -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);

Expand Down
51 changes: 31 additions & 20 deletions crhmcode/src/modules/ClassCRHMCanopy.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand All @@ -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);

Expand All @@ -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:
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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:

Expand Down Expand Up @@ -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:

Expand All @@ -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];
}
Expand Down
31 changes: 16 additions & 15 deletions crhmcode/src/modules/ClassCRHMCanopyClearing.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand All @@ -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);

Expand All @@ -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:
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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:

Expand Down Expand Up @@ -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:

Expand All @@ -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;
Expand Down
Loading