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

[samsrabin]

There are (at least) two different ways that per-age-class (_AP) outputs are be calculated in update_history_dyn().

When I look at such a variable, I expect the value in each age class to be the mean across all patches in that age class, weighted by patch area. This is how it's done for LAI (FATES_LAI_AP), requiring two steps. First, in patchloop:

hio_lai_si_age(io_si,cpatch%age_class) = hio_lai_si_age(io_si,cpatch%age_class) &
   + sum(cpatch%tlai_profile(:,:,:)) * cpatch%area

Then you need to divide by the total age class area. Since that's accumulated during patchloop, this happens after that loop:

! divide so-far-just-summed but to-be-averaged patch-age-class
! variables by patch-age-class area to get mean values
do ipa2 = 1, nlevage
   if (hio_area_si_age(io_si, ipa2) .gt. nearzero) then
         hio_lai_si_age(io_si, ipa2) = hio_lai_si_age(io_si, ipa2) / (hio_area_si_age(io_si, ipa2)*AREA)
         ...
   else
      hio_lai_si_age(io_si, ipa2) = 0._r8
      ...
   endif
end do

However, for FATES_FUEL_AMOUNT_AP, it's all done in one step within patchloop:

 ! Fuel sum [kg/m2]
 hio_fire_sum_fuel_si_age(io_si, cpatch%age_class) = hio_fire_sum_fuel_si_age(io_si, cpatch%age_class) +  &
    cpatch%sum_fuel * cpatch%area * AREA_INV

There's never a division by the total age class area, so this doesn't give the average fuel load in each patch age class (kg/m2) as the units of FATES_FUEL_AMOUNT_AP state. Instead, it gives each age class's contribution to the total (not per-m2) fuel load. So if you do an unweighted sum across all age classes, you'll get the total fuel load. And it turns out, FATES_FUEL_AMOUNT also gives the total fuel load, in contrast to its units of kg/m2. :

As expected, this isn't the case for LAI:

I don't know if this is intended behavior for FATES_FUEL_AMOUNT_AP, but as I mentioned it's certainly not what I would expect as a user of the output. Of outputs with fates_levage dimension (I haven't checked outputs with age class crossed by something else, like PFT), the following outputs (in addition to FATES_FUEL_AMOUNT_AP) appear to follow the same behavior:

• FATES_ZSTAR_AP
• FATES_SECONDAREA_ANTHRODIST_AP
• FATES_SECONDAREA_DIST_AP
• FATES_BURNFRAC_AP
• FATES_CANOPYAREA_AP
• FATES_VEGC_AP

[samsrabin]

I'm realizing now that my code doesn't have the fix to FATES_LAI_AP that Ryan made as part of #1174. Will add that, redo my run, and see if the LAI map looks any different.

[samsrabin]

Progress on this is blocked until #1207 is resolved.

[samsrabin]

Now that I know I can work around #1207, I can proceed on this.

[samsrabin]

This is still the case as of tag sci.1.76.4_api.35.1.0 (June 11).

[samsrabin]

Actually... I think #1174 didn't fully fix the LAI averaging problem. It looks fine after 12 months, which I think is why I told @rgknox that #1174 fixed things, but it turns out there's only one age class at that point. After 24 months (3 age classes):

That second map—the cross-ageclass mean of LAI with weights from FATES_PATCHAREA_AP—should be all zero. It looks similar with weights from FATES_CANOPYAREA_AP.

[samsrabin]

It turns out that FATES_FUEL_AMOUNT(_AP) is a bad example of the bad behavior, because it's always zero for the test runs I've been doing. FATES_VEGC(_AP) is a good example, though—the unweighted cross-ageclass sum of FATES_VEGC_AP is equal within 1e-6 to FATES_VEGC:

Note that this actually represents a third way (second wrong way) of doing things.

FATES_VEGC:

this%hvars(ih_totvegc_si)%r81d(io_si) =                      &
     this%hvars(ih_totvegc_si)%r81d(io_si)+ ccohort%n *        &
     total_m / m2_per_ha

FATES_VEGC_AP:

hio_biomass_si_age(io_si,cpatch%age_class) = hio_biomass_si_age(io_si,cpatch%age_class) &
     + total_m * ccohort%n * AREA_INV

/ m2_per_ha == * AREA_INV == * 1e-4, so it turns out that these are doing the exact same thing—no weights applied. I think here that means that FATES_VEGC is wrong but FATES_VEGC_AP is right.

[ckoven]

@samsrabin my hunch here is that this gets to the issue with averaging ratios rather than taking the ratio of averages? So properties that are conserved, like biomass, behave well when averaging over time, in the sense that the sum of the age-resolved variable equals the total-ecosystem variable. This isn't the case for variables like LAI, that have a ratio in them where both the numerator (leaf area) and denominator (patch area) change in time, and so the ratio isn't conserved. So perhaps a solution is to get rid of a patch-age-resolved LAI variable entirely, and replace it with a patch-age-resolved leaf area variable, which a user could use to calculate LAI by dividing it by FATES_PATCHAREA_AP?

[rgknox]

@samsrabin , sorry if this is a dumb question, but are you weighting you sum of the _AP variables by patch area? It seems to do this comparison, you would need that, looks like it would be FATES_PATCHAREA_AP.

[rgknox]

Actually... I think #1174 didn't fully fix the LAI averaging problem. It looks fine after 12 months, which I think is why I told @rgknox that #1174 fixed things, but it turns out there's only one age class at that point. After 24 months (3 age classes):

That second map—the cross-ageclass mean of LAI with weights from FATES_PATCHAREA_AP—should be all zero. It looks similar with weights from FATES_CANOPYAREA_AP.

oh, nevermind, I see you do use it here

[samsrabin]

@samsrabin my hunch here is that this gets to the issue with averaging ratios rather than taking the ratio of averages? So properties that are conserved, like biomass, behave well when averaging over time, in the sense that the sum of the age-resolved variable equals the total-ecosystem variable. This isn't the case for variables like LAI, that have a ratio in them where both the numerator (leaf area) and denominator (patch area) change in time, and so the ratio isn't conserved. So perhaps a solution is to get rid of a patch-age-resolved LAI variable entirely, and replace it with a patch-age-resolved leaf area variable, which a user could use to calculate LAI by dividing it by FATES_PATCHAREA_AP?

@ckoven I think you're right, thanks! And your solution sounds smart to me.