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

added support for 2m reference variables

parent 4648c7d0
......@@ -12,6 +12,7 @@ It is not the same as the displa calculated in CanopyTemperature
#include "clm_constants.h"
#include "landtype.h"
#include "frictionvelocity.h"
#include "qsat.h"
class BareGroundFluxes {
......@@ -25,8 +26,20 @@ private:
double um; // wind speed including the stablity effect [m/s]
double temp1; // relation for potential temperature profile
double temp2; // relation for specific humidity profile
double temp12m; // relation for potential temperature profile applied at 2-m
double temp22m; // relation for specific humidity profile applied at 2-m
double ustar; // friction velocity [m/s]
// local vars
double frac_veg_nosno;
double forc_q;
double forc_th;
double forc_hgt_u_patch;
double thm;
double thv;
double z0mg;
public:
/* BareGroundFluxes::InitializeFlux
......@@ -34,18 +47,18 @@ DESCRIPTION:
initialize fluxes and call MoninObukIni()
INPUTS:
Land [LandType] struct containing information about landtype
frac_veg_nosno [int] fraction of vegetation not covered by snow (0 OR 1) [-]
forc_u [double] atmospheric wind speed in east direction (m/s)
forc_v [double] atmospheric wind speed in north direction (m/s)
forc_q [double] atmospheric specific humidity (kg/kg)
forc_th [double] atmospheric potential temperature (Kelvin)
forc_hgt_u_patch [double] observational height of wind at pft level [m]
thm [double] intermediate variable (forc_t+0.0098*forc_hgt_t_patch)
thv [double] virtual potential temperature (kelvin)
t_grnd [double] ground temperature (Kelvin)
qg [double] ground specific humidity [kg/kg]
z0mg [double] roughness length over ground, momentum [m]
Land [LandType] struct containing information about landtype
frac_veg_nosno_in [int] fraction of vegetation not covered by snow (0 OR 1) [-]
forc_u [double] atmospheric wind speed in east direction (m/s)
forc_v [double] atmospheric wind speed in north direction (m/s)
forc_q_in [double] atmospheric specific humidity (kg/kg)
forc_th_in [double] atmospheric potential temperature (Kelvin)
forc_hgt_u_patch_in [double] observational height of wind at pft level [m]
thm_in [double] intermediate variable (forc_t+0.0098*forc_hgt_t_patch)
thv_in [double] virtual potential temperature (kelvin)
t_grnd [double] ground temperature (Kelvin)
qg [double] ground specific humidity [kg/kg]
z0mg_in [double] roughness length over ground, momentum [m]
OUTPUTS:
dlrad [double] downward longwave radiation below the canopy [W/m2]
......@@ -53,20 +66,28 @@ ulrad [double] upward longwave radiation above the canopy [W/m2]
*/
void InitializeFlux(
const LandType& Land,
const int& frac_veg_nosno,
const int& frac_veg_nosno_in,
const double& forc_u,
const double& forc_v,
const double& forc_q,
const double& forc_th,
const double& forc_hgt_u_patch,
const double& thm,
const double& thv,
const double& forc_q_in,
const double& forc_th_in,
const double& forc_hgt_u_patch_in,
const double& thm_in,
const double& thv_in,
const double& t_grnd,
const double& qg,
const double& z0mg,
const double& z0mg_in,
double& dlrad,
double& ulrad)
{
frac_veg_nosno = frac_veg_nosno_in;
forc_q = forc_q_in;
forc_th = forc_th_in;
forc_hgt_u_patch = forc_hgt_u_patch_in;
thm = thm_in;
thv = thv_in;
z0mg = z0mg_in;
if (!Land.lakpoi && !Land.urbpoi && frac_veg_nosno == 0) {
ur = std::max(1.0, std::sqrt(forc_u * forc_u + forc_v * forc_v));
dth = thm - t_grnd;
......@@ -131,6 +152,8 @@ z0qg [double] roughness length over ground, latent heat [m]
FrictionVelocityWind(forc_hgt_u_patch, displa, um, obu, z0mg, ustar);
FrictionVelocityTemperature(forc_hgt_t_patch, displa, obu, z0hg, temp1);
FrictionVelocityHumidity(forc_hgt_q_patch, forc_hgt_t_patch, displa, obu, z0hg, z0qg, temp1, temp2);
FrictionVelocityTemperature2m(obu, z0hg, temp12m);
FrictionVelocityHumidity2m(obu, z0hg, z0qg, temp12m, temp22m);
if (!Land.lakpoi && !Land.urbpoi && frac_veg_nosno == 0) {
tstar = temp1 * dth;
......@@ -174,6 +197,7 @@ qg_snow [double] specific humidity at snow surface [kg/kg]
qg_soil [double] specific humidity at soil surface [kg/kg]
qg_h2osfc [double] specific humidity at h2osfc surface [kg/kg]
t_soisno[nlevgrnd+nlevsno] [double] col soil temperature (Kelvin)
forc_pbot [double] atmospheric pressure (Pa)
OUTPUTS:
cgrnds [double] deriv, of soil sensible heat flux wrt soil temp [w/m2/k]
......@@ -189,6 +213,11 @@ qflx_evap_tot [double] qflx_evap_soi + qflx_evap_can + qflx_tran_ve
qflx_ev_snow [double] evaporation flux from snow (W/m**2) [+ to atm]
qflx_ev_soil [double] evaporation flux from soil (W/m**2) [+ to atm]
qflx_ev_h2osfc [double] evaporation flux from h2osfc (W/m**2) [+ to atm]
t_ref2m [double] 2 m height surface air temperature (Kelvin)
t_ref2m_r [double] Rural 2 m height surface air temperature (Kelvin)
q_ref2m [double] 2 m height surface specific humidity (kg/kg)
rh_ref2m_r [double] Rural 2 m height surface relative humidity (%)
rh_ref2m [double] 2 m height surface relative humidity (%)
*/
void ComputeFlux(
const LandType& Land,
......@@ -205,6 +234,7 @@ qflx_ev_h2osfc [double] evaporation flux from h2osfc (W/m**2) [+ to
const double& qg_soil,
const double& qg_h2osfc,
const double *t_soisno,
const double& forc_pbot,
double& cgrnds,
double& cgrndl,
double& cgrnd,
......@@ -217,11 +247,17 @@ qflx_ev_h2osfc [double] evaporation flux from h2osfc (W/m**2) [+ to
double& qflx_evap_tot,
double& qflx_ev_snow,
double& qflx_ev_soil,
double& qflx_ev_h2osfc) {
double& qflx_ev_h2osfc,
double& t_ref2m,
double& t_ref2m_r,
double& q_ref2m,
double& rh_ref2m,
double& rh_ref2m_r) {
if (!Land.lakpoi && !Land.urbpoi && frac_veg_nosno == 0) {
double rah, raw, raih, raiw;
double e_ref2m, de2mdT, qsat_ref2m, dqsat2mdT;
// Determine aerodynamic resistances
rah = 1.0 / (temp1 * ustar);
raw = 1.0 / (temp2 * ustar);
......@@ -259,6 +295,22 @@ qflx_ev_h2osfc [double] evaporation flux from h2osfc (W/m**2) [+ to
qflx_ev_snow = -raiw * (forc_q - qg_snow);
qflx_ev_soil = -raiw * (forc_q - qg_soil);
qflx_ev_h2osfc = -raiw * (forc_q - qg_h2osfc);
// 2 m height air temperature
t_ref2m = thm + temp1 * dth * (1.0 / temp12m - 1.0 / temp1);
// 2 m height specific humidity
q_ref2m = forc_q + temp2 * dqh * (1.0 / temp22m - 1.0 / temp2);
// 2 m height relative humidity
QSat(t_ref2m, forc_pbot, e_ref2m, de2mdT, qsat_ref2m, dqsat2mdT);
rh_ref2m = std::min(100.0, (q_ref2m / qsat_ref2m * 100.0));
if (Land.ltype == istsoil || Land.ltype == istcrop) {
rh_ref2m_r = rh_ref2m;
t_ref2m_r = t_ref2m;
}
}
} // ComputeFlux
......
......@@ -61,6 +61,8 @@ private:
double zldis; // reference height "minus" zero displacement height [m]
double temp1; // relation for potential temperature profile
double temp2; // relation for specific humidity profile
double temp12m; // relation for potential temperature profile applied at 2-m
double temp22m; // relation for specific humidity profile applied at 2-m
double tlbef; // leaf temperature from previous iteration [K]
double delq; // temporary
double dt_veg; // change in t_veg, last iteration (Kelvin)
......@@ -484,6 +486,8 @@ void StabilityIteration(
FrictionVelocityWind(forc_hgt_u_patch, displa, um, obu, z0mv, ustar);
FrictionVelocityTemperature(forc_hgt_t_patch, displa, obu, z0hv, temp1);
FrictionVelocityHumidity(forc_hgt_q_patch, forc_hgt_t_patch, displa, obu, z0hv, z0qv, temp1, temp2);
FrictionVelocityTemperature2m(obu, z0hv, temp12m);
FrictionVelocityHumidity2m(obu, z0hv, z0qv, temp12m, temp22m);
// save leaf temp and leaf temp delta from previous iteration
tlbef = t_veg;
......@@ -714,7 +718,12 @@ dlrad [double] downward longwave radiation below the canopy [W/m2]
ulrad [double] upward longwave radiation above the canopy [W/m2]
cgrnds [double] deriv, of soil sensible heat flux wrt soil temp [w/m2/k]
cgrndl [double] deriv of soil latent heat flux wrt soil temp [w/m**2/k]
cgrnd [double] deriv. of soil energy flux wrt to soil temp [w/m2/k]
cgrnd [double] deriv. of soil energy flux wrt to soil temp [w/m2/k]
t_ref2m [double] 2 m height surface air temperature (Kelvin)
t_ref2m_r [double] Rural 2 m height surface air temperature (Kelvin)
q_ref2m [double] 2 m height surface specific humidity (kg/kg)
rh_ref2m_r [double] Rural 2 m height surface relative humidity (%)
rh_ref2m [double] 2 m height surface relative humidity (%)
*/
void ComputeFlux(
const LandType& Land,
......@@ -747,22 +756,29 @@ void ComputeFlux(
double& ulrad,
double& cgrnds,
double& cgrndl,
double& cgrnd
)
double& cgrnd,
double& t_ref2m,
double& t_ref2m_r,
double& q_ref2m,
double& rh_ref2m,
double& rh_ref2m_r)
{
if (!Land.lakpoi && !Land.urbpoi && frac_veg_nosno != 0) {
t_veg_out = t_veg;
double e_ref2m, de2mdT, qsat_ref2m, dqsat2mdT;
// Energy balance check in canopy
double lw_grnd = (frac_sno * pow(t_soisno[nlevsno-snl], 4.0) + (1.0 - frac_sno - frac_h2osfc) *
pow(t_soisno[nlevsno], 4.0) + frac_h2osfc * pow(t_h2osfc, 4.0));
double err = sabv + air + bir * pow(tlbef, 3.0) * (tlbef + 4.0 * dt_veg) + cir * lw_grnd -
eflx_sh_veg - hvap * qflx_evap_veg;
// Fluxes from ground to canopy space
double delt = wtal * t_grnd - wtl0 * t_veg - wta0 * thm;
eflx_sh_grnd = cpair * forc_rho * wtg * delt;
// compute individual sensible heat fluxes
double delt_snow = wtal * t_soisno[nlevsno-snl] - wtl0 * t_veg - wta0 * thm;
eflx_sh_snow = cpair * forc_rho * wtg * delt_snow;
......@@ -771,6 +787,7 @@ void ComputeFlux(
double delt_h2osfc = wtal * t_h2osfc - wtl0 * t_veg - wta0 * thm;
eflx_sh_h2osfc = cpair * forc_rho * wtg * delt_h2osfc;
qflx_evap_soi = forc_rho * wtgq * delq;
// compute individual latent heat fluxes
double delq_snow = wtalq * qg_snow - wtlq0 * qsatl - wtaq0 * forc_q;
qflx_ev_snow = forc_rho * wtgq * delq_snow;
......@@ -778,6 +795,17 @@ void ComputeFlux(
qflx_ev_soil = forc_rho * wtgq * delq_soil;
double delq_h2osfc = wtalq * qg_h2osfc - wtlq0 * qsatl - wtaq0 * forc_q;
qflx_ev_h2osfc = forc_rho * wtgq * delq_h2osfc;
// 2 m height air temperature
t_ref2m = thm + temp1 * dth * (1.0 / temp12m - 1.0 / temp1);
t_ref2m_r = t_ref2m;
// 2 m height specific humidity
q_ref2m = forc_q + temp2 * dqh * (1.0 / temp22m - 1.0 / temp2);
// 2 m height relative humidity
QSat(t_ref2m, forc_pbot, e_ref2m, de2mdT, qsat_ref2m, dqsat2mdT);
rh_ref2m = std::min(100.0, q_ref2m / qsat_ref2m * 100.0);
rh_ref2m_r = rh_ref2m;
// Downward longwave radiation below the canopy
dlrad = (1.0 - emv) * emg * forc_lwrad + emv * emg * sb * pow(tlbef, 3.0) * (tlbef + 4.0 * dt_veg);
// Upward longwave radiation above the canopy
......@@ -797,8 +825,6 @@ void ComputeFlux(
// if (abs(err(p)) > 0.1_r8)
}
......
/* functions from FrictionVelocityMod.F90
FrictionVelocity:
Split into 3 functions - wind, temperature, humidity
!!DON'T NEED 2-meter variables, they are only used as feedbacks to ATM model and don't affect LSM and subsurface calculations!
Split into 5 functions - wind, temperature, humidity, 2-m temp, 2-m humidity
Also don't need u10 variables, they are used for dust model
The friction velocity scheme is based on the work of Zeng et al. (1998):
......@@ -43,7 +42,7 @@ void MoninObukIni(
double rib; // bulk Richardson number
double zeta; // dimensionless height used in Monin-Obukhov theory
// double ustar = 0.06; // friction velocity [m/s] -- in CLM, but not used
double wc = 0.5; // convective velocity [m/s]
const double wc = 0.5; // convective velocity [m/s]
// wind speed
if (dthv >= 0.0) {
......@@ -108,9 +107,9 @@ void FrictionVelocityWind(
const double& z0m,
double& ustar)
{
double zetam = 1.574; // transition point of flux-gradient relation (wind profile)
double zldis = forc_hgt_u_patch - displa; // reference height "minus" zero displacement heght [m]
double zeta = zldis / obu; // dimensionless height used in Monin-Obukhov theory
const double zetam = 1.574; // transition point of flux-gradient relation (wind profile)
const double zldis = forc_hgt_u_patch - displa; // reference height "minus" zero displacement heght [m]
const double zeta = zldis / obu; // dimensionless height used in Monin-Obukhov theory
if (zeta < (-zetam)) {
ustar = vkc * um / (std::log(-zetam * obu / z0m) - StabilityFunc1(-zetam) +
......@@ -145,9 +144,9 @@ void FrictionVelocityTemperature(
const double& z0h,
double& temp1)
{
double zetat = 0.465; // transition point of flux-gradient relation (temp. profile)
double zldis = forc_hgt_t_patch - displa; // reference height "minus" zero displacement heght [m]
double zeta = zldis / obu; // dimensionless height used in Monin-Obukhov theory
const double zetat = 0.465; // transition point of flux-gradient relation (temp. profile)
const double zldis = forc_hgt_t_patch - displa; // reference height "minus" zero displacement heght [m]
const double zeta = zldis / obu; // dimensionless height used in Monin-Obukhov theory
if (zeta < (-zetat)) {
temp1 = vkc / (std::log(-zetat * obu / z0h) - StabilityFunc2(-zetat) +
......@@ -187,7 +186,7 @@ void FrictionVelocityHumidity(
const double& temp1,
double& temp2)
{
double zetat = 0.465; // transition point of flux-gradient relation (temp. profile)
const double zetat = 0.465; // transition point of flux-gradient relation (temp. profile)
if (forc_hgt_q_patch == forc_hgt_t_patch && z0q == z0h) {
temp2 = temp1;
......@@ -207,3 +206,66 @@ void FrictionVelocityHumidity(
}
}
/* Temperature profile applied at 2-m
INPUTS:
obu [double] monin-obukhov length (m)
z0h [double] roughness length over vegetation, sensible heat [m]
OUTPUTS:
temp12m [double] relation for potential temperature profile applied at 2-m
*/
void FrictionVelocityTemperature2m(
const double& obu,
const double& z0h,
double& temp12m)
{
const double zldis = 2.0 + z0h;
const double zeta = zldis / obu;
const double zetat = 0.465; // transition point of flux-gradient relation (temp. profile)
if (zeta < -zetat) {
temp12m = vkc / (std::log(-zetat * obu / z0h) - StabilityFunc2(-zetat) + StabilityFunc2(z0h / obu)
+ 0.8 * (pow(zetat, -0.333) - pow(-zeta, -0.333)));
} else if (zeta < 0.0) {
temp12m = vkc / (std::log(zldis / z0h) - StabilityFunc2(zeta) + StabilityFunc2(z0h / obu));
} else if (zeta <= 1.0) {
temp12m = vkc / (std::log(zldis / z0h) + 5.0 * zeta - 5.0 * z0h / obu);
} else {
temp12m = vkc / (std::log(obu / z0h) + 5.0 - 5.0 * (z0h / obu) + (5.0 * std::log(zeta) + zeta - 1.0));
}
}
/* Humidity profile applied at 2-m
INPUTS:
obu [double] monin-obukhov length (m)
z0h [double] roughness length over vegetation, sensible heat [m]
z0q [double] roughness length over vegetation, latent heat [m]
OUTPUTS:
relation for specific humidity profile applied at 2-m
*/
void FrictionVelocityHumidity2m(
const double& obu,
const double& z0h,
const double& z0q,
const double& temp12m,
double& temp22m)
{
if (z0q == z0h) {
temp22m = temp12m;
} else {
const double zldis = 2.0 + z0q;
const double zeta = zldis / obu;
const double zetat = 0.465; // transition point of flux-gradient relation (temp. profile)
if (zeta < -zetat) {
temp22m = vkc / (std::log(-zetat * obu / z0q) - StabilityFunc2(-zetat) + StabilityFunc2(z0q / obu)
+ 0.8 * (pow(zetat, -0.333) - pow(-zeta, -0.333)));
} else if (zeta < 0.0) {
temp22m = vkc / (std::log(zldis / z0q) - StabilityFunc2(zeta) + StabilityFunc2(z0q / obu));
} else if (zeta <= 1.0) {
temp22m = vkc / (std::log(zldis / z0q) + 5.0 * zeta - 5.0 * z0q / obu);
} else {
temp22m = vkc / (std::log(obu / z0q) + 5.0 - 5.0 * z0q / obu + (5.0 * std::log(zeta) + zeta - 1.0));
}
}
}
\ No newline at end of file
......@@ -3,7 +3,48 @@
*/
#include "clm_constants.h"
#include "qsat.h"
// For water vapor (temperature range 0C-100C)
static const double a0 = 6.11213476;
static const double a1 = 0.444007856;
static const double a2 = 0.143064234e-01;
static const double a3 = 0.264461437e-03;
static const double a4 = 0.305903558e-05;
static const double a5 = 0.196237241e-07;
static const double a6 = 0.892344772e-10;
static const double a7 = -0.373208410e-12;
static const double a8 = 0.209339997e-15;
// For derivative:water vapor
static const double b0 = 0.444017302;
static const double b1 = 0.286064092e-01;
static const double b2 = 0.794683137e-03;
static const double b3 = 0.121211669e-04;
static const double b4 = 0.103354611e-06;
static const double b5 = 0.404125005e-09;
static const double b6 = -0.788037859e-12;
static const double b7 = -0.114596802e-13;
static const double b8 = 0.381294516e-16;
// For ice (temperature range -75C-0C)
static const double c0 = 6.11123516;
static const double c1 = 0.503109514;
static const double c2 = 0.188369801e-01;
static const double c3 = 0.420547422e-03;
static const double c4 = 0.614396778e-05;
static const double c5 = 0.602780717e-07;
static const double c6 = 0.387940929e-09;
static const double c7 = 0.149436277e-11;
static const double c8 = 0.262655803e-14;
// For derivative:ice
static const double d0 = 0.503277922;
static const double d1 = 0.377289173e-01;
static const double d2 = 0.126801703e-02;
static const double d3 = 0.249468427e-04;
static const double d4 = 0.313703411e-06;
static const double d5 = 0.257180651e-08;
static const double d6 = 0.133268878e-10;
static const double d7 = 0.394116744e-13;
static const double d8 = 0.498070196e-16;
void QSat(
const double& T,
......
......@@ -4,14 +4,18 @@
void MoninObukIni(const double& ur, const double& thv, const double& dthv, const double& zldis, const double& z0m, double& um,
double& obu);
void FrictionVelocityWind( const double& forc_hgt_u_patch, const double& displa, const double& um, const double& obu,
void FrictionVelocityWind(const double& forc_hgt_u_patch, const double& displa, const double& um, const double& obu,
const double& z0m, double& ustar);
void FrictionVelocityTemperature( const double& forc_hgt_t_patch, const double& displa, const double& obu, const double& z0h,
void FrictionVelocityTemperature(const double& forc_hgt_t_patch, const double& displa, const double& obu, const double& z0h,
double& temp1);
void FrictionVelocityHumidity( const double& forc_hgt_q_patch, const double& forc_hgt_t_patch, const double& displa,
void FrictionVelocityHumidity(const double& forc_hgt_q_patch, const double& forc_hgt_t_patch, const double& displa,
const double& obu, const double& z0h, const double& z0q, const double& temp1, double& temp2);
void FrictionVelocityTemperature2m(const double& obu, const double& z0h, double& temp12m);
void FrictionVelocityHumidity2m(const double& obu, const double& z0h, const double& z0q, const double& temp12m, double& temp22m);
#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