Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update to code for Esoil in sparse vegetation. #315

Merged
merged 5 commits into from
Dec 9, 2015
Merged

Update to code for Esoil in sparse vegetation. #315

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
649243c
Update to code for Esoil in sparse vegetation. I had intended to add…
tbohn Nov 17, 2015
d958240
Added comments.
tbohn Nov 30, 2015
c3b772b
Added comments.
tbohn Nov 30, 2015
19ce733
Replaced hard-coded numbers with constants defined in vicNl_def.h.
tbohn Dec 6, 2015
c48f178
Replaced variable rarc_soil with constant RARC_SOIL.
tbohn Dec 7, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion src/arno_evap.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,9 @@ double arno_evap(layer_data_struct *layer,
double max_infil;
double Evap;
double tmpsum;
/* best-fitting value of architectural resistance, determined from
evaluation of model performance at eddy covariane tower sites */
double rarc_soil = 100.0; // [s/m]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I still don't love this but if @bartnijssen is okay with it, we can keep it for now. Maybe we should #define this in vicNl_def.h?

When we port these changes to VIC.5, rarc_soil needs to go into the parameters_struct or into the soil_con_struct.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, I'm defining a constant RARC_SOIL in vicNl_def.h...


Evap = 0;

Expand All @@ -97,7 +100,7 @@ double arno_evap(layer_data_struct *layer,

/* Calculate the potential bare soil evaporation (mm/time step) */

Epot = penman(air_temp, elevation, rad, vpd, ra, 0.0, 0.0) * delta_t / SEC_PER_DAY;
Epot = penman(air_temp, elevation, rad, vpd, ra, 0.0, rarc_soil) * delta_t / SEC_PER_DAY;

/**********************************************************************/
/* Compute temporary infiltration rate based on given soil_moist. */
Expand Down
70 changes: 57 additions & 13 deletions src/func_surf_energy_bal.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -300,6 +300,9 @@ double func_surf_energy_bal(double Ts, va_list ap)
double tmp_displacement[3];
double tmp_roughness[3];
double tmp_ref_height[3];
double ga_veg;
double ga_bare;
double ga_average;

/************************************
Read variables from variable list
Expand Down Expand Up @@ -687,15 +690,65 @@ double func_surf_energy_bal(double Ts, va_list ap)
TMean, Tair, wind[UnderStory], roughness[UnderStory]);
else
Ra_used[0] = HUGE_RESIST;


/*************************************************
Compute aerodynamic resistance for the case of exposed soil between
plants (or gaps in canopy). Assume plants and exposed soil are well-
mixed, i.e., exposed soil is neither as disconnected from the
atmosphere as soil under veg or canopy, nor as exposed as a large
area of exposed soil. Rather, it is subject to some wind attenuation
from the surrounding plants. Thus, compute as the area-weighted
average of "pure" understory and "pure" exposed conditions.

NOTE: can't average the resistances; must convert to conductances,
average the conductances, and then convert the average conductance
to a resistance.
*************************************************/
if (veg_var->vegcover < 1) {
/** If Ra_used (which pertains to under vegetation at this point)
is non-zero, use it to compute area-weighted average **/
if (Ra_used[0] > 0) {
/** aerodynamic conductance under vegetation **/
ga_veg = 1/Ra_used[0];
/** compute aerodynamic resistance over exposed soil (Ra_bare) **/
tmp_wind[0] = wind[0];
tmp_wind[1] = -999.; // unused
tmp_wind[2] = -999.; // unused
tmp_height = soil_con->rough/0.123; // eqn from calc_veg_roughness()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No hardcoded values - if this is a constant, it needs to be set in one of the header files.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, I'm going to replace the -999s with ERROR (which is defined as -999 in vicNl_def.h), both here and in full_energy(), which is where I copied these lines from.

I'm defining a few other constants too:
RATIO_DH_HEIGHT_VEG 0.67 - used in calc_veg_displacement() and calc_veg_height()
RATIO_RL_HEIGHT_VEG 0.123 - used both here and in calc_veg_roughness()
WINDH_SOIL 10.0

tmp_displacement[0] = calc_veg_displacement(tmp_height);
tmp_roughness[0] = soil_con->rough;
tmp_ref_height[0] = 10; // wind measurement height over bare soil
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No hardcoded values (see above)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, but release 4.2 has a bunch of hardcoded values throughout the code. For example, the 0.123 is also used as a hardcoded value in calc_veg_roughness(); 0.67 is used in calc_veg_displacement() and calc_veg_height(); the -999s are used in full_energy() (I took these lines from full_energy()). I can replace all of these with constants defined in vicNl_def.h, both in this function and in the others that I just listed. But does this really need to be done in 4.2, when constants are overhauled in 5.0?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm just going ahead and defining these constants in vicNl_def.h and replacing the numbers here with those constants.

Error = CalcAerodynamic(0,0,0,soil_con->snow_rough,soil_con->rough,0,Ra_bare,tmp_wind,tmp_displacement,tmp_ref_height,tmp_roughness);
Ra_bare[0] /= StabilityCorrection(tmp_ref_height[0], tmp_displacement[0], TMean, Tair, tmp_wind[0], tmp_roughness[0]);
/** if Ra_bare is non-zero, compute area-weighted average
aerodynamic conductance **/
if (Ra_bare[0] > 0) {
/** aerodynamic conductance over exposed soil **/
ga_bare = 1/Ra_bare[0];
/** area-weighted average aerodynamic conductance **/
ga_average = veg_var->vegcover*ga_veg + (1-veg_var->vegcover)*ga_bare;
/** aerodynamic resistance is inverse of conductance **/
Ra_used[0] = 1/ga_average;
}
/** else aerodynamic resistance is zero **/
else {
Ra_used[0] = 0;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be zero or infinite (or does it not matter because it means that the value is not used?)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It should be zero. Ra of 0 means infinite conductance. Perhaps some of the confusion here is due to the naming of the variables. Before I added this feature, Ra_used was the Ra of the entire veg tile. Since vegcover was assumed == 1.0 in the past, Ra_used was also the Ra of the vegetated portion (which was 100% of the tile). When I added this block, I intended to adjust Ra_used for the presence of exposed soil. In this case, the old Ra_used now initially only pertains to the vegetated portion (which is < 100% of the tile), and gets reassigned to the value pertaining to the entire tile later on.

So, what we have is:
If vegetated Ra is 0, then vegetated conductance (ga) is infinite, so total tile ga is also infinite, so Ra_used = 1/ga is 0.
If exposed soil Ra is 0, Ra_used also = 0.

What if I add another intermediate variable to make this clearer: Ra_veg. Ra_veg will be assigned = Ra_used at the beginning of this block. Then, ga_veg = 1/Ra_veg. ga_avg = weighted average as currently stated. Then, Ra_used = 1/ga_avg. This will hopefully make it clearer that Ra_used = 0 if either Ra_veg or Ra_bare are 0.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes - that would be much easier to understand - please go ahead.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, done.

}
}
/** else aerodynamic resistance is zero **/
else {
Ra_used[0] = 0;
}
}

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jhamman commented: Can you add some comments here. There are a few hard coded literals and a lot of moving pieces without any explanation.
@tbohn's reply: OK, will do.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Still needs those comments. The hardcoded values need to be moved to header files.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But I added several lines of comments... It needs more? Which lines in particular need comments?

/*************************************************
Compute Evapotranspiration if not snow covered

Should evapotranspiration be active when the
ground is only partially covered with snow????

Use Arno Evap if LAI is set to zero (e.g. no
winter crop planted).
Use Arno Evap in the exposed soil portion, and/or
if LAI is zero.
*************************************************/
if ( VEG && !SNOWING && veg_var->vegcover > 0 ) {
Evap = canopy_evap(layer, veg_var, TRUE,
Expand All @@ -709,20 +762,11 @@ double func_surf_energy_bal(double Ts, va_list ap)
transp[i] = layer[i].evap;
layer[i].evap = 0;
}
tmp_wind[0] = wind[0];
tmp_wind[1] = -999.;
tmp_wind[2] = -999.;
tmp_height = soil_con->rough/0.123;
tmp_displacement[0] = calc_veg_displacement(tmp_height);
tmp_roughness[0] = soil_con->rough;
tmp_ref_height[0] = 10;
Error = CalcAerodynamic(0,0,0,soil_con->snow_rough,soil_con->rough,0,Ra_bare,tmp_wind,tmp_displacement,tmp_ref_height,tmp_roughness);
Ra_bare[0] /= StabilityCorrection(tmp_ref_height[0], tmp_displacement[0], TMean, Tair, tmp_wind[0], tmp_roughness[0]);
Evap *= veg_var->vegcover;
Evap += (1-veg_var->vegcover)
* arno_evap(layer, surf_atten*NetBareRad, Tair, vpd,
depth[0], max_moist * depth[0] * 1000.,
elevation, b_infilt, Ra_bare[0], delta_t,
elevation, b_infilt, Ra_used[0], delta_t,
resid_moist[0], frost_fract);
for (i=0; i<options.Nlayer; i++) {
layer[i].evap = veg_var->vegcover*transp[i] + (1-veg_var->vegcover)*layer[i].evap;
Expand Down