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

added SurfaceResistanceMod, QSatMod; started CanopyTemperatureMod; fixed minor...

added SurfaceResistanceMod, QSatMod; started CanopyTemperatureMod; fixed minor error in SurfaceRadiation
parent 6c948c38
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)
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)
target_include_directories (physics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
/* physics functions from CanopyTemperatureMod.F90
*/
#include <algorithm>
#include <cmath>
#include "clm_constants.hh"
/* SavePriorGroundTemp()
record t_h2osfc and t_soisno prior to updating
*/
void SavePriorGroundTemp(
const int& ctype,
const bool& lakpoi,
const int& nlevsno,
const int& nlevgrnd,
const double* t_soisno_,
double* tssbef_,
double* t_h2osfc_bef)
{
// reference variables with pointer offset to allow use of CLM/ELM index without modification
const double* tssbef = &tssbef_[nlevsno-1];
double* t_soisno = &t_soisno_[nlevsno-1];
if (!lakpoi) {
for (int i = -nlevsno+1; i <= nlevgrnd; i++) {
if ((ctype == icol_sunwall || ctype == icol_shadewall || ctype == icol_roof) && i > nlevurb) {
tssbef[i] = spval;
} else {
tssbef[i] = t_soisno[i];
}
// record t_h2osfc prior to updating
t_h2osfc_bef = t_h2osfc;
}
}
}
/* CalculateGroundTemp()
*/
void CalculateGroundTemp(
const bool& lakpoi,
const int& nlevsno,
const int& snl,
const double& frac_sno_eff,
const double& frac_h2osfc,
const double& t_h2osfc,
const double* t_soisno_,
double& t_grnd)
{
// reference variables with pointer offset to allow use of CLM/ELM index without modification
const double* t_soisno = &t_soisno_[nlevsno-1];
if (!lakpoi) {
if (snl < 0) {
t_grnd = frac_sno_eff * t_soisno[snl+1] + (1.0 - frac_sno_eff - frac_h2osfc) * t_soisno[1] + frac_h2osfc * t_h2osfc;
} else {
t_grnd = (1.0 - frac_h2osfc) * t_soisno[1] + frac_h2osfc * t_h2osfc;
}
}
}
/*
*/
void CalculateSoilAlpha(
)
{
if (!lakpoi) {
if (ctype == icol_road_perv) {
hr_road_perv = 0.0
}
double qred = 1.0;
if (ltype != istwet && ltype != istice && ltype != istice_mec) {
if (ltype == istsoil || ltype == istcrop) {
wx = (h2osoi_liq[1] / denh2o + h2osoi_ice[1] / denice) / dz[1];
fac = std::min(1.0, wx / watsat[1]);
fac = std::max(fac, 0.01);
psit = -sucsat[1] * pow(fac, (-bsw[1]));
psit = std::max(smpmin, psit);
// modify qred to account for h2osfc
hr = exp(psit / roverg / t_soisno[1]);
qred = (1.0 - frac_sno - frac_h2osfc) * hr + frac_sno + frac_h2osfc;
soilalpha = qred;
} else if (ctype == icol_road_perv) {
// Pervious road depends on water in total soil column
nlevbed = nlev2bed;
for (int j = 1; j <= nlevbed; j++) {
if (t_soisno[j] >= tfrz) {
vol_ice = std::min(watsat[j], h2osoi_ice[j] / (dz[j] * denice));
eff_porosity = watsat[j] - vol_ice;
vol_liq = std::min(eff_porosity, h2osoi_liq[j] / (dz[j] * denh2o));
fac = std::min(std::max(vol_liq - watdry[j], 0.0) / (watopt[j] - watdry[j]), 1.0);
} else {
fac = 0.0;
}
rootr_road_perv[j] = rootfr_road_perv[j] * fac;
hr_road_perv = hr_road_perv + rootr_road_perv[j];
}
// Allows for sublimation of snow or dew on snow
qred = (1.0 - frac_sno) * hr_road_perv + frac_sno;
// Normalize root resistances to get layer contribution to total ET
if (hr_road_perv > 0.0) {
for (int j = 1; j <= nlevsoi; j++) {
rootr_road_perv[j] = rootr_road_perv[j] / hr_road_perv;
}
}
soilalpha_u = qred;
} else if (ctype == icol_sunwall || ctype == icol_shadewall) {
qred = 0.0;
soilalpha_u = spval;
} else if (ctype == icol_roof || ctype == icol_road_imperv) {
qred = 1.0;
soilalpha_u = spval;
}
} else {
soilalpha = spval;
}
}
}
/*
*/
void CalculateHumidities(
)
{
if (!lakpoi) {
if (ltype == istsoil || ltype == istcrop) {
QSat(t_soisno[snl+1], forc_pbot, eg, degdT, qsatg, qsatgdT);
if (qsatg > forc_q && forc_q > qsatg) {
qsatg = forc_q;
qsatgdT = 0.0;
}
qg_snow = qsatg;
dqgdT = frac_sno * qsatgdT;
QSat(t_soisno[1] , forc_pbot, eg, degdT, qsatg, qsatgdT);
if (qsatg > forc_q && forc_q > hr * qsatg) {
qsatg = forc_q;
qsatgdT = 0.0;
}
qg_soil = hr * qsatg;
dqgdT = dqgdT + (1.0 - frac_sno - frac_h2osfc) * hr * qsatgdT;
// to be consistent with hs_top values in SoilTemp, set qg_snow to qg_soil for snl = 0 case
// this ensures hs_top_snow will equal hs_top_soil
if (snl >= 0) {
qg_snow; = qg_soil;
dqgdT = (1.0 - frac_h2osfc) * hr * dqgdT;
}
QSat(t_h2osfc, forc_pbot, eg, degdT, qsatg, qsatgdT);
if (qsatg > forc_q && forc_q > qsatg) {
qsatg = forc_q;
qsatgdT = 0.0;
}
qg_h2osfc = qsatg;
dqgdT = dqgdT + frac_h2osfc * qsatgdT;
qg = frac_sno_eff * qg_snow + (1.0 - frac_sno_eff - frac_h2osfc) * qg_soil + frac_h2osfc * qg_h2osfc;
} else {
QSat(t_grnd, forc_pbot, eg, degdT, qsatg, qsatgdT);
qg = qred * qsatg;
dqgdT = qred * qsatgdT;
if (qsatg > forc_q && forc_q > qred * qsatg) {
qg = forc_q;
dqgdT = 0.0;
}
qg_snow = qg;
qg_soil = qg;
qg_h2osfc = qg;
}
}
}
/*
*/
#include "clm_constants.hh"
void QSat(
const double& T,
const double& p,
double& es,
double& esdT,
double& qs,
double& qsdT)
{
// For water vapor (temperature range 0C-100C)
const double a0 = 6.11213476;
const double a1 = 0.444007856;
const double a2 = 0.143064234e-01;
const double a3 = 0.264461437e-03;
const double a4 = 0.305903558e-05;
const double a5 = 0.196237241e-07;
const double a6 = 0.892344772e-10;
const double a7 = -0.373208410e-12;
const double a8 = 0.209339997e-15;
// For derivative:water vapor
const double b0 = 0.444017302;
const double b1 = 0.286064092e-01;
const double b2 = 0.794683137e-03;
const double b3 = 0.121211669e-04;
const double b4 = 0.103354611e-06;
const double b5 = 0.404125005e-09;
const double b6 = -0.788037859e-12;
const double b7 = -0.114596802e-13;
const double b8 = 0.381294516e-16;
// For ice (temperature range -75C-0C)
const double c0 = 6.11123516;
const double c1 = 0.503109514;
const double c2 = 0.188369801e-01;
const double c3 = 0.420547422e-03;
const double c4 = 0.614396778e-05;
const double c5 = 0.602780717e-07;
const double c6 = 0.387940929e-09;
const double c7 = 0.149436277e-11;
const double c8 = 0.262655803e-14;
// For derivative:ice
const double d0 = 0.503277922;
const double d1 = 0.377289173e-01;
const double d2 = 0.126801703e-02;
const double d3 = 0.249468427e-04;
const double d4 = 0.313703411e-06;
const double d5 = 0.257180651e-08;
const double d6 = 0.133268878e-10;
const double d7 = 0.394116744e-13;
const double d8 = 0.498070196e-16;
double td,vp,vp1,vp2,T_limit;
T_limit = T - SHR_CONST_TKFRZ;
if (T_limit > 100.0) { T_limit=100.0; }
if (T_limit < -75.0) { T_limit=-75.0; }
td = T_limit;
if (td >= 0.0) {
es = a0 + td *(a1 + td*(a2 + td*(a3 + td*(a4
+ td*(a5 + td*(a6 + td*(a7 + td*a8)))))));
esdT = b0 + td*(b1 + td*(b2 + td*(b3 + td*(b4
+ td*(b5 + td*(b6 + td*(b7 + td*b8)))))));
} else {
es = c0 + td*(c1 + td*(c2 + td*(c3 + td*(c4
+ td*(c5 + td*(c6 + td*(c7 + td*c8)))))));
esdT = d0 + td*(d1 + td*(d2 + td*(d3 + td*(d4
+ td*(d5 + td*(d6 + td*(d7 + td*d8)))))));
}
es = es * 100.0; // pa
esdT = esdT * 100.0; // pa/K
vp = 1.0 / (p - 0.378 * es);
vp1 = 0.622 * vp;
vp2 = vp1 * vp;
qs = es * vp1; // kg/kg
qsdT = esdT * vp2 * p; // 1 / K
}
......@@ -400,7 +400,7 @@ void CanopySunShadeFractions(
double& fsun)
{
if (ltype < isturb_MIN) {
int ipar = 1; // The band index for PAR
int ipar = 0; // The band index for PAR
for (int iv = 0; iv < nrad; iv++) {
parsun_z[iv] = 0.0;
parsha_z[iv] = 0.0;
......
// functions from SurfaceResistanceMod.F90
#include <algorithm>
#include <cmath>
#include "clm_constants.hh"
/* calc_soilevap_stress()
DESCRIPTION:
Compute the Lee-Pielke beta factor to scale actual soil evaporation from potential evaporation.
This is the calc_beta_leepielke1992() function from CLM/ELM that gets called from calc_soilevap_stress().
code that gets called when use_vsfm == true is not currently included
INPUTS:
ltype [int] landunit type
lakpoi [bool] true => landunit is an urban point
frac_sno [double] fraction of ground covered by snow (0 to 1)
frac_h2osfc [double] fraction of ground covered by surface water (0 to 1)
watsat[nlevgrnd] [double] volumetric soil water at saturation (porosity)
watfc[nlevgrnd] [double] volumetric soil water at field capacity
h2osoi_liq[nlevgrnd+nlevsno] [double] ice lens (kg/m2)
h2osoi_ice[nlevgrnd+nlevsno] [double] liquid water (kg/m2)
dz[nlevgrnd+nlevsno] [double] layer thickness (m)
OUTPUTS:
soilbeta [double] factor that reduces ground evaporation
*/
void calc_soilevap_stress(
const int& ltype,
const bool& lakpoi,
const double& frac_sno,
const double& frac_h2osfc,
const double* watsat_,
const double* watfc_,
const double* h2osoi_liq_,
const double* h2osoi_ice_,
const double* dz_,
double& soilbeta)
{
// reference variables with pointer offset to allow use of CLM/ELM index without modification
const double* watsat = &watsat_[-1];
const double* watfc = &watfc_[-1];
const double* h2osoi_liq = &h2osoi_liq_[nlevsno-1];
const double* h2osoi_ice = &h2osoi_ice_[nlevsno-1];
const double* dz = &dz_[nlevsno-1];
// local variables
double fac, fac_fc, wx;
const double denice = SHR_CONST_RHOICE;
const double denh2o = SHR_CONST_RHOFW;
if (!lakpoi) {
if (ltype != istwet && ltype != istice && ltype != istice_mec) {
if (ltype == istsoil || ltype == istcrop) {
wx = (h2osoi_liq[1] / denh2o + h2osoi_ice[1] / denice) / dz[1];
fac = std::min(1.0, wx / watsat[1]);
fac = std::max(fac, 0.01);
// Lee and Pielke 1992 beta, added by K.Sakaguchi
if (wx < watfc[1] ) { //!when water content of ths top layer is less than that at F.C.
fac_fc = std::min(1.0, wx / watfc[1]); // eqn5.66 but divided by theta at field capacity
fac_fc = std::max(fac_fc, 0.01);
// modify soil beta by snow cover. soilbeta for snow surface is one
soilbeta = (1.0 - frac_sno - frac_h2osfc) * 0.25 * pow(1.0 - cos(SHR_CONST_PI * fac_fc), 2.0)
+ frac_sno + frac_h2osfc;
} else {
soilbeta = 1.0;
}
} else if (ltype == icol_road_perv) {
soilbeta = 0.0;
} else if (ltype == icol_sunwall || ltype == icol_shadewall) {
soilbeta = 0.0;
} else if (ltype == icol_roof || ltype == icol_road_imperv) {
soilbeta = 0.0;
}
} else {
soilbeta = 1.0;
}
}
}
......@@ -38,6 +38,9 @@ static const double zsno = 0.0024; // Roughness length for snow [m]
//from shr_const_mod.F90
static const double SHR_CONST_PI = 3.14159265358979323846;
static const double SHR_CONST_TKFRZ = 273.15;
static const double SHR_CONST_RHOFW = 1.000e3; // density of fresh water ~ kg/m^3
static const double SHR_CONST_RHOSW = 1.026e3; // density of sea water ~ kg/m^3
static const double SHR_CONST_RHOICE = 0.917e3; // density of ice ~ kg/m^3
static const double tfrz = SHR_CONST_TKFRZ;
#endif
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment