Skip to content

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

Open
sit23 wants to merge 12 commits into
ExeClim:masterfrom
sit23:polvani_kushner_2026
Open

Adding Polvani Kushner setup to Isca, including polar heating options by Regan Mudhar#308
sit23 wants to merge 12 commits into
ExeClim:masterfrom
sit23:polvani_kushner_2026

Conversation

@sit23

@sit23 sit23 commented Aug 25, 2026

Copy link
Copy Markdown
Contributor

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.

wseviour and others added 11 commits August 25, 2026 10:09
Adds namelist options A, B and P_opt controlling an additional
zonal-wavenumber jet-latitude forcing term applied to the equilibrium
temperature profile, following Garfinkel et al. (2013). Defaults to
A=B=0, which leaves teq unchanged from the existing Polvani_Kushner
behaviour.

Ported from rgnmudhar/Isca_rmudhar polvani_kushner_rm (e.g. commits
717f1bc, f9b2675), curated on top of the 2026 port of Will Seviour's
original PK commits.
…ut file

Ported from polar_heating()/heat_perturb()/combo_heat1() in
rgnmudhar/Isca_rmudhar polvani_kushner_rm:input/input_files.py. Rather than
committing the generated binary NetCDF (following the frierson_dry_heating
precedent of generating local heating input files at run time), this
generates the field from the repo's existing T42 grid file and ozone
climatology file. Verified to reproduce the original committed
input/asymmetry/w15a4p600f800g50_q6m2y45l800u200.nc from Isca_rmudhar to
within float32 rounding (max relative difference ~5e-7).
era_land_t42.nc was committed locally in Will Seviour's original PK commit,
but the test case script actually points at the shared input/land_masks/
copy that already exists on master, so the local copy was dead weight.
…r-2 heating)

Ported from rgnmudhar/Isca_rmudhar polvani_kushner_rm:exp/test_cases/polvani_kushner/polvani_kushner_test_case.py
(her PK_e0v4z13 configuration: eps=0, vtx_gamma=4, z_ozone=13, t_strat=216.65,
T42L60, dt_atmos=240s, sponge on), adapted to:
 - generate the combined polar + midlatitude wavenumber-2 heating input file
   on the fly via create_pk_heating_input_file.py rather than requiring a
   committed binary NetCDF (matching this repo's frierson_dry_heating
   precedent), using the w15a4p600f800g50_q6m2y45l800u200 combo she
   generated via combo_heat1()'s defaults
 - shorten the run loop from her 504-month production length to a 12-month
   illustrative example, consistent with the other test cases in this repo
Follows the existing pattern (e.g. frierson/frierson_dry_heating) where the
more specific name is checked after the general one so it can override
input_files/nml_out/codebase_to_use.
Matches the existing frierson_dry_heating/input/*.nc entry - the heating
input file is generated at test-case run time (see
create_pk_heating_input_file.py), not committed.
@sit23

sit23 commented Aug 25, 2026

Copy link
Copy Markdown
Contributor Author

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:

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.

The original Polvani & Kushner (2002) setup has no topography/land contrast
and no QBO relaxation - both were extras Will Seviour's original commits
added as options but the baseline test case shouldn't default to. Removing
the era_land_t42.nc topography input and spectral_init_cond_nml block lets
topography_option fall back to its 'flat' default; removing relax_to_qbo/
qbo_amp lets relax_to_qbo fall back to its default of False. Renamed the
experiment from polvani_kushner_qbo3 to polvani_kushner_default to match.

Smoke-tested for 2 days with no error.
@sit23

sit23 commented Aug 25, 2026

Copy link
Copy Markdown
Contributor Author

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

sit23 commented Aug 25, 2026

Copy link
Copy Markdown
Contributor Author

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

@sit23
sit23 marked this pull request as ready for review August 25, 2026 11:26
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants