Commit e5dda448 authored by Beisman, James Joseph's avatar Beisman, James Joseph
Browse files

organized functions from CanopyHydrology into a class and converted functions...

organized functions from CanopyHydrology into a class and converted functions to new indexing; added to indexing documentation
parent 5723a30d
......@@ -38,9 +38,9 @@ will be shifted by (nlevsno - 1), similar to the surface data.
ABOVE SURFACE DATA AND ABOVE AND BELOW SURFACE DATA:
snow or canopy bottom layer (adjacent to ground surface) - CLM: var(0), THIS PROJECT: var[nlevsno-1]
surface/top subsurface layer - CLM: var(1), THIS PROJECT: var[nlevsno]
bottom subsurface layer - CLM: var(nlevgrnd), THIS PROJECT: var[nlevgrnd+nlevsno-1]
top active snow layer - CLM: var(snl+1), THIS PROJECT: var[nlevsno-snl] -- snl is negative in CLM, positive here
second from top active snow layer - CLM: var(snl+2), THIS PROJECT: var[nlevsno-snl+1]
bottom subsurface layer - CLM: var(nlevgrnd), THIS PROJECT: var[nlevgrnd+nlevsno-1]
CONCURRENTLY WORKING WITH SUBSURFACE ONLY AND ABOVE SURFACE DATA:
subsurface data's index is offset from above surface data by nlevsno
......@@ -105,3 +105,23 @@ CLM variable indices new indices
.
-- ---
|nlevgrnd| bottom layer |nlevgrnd - 1|
Special zi mapping
CLM new indices
zi zi
grid grid
-- -5 0 ---
|-4| |0|
-- -4 1 ---
|-3| |1|
-- -3 2 ---
|-2| |2|
-- -2 3 ---
|-1| |3|
-- -1 4 --- [nlevsno - 1]
| 0| |4|
----- 0 ground interface 5 ----- [nlevsno]
| 1| |5|
-- 1 6 --- [nlevsno + 1]
\ No newline at end of file
add_library (physics physics_dummy.cc fracwet.cc canopy_interception.cc canopy_ground_flux.cc canopy_irrigation.cc coldstart_topography.cc coldstart_snow.cc snow_init.cc frac_h2osfc.cc SurfaceRadiationMod_functions.cc SurfaceResistanceMod_functions.cc QSatMod_functions.cc CanopyTemperatureMod_functions.cc FrictionVelocityMod_functions.cc)
add_library (physics physics_dummy.cc CanopyHydrology.cc coldstart_topography.cc coldstart_snow.cc SurfaceRadiationMod_functions.cc SurfaceResistanceMod_functions.cc QSatMod_functions.cc CanopyTemperatureMod_functions.cc FrictionVelocityMod_functions.cc)
target_include_directories (physics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
This diff is collapsed.
/*
DESCRIPTION:
Add liquid and solid water inputs to ground surface after interception and
canopy storage losses.
INPUTS:
ctype [int] urban column type
do_capsnow [bool] true => do snow capping
frac_veg_nosno [int] fraction of veg not covered by snow (0/1 now) [-]
forc_rain [double] rain rate (kg H2O/m**2/s, or mm liquid H2O/s)
forc_snow [double] snow rate (kg H2O/m**2/s, or mm liquid H2O/s)
qflx_candrip [double] rate of canopy runoff and snow falling off canopy [mm/s]
qflx_through_snow [double] direct rain throughfall [mm/s]
qflx_through_rain [double] direct snow throughfall [mm/s]
fracsnow [double] frac of precipitation that is snow [-]
fracrain [double] frac of precipitation that is rain [-]
qflx_irrig [double] irrigation amount (mm/s)
OUTPUTS:
qflx_prec_grnd [double] water onto ground including canopy runoff [kg/(m2 s)]
qflx_snwcp_liq [double] excess rainfall due to snow capping (mm H2O /s) [+]
qflx_snwcp_ice [double] excess snowfall due to snow capping (mm H2O /s) [+]
qflx_snow_grnd [double] snow on ground after interception (mm H2O/s) [+]
qflx_rain_grnd [double] rain on ground after interception (mm H2O/s) [+]
*/
#include <algorithm>
#include <cmath>
#include "clm_constants.hh"
void CanopyGroundFlux(const int& ctype,
const bool& do_capsnow,
const int& frac_veg_nosno,
const double& forc_rain,
const double& forc_snow,
const double& qflx_candrip,
const double& qflx_through_snow,
const double& qflx_through_rain,
const double& fracsnow,
const double& fracrain,
const double& qflx_irrig,
double& qflx_prec_grnd,
double& qflx_snwcp_liq,
double& qflx_snwcp_ice,
double& qflx_snow_grnd,
double& qflx_rain_grnd)
{
double qflx_prec_grnd_snow, qflx_prec_grnd_rain;
double qflx_dirct_rain, qflx_leafdrip; // these are not currently needed, but should be passed out if SedYield is desired
// Precipitation onto ground (kg/(m2 s))
if ((ctype != icol_sunwall) && (ctype != icol_shadewall)) {
if (frac_veg_nosno == 0) {
qflx_prec_grnd_snow = forc_snow;
qflx_prec_grnd_rain = forc_rain;
qflx_dirct_rain = forc_rain;
qflx_leafdrip = 0.0;
} else {
qflx_prec_grnd_snow = qflx_through_snow + (qflx_candrip * fracsnow);
qflx_prec_grnd_rain = qflx_through_rain + (qflx_candrip * fracrain);
qflx_dirct_rain = qflx_through_rain;
qflx_leafdrip = qflx_candrip * fracrain;
}
} else {
// Urban sunwall and shadewall have no intercepted precipitation
qflx_prec_grnd_snow = 0.0;
qflx_prec_grnd_rain = 0.0;
qflx_dirct_rain = 0.0;
qflx_leafdrip = 0.0;
}
// Add irrigation water directly onto ground (bypassing canopy interception)
qflx_prec_grnd_rain = qflx_prec_grnd_rain + qflx_irrig;
// Total water onto ground
qflx_prec_grnd = qflx_prec_grnd_snow + qflx_prec_grnd_rain;
if (do_capsnow) {
qflx_snwcp_liq = qflx_prec_grnd_rain;
qflx_snwcp_ice = qflx_prec_grnd_snow;
qflx_snow_grnd = 0.0;
qflx_rain_grnd = 0.0;
} else {
qflx_snwcp_liq = 0.0;
qflx_snwcp_ice = 0.0;
qflx_snow_grnd = qflx_prec_grnd_snow; // ice onto ground (mm/s)
qflx_rain_grnd = qflx_prec_grnd_rain; // liquid water onto ground (mm/s)
}
}
/*
DESCRIPTION:
Calculate interception and partition incoming precipitation
into canopy storage, canopy runoff, rain and snow throughfall.
Also calculates fraction of precipitation that is rain/snow.
INPUTS:
ltype [int] landunit type
urbpoi [bool] true => landunit is an urban point
ctype [int] urban column type
frac_veg_nosno [int] fraction of veg not covered by snow (0/1 now) [-]
forc_rain [double] rain rate (kg H2O/m**2/s, or mm liquid H2O/s)
forc_snow [double] snow rate (kg H2O/m**2/s, or mm liquid H2O/s)
dewmx [double] Maximum allowed dew [mm]
elai [double] one-sided leaf area index with burying by snow
esai [double] one-sided stem area index with burying by snow
dtime [double] time step length (sec)
OUTPUTS:
h2ocan [double] total canopy water (mm H2O)
qflx_candrip [double] rate of canopy runoff and snow falling off canopy [mm/s]
qflx_through_snow [double] direct rain throughfall [mm/s]
qflx_through_rain [double] direct snow throughfall [mm/s]
qflx_prec_intr [double] interception of precipitation [mm/s]
fracsnow [double] frac of precipitation that is snow [-]
fracrain [double] frac of precipitation that is rain [-]
*/
#include <algorithm>
#include <cmath>
#include "clm_constants.hh"
void CanopyInterception(const int& ltype,
const bool& urbpoi,
const int& ctype,
const int& frac_veg_nosno,
const double& forc_rain,
const double& forc_snow,
const double& dewmx,
const double& elai,
const double& esai,
const double& dtime,
double& h2ocan,
double& qflx_candrip,
double& qflx_through_snow,
double& qflx_through_rain,
double& qflx_prec_intr,
double& fracsnow,
double& fracrain)
{
//Canopy interception/storage and throughfall
//Add precipitation to leaf water
if (ltype == istsoil || ltype == istwet || urbpoi || ltype == istcrop) {
qflx_candrip = 0.0; // rate of canopy runoff
qflx_through_snow = 0.0; // rain precipitation direct through canopy
qflx_through_rain = 0.0; // snow precipitation direct through canopy
qflx_prec_intr = 0.0; // total intercepted precipitation
fracsnow = 0.0; // fraction of input precip that is snow
fracrain = 0.0; // fraction of input precip that is rain
if (ctype != icol_sunwall && ctype != icol_shadewall) {
if (frac_veg_nosno == 1 && (forc_rain + forc_snow) > 0.0) {
// determine fraction of input precipitation that is snow and rain
fracsnow = forc_snow/(forc_snow + forc_rain);
fracrain = forc_rain/(forc_snow + forc_rain);
// The leaf water capacities for solid and liquid are different,
// generally double for snow, but these are of somewhat less
// significance for the water budget because of lower evap. rate at
// lower temperature. Hence, it is reasonable to assume that
// vegetation storage of solid water is the same as liquid water.
double h2ocanmx = dewmx * (elai + esai);
// Coefficient of interception
// set fraction of potential interception to max 0.25
double fpi = 0.25 * (1.0 - exp(-0.5*(elai + esai)));
// Direct throughfall
qflx_through_snow = forc_snow * (1.0-fpi);
qflx_through_rain = forc_rain * (1.0-fpi);
// Intercepted precipitation [mm/s]
qflx_prec_intr = (forc_snow + forc_rain) * fpi;
// Water storage of intercepted precipitation and dew
h2ocan = std::max(0.0, (h2ocan + dtime*qflx_prec_intr));
// Initialize rate of canopy runoff and snow falling off canopy
qflx_candrip = 0.0;
// Excess water that exceeds the leaf capacity
double xrun = (h2ocan - h2ocanmx)/dtime;
// Test on maximum dew on leaf
// Note if xrun > 0 then h2ocan must be at least h2ocanmx
if (xrun > 0.0) {
qflx_candrip = xrun;
h2ocan = h2ocanmx;
}
}
}
} else if (ltype == istice||ltype == istice_mec) {
h2ocan = 0.0;
qflx_candrip = 0.0;
qflx_through_snow = 0.0;
qflx_through_rain = 0.0;
qflx_prec_intr = 0.0;
fracsnow = 0.0;
fracrain = 0.0;
}
}
/*
DESCRIPTION:
Determine whether we're irrigating here; set qflx_irrig appropriately.
INPUTS:
irrig_rate [double] current irrigation rate (applied if n_irrig_steps_left > 0) [mm/s]
INOUT:
n_irrig_steps_left [int] number of time steps for which we still need to irrigate today
OUTPUTS:
qflx_irrig [double] irrigation amount (mm/s)
*/
void CanopyIrrigation(const double& irrig_rate,
int& n_irrig_steps_left,
double& qflx_irrig)
{
if (n_irrig_steps_left > 0) {
qflx_irrig = irrig_rate;
n_irrig_steps_left -= 1;
} else {
qflx_irrig = 0.0;
}
}
/*
DESCRIPTION:
Determine fraction of land surfaces which are submerged
based on surface microtopography and surface water storage.
!!h2osoi_liq - indexing in CLM, (-4,..,15) (top of snow,..,bottom of bedrock) -- for now, I decouple
the snow layers CLM(-4,..,0), this code (5,..,0) from the subsurface CLM(1,..15) this code (0,..14)
INPUTS:
dtime [double] land model time step (sec)
ltype [int] landunit type
micro_sigma [double] microtopography pdf sigma (m)
h2osno [double] snow water (mm H2O)
no_update [bool] flag to make calculation w/o updating variables
OUTPUTS:
h2osfc [double] surface water (mm)
h2osoi_liq [double] liquid water (kg/m2) -- just the top subsurface layer for this function
frac_sno [double] fraction of ground covered by snow (0 to 1)
frac_sno_eff [double] effective fraction of ground covered by snow (0 to 1)
qflx_h2osfc2topsoi [double] liquid water coming from surface standing water top soil (mm H2O/s)
frac_h2osfc [double] fractional area with surface water greater than zero (0 to 1)
*/
#include <algorithm>
#include <cmath>
#include "clm_constants.hh"
void FracH2OSfc(
const double& dtime,
const int& ltype,
const double& micro_sigma,
const double& h2osno,
const bool& no_update,
double& h2osfc,
double& h2osoi_liq,
double& frac_sno,
double& frac_sno_eff,
double& qflx_h2osfc2topsoi,
double& frac_h2osfc)
{
double d,fd,dfdd,sigma;
// arbitrary lower limit on h2osfc for safer numerics...
double min_h2osfc = 1.e-8;
qflx_h2osfc2topsoi = 0.0;
// h2osfc only calculated for soil vegetated land units
if (ltype == istsoil || ltype == istcrop) {
// Use newton-raphson method to iteratively determine frac_h20sfc
// based on amount of surface water storage (h2osfc) and
// microtopography variability (micro_sigma)
if (h2osfc > min_h2osfc) {
// a cutoff is needed for numerical reasons...(nonconvergence after 5 iterations)
d = 0.0;
sigma = 1.0e3 * micro_sigma; // convert to mm
for (int l = 0; l < 10; l++) {
fd = 0.5 * d * (1.0 + erf(d / (sigma * sqrt(2.0)))) + sigma / sqrt(2.0 * SHR_CONST_PI) * exp(-pow(d, 2) / (2.0 * pow(sigma, 2))) - h2osfc;
dfdd = 0.5 * (1.0 + erf(d / (sigma * sqrt(2.0))));
d = d - fd/dfdd;
}
// update the submerged areal fraction using the new d value
frac_h2osfc = 0.5 * (1.0 + erf(d / (sigma * sqrt(2.0))));
} else {
frac_h2osfc = 0.0;
h2osoi_liq = h2osoi_liq + h2osfc;
qflx_h2osfc2topsoi = h2osfc / dtime;
h2osfc = 0.0;
}
if (!no_update) {
// adjust fh2o, fsno when sum is greater than zero
if (frac_sno > (1.0 - frac_h2osfc) && h2osno > 0.0) {
if (frac_h2osfc > 0.01) {
frac_h2osfc = std::max((1.0 - frac_sno), 0.01);
frac_sno = 1.0 - frac_h2osfc;
} else {
frac_sno = 1.0 - frac_h2osfc;
}
frac_sno_eff = frac_sno;
}
}
} else {
frac_h2osfc = 0.0;
}
}
/*
DESCRIPTION:
Determine fraction of vegetated surfaces which are wet and
fraction of elai which is dry. The variable ``fwet'' is the
fraction of all vegetation surfaces which are wet including
stem area which contribute to evaporation. The variable ``fdry''
is the fraction of elai which is dry because only leaves
can transpire. Adjusted for stem area which does not transpire.
INPUTS:
frac_veg_nosno [int] fraction of veg not covered by snow (0/1 now) [-]
dewmx [double] Maximum allowed dew [mm]
elai [double] one-sided leaf area index with burying by snow
esai [double] one-sided stem area index with burying by snow
h2ocan [double] total canopy water (mm H2O)
OUTPUTS:
fwet [double] fraction of canopy that is wet (0 to 1)
fdry [double] fraction of foliage that is green and dry [-] (new)
*/
#include <algorithm>
#include <cmath>
void FracWet(const int& frac_veg_nosno,
const double& dewmx,
const double& elai,
const double& esai,
const double& h2ocan,
double& fwet,
double& fdry)
{
if (frac_veg_nosno == 1) {
if (h2ocan > 0.0) {
double vegt = frac_veg_nosno * (elai + esai);
double dewmxi = 1.0 / dewmx;
fwet = pow(((dewmxi/vegt)*h2ocan), 2.0/3.0);
fwet = std::min(fwet, 1.0); //Check for maximum limit of fwet
} else {
fwet = 0.0;
}
fdry = (1.0 - fwet) * elai / (elai + esai);
} else {
fwet = 0.0;
fdry = 0.0;
}
}
#ifndef LANDTYPE_H_
#define LANDTYPE_H_
struct LandType {
int ltype;
int ctype;
bool urbpoi;
bool lakpoi;
};
#endif
/*
DESCRIPTION:
Initialization snow layer(s) if the snow accumulation exceeds 10 mm.
INPUTS:
OUTPUTS:
*/
#include <algorithm>
#include <cmath>
#include "clm_constants.hh"
void SnowInit(
const double& dtime,
const int& ltype,
const bool& urbpoi,
const bool& do_capsnow,
const int& oldfflag,
const double& forc_t,
const double& t_grnd,
const double& qflx_snow_grnd,
const double& qflx_snow_melt,
const double& n_melt,
double& snow_depth,
double& h2osno,
double& int_snow,
double* swe_old,
double* h2osoi_liq,
double* h2osoi_ice,
double* t_soisno,
double* frac_iceold,
int& snl,
double* dz,
double* z,
double* zi,
double* snw_rds,
int& newnode,
double& qflx_snow_h2osfc,
double& frac_sno_eff,
double& frac_sno)
{
double rpi = SHR_CONST_PI;
double temp_snow_depth, dz_snowf, newsnow, bifall, snowmelt, accum_factor, temp_intsnow, z_avg;
// Determine snow height and snow water
// Use Alta relationship, Anderson(1976); LaChapelle(1961),
// U.S.Department of Agriculture Forest Service, Project F,
// Progress Rep. 1, Alta Avalanche Study Center:Snow Layer Densification.
qflx_snow_h2osfc = 0.0;
// set temporary variables prior to updating
temp_snow_depth = snow_depth;
// save initial snow content
for(int j = nlevsno - 1; j >= snl; j--) {
swe_old[j] = 0.0;
}
for (int j = 0; j < snl; j++) {
swe_old[j] = h2osoi_liq[j] + h2osoi_ice[j];
}
if (do_capsnow) {
dz_snowf = 0.0;
newsnow = qflx_snow_grnd * dtime;
frac_sno = 1.0;
int_snow = 5.e2;
} else {
if (forc_t > tfrz + 2.0) {
bifall = 50.0 + 1.7 * pow(17.0, 1.5);
} else if (forc_t > tfrz - 15.0) {
bifall = 50.0 + 1.7 * pow((forc_t - tfrz + 15.0), 1.5);
} else {
bifall = 50.0;
}
// all snow falls on ground, no snow on h2osfc
newsnow = qflx_snow_grnd * dtime;
// update int_snow
int_snow = std::max(int_snow, h2osno); //h2osno could be larger due to frost
// snowmelt from previous time step * dtime
snowmelt = qflx_snow_melt * dtime;
// set shape factor for accumulation of snow
accum_factor = 0.1;
/*====================== FSCA PARAMETERIZATIONS ======================
fsca parameterization based on *changes* in swe
first compute change from melt during previous time step */
if (h2osno > 0.0) {
if(snowmelt > 0.0) {
double smr = std::min(1.0, (h2osno / int_snow));
frac_sno = 1.0 - pow((acos(std::min(1.0, (2.0 * smr - 1.0))) / rpi), n_melt);
}
// update fsca by new snow event, add to previous fsca
if (newsnow > 0.0) {
double fsno_new = 1.0 - (1.0 - tanh(accum_factor * newsnow)) * (1.0 - frac_sno);
frac_sno = fsno_new;
// reset int_snow after accumulation events
temp_intsnow = (h2osno + newsnow) / (0.5 * (cos(rpi * pow( (1.0 - std::max(frac_sno, 1.e-6)), (1.0 / n_melt))) + 1.0));
int_snow = std::min(1.e8, temp_intsnow);
}
/*====================================================================*/
// for subgrid fluxes
if (subgridflag == 1 && !urbpoi) {
if (frac_sno > 0.0) {
snow_depth = snow_depth + newsnow / (bifall * frac_sno);
} else {
snow_depth = 0.0;
}
} else {
// for uniform snow cover
snow_depth = snow_depth + newsnow / bifall;
}
if (oldfflag == 1) {
// snow cover fraction in Niu et al. 2007
if (snow_depth > 0.0) {
frac_sno = tanh(snow_depth / (2.5 * zlnd * pow((std::min(800.0, ((h2osno+ newsnow) / snow_depth / 100.0))), 1.0))); // why to the power of 1.0??
}
if (h2osno < 1.0) {
frac_sno = std::min(frac_sno, h2osno);
}
}
} else { //h2osno == 0
// initialize frac_sno and snow_depth when no snow present initially
if (newsnow > 0.0) {
z_avg = newsnow / bifall;
frac_sno = tanh(accum_factor * newsnow);
// make int_snow consistent w/ new fsno, h2osno
int_snow = 0.0; //reset prior to adding newsnow below
temp_intsnow = (h2osno + newsnow) / (0.5 * (cos(rpi * pow( (1.0 - std::max(frac_sno, 1.e-6)), (1.0 / n_melt))) + 1.0));
int_snow = std::min(1.e8, temp_intsnow);
// update snow_depth and h2osno to be consistent with frac_sno, z_avg
if (subgridflag == 1 && !urbpoi) {
snow_depth = z_avg / frac_sno;
} else {
snow_depth = newsnow / bifall;
}
// use n&y07 formulation
if (oldfflag == 1) {
// snow cover fraction in Niu et al. 2007
if(snow_depth > 0.0) {
frac_sno = tanh(snow_depth / (2.5 * zlnd * pow((std::min(800.0, ((h2osno+ newsnow) / snow_depth / 100.0))), 1.0))); // why to the power of 1.0??
}
}
} else {
z_avg = 0.0;
snow_depth = 0.0;
frac_sno = 0.0;
}
}
// no snow on surface water
qflx_snow_h2osfc = 0.0;
// update h2osno for new snow
h2osno = h2osno + newsnow;
int_snow = int_snow + newsnow;
// update change in snow depth
dz_snowf = (snow_depth - temp_snow_depth) / dtime;
} // end else do_capsnow
// set frac_sno_eff variable
if (ltype == istsoil || ltype == istcrop) {
if (subgridflag == 1) {
frac_sno_eff = frac_sno;
} else {
frac_sno_eff = 1.0;
}
} else {
frac_sno_eff = 1.0;
}
if (ltype == istwet && t_grnd > tfrz) {