EAMxx: Fixed temporal interpolation for external forcings

[TaufiqHassan]

Fixes temporal interpolation for external forcings in MAMxx to match EAM behavior (as discussed in #7315).

Previously, prescribed vs EAMxx elevated emissions (bc_a4 for example) showed significantly large discrepancies (see also in https://acme-climate.atlassian.net/wiki/spaces/EAMXX/pages/5240389767/Emission+treatment):

   months     prescribed_emis          EAMxx    rel_diff (%)
0     jan  50.720847076634016        41.87359    21.128494
1     feb  33.052801575716245       33.808193    -2.234345
2     mar  34.564076987260805        34.10844     1.335847
3     apr    33.6571496550456        34.85354    -3.432618
4     may   36.04996553155214        46.04442   -21.706112
5     jun   56.01655671895189        71.36652   -21.508638
6     jul    86.6789942679998        98.19908   -11.731359
7     aug  109.69631225550144      105.578384     3.900351
8     sep   101.4825383463081         72.2314    40.496431
9     oct   43.06931282750091        37.60891    14.518909
10    nov   32.16846281675161       40.085526   -19.750428
11    dec  47.985367140309215       49.352615    -2.770366

As opposed to EAM-

   months     prescribed_emis            EAM        rel_diff (%)
0     jan   65.53070785855884    58.4741778369735    12.067771
1     feb  32.950843360402935   36.67505759330244   -10.154624
2     mar  33.296972174272504   32.59856770545893     2.142439
3     apr  28.873338098498838    31.1443850952515    -7.291995
4     may  41.523370000927116   41.45759393878491     0.158659
5     jun   49.88711151401323  54.040493434934405    -7.685685
6     jul   85.16360864281151    82.8069247826244     2.845999
7     aug   94.06291468483843   90.19021139784672     4.293929
8     sep   74.23516764683626   70.60245248538006      5.14531
9     oct   33.96152809508068  37.868317940064216   -10.316777
10    nov  31.705675644882124   33.74076635265753    -6.031549
11    dec   44.28931694279417   46.16966256252293    -4.072687

With this fix, prescribed vs EAMxx show similar behaviour as EAM
(Using identical emissions in the reference EAM run above)

   months     prescribed_emis      EAMxx   rel_diff (%)
0     jan   65.53070828175305   58.45631    12.102026
1     feb    32.9508434047515  36.658974   -10.115205
2     mar   33.29697212313313  32.588203     2.174924
3     apr   28.87333817845113  31.138515    -7.274519
4     may   41.52337000418099  41.451073     0.174416
5     jun   49.88711150420784  54.034306    -7.675113
6     jul   85.16360885120517  82.792755       2.8636
7     aug   94.06291452961088   90.16316     4.325217
8     sep   74.23516786718005   70.56994     5.193754
9     oct  33.961528086365256   37.84876   -10.270431
10    nov   31.70567563396028  33.733593    -6.011566
11    dec  44.289317262530815  46.165443    -4.063919

Note:

• These changes apply only to external forcings. Other tracer types (e.g., oxidants, LINOZ) still follow the existing/default method.
• The current implementation is only compatible with emission datasets that follow a periodic annual cycle.

[NBFB]

[odiazib]

@TaufiqHassan Thank you for doing all this work. I recommend that we start using the DataInterpolation class that @bartgol implemented. By using this class, we will not need the time interpolation scheme implemented in this PR. I believe using this class is a better option because we will avoid duplicating code and maintain consistency with the rest of EAMxx.

I have started looking at the SPA code (see here: SPA Code), where this DataInterpolation is used. The only issue we need to address before starting to refactor the MAM4xx code is that DataInterpolation does not have the vertical interpolation needed by MAM4xx.

@singhbalwinder @TaufiqHassan @kaizhangpnl @mahf708
One option would be to add the MAM4xx vertical interpolation to DataInterpolation. The other option would be to review the current vertical interpolations in DataInterpolation and choose one option (see here).

I can work on refactoring the MAM4xx reader to use DataInterpolation. However, we must decide what our next steps are for vertical interpolation.

Finally, I am also okay with merging this PR and updating the code with DataInterpolation in the coming week.

[TaufiqHassan]

@TaufiqHassan Thank you for doing all this work. I recommend that we start using the DataInterpolation class that @bartgol implemented. By using this class, we will not need the time interpolation scheme implemented in this PR. I believe using this class is a better option because we will avoid duplicating code and maintain consistency with the rest of EAMxx.

I have started looking at the SPA code (see here: SPA Code), where this DataInterpolation is used. The only issue we need to address before starting to refactor the MAM4xx code is that DataInterpolation does not have the vertical interpolation needed by MAM4xx.

@singhbalwinder @TaufiqHassan @kaizhangpnl @mahf708 One option would be to add the MAM4xx vertical interpolation to DataInterpolation. The other option would be to review the current vertical interpolations in DataInterpolation and choose one option (see here).

I can work on refactoring the MAM4xx reader to use DataInterpolation. However, we must decide what our next steps are for vertical interpolation.

Finally, I am also okay with merging this PR and updating the code with DataInterpolation in the coming week.

@odiazib, that sounds great. This PR follows the same approach used in DataInterpolation, though in a more lightweight form. I agree that optimizing the code by removing any redundant or duplicate logic is a good idea. Please feel free to merge it whenever you think it’s appropriate, or wait until the DataInterpolation integration is complete. I’m happy with either option.

[bartgol]

What special vertical interpolation algorithm does MAM need?

[mahf708]

SPA uses vertical interpolation too; am I misunderstanding your point?

My vote would to be merge this for now, and then you can take time reworking the infra on the mam side. If you prefer to get input from me and Luca, we can help you. I happen to have strong opinions on design, etc., and how things can be simplified in mam, and I think Luca will have the final say since he knows all the details better than me. Though you should feel free to go your way as well

[singhbalwinder]

@kaizhangpnl and @TaufiqHassan: If this functionality is a blocker for you, we can merge this so that it is available for you to use; otherwise, it will be better to use the DataInterpolation class.

@bartgol: @odiazib implemented vertical implementations that work for the pressure level and height (in meters?). Here is the file. I may be wrong, but I think SPA vertical interpolation works only for the pressure levels.

[kaizhangpnl]

@kaizhangpnl and @TaufiqHassan: If this functionality is a blocker for you, we can merge this so that it is available for you to use; otherwise, it will be better to use the DataInterpolation class.

@bartgol: @odiazib implemented vertical implementations that work for the pressure level and height (in meters?). Here is the file. I may be wrong, but I think SPA vertical interpolation works only for the pressure levels.

I can use Taufiq's branch for my upcoming testing.

[TaufiqHassan]

@kaizhangpnl and @TaufiqHassan: If this functionality is a blocker for you, we can merge this so that it is available for you to use; otherwise, it will be better to use the DataInterpolation class.

This isn’t a blocker for my work. I can always use this branch if needed.

[odiazib]

@bartgol

What special vertical interpolation algorithm does MAM need?
@mahf708
SPA uses vertical interpolation too; am I misunderstanding your point?

In MAM4xx, we have at least two types of vertical interpolation.

Case One: It uses pressure.
See the function call here:
Function Call
See the function implementation here:
Function Implementation

Case Two: It uses altitude.
See the function call here:
Function Call
See the function implementation here:
Function Implementation

When I implemented this reader, I wanted to use SPA vertical interpolation, but at that time, I thought it would be easier to keep the vertical interpolation as it is implemented in EAM.

[mahf708]

Looks good to me.

I think this introduces some different logic, but given that this will be rewritten, I think it is ok to integrate it for now, and worry about other details later.

I think integrating this now is likely the better option than waiting for the larger refactor of all DataInterpolation business in MAM, and it would likely be less frustrating than maintaining a branch especially with other fixes coming in, etc.

[mahf708]

@odiazib thanks for clarification --- I agree with your inclination to go for pressure-based interpolation, and that should likely be largely plug-n-play, right?

[odiazib]

@odiazib thanks for clarification --- I agree with your inclination to go for pressure-based interpolation, and that should likely be largely plug-n-play, right?

For the pressure-based interpolation, there were differences between MAM4 and SPA in how extrapolation was handled. For the interpolation based on altitude, I believe the option would be to change the approach to how the data is generated in order to switch from altitude-based to pressure-based interpolation. Regarding vertical emission, which is what @TaufiqHassan is addressing with the time interpolation, vertical interpolation is altitude-based (is this correct, @TaufiqHassan?). For this, I would like to hear from @kaizhangpnl and @TaufiqHassan, who have a better understanding of these details.

[TaufiqHassan]

@odiazib Yes, these are attitude-based interpolations. iirc MAM uses an interval-based conservative rebinning between the source altitudes (km) and model interface heights (zint after reversing and converting to km).

[kaizhangpnl]

@odiazib thanks for clarification --- I agree with your inclination to go for pressure-based interpolation, and that should likely be largely plug-n-play, right?

For the pressure-based interpolation, there were differences between MAM4 and SPA in how extrapolation was handled. For the interpolation based on altitude, I believe the option would be to change the approach to how the data is generated in order to switch from altitude-based to pressure-based interpolation. Regarding vertical emission, which is what @TaufiqHassan is addressing with the time interpolation, vertical interpolation is altitude-based (is this correct, @TaufiqHassan?). For this, I would like to hear from @kaizhangpnl and @TaufiqHassan, who have a better understanding of these details.

The original MAM implementation follows the AeroCom protocol (https://acp.copernicus.org/articles/6/4321/2006/acp-6-4321-2006.html) by prescribing elevated emissions using specified injection heights. These emissions should be distributed across different pressure levels based on the instantaneous atmospheric state. Therefore, it's more straightforward to interpolate the emissions from specified injection altitudes to the model pressure levels using predicted geopotential height fields.

[bartgol]

DataInterpolation uses VerticalRemapper for the vertical interpolation, and knows very little about it. VerticalRemapper uses "pressure", but it really doesn't know what "pressure" means, so it could be fed other fields. However, it DOES assume that the "coordinate" (which it calls "pressure") is increasing with the view level index. That works for pressure, but not for elevation.

Supporting altitude is really just about supporting an x coordianate that is decreasing with index in ekat::LinInterp (although I'm not sure if the src altitude is in the input file or computed on the fly from state vars). A while ago I started tinkering with a "BackwardView" class that wrapped a 1d view and accessed it "from the end". If altitude interpolation is what you guys need, we can try to revive that effort (or maybe find a more hard-coded solution).

In short, for altitude vremap, 90% of the work is supporting a "bwd" x in LinInterp, while the remaining 10% is changing the var names for the "vertical coordinate" in VerticalRemapper and DataInterpolation to avoid confusion (e.g., pressure->vert_coordinate).

I agree with Naser: merge this if you need it. If you are ok with pressure-based interpolation, switching to DataInterpolation should be "quick". If you want to keep altitude-based interpolation, it will take some time...

[kaizhangpnl]

Thank you all for the helpful comments. Based on the discussions, I suggest we merge this and take it as a reference for verifying future implementations, regardless of which option we ultimately adopt.

[kaizhangpnl]

Please include a note in the PR description indicating that the current implementation is only compatible with emission datasets that follow a periodic annual cycle.

[TaufiqHassan]

Added in description.

[singhbalwinder]

Thank you all. All the fails are expected. I will start merging this PR.