Skip to content

GitLab

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

minor formating changes; added numrad to header

parent 8ebaa41a
Hide whitespace changes
Inline Side-by-side
INDEX_MAP_README.txt
View file @ bf864188
......@@ -124,4 +124,5 @@ grid grid
| 0| |4|
----- 0 ground interface 5 ----- [nlevsno]
| 1| |5|
-- 1 6 --- [nlevsno + 1]
\ No newline at end of file
-- 1 6 --- [nlevsno + 1]
\ No newline at end of file
src/CMakeLists.txt
View file @ bf864188
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)
add_library (physics physics_dummy.cc CanopyHydrology.cc coldstart_topography.cc coldstart_snow.cc SurfaceRadiation.cc SurfaceResistanceMod_functions.cc QSatMod_functions.cc CanopyTemperatureMod_functions.cc FrictionVelocityMod_functions.cc)
target_include_directories (physics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
src/FrictionVelocityMod_functions.cc
View file @ bf864188
......@@ -97,7 +97,7 @@ forc_hgt_u_patch [double] observational height of wind at pft level [m]
displa [double] displacement height (m)
um [double] wind speed including the stability effect [m/s]
obu [double] monin-obukhov length (m)
z0m [double] roughness length over vegetation, momentum [m] [lbn:ubn]
z0m [double] roughness length over vegetation, momentum [m]
OUTPUTS:
ustar [double] friction velocity [m/s]
......
src/SurfaceRadiation.cc
View file @ bf864188
......@@ -39,35 +39,35 @@ fsa_r [double] rural solar radiation absorbed (total) (W/m**2)
fsun [double] sunlit fraction of canopy
sabg_lyr[nlevsno+1] [double] absorbed radiative flux (pft,lyr) [W/m2]
*/
void SurfRadZeroFluxes (
const LandType& Land,
double& sabg_soil,
double& sabg_snow,
double& sabg,
double& sabv,
double& fsa,
double& fsa_r,
double& fsun,
double* sabg_lyr)
{
// Initialize fluxes
if (!Land.urbpoi) {
sabg_soil = 0.0;
sabg_snow = 0.0;
sabg = 0.0;
sabv = 0.0;
fsa = 0.0;
if (Land.ltype == istsoil || Land.ltype == istcrop) { fsa_r = 0.0; }
for (int j = 0; j < nlevsno + 1; j++) { sabg_lyr[j] = 0.0; }
void SurfRadZeroFluxes (
const LandType& Land,
double& sabg_soil,
double& sabg_snow,
double& sabg,
double& sabv,
double& fsa,
double& fsa_r,
double& fsun,
double* sabg_lyr) {
// Initialize fluxes
if (!Land.urbpoi) {
sabg_soil = 0.0;
sabg_snow = 0.0;
sabg = 0.0;
sabv = 0.0;
fsa = 0.0;
if (Land.ltype == istsoil || Land.ltype == istcrop) { fsa_r = 0.0; }
for (int j = 0; j < nlevsno + 1; j++) { sabg_lyr[j] = 0.0; }
}
// 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 (Land.urbpoi) { fsun = 0.0; }
}
// 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 (Land.urbpoi) { fsun = 0.0; }
}
/* SurfRadAbsorbed()
/* SurfaceRadiation::SurfRadAbsorbed()
DESCRIPTION: calculate solar flux absorbed by canopy, soil, snow, and ground
INPUTS:
......@@ -96,71 +96,70 @@ sabg_soil [double] solar radiation absorbed by soil (W/m**2)
sabg_snow [double] solar radiation absorbed by snow (W/m**2))
*/
void SurfRadAbsorbed(
const LandType& Land,
const int& snl,
const double* ftdd,
const double* ftid,
const double* ftii,
const double* forc_solad,
const double* forc_solai,
const double* fabd,
const double* fabi,
const double* albsod,
const double* albsoi,
const double* albsnd_hst,
const double* albsni_hst,
const double* albgrd,
const double* albgri,
double& fsa_r,
double& sabv,
double& fsa,
double& sabg,
double& sabg_soil,
double& sabg_snow)
{
double absrad, cad[numrad], cai[numrad];
if (!Land.urbpoi) {
for (int ib = 0; ib < numrad; ib++) {
cad[ib] = forc_solad[ib] * fabd[ib];
cai[ib] = forc_solai[ib] * fabi[ib];
sabv += cad[ib] + cai[ib];
fsa += cad[ib] + cai[ib];
if (ib == 0) { parveg = cad[ib] + cai[ib]; }
if (Land.ltype == istsoil || Land.ltype == istcrop) { fsa_r += cad[ib] + cai[ib]; }
// Transmitted = solar fluxes incident on ground
trd[ib] = forc_solad[ib] * ftdd[ib];
tri[ib] = forc_solad[ib] * ftid[ib] + forc_solai[ib]*ftii[ib];
// Solar radiation absorbed by ground surface
// calculate absorbed solar by soil/snow separately
absrad = trd[ib] * (1.0 - albsod[ib]) + tri[ib] * (1.0 - albsoi[ib]);
sabg_soil += absrad;
absrad = trd[ib] * (1.0 - albsnd_hst[ib]) + tri[ib] * (1.0 - albsni_hst[ib]);
sabg_snow += absrad;
absrad = trd[ib]*(1.0 - albgrd[ib]) + tri[ib] * (1.0 - albgri[ib]);
sabg += absrad;
fsa += absrad;
if (Land.ltype == istsoil || Land.ltype == istcrop) { fsa_r += absrad; }
if (snl == 0) {
sabg_snow = sabg;
sabg_soil = sabg;
}
if (subgridflag == 0) {
sabg_snow = sabg;
sabg_soil = sabg;
}
} // end of numrad
} // end if not urban
}
void SurfRadAbsorbed(
const LandType& Land,
const int& snl,
const double* ftdd,
const double* ftid,
const double* ftii,
const double* forc_solad,
const double* forc_solai,
const double* fabd,
const double* fabi,
const double* albsod,
const double* albsoi,
const double* albsnd_hst,
const double* albsni_hst,
const double* albgrd,
const double* albgri,
double& fsa_r,
double& sabv,
double& fsa,
double& sabg,
double& sabg_soil,
double& sabg_snow) {
double absrad, cad[numrad], cai[numrad];
if (!Land.urbpoi) {
for (int ib = 0; ib < numrad; ib++) {
cad[ib] = forc_solad[ib] * fabd[ib];
cai[ib] = forc_solai[ib] * fabi[ib];
sabv += cad[ib] + cai[ib];
fsa += cad[ib] + cai[ib];
if (ib == 0) { parveg = cad[ib] + cai[ib]; }
if (Land.ltype == istsoil || Land.ltype == istcrop) { fsa_r += cad[ib] + cai[ib]; }
// Transmitted = solar fluxes incident on ground
trd[ib] = forc_solad[ib] * ftdd[ib];
tri[ib] = forc_solad[ib] * ftid[ib] + forc_solai[ib]*ftii[ib];
// Solar radiation absorbed by ground surface
// calculate absorbed solar by soil/snow separately
absrad = trd[ib] * (1.0 - albsod[ib]) + tri[ib] * (1.0 - albsoi[ib]);
sabg_soil += absrad;
absrad = trd[ib] * (1.0 - albsnd_hst[ib]) + tri[ib] * (1.0 - albsni_hst[ib]);
sabg_snow += absrad;
absrad = trd[ib]*(1.0 - albgrd[ib]) + tri[ib] * (1.0 - albgri[ib]);
sabg += absrad;
fsa += absrad;
if (Land.ltype == istsoil || Land.ltype == istcrop) { fsa_r += absrad; }
if (snl == 0) {
sabg_snow = sabg;
sabg_soil = sabg;
}
if (subgridflag == 0) {
sabg_snow = sabg;
sabg_soil = sabg;
}
} // end of numrad
} // end if not urban
}
/* SurfRadLayers()
/* SurfaceRadiation::SurfRadLayers()
DESCRIPTION: compute absorbed flux in each snow layer and top soil layer
INPUTS:
......@@ -180,111 +179,109 @@ sub_surf_abs_SW [double] percent of solar radiation absorbed below first sn
sabg_pen [double] solar (rural) radiation penetrating top soisno layer (W/m**2)
sabg_lyr[nlevsno+1] [double] absorbed radiative flux (pft,lyr) [W/m2]
*/
void SurfRadLayers(
const LandType& Land,
const int& snl,
const int& subgridflag,
const double& sabg,
const double& sabg_snow,
const double& snow_depth,
const double* flx_absdv,
const double* flx_absdn,
const double* flx_absiv,
const double* flx_absin,
double& sub_surf_abs_SW,
double& sabg_pen,
double* sabg_lyr)
{
double err_sum = 0.0;
double sabg_snl_sum;
// compute absorbed flux in each snow layer and top soil layer,
// based on flux factors computed in the radiative transfer portion of SNICAR.
if (!Land.urbpoi) {
sabg_snl_sum = 0.0;
sub_surf_abs_SW = 0.0;
// CASE1: No snow layers: all energy is absorbed in top soil layer
if (snl == 0) {
for (int i = 0; i <= nlevsno; i++) { sabg_lyr[i] = 0.0; }
sabg_lyr[nlevsno] = sabg;
sabg_snl_sum = sabg_lyr[nlevsno];
} else { // CASE 2: Snow layers present: absorbed radiation is scaled according to flux factors computed by SNICAR
for (int i = 0; i < nlevsno + 1; i++) {
sabg_lyr[i] = flx_absdv[i] * trd[0] + flx_absdn[i] * trd[1] + flx_absiv[i] * tri[0] + flx_absin[i] * tri[1];
// summed radiation in active snow layers:
// if snow layer is at or below snow surface
if (i >= nlevsno - snl) { sabg_snl_sum += sabg_lyr[i]; }
// if snow layer is below surface snow layer accumulate subsurface flux as a diagnostic
if (i > nlevsno - snl) { sub_surf_abs_SW += + sabg_lyr[i]; }
}
// Divide absorbed by total, to get % absorbed in subsurface
if (sabg_snl_sum != 0.0) { // inequality with double?
sub_surf_abs_SW = sub_surf_abs_SW / sabg_snl_sum;
} else {
sub_surf_abs_SW = 0.0;
}
// Error handling: The situation below can occur when solar radiation is
// NOT computed every timestep.
// When the number of snow layers has changed in between computations of the
// absorbed solar energy in each layer, we must redistribute the absorbed energy
// to avoid physically unrealistic conditions. The assumptions made below are
// somewhat arbitrary, but this situation does not arise very frequently.
// This error handling is implemented to accomodate any value of the
// radiation frequency.
// change condition to match sabg_snow isntead of sabg
if (std::abs(sabg_snl_sum - sabg_snow) > 0.00001) {
if (snl == 0) {
for (int j = 0; j < nlevsno; j++) { sabg_lyr[j] = 0.0; }
sabg_lyr[nlevsno] = sabg;
} else if (snl == 1) {
for (int j = 0; j < nlevsno - 1; j++) { sabg_lyr[j] = 0.0; }
sabg_lyr[nlevsno - 1] = sabg_snow * 0.6;
sabg_lyr[nlevsno] = sabg_snow * 0.4;
void SurfRadLayers(
const LandType& Land,
const int& snl,
const int& subgridflag,
const double& sabg,
const double& sabg_snow,
const double& snow_depth,
const double* flx_absdv,
const double* flx_absdn,
const double* flx_absiv,
const double* flx_absin,
double& sub_surf_abs_SW,
double& sabg_pen,
double* sabg_lyr) {
double err_sum = 0.0;
double sabg_snl_sum;
// compute absorbed flux in each snow layer and top soil layer,
// based on flux factors computed in the radiative transfer portion of SNICAR.
if (!Land.urbpoi) {
sabg_snl_sum = 0.0;
sub_surf_abs_SW = 0.0;
// CASE1: No snow layers: all energy is absorbed in top soil layer
if (snl == 0) {
for (int i = 0; i <= nlevsno; i++) { sabg_lyr[i] = 0.0; }
sabg_lyr[nlevsno] = sabg;
sabg_snl_sum = sabg_lyr[nlevsno];
} else { // CASE 2: Snow layers present: absorbed radiation is scaled according to flux factors computed by SNICAR
for (int i = 0; i < nlevsno + 1; i++) {
sabg_lyr[i] = flx_absdv[i] * trd[0] + flx_absdn[i] * trd[1] + flx_absiv[i] * tri[0] + flx_absin[i] * tri[1];
// summed radiation in active snow layers:
// if snow layer is at or below snow surface
if (i >= nlevsno - snl) { sabg_snl_sum += sabg_lyr[i]; }
// if snow layer is below surface snow layer accumulate subsurface flux as a diagnostic
if (i > nlevsno - snl) { sub_surf_abs_SW += + sabg_lyr[i]; }
}
// Divide absorbed by total, to get % absorbed in subsurface
if (sabg_snl_sum != 0.0) { // inequality with double?
sub_surf_abs_SW = sub_surf_abs_SW / sabg_snl_sum;
} else {
for (int j = 0; j <= nlevsno; j++) { sabg_lyr[j] = 0.0; }
sabg_lyr[nlevsno - snl] = sabg_snow * 0.75;
sabg_lyr[nlevsno - snl + 1] = sabg_snow * 0.25;
sub_surf_abs_SW = 0.0;
}
}
// If shallow snow depth, all solar radiation absorbed in top or top two snow layers
// to prevent unrealistic timestep soil warming
if (subgridflag == 0) {
if (snow_depth < 0.1) {
// Error handling: The situation below can occur when solar radiation is
// NOT computed every timestep.
// When the number of snow layers has changed in between computations of the
// absorbed solar energy in each layer, we must redistribute the absorbed energy
// to avoid physically unrealistic conditions. The assumptions made below are
// somewhat arbitrary, but this situation does not arise very frequently.
// This error handling is implemented to accomodate any value of the
// radiation frequency.
// change condition to match sabg_snow isntead of sabg
if (std::abs(sabg_snl_sum - sabg_snow) > 0.00001) {
if (snl == 0) {
for (int j = 0; j < nlevsno; j++) { sabg_lyr[j] = 0.0; }
sabg_lyr[nlevsno] = sabg;
} else if (snl == 1) {
for (int j = 0; j < nlevsno - 1; j++) { sabg_lyr[j] = 0.0; }
sabg_lyr[nlevsno - 1] = sabg;
sabg_lyr[nlevsno] = 0.0;
sabg_lyr[nlevsno - 1] = sabg_snow * 0.6;
sabg_lyr[nlevsno] = sabg_snow * 0.4;
} else {
for (int j = 0; j <= nlevsno; j++) { sabg_lyr[j] = 0.0; }
sabg_lyr[nlevsno - snl] = sabg_snow * 0.75;
sabg_lyr[nlevsno - snl + 1] = sabg_snow * 0.25;
}
}
// If shallow snow depth, all solar radiation absorbed in top or top two snow layers
// to prevent unrealistic timestep soil warming
if (subgridflag == 0) {
if (snow_depth < 0.1) {
if (snl == 0) {
for (int j = 0; j < nlevsno; j++) { sabg_lyr[j] = 0.0; }
sabg_lyr[nlevsno] = sabg;
} else if (snl == 1) {
for (int j = 0; j < nlevsno - 1; j++) { sabg_lyr[j] = 0.0; }
sabg_lyr[nlevsno - 1] = sabg;
sabg_lyr[nlevsno] = 0.0;
} else {
for (int j = 0; j <= nlevsno; j++) { sabg_lyr[j] = 0.0; }
sabg_lyr[nlevsno - snl] = sabg_snow * 0.75;
sabg_lyr[nlevsno - snl + 1] = sabg_snow * 0.25;
}
}
}
}
// Error check - This situation should not happen:
for (int j = 0; j <= nlevsno; j++) { err_sum += sabg_lyr[j]; }
assert(!(std::abs(err_sum - sabg_snow) > 0.00001));
// Diagnostic: shortwave penetrating ground (e.g. top layer)
if (Land.ltype == istsoil || Land.ltype == istcrop) { sabg_pen = sabg - sabg_lyr[snl]; }
}
// Error check - This situation should not happen:
for (int j = 0; j <= nlevsno; j++) { err_sum += sabg_lyr[j]; }
assert(!(std::abs(err_sum - sabg_snow) > 0.00001));
// Diagnostic: shortwave penetrating ground (e.g. top layer)
if (Land.ltype == istsoil || Land.ltype == istcrop) { sabg_pen = sabg - sabg_lyr[snl]; }
}
}
/* SurfRadReflected()
/* SurfaceRadiation::SurfRadReflected()
DESCRIPTION: calculate reflected solar radiation
INPUTS:
......@@ -297,35 +294,33 @@ forc_solai[numrad] [double] diffuse radiation (W/m**2)
OUTPUTS:
fsr [double] solar radiation reflected (W/m**2)
*/
void SurfRadReflected(
const LandType& Land,
const double* albd,
const double* albi,
const double* forc_solad,
const double* forc_solai,
double& fsr)
{
double fsr_vis_d, fsr_nir_d, fsr_vis_i, fsr_nir_i, rvis, rnir;
// Radiation diagnostics
if (!Land.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;
} 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];
fsr_vis_i = albi[0] * forc_solai[0];
fsr_nir_i = albi[1] * forc_solai[1];
fsr = fsr_vis_d + fsr_nir_d + fsr_vis_i + fsr_nir_i;
void SurfRadReflected(
const LandType& Land,
const double* albd,
const double* albi,
const double* forc_solad,
const double* forc_solai,
double& fsr) {
double fsr_vis_d, fsr_nir_d, fsr_vis_i, fsr_nir_i, rvis, rnir;
// Radiation diagnostics
if (!Land.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;
} 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];
fsr_vis_i = albi[0] * forc_solai[0];
fsr_nir_i = albi[1] * forc_solai[1];
fsr = fsr_vis_d + fsr_nir_d + fsr_vis_i + fsr_nir_i;
}
}
}
/* CanopySunShadeFractions()
/* SurfaceRadiation::CanopySunShadeFractions()
DESCRIPTION:
This subroutine calculates and returns:
1) absorbed PAR for sunlit leaves in canopy layer
......@@ -358,68 +353,61 @@ laisun [double] sunlit leaf area
laisha [double] shaded leaf area
fsun [double] sunlit fraction of canopy
*/
void CanopySunShadeFractions(
const LandType& Land,
const int& nrad,
const double& elai,
const double* tlai_z,
const double* fsun_z,
const double* forc_solad,
const double* forc_solai,
const double* fabd_sun_z,
const double* fabd_sha_z,
const double* fabi_sun_z,
const double* fabi_sha_z,
double* parsun_z,
double* parsha_z,
double* laisun_z,
double* laisha_z,
double& laisun,
double& laisha,
double& fsun)
{
if (!Land.urbpoi) {
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;
laisun_z[iv] = 0.0;
laisha_z[iv] = 0.0;
}
// Loop over patches to calculate laisun_z and laisha_z for each layer.
// Derive canopy laisun, laisha, and fsun from layer sums.
// If sun/shade big leaf code, nrad=1 and fsun_z(p,1) and tlai_z(p,1) from
// SurfaceAlbedo is canopy integrated so that layer value equals canopy value.
laisun = 0.0;
laisha = 0.0;
for (int iv = 0; iv < nrad; iv++) {
laisun_z[iv] = tlai_z[iv] * fsun_z[iv];
laisha_z[iv] = tlai_z[iv] * (1.0 - fsun_z[iv]);
laisun += laisun_z[iv];
laisha += laisha_z[iv];
}
if (elai > 0.0) {
fsun = laisun / elai;
} else {
fsun = 0.0;
}
// Absorbed PAR profile through canopy
// If sun/shade big leaf code, nrad=1 and fluxes from SurfaceAlbedo
// are canopy integrated so that layer values equal big leaf values.
for (int iv = 0; iv < nrad; iv++) {
parsun_z[iv] = forc_solad[ipar] * fabd_sun_z[iv] + forc_solai[ipar] * fabi_sun_z[iv];
parsha_z[iv] = forc_solad[ipar] * fabd_sha_z[iv] + forc_solai[ipar] * fabi_sha_z[iv];
void CanopySunShadeFractions(
const LandType& Land,
const int& nrad,
const double& elai,
const double* tlai_z,
const double* fsun_z,
const double* forc_solad,
const double* forc_solai,
const double* fabd_sun_z,
const double* fabd_sha_z,
const double* fabi_sun_z,
const double* fabi_sha_z,
double* parsun_z,
double* parsha_z,
double* laisun_z,
double* laisha_z,
double& laisun,
double& laisha,
double& fsun) {
if (!Land.urbpoi) {
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;
laisun_z[iv] = 0.0;
laisha_z[iv] = 0.0;
}
// Loop over patches to calculate laisun_z and laisha_z for each layer.
// Derive canopy laisun, laisha, and fsun from layer sums.
// If sun/shade big leaf code, nrad=1 and fsun_z(p,1) and tlai_z(p,1) from
// SurfaceAlbedo is canopy integrated so that layer value equals canopy value.
laisun = 0.0;
laisha = 0.0;
for (int iv = 0; iv < nrad; iv++) {
laisun_z[iv] = tlai_z[iv] * fsun_z[iv];
laisha_z[iv] = tlai_z[iv] * (1.0 - fsun_z[iv]);
laisun += laisun_z[iv];
laisha += laisha_z[iv];
}
if (elai > 0.0) {
fsun = laisun / elai;
} else {
fsun = 0.0;
}
// Absorbed PAR profile through canopy
// If sun/shade big leaf code, nrad=1 and fluxes from SurfaceAlbedo
// are canopy integrated so that layer values equal big leaf values.
for (int iv = 0; iv < nrad; iv++) {
parsun_z[iv] = forc_solad[ipar] * fabd_sun_z[iv] + forc_solai[ipar] * fabi_sun_z[iv];
parsha_z[iv] = forc_solad[ipar] * fabd_sha_z[iv] + forc_solai[ipar] * fabi_sha_z[iv];
}
}
}
}
};
\ No newline at end of file
};
src/clm_constants.hh
View file @ bf864188
......@@ -22,7 +22,8 @@ static const int subgridflag = 0;
//from clm_varpar.F90