Skip to content

Conversation

@samsrabin
Copy link
Contributor

@samsrabin samsrabin commented Jun 30, 2025

Description:

This PR adds the ability to track how much of a site's forest area is in each "bin" of distance from a forest edge. This will eventually be used for the escaped fire work (forest closer to edge will have a higher fire danger index), but for now it's purely diagnostic. I will also hide all this behind a feature flag, so none of the new variables will be allocated and no calculations will happen unless the user explicitly requests it. 2025-09-10: I ended up running into problems with not allocating the variables, so I reverted that bit in 4d3b7d4.

Includes changes from #1426; will look bigger than it really is until that is merged.

Corresponding CTSM PR is ESCOMP/CTSM#3340.

Detailed description

New parameter file parameters

The first step in this work was empirical analysis to develop a function for assigning forest to different edge distance bins (including what set of bins to use). I made a Python project to explore the possibilities and published it here, with lin_edgeareas_explore.py being the script to run if you want to generate the parameter sets and figures yourself. If you look through that script you can see various commented-out options I went through during the process; the v0.1.1 tag shows the exact setup needed to reproduce my final parameter set. (I haven't tested the different options with v0.1.1, so some or all of them may be broken.) Briefly, given land cover and forest edge bin area from Lin's remote sensing work, it fits (for each bin) the fraction of the gridcell's forest area in the bin (Y variable) to the total forested fraction of the gridcell. It tries the following fits, choosing the best for each bin based on lowest AIC value (bold indicates fit types that actually get chosen in v0.1.1):

  • Gaussian
  • Skewed Gaussian
  • Lognormal
  • Logistic
  • Exponential
  • Linear
  • Quadratic

It then saves a .cdl file skeleton with the resulting parameters, in addition to making a couple of figures illustrating the results.

This PR adds to the FATES parameter file:

  • Variable fates_forest_tree_fraction_threshold: Tree fraction above which patch is considered "forest." Set to 0.5; need to check whether this makes sense based on Lin's landcover data.
  • Dimension fates_edgeforest_bins: Number of edge forest bins, including "deep forest." All the following parameter file variables have this dimension.
  • Variable fates_edgeforest_bin_edges: Boundaries of forest edge bins (for each bin, include value closest to zero): 0, 30, 60, 120, 150, and 300.
  • Variables fates_edgeforest_gaussian_* (3): Parameters describing fits for bins with Gaussian fits.
  • Variables fates_edgeforest_lognormal_* (3): Parameters describing fits for bins with lognormal fits.
  • Variables fates_edgeforest_quadratic_* (3): Parameters describing fits for bins with quadratic fits.

Each edge forest bin has only one fit type, so the .cdl file has _ for that bin in parameters of the other types.

New code

This PR adds three modules to the FATES code. The first, main/FatesEcotypesMod.F90, is used to determine whether a patch is "forest" or not. I named it FatesEcotypesMod instead of something about forest because it could be extended in the future for other land cover classes. The is_patch_forest() function checks whether a patch is forest based on the fates_forest_tree_fraction_threshold parameter (see above) as well as an optional grass biomass threshold (above which a patch can't be forest). The latter is not currently used for anything except experimental outputs, as there is no corresponding variable on the parameter file.

The other two have to do with edge forest specifically. The edge forest parameters are read in main/FatesEdgeForestParamsMod.F90 to avoid junking up EdParamsMod. Those parameters are then used in main/FatesEdgeForestMod.F90, the main purpose of which is get_fraction_of_edgeforest_in_each_bin(). That subroutine loops through a site's patches from nearest to farthest from edge, assigning their area to the various edge forest bins. It is called via the public subroutine calculate_edgeforest_area() in CalculateTreeGrassAreaSite()(moved to FatesEdgeForestMod from main/EdTypesMod.F90, I think because of circular dependencies). I also added a call of it to canopy_summarization() (biogeochem/EDCanopyStructureMod.F90), but that's not actually called anywhere.

New outputs

  • FATES_AREA_PLANTS_AP and FATES_AREA_TREES_AP: Age-structured versions of existing variables; these might be useful for experimentation as this work continues.
  • FATES_IS_FOREST and FATES_IS_FOREST_AP: Fraction of the site (or each age bin in the site) that is forest according to the fates_forest_tree_fraction_threshold parameter (see above).
  • FATES_IS_FOREST_PCTX: Fraction of the site that is forest according to a tree cover threshold of X% instead of the fates_forest_tree_fraction_threshold parameter.
  • FATES_IS_FOREST_PCTX_0GRASS: Fraction of the site that is forest according to a tree cover threshold of X% and a grass biomass maximum of zero (although as noted above I currently have the grass biomass bit backwards).
  • FATES_FOREST_AREA_EB: Fraction of the site's forest area that is in each edge bin. Has new fates_levedge dimension.

New namelist parameters

If the HLM's use_edge_forest parameter is .false.:

  • Arrays are allocated but only for one bin (modeled after hlm_use_tree_damage).
  • Edge forest area is not calculated.
  • FATES_FOREST_AREA_EB output can be saved but will just have the initialization value.

Checking results

Bins are labeled by their minimum distance to edge (m). Adjusted/normalized Y values show the results after adjusting the raw fit values so that the sum of all curves at any X value is equal to 1.

Expectation from my Python notebook based on Lin's work (lin_edgeareas_explore.py here):
screenshot_4115

✅ Functional testing:
edge_forest_plot_frac_in_every_bin 1

✅ Results at end of a four-year test (comparing to dashed lines from above):
screenshot_4117

Collaborators:

@mpaiao, @ckoven, @menglinmet

Expectation of Answer Changes:

None

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

Integrator

  • FATES PASS/FAIL regression tests were run
  • Evaluation of test results for answer changes was performed and results provided
  • FATES-CLM6 Code Freeze: satellite phenology regression tests are b4b

If satellite phenology regressions are not b4b, please hold merge and notify the FATES development team.

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 and others added 30 commits October 16, 2024 11:58
Add FatesEcotypesMod with is_patch_forest() function.

Updates new Patch logical is_forest based on tree fraction, immediately after total_tree_area is updated.

Set up for adding IS_FOREST outputs based on grass biomass.

(cherry picked from commit 7ee78b4dc44d31eecbd9e98cde264392131bf93d)
Add parameter file fates_params_shuman24.cdl.

As fates_params_default.cdl, but with only:
- broadleaf_evergreen_tropical_tree
- broadleaf_hydrodecid_tropical_tree
- c4_grass

Change shuman24 PFT names and some values.

More changes based on Shuman et al. (2024) Table 2.

Add params from fates_params_sci.1.45.0_api.15.0.0.cdl, for ref.

Add Jackie's parameter file, for reference.

Add fates_params_shuman24.20240819.cdl (just extracted PFTs).

Apply Jackie's parameter file changes.

To fates_params_shuman24.20240819.cdl. Minimal changes that match her CHANGES, not her VALUES.

(cherry picked from commit 6151cd9bd255e8c82ea572bde39d6742b9ae09c7)
- FATES_AREA_PLANTS_AP
- FATES_AREA_TREES_AP
- FATES_IS_FOREST
- FATES_IS_FOREST_AP
- FATES_IS_FOREST_PCT10
- FATES_IS_FOREST_PCT10_0GRASS
- FATES_IS_FOREST_PCT25
- FATES_IS_FOREST_PCT25_0GRASS
- FATES_IS_FOREST_PCT50
- FATES_IS_FOREST_PCT50_0GRASS
- FATES_IS_FOREST_PCT75
- FATES_IS_FOREST_PCT75_0GRASS
- FATES_IS_FOREST_PCT90
- FATES_IS_FOREST_PCT90_0GRASS
Uses new dimension, site x fates_levedgebin.
Global 4x5 test of this change (on a different branch): Results in absolute differences in bin area of < 0.0001 m2 in any given gridcell. Total sum of abs. differences < 0.22 m2.
It's actually now defined in EDParamsMod so it can be read from the parameter file. This instance was left over from the move to the parameter file.
This should be impossible because they all use the same netCDF dimension.
New sapwood, agb and leaf allometries for Grasses, by Xiulin Gao.
Patch level memory structures were refactor to be dynamic, paving the way for higher nclmax.

# Conflicts:
#	biogeochem/FatesPatchMod.F90
Fire-weather refactor (non-bfb)

This moves the effective wind speed and associated methods into the
site-level fire weather class.

When merging into edge-area-202410: Had added "currentPatch%is_forest = is_patch_forest()" and "call calculate_edge_area()" to subroutine wind_effect. These are now in subroutine CalculateTreeGrassAreaSite.
Changes the patch insertion method

This update changes the patch ordering method to account for landuse and
no comp pft label as well as patch age.
@samsrabin samsrabin requested a review from adrifoster July 23, 2025 17:07
@samsrabin
Copy link
Contributor Author

@mpaiao @ckoven @adrifoster This is ready to review!

Adrianna, I added you because I ended up needing to make some changes to the test infrastructure. Let me know if you'd like me to walk you through them.

@rgknox rgknox self-requested a review July 28, 2025 19:08
@samsrabin samsrabin moved this from Finding Reviewers to Under Review in FATES Pull Request Planning and Status Aug 25, 2025
Copy link
Contributor

@adrifoster adrifoster left a comment

Choose a reason for hiding this comment

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

Looks nice! I like all the unit and functional tests but I'm not sure about doing I/O for the unit tests. Generally we have maintained that the unit tests should do no I/O. So if you need access to a parameter somehow then maybe something needs to get refactored so that we don't have to use the actual parameter for the unit test?

Copy link
Contributor

@ckoven ckoven left a comment

Choose a reason for hiding this comment

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

@samsrabin this looks great to me, thanks for doing it! I have a couple comments, all minor stuff.

Comment on lines 74 to 76
if (is_patch_forest .and. present(grass_biomass_threshold)) then
is_patch_forest = .not. does_patch_have_grass_bmthresh(patchptr, grass_biomass_threshold)
end if
Copy link
Contributor

Choose a reason for hiding this comment

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

Would be useful to explain the logic here around the use of the grass biomass threshold for forest definition. Is there a reference to cite for this definition of forest?

Copy link
Contributor

Choose a reason for hiding this comment

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

Though I am also noticing that the grass_biomass_threshold isn't actually used in the current logic, other than in some diagnostic history variables and unit tests.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yeah, grass biomass was just an idea I had for another constraint. If you still want me to explain the logic, let me know; otherwise please resolve this comment.

f = 0
p = 0
index_forestpatches_to_allpatches(:) = 0
currentPatch => site%oldest_patch
Copy link
Contributor

Choose a reason for hiding this comment

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

there should be probably be a pair of warnings both here and in assign_patches_to_bins that the order of patch looping (oldest to youngest) needs to match in the two cases.

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 don't think that's true, actually. The ordering of the loop in AssignPatchesToBins() shouldn't matter, because each forest patch should have a unique value in index_forestpatches_to_allpatches (as tested in indexx_lowTie(), indexx_highTie(), and indexx_all_equal().

@samsrabin
Copy link
Contributor Author

@adrifoster I've changed this to hard-code the necessary parameter values, so unit tests no longer read the parameter file.

subroutine FatesParamsInitForFactory()
! Initialize some parameters that are needed for unit-testing factories

allocate(ED_val_history_ageclass_bin_edges(7))
Copy link
Contributor

Choose a reason for hiding this comment

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

@samsrabin these are also still allocated here: https://github.com/samsrabin/fates/blob/edge-area-202410/main/EDParamsMod.F90#L119

Do we need to remove those statements?

double fates_dev_arbitrary_pft(fates_pft) ;
fates_dev_arbitrary_pft:units = "unknown" ;
fates_dev_arbitrary_pft:long_name = "Unassociated pft dimensioned free parameter that developers can use for testing arbitrary new hypotheses" ;
double fates_edgeforest_bin_edges(fates_edgeforest_bins) ;
Copy link
Contributor

Choose a reason for hiding this comment

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

group name-spacing (edgeforest), nice

Copy link
Contributor

Choose a reason for hiding this comment

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

extra points for inserting it in alphabetical order!


n_forest_patches = 0
area_forest_patches = 0._r8
area_site = 0._r8
Copy link
Contributor

Choose a reason for hiding this comment

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

@samsrabin , since there is no filtering of the the patch loop (ie it traverses all patches from youngest to oldest), and area_site is not inside any logic filters, shouldn't area_site always equal 10k (ie "AREA")? Was this not expected or was this a check?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

output: diagnostics parameters: new Pertaining to adding new parameters to the parameter file software: API Pertaining to specific API updates with any host land model status: Needs Review type: enhancement

Projects

Status: Under Review

Development

Successfully merging this pull request may close these issues.

5 participants