EAMxx: Data interpolation for MAM4 microphysics.

[odiazib]

Our current implementation of the reader for Linoz, invariants, and vertical emission files uses an approach based on the former SPA reader. @bartgol implemented the DataInterpolation class, which provides more flexibility to handle different types of files. Here, we are using this DataInterpolation class to read the necessary files for the MAM4 microphysics interface. Because MAM4 still uses the vertical interpolation from EAM (the ported version), I have added three new vertical interpolation types to this class and created a VerticalRemapperMAM4 class for MAM4 vertical interpolation.

This PR also fixes the time interpolation issue for Linoz, invariant, and vertical emission fields.

Because we will use single-year NetCDF files for elevated emissions, invariants, and Linoz, the following namelist parameters were deleted: mam4_linoz_ymd, mam4_oxid_ymd, and elevated_emiss_ymd.

I implemented the VerticalRemapperMAM4 class, which inherits from VerticalRemapper, to invoke MAM4xx vertical interpolation routines that were ported from EAM. The DataInterpolation class handles VerticalRemapperMAM4 as a custom vertical remapper.

In particular, VerticalRemapperMAM4 supports three types of vertical interpolation:

1. `VerticalRemapperMAM4::MAM4_PSRef` This is equivalent to `DataInterpolation::VRemapType::Dynamic3DRef`, where the source pressure field is computed as:

psrc​(icol,k)=ps_v(icol)⋅hybm(k)+P0​⋅hyam(k)

The main difference is that MAM4_PSRef invokes the mam4::vertical_interpolation::vert_interp routine.

I added a boolean parameter, mam4_use_mam4xx_oxi_vert_remap, to allow users to switch between MAM4_PSRef and DataInterpolation::VRemapType::Dynamic3DRef. Currently, this parameter is hardcoded to true, meaning oxidant vertical interpolation will be performed using the MAM4xx routine. However, if the evaluation team is interested in testing the use of DataInterpolation::VRemapType::Dynamic3DRef, this parameter can be moved to the namelist.

2. VerticalRemapperMAM4::MAM4_ZONAL

This vertical interpolation type is similar to DataInterpolation::Static1D, where a 1D variable is used as the independent variable. In the case of MAM4_ZONAL, the independent variable is lev with units of mbar. The main difference is that MAM4_ZONAL invokes the mam4::vertical_interpolation::vert_interp routine. Note that in this case, unit conversion from mbar to Pascals is handled by mam4::vertical_interpolation::vert_interp.

3. VerticalRemapperMAM4::MAM4_ELEVATED_EMISSIONS
This vertical interpolation type uses altitude as the independent variable. There is no equivalent in the VerticalRemapper for this case.

Tests (DataInterpolation for elevated emissions)

I tested this PR using only the DataInterpolation for elevated emissions. In this case, we expect the tests to pass, including baselines (BFB). However, I noted that the single-precision tests fail for gcc-openmp and gcc-cuda during the baseline comparison. Additionally, the test REP_D_Ln5.ne4pg2_oQU480.F2010-EAMxx-MAM4xx.ghci-snl-cpu_gnu is also failing in the baseline comparison.

gcc-cuda/sp and gcc-openmp/sp

Click to view detailed comparison

Most of the standalone tests are passing, except for gcc-cuda/sp and gcc-openmp/sp. It is important to note that the double-precision tests are passing, which means this implementation using the DataInterpolation class is producing BFB results consistent with the former reader (using @TaufiqHassan's correction for time interpolation).

REP_D_Ln5.ne4pg2_oQU480.F2010-EAMxx-MAM4xx.ghci-snl-cpu_gnu

Click to view detailed comparison

This test is failing with DIFFs, and these DIFFs increase as the number of time steps increases. I carefully reviewed this test and concluded that these DIFFs are very small in the first iteration. However, as MAM4 affects host model variables (in particular, temperature), these DIFFs increase over time. I conducted a couple of experiments where I ran both the old and new readers, and I found that the DIFFs remain small. For example:

Sector	BFB Hash	Value
forc_new_8_sector_0	12e765df7450c4b7	243788.7310767502
forc_old_8_sector_0	12e765df7450e328	243788.7310767502
forc_new_8_sector_1	4805cb37e8068d23	151464.7675245383
forc_old_8_sector_1	4805cb37e80691e7	151464.7675245383
forc_new_8_sector_2	3e98c5625d38c1fd	262341.6782332305
forc_old_8_sector_2	3e98c5625d38caa5	262341.6782332306

Using preprocessor directives (USE_OLD_LINOZ_FILE_READ and USE_OLD_VERTICAL_FILE_READ) was causing confusion for the evaluation team. Therefore, I decided to delete the former reader.

Tests (DataInterpolation for both invariants and linoz)

Before this PR, we knew there was a time interpolation issue inherited from the former SPA reader. Therefore, using the DataInterpolation class for invariants and Linoz will produce NBFB results in all tests where microphysics processes are involved. Thus, baseline tests are expected to fail for gcc-cuda, gcc-openmp, and REP_D_Ln5.ne4pg2_oQU480.F2010-EAMxx-MAM4xx.ghci-snl-cpu_gnu.

During the testing of DataInterpolation for invariants and Linoz, we found a bug in the former reader and noted that the DataInterpolation class only works for files containing single-year data. The bug in the former reader involved selecting the wrong year for interpolation during a restart, which caused the restart to be NBFB. This bug was fixed by replacing the NetCDF files for oxidants and Linoz. The new files contain only single-year data to accommodate the current features of the DataInterpolation class.

NC files

Elevated emissions

The DataInterpolation class relies on the time variable for time interpolation. In the elevated emission files currently used in the master branch, the time variable is not present. Instead, these files have the date variable. @meng630 helped add the time variable to these files. Note that we do not expect any changes due to these file updates.

Invariants

We are using new invariant files for the year 2015. We note that the current files do not contain data from 2010.

Linoz

We are using a single-year file for Linoz for the year 2010. In addition, we replaced the time variable and created a new time variable that does not assume years of 365 days and months of 30 days, which was causing interpolation issues in the DataInterpolation class. Note that the DataInterpolation class uses time instead of date for time interpolation.

NBFB

This PR introduces modifications that, in theory, should not produce significant changes and is expected to be a BFB PR. However, we have identified a few tests that are failing due to baseline comparisons. As a result, this PR will be classified as NBFB.

These tests are failing because of the baselines:

gcc-openmp/sp and gcc-cuda/dbg
mam4_aero_microphys_standalone_baseline_cmp

DIFF REP_D_Ln5.ne4pg2_oQU480.F2010-EAMxx-MAM4xx.ghci-snl-cpu_gnu.eamxx-L72 (phase BASELINE)

Key Findings

@singhbalwinder and @odiazib investigated the cause of the BFB behavior and found that the issues primarily arise from the alpha parameter involved in the linear interpolation. In the `DataInterpolation` class, our current implementation for reading tracers in `mam4xx` utilizes the following linear interpolation formula:

x_interpolated = x1 * alpha + x2 * (1 - alpha)

where:

• x1 is the value at time 1,
• x2 is the value at time 2,
• x_interpolated is the value at the current time, and
• alpha = (t_now - t1) / (t2 - t1).

1. DIFFs in Standalone Tests: We noted that DIFFs in standalone tests are only present in the single precision cases (CUDA and OpenMP). In these instances, the alpha parameter is hard-coded to double in the DataInterpolation class, leading to the observed DIFFs. For reference, see the implementation here.

2. Impact on t_now: In the case of REP_D_Ln5.ne4pg2_oQU480.F2010-EAMxx-MAM4xx.ghci-snl-cpu_gnu.eamxx-L72, we found that the changes in this PR are producing small variations in the value of t_now. These differences are in the last digits of t_now and propagate to alpha.

[NBFB]

[bartgol]

I am not entirely sure if this should be an atm proc, or if you should put the data read/interpolation feature inside the already existing eamxx_mam_generic_process_interface.*pp files.

That said, I have some thoughts, and some cleanup suggestions. I did not get into the details of the mam processes or in the reader utils.

[bartgol]

This is fine, but FYI: if you have multiple copies of a parameter (like mam4_linoz_file_name here), the metadata type and doc only need to be set once in the first occurrence. It's ok to put it in all of them, but for long doc strings, it may help to just add it once.

Again, this is fine.

[odiazib]

Since I plan to delete the new process mam4_interpolation_aero_microphys, all duplicate parameters will be removed. I will move the reader to the aero_microphys atm process.

[bartgol]

A few thoughts on this:

1. You only need to compute vertical layer heights if you have an elevation-based interpolation no? Wouldn't it make sense to put the calc of zmid in an if statement?

2. There seem to be quite a bit of code related to computing this quantity. Wild idea: why not storing a VerticalLayerDiagnostic obj in this class, and running the diag to compute the zmid field? Wouldn't you be able to get rid of a bunch of code in this class (e.g., buffers, and maybe even wet/dry stuff)?

[odiazib]

1. Yes, it does. I copied and pasted this block of code from the microphysics interface, where zmid is needed. I will delete this interface.
2. The wet/dry block of code is from the microphysics interface. I will explore the VerticalLayerDiagnostic approach.

[odiazib]

@bartgol Can we review the VerticalLayerDiagnostic in a follow-up pull request (PR)?

[bartgol]

There's a few unused typedefs here. Let's keep the class header lean (and easier to read).

[odiazib]

This interface will be removed.

[bartgol]

I think you can get rid of the set_xyz_pressure impl, since the base class impl should now be enough.

[odiazib]

I removed set_source_pressure(const Field& p). For the elevated emission case, which is Custom type in DataInterpolation, I am still using set_target_pressure(const Field& p) and set_source_pressure(const std::string& file_name).

[bartgol]

These typedefs are also unused and can be removed. Even ThreadTeam is not strictly needed, since you can just use auto in the lambda (but at least you could mv the using line to the cpp, where it's used, so you can rm the mam4.hpp include, and keep the header leaner).

[singhbalwinder]

These comments are likely for old version of this PR (likely on the deleted files for which the code has been moved). Ignore if not relevant.

[singhbalwinder]

Thanks, Oscar for taking on this major task! I have left some minor comments and suggestions. Overall looks great!

Could you please add some comments to this code? It would be good to add comments to at least separate out sections for different types of files (e.g., oxidants, LINOZ, elevated emissions).
Also adding comments about what kind of interpolation they need (time, horizontal, vertical...and what kind of vertical). It will also be good to mention the zonal LINOZ file (lat,lev) in the comments here.

[singhbalwinder]

We already have all these names (linoz and oxids) in a vector. Can we use a for-loop instead of explicitly using names?

[TaufiqHassan]

Hi @odiazib, this looks great! I've looked at the external forcings from a ne4pg2 simulation and compared the results against an older simulation using PR#7381

The results look consistent as expected from the description. Testing it on ne30pg2 grid against EAM.

[TaufiqHassan]

The ne30pg2 comparison against EAM looks consistent with different 2016 emissions data.

[singhbalwinder]

@odiazib has identified reasons for the NBFB behavior we are seeing for some tests (see PR description). I am planning to merge this PR today.