[Bug] Interpolation does not pick correct time bounds for LAI

nicholasbalasus commented 6 months ago

Name and Institution (Required)

Name: Nick Balasus Institution: Harvard University

Confirm you have reviewed the following documentation

[x] Contributing guidelines

New HEMCO feature or discussion

This is an issue found while looking at the mercury simulation with @bgeyman.

In GEOS-Chem, leaf area index comes from files with a pattern like this: $ROOT/Yuan_XLAI/v2021-06/Yuan_proc_MODIS_XLAI.025x025.$YYYY.nc. These values are at 0.25° x 0.25° and are spaced 8 days apart.

Times for Yuan_proc_MODIS_XLAI.025x025.2019.nc:
[2019-01-02, 2019-01-10, 2019-01-18, 2019-01-26, 2019-02-03, 2019-02-11, 2019-02-19, 2019-02-27, 2019-03-07, 
 2019-03-15, 2019-03-23, 2019-03-31, 2019-04-08, 2019-04-16, 2019-04-24, 2019-05-02, ..., 2019-12-28]

Leaf area index uses the I time flag for interpolation to get daily values from the values spaced 8 days apart. As far as I can tell, the 73 XLAI entires are the only entries in any of the HEMCO config files that use the I flag. It something looks like this: * XLAI00 $ROOT/Yuan_XLAI/v2021-06/Yuan_proc_MODIS_XLAI.025x025.$YYYY.nc XLAI00 2000-2020/1-12/1-31/0 I ...

Using the attached run.sh, a 4° x 5° global methane simulation is run for 2019-2021, archiving only daily/instantaneous Met_MODISLAI. This is compared to the input leaf area files (Yuan_proc_MODIS_XLAI.025x025.YYYY.nc). Using the summed global leaf area as a diagnostic, it can be seen that the correct time interpolation bounds are not being chosen: lai_issue

For example, when the simulation date was 2019-01-03, the bounding dates for interpolation should've been [2019-01-02, 2019-01-10]. However, making HEMCO verbose, it can be seen that the bounding dates are actually [2019-01-02, 2019-05-02].

This issue can be traced back to the GetIndex2Interp subroutine in hcoio_util_mod.F90. The lower time bound (indicated by tidx1) is being selected correctly. To select the upper time bound (indicated by tidx2), this code is used: https://github.com/geoschem/HEMCO/blob/a5d5169aee598fc62fb0907d5aef5d7a5353ea41/src/Core/hcoio_util_mod.F90#L1160-L1219

The logic of this code is to start with the lower time bound (201901020000) and add first years (1e8) then months (1e6) then days (1e4). The value is added until the time exceeds the max time available for interpolation (201912280000) or an available time is landed upon. This is done instead of taking the next available time to accommodate some interpolation scenarios as described in this comment and mentioned by Bob/Christoph elsewhere. https://github.com/geoschem/HEMCO/blob/a5d5169aee598fc62fb0907d5aef5d7a5353ea41/src/Core/hcoio_util_mod.F90#L1147-L1154

Unfortunately, for the leaf area index files, when the loop gets to adding months after the year loop is not successful, it lands on the matching date of 201905020000 and thus chooses this instead of 20190110000. To avoid this issue, I suggest we remove the month loop. This should keep the functionality originally intended (as described in the wiki link above) while fixing the interpolation for leaf area index files. This simple fix is in my fork and can be run with the attached run-with-fix.sh. The plot below shows the fix is working: lai_issue_with_fix

I'm just unsure if this change will break something else, so I would appreciate any thoughts before opening the PR! Thank you.

run.txt run-with-fix.txt make-plots.txt

yantosca commented 6 months ago

Thanks @nicholasbalasus and @bgeyman. I thought we had fixed this but apparently not. You can go ahead and open the PR. We will have to add this into 14.4.0, as this will probably change fullchem simulation results (and thus will have to be benchmarked).

nicholasbalasus commented 6 months ago

Thanks, @yantosca!

nicholasbalasus commented 5 months ago

Closed by #254

geoschem / HEMCO