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

converted SurfaceRadiation and SurfaceResistance to new indexing

parent 27d38a64
......@@ -32,6 +32,7 @@ DESCRIPTION: zero out fluxes before surface radiation calculations
INPUTS:
ltype [int] landunit type
urbpoi [bool] true => landunit is an urban point
OUTPUTS:
sabg_soil [double] solar radiation absorbed by soil (W/m**2)
......@@ -45,6 +46,7 @@ sabg_lyr[nlevsno+1] [double] absorbed radiative flux (pft,lyr) [W/m2]
*/
void SurfRadZeroFluxes (
const int& ltype,
const bool& urbpoi,
double& sabg_soil,
double& sabg_snow,
......@@ -56,7 +58,7 @@ void SurfRadZeroFluxes (
double* sabg_lyr)
{
// Initialize fluxes
if (ltype < isturb_MIN) {
if (!urbpoi) {
sabg_soil = 0.0;
sabg_snow = 0.0;
sabg = 0.0;
......@@ -69,7 +71,7 @@ void SurfRadZeroFluxes (
// zero-out fsun for the urban patches
// the non-urban patches were set prior to this call
// and split into ed and non-ed specific functions
if (ltype >= isturb_MIN && ltype <= isturb_MAX) { fsun = 0.0; }
if (urbpoi) { fsun = 0.0; }
}
/* SurfRadAbsorbed()
......@@ -77,6 +79,7 @@ DESCRIPTION: calculate solar flux absorbed by canopy, soil, snow, and ground
INPUTS:
ltype [int] landunit type
urbpoi [bool] true => landunit is an urban point
nband [int] number of solar radiation waveband classes (equal to numrad)
snl [int] number of snow layers
ftdd[numrad] [double] down direct flux below canopy per unit direct flux
......@@ -107,6 +110,7 @@ tri[numrad] [double] transmitted solar radiation: diffuse (W/m**2)
*/
void SurfRadAbsorbed(
const int& ltype,
const bool& urbpoi,
const int& nband,
const int& snl,
const double* ftdd,
......@@ -135,7 +139,7 @@ void SurfRadAbsorbed(
{
double absrad, cad[nband], cai[nband];
if (ltype < isturb_MIN) {
if (!urbpoi) {
for (int ib = 0; ib < nband; ib++) {
cad[ib] = forc_solad[ib] * fabd[ib];
......@@ -178,6 +182,7 @@ DESCRIPTION: compute absorbed flux in each snow layer and top soil layer
INPUTS:
ltype [int] landunit type
urbpoi [bool] true => landunit is an urban point
snl [int] number of snow layers
subgridflag [int] use subgrid fluxes
sabg [double] solar radiation absorbed by ground (W/m**2)
......@@ -197,6 +202,7 @@ sabg_lyr[nlevsno+1] [double] absorbed radiative flux (pft,lyr) [W/m2]
*/
void SurfRadLayers(
const int& ltype,
const bool& urbpoi,
const int& snl,
const int& subgridflag,
const double& sabg,
......@@ -218,7 +224,7 @@ void SurfRadLayers(
// compute absorbed flux in each snow layer and top soil layer,
// based on flux factors computed in the radiative transfer portion of SNICAR.
if (ltype < isturb_MIN) {
if (!urbpoi) {
sabg_snl_sum = 0.0;
sub_surf_abs_SW = 0.0;
......@@ -305,7 +311,7 @@ void SurfRadLayers(
DESCRIPTION: calculate reflected solar radiation
INPUTS:
ltype [int] landunit type
urbpoi [bool] true => landunit is an urban point
albd[numrad] [double] surface albedo (direct)
albi[numrad] [double] surface albedo (diffuse)
forc_solad[numrad] [double] direct beam radiation (W/m**2)
......@@ -315,7 +321,7 @@ OUTPUTS:
fsr [double] solar radiation reflected (W/m**2)
*/
void SurfRadReflected(
const int& ltype,
const bool& urbpoi,
const double* albd,
const double* albi,
const double* forc_solad,
......@@ -326,14 +332,12 @@ void SurfRadReflected(
double fsr_vis_d, fsr_nir_d, fsr_vis_i, fsr_nir_i, rvis, rnir;
// Radiation diagnostics
if (ltype < isturb_MIN) {
if (!urbpoi) {
// NDVI and reflected solar radiation
rvis = albd[0] * forc_solad[0] + albi[0] * forc_solai[0];
rnir = albd[1] * forc_solad[1] + albi[1] * forc_solai[1];
fsr = rvis + rnir;
}
if (ltype >= isturb_MIN && ltype <= isturb_MAX) {
} else {
// Solar reflected per unit ground area (roof, road) and per unit wall area (sunwall, shadewall)
fsr_vis_d = albd[0] * forc_solad[0];
fsr_nir_d = albd[1] * forc_solad[1];
......@@ -356,7 +360,7 @@ This subroutine calculates and returns:
7) sunlit fraction of canopy
INPUTS:
ltype [int] landunit type
urbpoi [bool] true => landunit is an urban point
nrad [int] number of canopy layers
elai [double] one-sided leaf area index
tlai_z[nlevcan] [double] tlai increment for canopy layer
......@@ -379,7 +383,7 @@ fsun [double] sunlit fraction of canopy
*/
void CanopySunShadeFractions(
const int& ltype,
const bool& urbpoi,
const int& nrad,
const double& elai,
const double* tlai_z,
......@@ -399,7 +403,7 @@ void CanopySunShadeFractions(
double& laisha,
double& fsun)
{
if (ltype < isturb_MIN) {
if (!urbpoi) {
int ipar = 0; // The band index for PAR
for (int iv = 0; iv < nrad; iv++) {
parsun_z[iv] = 0.0;
......
......@@ -13,7 +13,7 @@ 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
lakpoi [bool] true => landunit is a lake 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)
......@@ -30,36 +30,27 @@ void calc_soilevap_stress(
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_,
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) {
// local variables
double fac, fac_fc, wx;
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]);
wx = (h2osoi_liq[nlevsno] / denh2o + h2osoi_ice[nlevsno] / denice) / dz[nlevsno];
fac = std::min(1.0, wx / watsat[0]);
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
if (wx < watfc[0] ) { //!when water content of ths top layer is less than that at F.C.
fac_fc = std::min(1.0, wx / watfc[0]); // 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)
......
......@@ -30,17 +30,32 @@ static const int icol_shadewall = isturb_MIN*10 + 3;
static const int icol_road_imperv = isturb_MIN*10 + 4;
static const int icol_road_perv = isturb_MIN*10 + 5;
//from clm_varcon.F90
static const double snw_rds_min = 54.526; //minimum allowed snow effective radius (also "fresh snow" value) [microns]
static const double zlnd = 0.01; // Roughness length for soil [m]
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;
static const double SHR_CONST_BOLTZ = 1.38065e-23; // Boltzmann's constant ~ J/K/molecule
static const double SHR_CONST_AVOGAD = 6.02214e26; // Avogadro's number ~ molecules/kmole
static const double SHR_CONST_MWWV = 18.016; // molecular weight water vapor
static const double SHR_CONST_RGAS = SHR_CONST_AVOGAD*SHR_CONST_BOLTZ; // Universal gas constant ~ J/K/kmole
static const double SHR_CONST_RWV = SHR_CONST_RGAS/SHR_CONST_MWWV; // Water vapor gas constant ~ J/K/kg
static const double SHR_CONST_G = 9.80616; // acceleration of gravity ~ m/s^2
static const double SHR_CONST_LATICE = 3.337e5; // latent heat of fusion ~ J/kg
static const double SHR_CONST_LATVAP = 2.501e6; // latent heat of evaporation ~ J/kg
static const double SHR_CONST_LATSUB = SHR_CONST_LATICE + SHR_CONST_LATVAP; // latent heat of sublimation ~ J/kg
//from clm_varcon.F90
static const double snw_rds_min = 54.526; //minimum allowed snow effective radius (also "fresh snow" value) [microns]
static const double zlnd = 0.01; // Roughness length for soil [m]
static const double zsno = 0.0024; // Roughness length for snow [m]
static const double tfrz = SHR_CONST_TKFRZ; // freezing temperature [K]
static const double roverg = SHR_CONST_RWV/SHR_CONST_G*1000.0; // Rw/g constant = (8.3144/0.018)/(9.80616)*1000. mm/K
static const double denice = SHR_CONST_RHOICE;
static const double denh2o = SHR_CONST_RHOFW;
static const double hvap = SHR_CONST_LATVAP; // Latent heat of evap for water [J/kg]
static const double hsub = SHR_CONST_LATSUB; // Latent heat of sublimation [J/kg]
static const double hfus = SHR_CONST_LATICE; // Latent heat of fusion for ice [J/kg]
#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