-
Notifications
You must be signed in to change notification settings - Fork 92
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.
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
Refactor and standardize per-ageclass history #1252
base: main
Are you sure you want to change the base?
Conversation
fa87c89
to
1d1c44f
Compare
The goal here is to simplify things and improve my understanding before more substantive refactoring and fixing. Changes: - Do everything within one patch loop, rather than looping through to get numerators and denominators separately, then looping again to do the division. - Use new "weight" variable instead of repeating (e.g.) "cpatch*area * AREA_INV". - Replace "hio_area_si_age(io_si, ipa2)*AREA" with "sites(s)%area_by_age(ipa2)" - Replace "sites(s)%area_by_age(cpatch%age_class)" with new "age_class_area" variable. Should be identical within roundoff precision. Affected variables: - FATES_BURNFRAC_AP - FATES_CANOPYAREA_AP - FATES_FIRE_INTENSITY_BURNFRAC_AP - FATES_FUEL_AMOUNT_AP - FATES_GPP_AP - FATES_LAI_AP - FATES_LBLAYER_COND_AP - FATES_NCL_AP - FATES_NPP_AP - FATES_SCORCH_HEIGHT_APPF - FATES_SECONDAREA_ANTHRODIST_AP - FATES_SECONDAREA_DIST_AP - FATES_STOMATAL_COND_AP - FATES_VEGC_AP - FATES_ZSTAR_AP Also renames various internal FatesHistoryInterfaceMod variables from "area" to "fracarea" for clarity (units are m2/m2, not m2).
To avoid similar mistakes in the future, adds function SumMortForHistory to fates_cohort_type.
These are useful for checking whether the per-ageclass versions are normalized correctly. * FATES_CANOPYAREA * FATES_NCL * FATES_PATCHAREA * FATES_SCORCH_HEIGHT_PF * FATES_SECONDAREA_ANTHRODIST * FATES_SECONDAREA_DIST * FATES_ZSTAR
Now dividing by total site area instead of age-class area: - FATES_GPP_AP - FATES_NPP_AP - FATES_LAI_AP - FATES_NCL_AP - FATES_SCORCH_HEIGHT_APPF Now dividing by site-wide canopy area instead of age-class canopy area: - FATES_LBLAYER_COND_AP - FATES_STOMATAL_COND_AP
1d1c44f
to
7e1bf60
Compare
Tell the user how to get the value they probably want---the value on each age-class---rather than what each age-class contributes to the cross-ageclass area-weighted mean.
- Use descriptive weight names - Remove unused variable ipa - Organize into sections - Add a TODO
- Use descriptive weight names - Remove unused variable ipa - Fix indentation
d249965
to
8356d40
Compare
@@ -358,7 +358,7 @@ module FatesHistoryInterfaceMod | |||
|
|||
integer :: ih_primaryland_fusion_error_si | |||
|
|||
integer :: ih_area_si_landuse | |||
integer :: ih_fracarea_si_landuse |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I renamed a bunch of area
variables to be fracarea
instead, as "area" implies m2 but most of these are fractional. There are still some area
variables in other subroutines, like ih_area_plant_si
, that should probably also be changed for consistency.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
thanks
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I've now changed other such variables throughout this module to match.
|
real(r8) :: mort_natural | ||
real(r8) :: mort_logging | ||
real(r8) :: mort_sum | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Lets double check the units on the natural versus the logging mortality, I seem to remember they were different, such as their usage here:
https://github.com/NGEET/fates/blob/main/biogeochem/EDMortalityFunctionsMod.F90#L347-L354
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yep, that's what the per_year
option is for—depending on its value, one will be converted into units of the other. I'll add units to comments in revision.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No, the problem is that one is in # individuals and the other is fraction.
- Fix this.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@samsrabin here is an example of how to get total carbon biomass from the cohort, if you want to add it into this routine. This will help you convert from number of individuals to kgC:
https://github.com/NGEET/fates/blob/main/main/FatesHistoryInterfaceMod.F90#L3895-L3903
I looked at how this sum mortality was being used, and I had the units wrong, it looks like the function is being called with the expectation that the output units are not in carbon but the abundance rate.
Part of the confusion is related to our incorrect labeling of the rates on the cohort structures. Here is where we apply mortality to actually change number density, dndt is a fraction per time (because that fraction is multiplying by cohort%n:
https://github.com/NGEET/fates/blob/main/main/EDMainMod.F90#L804
https://github.com/NGEET/fates/blob/main/biogeochem/EDMortalityFunctionsMod.F90#L352-L354
@rgknox @ckoven This is ready for re-review. The only remaining thing I need to ask about (marked by TODOs in the code) is that the reproductive biomass seems to not be included in total biomass in some parts of the history code, whereas it is included in others. Do you think there's a reason for that? |
Used in history outputs to get biomass of each organ, plus alive and total biomass.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Review part 1/2. Really love the work to break things up more here. There was one minor area of potential clean up that I wanted to have double checked.
iscag_anthrodist = get_age_class_index(cpatch%age_since_anthro_disturbance) | ||
|
||
hio_agesince_anthrodist_si_age(io_si,ageclass_since_anthrodist) = & | ||
hio_agesince_anthrodist_si_age(io_si,ageclass_since_anthrodist) & | ||
hio_agesince_anthrodist_si(io_si) = & | ||
hio_agesince_anthrodist_si(io_si) & |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Since we're removing the second dimension here can we remove this call to get_age_class_index
?
real(r8) :: repro_m ! Total reproductive mass (on plant) "" | ||
real(r8) :: patch_area_div_site_area ! Weighting based on patch area relative to site area | ||
real(r8) :: patch_canarea_div_site_area ! Weighting based on patch canopy area relative to site area | ||
real(r8) :: cohort_n_div_site_area ! Weighting based on cohort density relative to site area |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
FYI, I created a tasklist to track these sort of simple improvements for future PRs: #1282
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Review 2/2. I have one question on the hio_c_stomata_si_age
calculation. I also have one minor suggested change to condense loops, but it can be delayed if necessary to move forward. Great work.
do ft = 1,numpft | ||
hio_scorch_height_si_pft(io_si,ft) = hio_scorch_height_si_pft(io_si,ft) + & | ||
cpatch%Scorch_ht(ft) * patch_area_div_site_area | ||
end do | ||
|
||
hio_ncl_si_age(io_si,cpatch%age_class) = hio_ncl_si_age(io_si,cpatch%age_class) & | ||
+ cpatch%ncl_p * patch_area_div_site_area | ||
do ft = 1,numpft | ||
iagepft = get_agepft_class_index(cpatch%age,ft) | ||
hio_scorch_height_si_agepft(io_si,iagepft) = hio_scorch_height_si_agepft(io_si,iagepft) + & | ||
cpatch%Scorch_ht(ft) * patch_area_div_site_area | ||
end do |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
do ft = 1,numpft | |
hio_scorch_height_si_pft(io_si,ft) = hio_scorch_height_si_pft(io_si,ft) + & | |
cpatch%Scorch_ht(ft) * patch_area_div_site_area | |
end do | |
hio_ncl_si_age(io_si,cpatch%age_class) = hio_ncl_si_age(io_si,cpatch%age_class) & | |
+ cpatch%ncl_p * patch_area_div_site_area | |
do ft = 1,numpft | |
iagepft = get_agepft_class_index(cpatch%age,ft) | |
hio_scorch_height_si_agepft(io_si,iagepft) = hio_scorch_height_si_agepft(io_si,iagepft) + & | |
cpatch%Scorch_ht(ft) * patch_area_div_site_area | |
end do | |
do ft = 1,numpft | |
iagepft = get_agepft_class_index(cpatch%age,ft) | |
hio_scorch_height_si_pft(io_si,ft) = hio_scorch_height_si_pft(io_si,ft) + & | |
cpatch%Scorch_ht(ft) * patch_area_div_site_area | |
hio_scorch_height_si_agepft(io_si,iagepft) = hio_scorch_height_si_agepft(io_si,iagepft) + & | |
cpatch%Scorch_ht(ft) * patch_area_div_site_area | |
end do | |
hio_ncl_si_age(io_si,cpatch%age_class) = hio_ncl_si_age(io_si,cpatch%age_class) & | |
+ cpatch%ncl_p * patch_area_div_site_area |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Suggestion to condense numpft
loops. Not critical.
if (site_canopy_area .gt. nearzero) then | ||
patch_canarea_div_site_canarea = cpatch%total_canopy_area / site_canopy_area | ||
hio_c_stomata_si_age(io_si,cpatch%age_class) = & | ||
hio_c_stomata_si_age(io_si,cpatch%age_class) + & | ||
cpatch%c_stomata * mol_per_umol & | ||
* patch_canarea_div_site_canarea |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Apologies if this was discussed elsewhere and I missed it. The previous version of this calculation was dividing through the sum of the total_canopy_area
just for the given age class:
fates/main/FatesHistoryInterfaceMod.F90
Lines 5603 to 5606 in 1be3962
! Normalize resistance diagnostics | |
if (canopy_area_by_age(ipa2) .gt. nearzero) then | |
hio_c_stomata_si_age(io_si,ipa2) = & | |
hio_c_stomata_si_age(io_si,ipa2) / canopy_area_by_age(ipa2) |
fates/main/FatesHistoryInterfaceMod.F90
Lines 5297 to 5298 in 1be3962
canopy_area_by_age(cpatch%age_class) = & | |
canopy_area_by_age(cpatch%age_class) + cpatch%total_canopy_area |
This change appears to introduce dividing the by the site_canopy_area
. Do I have that correct and is that the intent?
Description:
As noted in Issue #1205, some per-ageclass history variables are weighted based on a patch's area (or a cohort's density) relative to age class area, whereas others are weighted based on that relative to total site area. To avoid confusion, this should be standardized. Discussion led us to decide on standardizing to the latter to avoid order-of-operations issues.
During development, I discovered and fixed two other bugs:
FATES_VEGC_APPF
was incorrectly not being updated for new cohorts (contributes to Some biomass outputs only updated for non-new cohorts #1261).FATES_MORTALITY_CANOPY_SZAP
andFATES_MORTALITY_USTORY_SZAP
were incorrectly excluding termination and logging mortality.Supersedes PR #1225, where I had made the opposite change (standardizing all to use per-ageclass area).
Collaborators:
@ckoven, @rgknox
Expectation of Answer Changes:
These variables are changed to be weighted relative to site area (rather than age-class area):
FATES_GPP_AP
FATES_LAI_AP
FATES_NCL_AP
FATES_NPP_AP
FATES_SCORCH_HEIGHT_APPF
These variables are changed to be weighted relative to site-wide canopy area (rather than age-class canopy area):
FATES_LBLAYER_COND_AP
FATES_STOMATAL_COND_AP
These variables are affected by bug fixes:
FATES_MORTALITY_CANOPY_SZAP
(logging and termination mortality now included)FATES_MORTALITY_USTORY_SZAP
(logging and termination mortality now included)FATES_VEGC_APPF
(see Some biomass outputs only updated for non-new cohorts #1261)Other per-ageclass variables see roundoff-level changes as a result of refactoring (absolute RMSE < 1e-12, relative RMSE < 1e-13). In my testing (
SMS_Lm49.f10_f10_mg37.I2000Clm60Fates.derecho_intel.clm-FatesColdAllVarsMonthly
), the following were affected:FATES_BURNFRAC_AP
FATES_DDBH_CANOPY_SZAP
FATES_DDBH_USTORY_SZAP
FATES_FIRE_INTENSITY_BURNFRAC_AP
FATES_FUEL_AMOUNT_AP
FATES_FUEL_AMOUNT_APFC
FATES_MORTALITY_CFLUX_CANOPY
FATES_MORTALITY_CFLUX_PF
FATES_MORTALITY_CFLUX_USTORY
FATES_NPLANT_CANOPY_SZAP
FATES_NPLANT_SZAP
FATES_NPLANT_SZAPPF
FATES_NPLANT_USTORY_SZAP
FATES_NPP_APPF
FATES_SCORCH_HEIGHT_PF
FATES_VEGC_AP
Answers from that test (at 1aff666) compared to baseline can be found here. Note that
FATES_NPLANT_SZAPPF
was not tested to see if it sums correctly across the age-class axis because I don't have a way to handle triplexed variables; however, it was only affected at roundoff level by this work.Checklist
All checklist items must be checked to enable merging this pull request:
Contributor
Integrator
Documentation
Test Results:
CTSM (or) E3SM (specify which) test hash-tag:
CTSM (or) E3SM (specify which) baseline hash-tag:
FATES baseline hash-tag:
Test Output: