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.

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

Open
wants to merge 24 commits into
base: main
Choose a base branch
from

Conversation

samsrabin
Copy link
Contributor

@samsrabin samsrabin commented Sep 16, 2024

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:

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:

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

  • The in-code documentation has been updated with descriptive comments
  • The documentation has been assessed to determine if updates are necessary: Yes

Integrator

  • FATES PASS/FAIL regression tests were run
  • Evaluation of test results for answer changes was performed and results provided

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:

@samsrabin samsrabin added bug - software engineering PR status: Not Ready The author is signaling that this PR is a work in progress and not ready for integration. history output Pertaining to FATES history output variables refactor Restructures code without changing functionality labels Sep 16, 2024
@samsrabin samsrabin self-assigned this Sep 16, 2024
@samsrabin samsrabin linked an issue Sep 16, 2024 that may be closed by this pull request
@samsrabin samsrabin force-pushed the refactor-history-2 branch 2 times, most recently from fa87c89 to 1d1c44f Compare October 10, 2024 22:50
@samsrabin samsrabin removed the refactor Restructures code without changing functionality label Oct 11, 2024
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
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
@samsrabin samsrabin marked this pull request as ready for review October 13, 2024 21:10
@samsrabin samsrabin changed the title [WIP] Refactor and standardize per-ageclass history Refactor and standardize per-ageclass history Oct 13, 2024
@samsrabin samsrabin added PR status: Needs Review PR status: Not Ready The author is signaling that this PR is a work in progress and not ready for integration. and removed PR status: Not Ready The author is signaling that this PR is a work in progress and not ready for integration. PR status: Needs Review labels Oct 13, 2024
@@ -358,7 +358,7 @@ module FatesHistoryInterfaceMod

integer :: ih_primaryland_fusion_error_si

integer :: ih_area_si_landuse
integer :: ih_fracarea_si_landuse
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 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.

Copy link
Contributor

Choose a reason for hiding this comment

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

thanks

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've now changed other such variables throughout this module to match.

@samsrabin
Copy link
Contributor Author

samsrabin commented Oct 21, 2024

  • Rename and add comments to update_history_dyn1() and update_history_dyn2() to indicate they're intended for site-level vs. multiplexed outputs, respectively.

real(r8) :: mort_natural
real(r8) :: mort_logging
real(r8) :: mort_sum

Copy link
Contributor

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

Copy link
Contributor Author

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.

Copy link
Contributor Author

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.

Copy link
Contributor

@rgknox rgknox Nov 18, 2024

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

@samsrabin
Copy link
Contributor Author

@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?

@samsrabin samsrabin requested a review from rgknox October 29, 2024 17:13
Used in history outputs to get biomass of each organ, plus alive and total biomass.
Copy link
Contributor

@glemieux glemieux left a 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.

Comment on lines +3388 to +3391
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) &
Copy link
Contributor

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
Copy link
Contributor

@glemieux glemieux Nov 9, 2024

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

Copy link
Contributor

@glemieux glemieux left a 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.

Comment on lines +4748 to +4759
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
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
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

Copy link
Contributor

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.

Comment on lines +5709 to +5714
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
Copy link
Contributor

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:

! 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)

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?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
Status: Under Review
Development

Successfully merging this pull request may close these issues.

Incorrect (?) calculation of per-age-class outputs
4 participants