Refactor 20-year running means of crop GDD accumulation

[samsrabin]

Description of changes

Three variables track the 20-year running mean of GDD accumulation (base temperatures 0, 8, and 10°C) during the "growing season" (April through September in the Northern Hemisphere, October through March in the Southern Hemisphere). This PR refactors those to use accumulMod, resolving overly-strong weighting of the first few years after a crop becomes active.

Specific notes

Contributors other than yourself, if any: @billsacks

CTSM Issues Fixed (include github issue #):

• Fixes Crop GDD accumulation gives large weight to the first years #75.

Are answers expected to change (and if so in what way)? Yes. GDD020, GDD820, and GDD1020 will differ, most strongly in the first few years after a crop becomes active. This will have downstream effects on lots of variables, since those are used in determining sowing date and maturity requirements.

Any User Interface Changes (namelist or namelist defaults changes)? No.

Testing performed, if any: A long crop test, ERS_Ly5_P144x1.f10_f10_mg37.IHistClm51BgcCrop.cheyenne_intel.clm-cropMonthOutput, runs completely but—as expected—shows answer changes.

[samsrabin]

I'd like to do a long-ish test run—say, 1960 through 2014—to test my hypothesis that this change will matter less the longer a stand is active. Because this would be easier with all the changes from #1863, I'm marking this as blocked:dependency. However, I will also mark it as "ready for review" because I want the code to be in good shape before spending the CPU-hours on that test.

[samsrabin]

Testing confirms that the most extreme differences do indeed decrease over time:

[billsacks]

Thanks a lot for working on this @samsrabin ! I haven't given this a careful review at this point - just skimmed through it. And I probably won't have time to get my head into it deeply enough to give it a careful review in the next couple of weeks. So if you want faster turnaround of this, it should probably be assigned to someone else, like @slevis-lmwg . For now I just have one comment / question - see below.

[billsacks]

We normally initialize history fields to spval so that grid cell averages will come out correctly for fields that don't have valid values over some patches - e.g., non-vegetated (or non-crop) patches in this case. Can you comment on why you made this change? And have you confirmed that it works to output grid cell averages of these fields with this change to NaN?

[samsrabin]

Good point! I think I was trying to stick as closely as possible to how things were done for some other variable that used accumulMod, but of course the history field is different from the accumulation field, so I don't need to worry about matching the former.

• Change back to spval.

[slevis-lmwg]

@samsrabin and I discussed this PR today.
@samsrabin will post a summary of our conversation.

[samsrabin]

@slevis-lmwg and I spoke about this today, and we realized that accumulMod's update_accum_field_runaccum() actually uses almost the same algorithm for "running mean." This prompted two questions.

Why did I get non-roundoff answer differences?

Possible, non-mutually-exclusive factors we discussed:

• Something funky happening when restarting from a version that didn't have the accumulMod information. A cold-start experiment showing roundoff-level differences would point to this. But I doubt it, because in my experiment I artificially reset the variables upon restart.
• Something different about when the old vs. new versions get reset. E.g., if some crop's area goes to zero, one is reset while the other isn't. I would think this to be a small fraction of gridcells in any given year, though.

Something I realized after our conversation, which is definitely contributing:

• The old method assigns a weight of 19/20 to the previous value regardless of how long the patch has been "active." In contrast, accumulMod correctly uses ½, ⅔, ¾, …, 18/19, 19/20, 19/20, …
• Test this by amending the old method to do the same.

Is this work worth it?

Thinking that there was no difference between the old and new algorithm, @slevis-lmwg and I had thought it might be more effort than it was worth to get down to roundoff-level differences, because this would essentially just be a code cleanup PR. A true running mean would require that we retain the last 20 years' values for every crop patch, which would presumably be too memory-intensive.

However, since it looks like there is actually an error in the old method that accumulMod fixes, I think I'd like to proceed. At the very least, I want to do the test mentioned above (fix the bug in the old version) to see if that alone gives us roundoff-level differences. If so, this PR would (I think) be ready to come in. If not, maybe we'd decide to put it on ice in favor of other priorities.

Alternatively, I could do a smaller PR just fixing the bug while still using the old code.

[samsrabin]

I should add, @slevis-lmwg and I also realized there's another bug that the PR doesn't yet address.

As I mentioned in the PR description, these variables

track the 20-year running mean of GDD accumulation (base temperatures 0, 8, and 10°C) during the "growing season" (April through September in the Northern Hemisphere, October through March in the Southern Hemisphere).

Since they're updated (and each year's tracking variable is reset to 0) at the end of the calendar year, the values for a gridcell in the Southern Hemisphere are approximately half what they should be in the first year! Subsequent years are a bit weird as well, since they're the sum of "January through March in the previous growing season" and "October through December in the current growing season."

Filed as issue #2121.

[samsrabin]

Testing a pure fix of the 19/20 bug unfortunately did give differences from 255649b:

So there is actually a difference, beyond just the 19/20 thing, in what's happening with the old method vs. the accumulMod method. Not sure if it's worth trying to track those down, especially considering how small they are relative to the impact of the bugfix itself:

For an illustration of what's happening with a representative patch:

As you can see, the "fixed" method (and thus also accumulMod) has the advantage of the GDD020 value being reasonable soon after patch establishment. The "buggy" method only approaches a reasonable value over time.

[samsrabin]

I don't think the Southern Hemisphere issue is actually worth fixing. See this comment on #2121, which I've now closed.

[samsrabin]

@slevis-lmwg, would it be okay to add some discussion of this—specifically, whether it's worth trying to figure out the causes of the differences shown in the top figure here—to our tillage re-review tomorrow?

[samsrabin]

@slevis-lmwg, I don't think we ended up discussing this as I suggested in my last comment. Could we schedule a call?

[samsrabin]

@slevis-lmwg raised the idea that this might be a roundoff-level difference that accumulates into something seemingly significant. I'll plan to do a short test with per-timestep outputs to check.

[slevis-lmwg]

@samsrabin and I also agreed that this PR's changes are worth coming in regardless:

• roundoff is even better
• if not roundoff, the difference is still small and it's worth moving these calculations to our existing accumulMod tools rather than recreating accumulation algorithms locally

[samsrabin]

I'm not sure if this qualifies as roundoff, but:

nccmp -dfs -N -c 1 -x date_written,time_written /glade/scratch/samrabin/archive/refactor-crop-gdd-runmeans-dev.20231108v2/lnd/hist/refactor-crop-gdd-runmeans-dev.20231108v2.clm2.h0.1902-01-01-00000.nc /glade/scratch/samrabin/archive/fix-gdd20-weighting-test.20231108v2/lnd/hist/fix-gdd20-weighting-test.20231108v2.clm2.h0.1902-01-01-00000.nc

(i.e., looking at the first pft-timestep that differs between the old and new output files) gives

DIFFER : VARIABLE : GDD0 : POSITION : [2,66518] : VALUES : 0.132446 <> 0.132196
DIFFER : VARIABLE : GDD10 : POSITION : [2,66597] : VALUES : 0.00671405 <> 0.00662143
DIFFER : VARIABLE : GDD8 : POSITION : [2,66597] : VALUES : 0.0483807 <> 0.0482881

Those are differences of about -0.19%, -1.38%, and -0.19%, respectively.

What do you think, @slevis-lmwg?

(Outputs preserved for posterity in /glade/work/samrabin/fix-gdd20-weighting-test/.)

[slevis-lmwg]

You mentioned "first pft-timestep that differs" which makes me wonder how long it takes for this to occur. Otherwise, these diffs do seem small, though I'm more used to looking at the cprnc files for absolute and normalized diffs from a baseline.

[samsrabin]

This is the second timestep of the first year I expected any differences.

So is this ready to come in, then?

[slevis-lmwg]

Based on this earlier post, I think so.

[ekluzek]

@slevis-lmwg this is something that @samsrabin would like to bring in right away. Could you do a code review here?

[ekluzek]

There's actually not a lot here in terms of code. But, the code is improved here by removing a bunch of more complex code with more straightforward code that's reusing standard accumulation calls. So it's much more readable and maintainable. So thanks for bringing in this nice refactoring @samsrabin!

[samsrabin]

Thanks, @ekluzek! @slevis-lmwg actually already reviewed the latest version of the code; he just wanted me to do that test. So this should be good to go now.

[slevis-lmwg]

@samsrabin and I went over this before. I took a quick look again and will approve. Thank you for this work @samsrabin!