Derive the EAMxx source vertical grid file at run time

[chengzhuzhang]

Fixes #863.

Problem

zppy hardcoded an L128v1 vertical coordinate file (e3sm_to_cmip_data/grids/vert_L128.nc) as ncremap --vrt_in whenever prc_typ == 'eamxx'. E3SM#8692 makes L128v4 the default EAMxx grid, and because the level count is still 128, ncremap does not error on v4 input — it interpolates with v1 coefficients and emits silently wrong pressure-level data into ts_vrt_remap/ and cmip_ts/.

Setting vrt_in_file explicitly already worked as a user override, so this is a bad default rather than a missing capability. The problem is that a user who does not know to set it gets no error.

Change

Derive the source vertical grid from the run's own output rather than shipping a grid-version-specific data file:

run_nco ncks -O -v hyai,hybi,hyam,hybm,P0 ${file} vrt_in.nc 2> /dev/null || \
  { run_nco ncks -O -v hyai,hybi,hyam,hybm ${file} vrt_in.nc && \
    run_nco ncap2 -A -s 'P0=100000.0' vrt_in.nc ; }

(the fallback is the operation recorded in vert_L128.nc's own history attribute). This is grid-agnostic — L72, L128v1, L128v4 and any future EAMXX_VGRID work with no staged file and no zppy release.

• ts.bash derives from the raw history file (head -n 1 input.txt), which always carries hy*, so the [ts] path does not depend on extra_vars. Only ps remains load-bearing there, for --ps_nm. Failure gets its own status code 7.
• e3sm_to_cmip.bash derives lazily, inside the loop, from the file it is about to remap, guarded so it runs once. Eager derivation would have been a regression: interp_vars defaults to EAM names, so that loop usually remaps nothing for EAMxx, and a job that never calls ncremap must not fail on a missing derivation. Status code 5, with a message naming [ts] extra_vars as the fix.
• Both templates set a shell vrt_in_file, so an explicit vrt_in_file still takes precedence and short-circuits the derivation entirely.
• vrt_remap_file / cmip_plevdata (the --vrt_out plev19 target) are untouched — the target grid is unaffected by this.

Docs, default.ini comments, and the examples/post.v3.eamxx.cfg header no longer claim a vert_L128.nc default. No references to that file remain anywhere in zppy/, docs/, or examples/.

Tests

New tests/test_vertical_remap.py renders both templates and asserts three cases: derived (no vert_L128, derivation emitted), user-supplied vrt_in_file (path wins, no derivation), and prc_typ = eam (no --vrt_in / --ps_nm at all). 96 unit tests pass; pre-commit clean.

Verification

Checked against the L128v1 EAMxx case from examples/post.v3.eamxx.cfg (ne256pg2 ... F20TR-SCREAMv1, ne30pg2 monthly output, 1995-1999), whose ts/ and ts_vrt_remap/ output from the old staged-file code is still on scratch.

• The derived grid file reproduces vert_L128.nc bit-for-bit. hyai, hyam, hybi, hybm and P0 (double, 100000.0) are bit-identical, whether derived from the raw history file (the [ts] path) or from the split per-variable ts file (the [e3sm_to_cmip] path). lev top/bottom match at 2.5802608 / 998.49646 hPa. The derived file additionally carries ilev, which ncks pulls in as a coordinate; ncremap does not use it.
• ncremap output is bit-identical. Same ncremap line on a two-month subset of T_mid_199501_199912.nc, derived vs staged --vrt_in: maxabsdiff = 0.0 over all (2, 19, 180, 360) values, with no change in missing-value count. The same holds against the ts_vrt_remap/T_mid_199501_199912.nc that the old code wrote in July, so this is a real before/after comparison and not just self-consistency within one session.
• P0 is genuinely required, not assumed by NCO. With P0 stripped from an otherwise identical vrt_in file, ncremap fails: ERROR Failed to vertically interpolate. cmd_rgr[0] failed. The ncap2 step is load-bearing, hence the status check around it.
• P0 handling is future-proof. EAMxx writes no P0 today, but later versions will, so the derivation asks for it first and synthesizes 100000.0 only when that extraction fails. Checked against real files: a source without P0 gives the synthesized 100000; a source carrying P0 = 99999 keeps 99999 rather than being overwritten; a source missing the hybrid coefficients still fails loudly with the NCO error and trips the status code.
• L128v4 is picked up correctly, and the old default was badly wrong on it. Checked against an L128v4 EAMxx run (ne256pg2_ne256pg2.F2010xx-ZM-CICE.260306.splitform_TMSoff_Recipe2_updated, 1ma_ne30pg2 monthly output, 2010-01). The derivation succeeds on the raw file. That file has no P0, so the fallback supplies 100000.0. All four hybrid coefficients differ from vert_L128.nc (max |Δhybm| = 0.29), and the derived lev top is 0.593 hPa vs v1's 2.58 hPa (ilev top 0.499 vs 2.26 hPa).
  • Same ncremap line on T_mid to plev19, derived vs staged --vrt_in: the staged v1 file puts T wrong by up to 25.6 K (mean −21 K at 300–400 hPa), gives errors of several K throughout the stratosphere, and returns all-missing at 1 hPa, which is above v1's top.
  • Independent check at 500 hPa against the native level nearest 500 hPa, computed from the file's own coefficients: the derived output has a mean T500 of 257.97 K, vs 257.71 K for the native-nearest value (mean |Δ| 0.56 K, i.e. interpolation spacing). The staged output gives 274.43 K (mean |Δ| 16.7 K). This is the silent corruption described in Support new EAMxx L128v4 vertical grid #863. Nothing errors out.

🤖 Generated with Claude Code

[chengzhuzhang]

End-to-end test: L128v4 EAMxx, [ts] + [e3sm_to_cmip]

I ran zppy from this branch (e6b2f380) on Perlmutter with the L128v4 EAMxx run: ne256pg2_ne256pg2.F2010xx-ZM-CICE.260306.splitform_TMSoff_Recipe2_updated, using its 1ma_ne30pg2 monthly output for 2010–2012. This config doesn't set vrt_in_file, so both jobs used the new run-time derivation.

Config (post.eamxx.L128v4.vrtin.cfg; the same as the provenance copy zppy saved)

[default]
input = /pscratch/sd/t/terai/e3sm_scratch/pm-gpu/ne256pg2_ne256pg2.F2010xx-ZM-CICE.260306.splitform_TMSoff_Recipe2_updated
output = /pscratch/sd/c/chengzhu/tests/eamxx_L128v4_vrtin_09292026
case = ne256pg2_ne256pg2.F2010xx-ZM-CICE.260306.splitform_TMSoff_Recipe2_updated
www = /global/cfs/cdirs/e3sm/www/chengzhu/tests/eamxx_L128v4_vrtin_09292026
partition = "debug"
account = "e3sm"
campaign = "water_cycle"
debug = False

[ts]
active = True
walltime = "00:30:00"
years = "2010:2012:3"
ts_num_years = 3

  [[atm_monthly_180x360_aave]]
  input_component = "eamxx"
  input_subdir = "run/"
  case = "1ma_ne30pg2"
  input_files = "AVERAGE.nmonths_x1"
  frequency = "monthly"
  mapping_file = map_ne30pg2_to_cmip6_180x360_aave.20200201.nc
  vars="ps,surf_radiative_T,SeaLevelPressure,IceWaterPath,qv_2m,precip_liq_surf_mass_flux,precip_ice_surf_mass_flux,omega_at_500hPa,omega_at_700hPa,omega_at_850hPa,T_mid_at_700hPa,T_2m,surface_upward_latent_heat_flux,surf_sens_flux,z_mid_at_700hPa,wind_speed_10m,surf_evap,U_at_10m_above_surface,V_at_10m_above_surface,LW_clrsky_flux_dn_at_model_bot,LW_clrsky_flux_up_at_model_top,LW_flux_dn_at_model_bot,LW_flux_up_at_model_bot,LW_flux_up_at_model_top,SW_clrsky_flux_dn_at_model_bot,SW_clrsky_flux_dn_at_model_top,SW_clrsky_flux_up_at_model_bot,SW_clrsky_flux_up_at_model_top,SW_flux_dn_at_model_bot,SW_flux_dn_at_model_top,SW_flux_up_at_model_bot,SW_flux_up_at_model_top,ShortwaveCloudForcing,LongwaveCloudForcing,isccp_cldtot,U,V,T_mid,z_mid,omega,RelativeHumidity,p_mid,qv"
  extra_vars= "ps,hyai,hyam,hybi,hybm,area,landfrac,ocnfrac"
  vrt_remap_vars = "U,V,T_mid,z_mid,omega,RelativeHumidity,p_mid,qv"

[e3sm_to_cmip]
active = True
frequency = "monthly"
ts_grid = "180x360_aave"
ts_num_years = 3
walltime = "00:30:00"
years = "2010:2012:3"
qos = "debug"

  [[atm_monthly_180x360_aave]]
  input_component = "eamxx"
  case = "1ma_ne30pg2"
  input_files = "AVERAGE.nmonths_x1"

Results

Job	Status	Notes
ts_atm_monthly_180x360_aave_2010-2012-0003	OK (2m04s)	All 43 variables written with 36 time steps each. The 8 vrt_remap_vars were remapped to plev19 in ts_vrt_remap/.
e3sm_to_cmip_atm_monthly_180x360_aave_2010-2012-0003	OK (44s)	28/28 handlers converted with 0 failures. The 13 variables marked "non-derivable" are expected, because their EAMxx inputs weren't in vars.

• Derivation used: neither generated script refers to vert_L128. The [ts] job derived vrt_in.nc and remapped the 8 vrt_remap_vars. The [e3sm_to_cmip] job made no ncremap calls, because all of its 3D inputs were symlinks to files [ts] had already remapped (ts_vrt_remap/), so its derivation block wasn't run in this test. That block is covered by the unit tests and was checked by hand on a split ts file.
• plev19 output is physically reasonable (cmip_ts, shape (36, 19, 180, 360)):
  • ta ranges from 160 to 312 K, with a mean of 257 K at 1 hPa. 1 hPa is below the L128v4 top (0.593 hPa), so the top level holds real data. With the old L128v1 file, the 1 hPa level was all missing.
  • zg goes up to about 50 km, with a mean of 47.4 km at 1 hPa.
  • ua ranges from −96 to 184 m/s.
  • About 48% of the 1000 hPa level is masked where the surface pressure is below 1000 hPa. That's the expected below-ground masking.

Output: /pscratch/sd/c/chengzhu/tests/eamxx_L128v4_vrtin_09292026/post/

Unrelated issue found

The CMORized hur is a 0–1 fraction but is labeled %, because the EAMxx hur handler in e3sm_to_cmip has no unit conversion. This comes from e3sm_to_cmip, not from this PR. It's tracked in E3SM-Project/e3sm_to_cmip#350, with a fix in E3SM-Project/e3sm_to_cmip#351.

[chengzhuzhang]

@czender hey Charlie, this PR is made to support generating vrt_in file on the fly in zppy to replace the use of a static file, so that when eamxx vertical grid gets changed, zppy workflow remain works. Could you help review this PR, especially this block that derives vrt_in? Thank you.