Adding Polvani Kushner setup to Isca, including polar heating options by Regan Mudhar

[sit23]

The Polvani-Kushner setup for a Newtonian-relaxation representation of the stratosphere has been a feature in Isca's branches for many years, but has yet to be incorporated into the master branch. The original code for PK was written by @wseviour, and was subsequently adapted by @rgnmudhar for her paper here:
https://doi.org/10.1029/2023JD040416

Getting this combination of commits into the master is slightly complicated, as I wanted to preserve as much of the history as possible, despite the original PK modifications being on Will's fork and the subsequent additions being on Regan's I've therefore asked Claude to pull these things together and create two new test cases for Isca:

1. Standard PK setup
2. Regan's PK setup with midlatitude and polar heating options

Claude has verified that setup 2 produces the same results as the same experiments on Regan's fork - https://github.com/rgnmudhar/Isca_rmudhar/tree/polvani_kushner_rm.

Now going to run the trip tests to check this hasn't broken anything else, then will merge in.

[sit23]

Here is Claude's description of the P/R:

Polvani-Kushner forcing: 2026 port

This branch (polvani_kushner_2026, off master) brings the Polvani & Kushner
(2002) stratospheric-vortex extension to hs_forcing.F90 into Isca proper,
combining Will Seviour's original implementation with Regan Mudhar's later
extensions, and adds two example test cases.

Both original branches predate this:

• Will Seviour: sit23/Isca:polvani_kushner_st
• Regan Mudhar: rgnmudhar/Isca_rmudhar:polvani_kushner_rm

Regan's branch had manually copied Will's original modifications rather than
building on his commits, so this PR reconstructs proper authorship for both.

What's on this branch

Will Seviour's PK core (5 commits, cherry-picked, authorship preserved)

81b0702c, e30d1c8b, 20f020f7, f3da7b20, 76b0d5a0 — cherry-picked
individually from polvani_kushner_st rather than applied as a squashed diff,
so each commit could 3-way-merge against current master and pick up the
five years of unrelated master drift (equinox_day, the COLUMN_MODEL ifdef,
etc.) cleanly instead of reverting it. Together these add:

• equilibrium_t_option = 'Polvani_Kushner': a US-standard-atmosphere-based
  summer/winter temperature profile with a polar vortex, following the P-K02
  formulation (tstd_summer/tstd_winter, vtx_edge/vtx_width/vtx_gamma,
  z_ozone, strat_vtx)
• sponge_flag: simple upper-level sponge-layer damping
• local_heating_option = 'Polar': an analytic Gaussian polar heating option
• relax_to_qbo: a simple analytical QBO-like equatorial wind relaxation

Two of Will's original commits were excluded as out of scope: one HPC
environment workaround (edits necessary to run on BluePebble, which also
pulled in an unrelated 3961-line comp.F90) and one related NetCDF-deflate
workaround for the same machine. Both are unrelated to the PK feature itself.

Regan Mudhar's extensions (2 commits, curated, authorship preserved)

Her branch has ~100 commits, but it forked from a master that predates
Will's own commits and diverged well before today's master too — so a
literal replay would have reintroduced stale code and fought both Will's
commits and years of subsequent master changes it never saw. Instead, her
intentional changes were isolated by diffing her branch's tip directly
against Will's tip (both from the same era, so this cancels out the
master-drift noise) and re-applied as clean commits on top of the ported
history, attributed to her:

• a1d179f0 — a jet-latitude control option (A, B, P_opt), following
  Garfinkel et al. (2013), added as a small perturbation to the equilibrium
  temperature's latitudinal profile. Defaults to A=B=0, i.e. off.
• 35c5983b — a generator for her combined polar-cap + midlatitude
  wavenumber-heating field (polar_heating() + heat_perturb() from her
  input/input_files.py, matching her combo_heat1() defaults). This
  reproduces her originally-committed
  input/asymmetry/w15a4p600f800g50_q6m2y45l800u200.nc to within float32
  rounding (max relative difference ~5e-7) — see Verification below.

Not ported: her postprocessing/plevel_interpolation workflow scripts and
machine-specific paths, which are personal-workflow tooling rather than a
generalisable Isca feature.

One thing worth noting for review: her tip-of-branch
polvani_kushner_test_case.py was actually pointed at her midlat-only
heating file, not the combined polar+midlat one — the combined file exists
in her repo as an artifact of an earlier, different run. The new
polvani_kushner_polar_heating_test_case.py test case here uses the
combined polar+midlat field, since that's the case actually being added.

Also note: her version of hs_forcing.F90 had dropped Will's QBO
relaxation and analytic Polar heating option (she never saw those commits,
and replaced polar heating with the more general file-driven approach
instead). Both are kept on this branch as available options, since they're
part of what Will's commits add and don't conflict with anything.

Test cases (exp/test_cases/polvani_kushner/)

• polvani_kushner_test_case.py — the baseline PK setup (Will's
  configuration: sponge on, QBO relaxation on, T42L40).
• polvani_kushner_polar_heating_test_case.py — PK plus the combined polar +
  midlatitude wavenumber-2 heating perturbation (Regan's PK_e0v4z13
  configuration: eps=0, vtx_gamma=4, z_ozone=13, t_strat=216.65,
  T42L60, dt_atmos=240s), generating its input file at run time via
  create_pk_heating_input_file.py rather than committing a binary NetCDF —
  following the existing frierson_dry_heating precedent. Run length was
  trimmed from her 504-month production setting to a 12-month illustrative
  example, matching the other test cases in the repo.

Both are registered in exp/test_cases/trip_test/trip_test_functions.py,
following the existing frierson/frierson_dry_heating pattern where the
more specific name overrides the general one.

Small cleanups made along the way

• Dropped a duplicate era_land_t42.nc that Will's original commit had
  committed locally under exp/test_cases/polvani_kushner/input/ — the test
  case actually points at the shared input/land_masks/ copy that already
  exists on master, so the local one was dead weight.
• Fixed a latent bug uncovered while merging Will's commits: one of his
  hunks (in an unrelated stratosphere_t_option == 'extend_tp' branch, not
  used by PK) changed a loop bound from size(t,2) to size(t,1). This
  turned out to be an artifact of the ~2020-era hs_forcing.F90 both his and
  Regan's branches forked from, since fixed on current master — the
  cherry-pick conflict was resolved to keep master's correct version rather
  than reintroducing it.
• Added exp/test_cases/polvani_kushner/input/*.nc to .gitignore, matching
  the existing frierson_dry_heating entry, since the heating input file is
  generated at test-case run time.

Verification

To confirm the ported branch reproduces Regan's original polar+midlat
heating experiment, I built both branches independently (this branch, and
her original Isca_rmudhar:polvani_kushner_rm as a reference) and ran her
exact PK_e0v4z13 configuration (pointed at the combined heating field) on
both for a short 2-day integration, from identical initial conditions.

Building her ~2021-era branch against today's compiler surfaced three
environment-compatibility issues, all already fixed on current master and
unrelated to any PK physics — patched locally in the reference checkout only
(not part of this branch) purely so it would compile:

• an sh Python library API change (sh>=2 returns str instead of a
  RunningCommand with .stdout) in several places in isca/codebase.py
  and isca/helpers.py
• modern gfortran's stricter handling of BOZ literal constants
  (-fallow-invalid-boz -fallow-argument-mismatch, missing from her
  mkmf.template.ubuntu_conda)
• a gettid() redeclaration conflict with glibc in affinity.c

With those purely-toolchain issues out of the way, both branches compiled
and ran cleanly. Comparing atmos_daily.nc output after 2 days:

field	max abs diff	max relative diff
ucomp	9.5e-7	1.1e-7
vcomp	9.5e-7	2.6e-7
temp	3.1e-5	1.1e-7
teq	3.1e-5	9.8e-8
local_heating	3.6e-11	5.3e-7
ps	7.8e-3	7.8e-8
height	7.8e-3	8.7e-8

All differences are at or below float32 machine epsilon (~1.2e-7) — i.e. the
port reproduces her original results to within floating-point rounding, not
a meaningful physical difference.

The baseline polvani_kushner_test_case.py (Will's QBO-relaxation
configuration) was also smoke-tested for 2 days and runs without error.

[sit23]

Have made some modifications to the PK test script as it included the QBO relaxation and land, neither of which we want for a more vanilla PK test case. Also asked Claude to run trip tests with the following results:

Update: baseline test case now matches Polvani & Kushner (2002) exactly

The original P-K02 setup has no land/topography and no QBO relaxation -
these were extras Will's commits added as options, but the baseline test
case shouldn't default to them. polvani_kushner_test_case.py (renamed
polvani_kushner_default) now:

• drops the era_land_t42.nc topography input, so topography_option falls
  back to its 'flat' default (no land/topography)
• drops relax_to_qbo/qbo_amp, so relax_to_qbo falls back to its
  default of False (no QBO relaxation)

Both options remain available in hs_forcing_nml for anyone who wants them -
this just changes what the example test case defaults to. Re-smoke-tested
for 2 days with no error.

Update: full trip_test run

Ran the full trip_test suite (master vs. this branch's tip, 3-day runs),
excluding Socrates-dependent cases (socrates_aquaplanet,
socrates_aquaplanet_cloud, ape_aquaplanet - GFDL_SOC isn't set up on this
host) and the two new PK test cases themselves (they use
equilibrium_t_option='Polvani_Kushner', which doesn't exist on master, so
there's nothing to bit-compare against there - they were already verified
separately, see above).

12 of 16 passed bit-identical. Of the 4 "failures":

• axisymmetric, grey_mars, radiative_eq_mars fail identically on
  master itself (a pre-existing interpolator_init ordering bug, and a
  pre-existing "run length must be multiple of atmosphere time step" issue
  with trip_test's hardcoded 3-day comparison window for the Mars cases) -
  confirmed unrelated to this PR.
• giant_planet showed a real ~2% difference from day 1. hs_forcing.F90's
  actual forcing subroutines are never called for this test case (it uses
  the moist idealized_moist_phys physics path, which only calls
  hs_forcing_init for diagnostic setup - confirmed via the call graph), so
  I ran giant_planet twice with the identical, unmodified master
  binary to check whether the model is bit-reproducible at all on this
  host. It isn't: two runs of the exact same code differed by a comparable
  magnitude (1.7e-2) to the master-vs-branch difference (7.7e-2). This
  clears the PR - giant_planet simply isn't run-to-run reproducible on
  this machine (most likely MPI reduction-order non-determinism, on a
  heavily-loaded shared host), independent of any code change.

[sit23]

So, this looks like it's pretty ready to be merged, with those edits pushed, which I'll do now.

[rgnmudhar]

Hi Stephen (and Claude) - thanks a lot for your work to incorporate this into the main Isca branch!

Overall it looks great! I mostly made comments relating to what we do or don't want to include.

In general, I think having a test case that is exactly like the Polvani & Kushner (2002) set-up makes sense, and then also another, separate one with the heating options. Crucially the Lindgren et al. 2019 midlatitude heating should be there (should probably rename the script to generic heating_test_case). For the polar heating, I am less certain as it's not really P-K related! I think it's nice not to lose it, so, if you agree with keeping it, then I think there are 2 key things:

1. I think you could move the polar_heating_XXX variables (L93-98 of the polvani_kushner_test_case, i.e. Will's heating) into the heating_test_case
2. for the version that takes an input file, optionalise it so you can either choose no heating (/ Will's heating) / Orlanski & Solman 2009 polar heating / Lingren et al. 2018 midlatitude heating / polar + midlatitude heating - this would require editing both the python script to create the input files and the test_case script

I also wondered if having a separate test case again for QBO would be good... could be worth checking a few points with @wseviour on this.

I hope everything is clear and I didn't miss anything! Let me know if you need a chat.

Thanks,
Regan

[rgnmudhar]

I personally feel that the jet-latitude control stuff is not part of the "base" P-K set-up so would not be crucial to include at this point if you really wanted to keep it clean. The original paper actually does it in Held-Suarez not Polvani-Kushner so could be worth putting in H-S instead/too (if you decide to keep it here). If keeping, in L101 I would add the link to the paper https://doi.org/10.1175/JCLI-D-12-00301.1

[rgnmudhar]

I would again add the link to the paper in the comment here on L106 https://doi.org/10.1029/2001GL014284

[rgnmudhar]

Change the comment here to "height of stratospheric heating source (default 13km (~200 hPa), 100hPa (~16km) used in in the P-K paper)"

[rgnmudhar]

Personally feel the polar cap heating is a separate thing to P-K (purely part of my AA experiments) so again could either remove the heating option stuff entirely, or at least keep references to the midlatitude wave-2 heating perturbation, or keep both but make it optional (see my idea below). Though the midlatitude heating is also not part of the base P-K set-up, it is what was necessary to introduce vortex variability, so could definitely see more of an argument for including.

One thought is to have a kind of "if/else" statement for the local_heating_option to give more flexibility. Something like:

# Generate the polar and/or midlatitude wavenumber-2 heating input file
# (all zeros unless the two functions below are given non-default arguments)
# For the polar heating see: https://doi.org/10.1175/2010JAS3267.1 and https://doi.org/10.1029/2023GL105132
# For the midlatitude heating see: https://doi.org/10.1029/2018JD028537
# rather than committing a large binary to git - see
# src/extra/python/scripts/create_pk_heating_input_file.py.
input_dir = os.path.join(GFDL_BASE, 'exp/test_cases/polvani_kushner/input')
os.makedirs(input_dir, exist_ok=True)
heat_type = 'midlat' # default midlat only, other options 'polar', 'midlat_polar', None
if heat_type == 'midlat':
          heating_file_path, heating_var_name = create_midlat_heating_file(input_dir)
elif heat_type == 'polar':
          heating_file_path, heating_var_name = create_polar_heating_file(input_dir)
elif heat_type == 'midlat_polar':
          heating_file_path, heating_var_name = create_polar_and_midlat_heating_file(input_dir)

But then the script creating the input file needs to be edited too and you would probably need to do something with subsequent lines in the present script related to the heating, e.g. L40, L57, and L115-118.

[rgnmudhar]

Also could be worth folding in the other kind of "polar heating" from the other test_case script by Will to this script, so the other is a completely clean set-up that follows the P-K paper exactly

[rgnmudhar]

change to 16. and the comment to be consistent with the "polar_heating_test_case" version i.e. "height of stratospheric heating source (default 16km (~100 hPa) used in in the P-K paper)"

[rgnmudhar]

Don't need qbofactr if not keeping qbo nudging

[rgnmudhar]

As per my earlier comment, could this be rewritten (perhaps two new separate functions?) to create a polar heating only and midlatitude heating only file, as well as this current combined one?

[rgnmudhar]

Actually I think the default number of levels should be 40, can keep 60 in the comment maybe?

[rgnmudhar]

Could remove L930-940 if decide not to keep this alternative type of polar heating from Will

[rgnmudhar]

Would need to edit if changing to generic heating and also having a separate test case for QBO