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

[nicholasbalasus]

Name and Institution (Required)

Name: Nick Balasus
Institution: Harvard University

Confirm you have reviewed the following documentation

• 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:

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:

HEMCO/src/Core/hcoio_util_mod.F90

Lines 1160 to 1219 in a5d5169

	! Search for a time slice in the future that has the same
	! month/day/hour as currently selected time slice.
	tmpYMDhm = availYMDhm(tidx1)
	DO
	! Increase by one year
	tmpYMDhm = tmpYMDhm + 1.0e8_dp
	! Exit if we are beyond available dates
	IF ( tmpYMDhm > availYMDhm(nTime) ) EXIT
	! Check if there is a time slice with that date
	DO I = tidx1,nTime
	IF ( tmpYMDhm == availYMDhm(I) ) THEN
	tidx2 = I
	EXIT
	ENDIF
	ENDDO
	IF ( tidx2 > 0 ) EXIT
	ENDDO
	! Repeat above but now only modify month.
	IF ( tidx2 < 0 ) THEN
	tmpYMDhm = availYMDhm(tidx1)
	DO
	! Increase by one month
	tmpYMDhm = tmpYMDhm + 1.0e6_dp
	! Exit if we are beyond available dates
	IF ( tmpYMDhm > availYMDhm(nTime) ) EXIT
	! Check if there is a time slice with that date
	DO I = tidx1,nTime
	IF ( ABS( tmpYMDhm - availYMDhm(I) ) < EPSILON ) THEN
	tidx2 = I
	EXIT
	ENDIF
	ENDDO
	IF ( tidx2 > 0 ) EXIT
	ENDDO
	ENDIF
	! Repeat above but now only modify day
	IF ( tidx2 < 0 ) THEN
	tmpYMDhm = availYMDhm(tidx1)
	DO
	! Increase by one day
	tmpYMDhm = tmpYMDhm + 1.0e4_dp
	! Exit if we are beyond available dates
	IF ( tmpYMDhm > availYMDhm(nTime) ) EXIT
	! Check if there is a time slice with that date
	DO I = tidx1,nTime
	IF ( tmpYMDhm == availYMDhm(I) ) THEN
	tidx2 = I
	EXIT
	ENDIF
	ENDDO
	IF ( tidx2 > 0 ) EXIT
	ENDDO

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.

HEMCO/src/Core/hcoio_util_mod.F90

Lines 1147 to 1154 in a5d5169

	! If we are inside the data range but none of the time slices
	! matches the preferred datetime, get the second time slices that
	! shall be used for data interpolation. This not necessarily needs
	! to be the consecutive time slice. For instance, imagine a data
	! set that contains montlhly data for years 2005 and 2010. For
	! Feb 2007, we would want to interpolate between Feb 2005 and Feb
	! 2010 data. The index tidx1 already points to Feb 2005, but the
	! upper index tidx2 needs to be set accordingly.

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:

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]

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]

Thanks, @yantosca!

[nicholasbalasus]

Closed by #254