Cg_drag and mg_drag integration from Ross' fork into main version

[sit23]

@rosscastle did a lot of work with cg_drag and mg_drag in his PhD, and his versions of Isca have not been merged from his fork. Here we attempt (with Claude's help) to do that merging procedure.

Merge orographic and non-orographic gravity wave drag schemes from Ross Castle's fork

Summary

This PR merges two gravity wave drag (GWD) schemes that had been developed independently on separate, never-merged branches of Ross Castle's fork, into a single working branch on sit23/Isca:

• Non-orographic / convective GWD (cg_drag), from ross/gwdp_chaim_steven
• Orographic / mountain GWD (mg_drag), from ross/realistic_mgdrag

Both branches independently modified idealized_moist_phys.F90 and damping_driver.f90, so these were hand-reconciled rather than merged blindly. Two pre-existing bugs were found and fixed while building short on/off comparison test cases to verify both schemes (see below). One suggested fix from an open upstream issue was investigated and confirmed unnecessary for this branch, with empirical evidence.

Files changed

src/atmos_param/cg_drag/cg_drag.f90                 | 642 ++++++++++++---------
src/atmos_param/damping_driver/damping_driver.f90   |  31 +-
src/atmos_param/mg_drag/mg_drag.f90                 | 257 ++++++---
src/atmos_spectral/driver/solo/idealized_moist_phys.F90 |   2 +-
4 files changed, 583 insertions(+), 349 deletions(-)

1. Non-orographic gravity wave drag — cg_drag

Based on Alexander & Dunkerton (JAS, 1999). This version descends from an implementation by Chaim Garfinkel, with edits from Stephen Thomson, further developed on Ross Castle's fork.

Chaim Garfinkel's contribution (Jan 2017)

Four parameters were added to the source-level wave spectrum calculation (!add by chaim jan 2017, cg_drag.f90):

• kelvin_kludge (default 1.0) — an asymmetric tuning factor applied only to the westward-propagating branch of the source-level wave momentum flux spectrum (B0(n) = B0(n)*kelvin_kludge where c0mu0(n) < 0). This is a pragmatic correction for the eastward/westward wave-stress asymmetry needed to sustain a realistic QBO without resolving actual equatorial Kelvin waves — hence the name.
• weightminus1, weightminus2, weighttop — vertical-weighting parameters intended for smoothing the wave forcing profile near the model top. These are declared and exposed in the namelist, but their application in the vertical-redistribution logic is currently commented out in the source (marked cig commented out) — left in for reference rather than active in this version.

Parallelization fixes (Stephen Thomson, gwd_attempt branch, inherited via Ross's fork)

The version of cg_drag.f90 on master has two independent bugs that break multi-core runs — see ExeClim/Isca#293:

1. source_level/source_amp are allocated with local per-PE array bounds (allocate(source_level(idf,jdf))) but were accessed with an incorrect global offset added on top (source_level(i+is-1, j+js-1)), causing an out-of-bounds crash on master for NCORES > 1. Fixed by indexing with the local (i,j) directly (the buggy lines are left in as commented-out dead code for reference).
2. gwd_u/gwd_v (which cache the computed forcing between cg_drag_freq recalculation cycles) have the same kind of inconsistency on master: allocated with local bounds but accessed via a global-position slice, gwd_u(is:ie,js:je,:). This branch resolves it from the other end — gwd_u/gwd_v are instead allocated with global-position-labelled bounds (allocate(gwd_u(is:ie,js:je,kmax))), making the existing gwd_u(is:ie,js:je,:) access correct rather than needing to change the access pattern itself.

The #293 issue thread suggested a third code change (removing the is:ie,js:je slice from the gwd_u/gwd_v assignment) that the reporter needed for their own unfixed baseline, but wasn't sure was the right fix. We verified empirically that this branch's approach doesn't need it: ran cg_drag "on" at NCORES=1 vs NCORES=4 for a 1-day integration spanning multiple full cg_drag_freq (6h) recalculation cycles — udt_cgwd and ucomp output were bit-identical (max abs diff 0.0, every field, every one of 24 hourly snapshots) between core counts.

2. Orographic gravity wave drag — mg_drag

From Ross Castle's realistic_mgdrag branch. Reads subgrid-orography standard deviation via source_of_sgsmtn='input' — a field regridded from ERA5's sdor parameter — rather than source_of_sgsmtn='computed' (which needs raw high-resolution elevation data not available in this environment). Verified against ERA5 ground truth: peak subgrid variance correctly located over the Andes/Altiplano (~662 m), Himalaya (~342 m), zero over open ocean.

3. Bugs found and fixed while verifying both schemes

a) mg_drag.f90 restart I/O domain bug (this PR, commit a1d4b05f)

mg_drag's read_data/write_data calls for its Ghprime (subgrid-orography variance) restart field omitted the grid_domain argument that every other domain-decomposed field in Isca passes (c.f. hs_forcing.F90's tg_prev). Without it, write_data at the end of each month doesn't gather all PEs' local latitude bands into the true global field — it writes only one PE's local tile. The following month's read_data then loads that malformed file into every PE's local array (the shapes happen to coincidentally match under a uniform domain decomposition), so the "global" field becomes one PE's tile repeated across the latitude axis instead of the real topography.

This only manifests for source_of_sgsmtn='input' runs chained across multiple restart months — 'computed' mode recomputes fresh from raw elevation every month rather than reading back its own restart, so this code path was never previously exercised. Since Ghprime is the literal array used in the drag calculation (zsvar = Ghprime**2), this was a real physics bug affecting 5 of 6 months of any multi-month 'input'-mode integration, not just a diagnostic display issue.

Fixed by importing grid_domain from transforms_mod and adding call set_domain(grid_domain) plus grid_domain as the 4th argument to both the read_data and write_data calls, matching hs_forcing.F90's established pattern exactly. Verified with a 2-month smoke test before the full rerun: month 2 is now bit-identical to month 1 (previously diverged almost completely — spatial correlation ~0.02).

b) ERA5 topography regridding — latitude order (not in this PR's diff; fixed in isca-experiments)

The script that built the ERA5-derived mg_drag.res.nc input file (isca-experiments/gwd_merge_comparison/regrid_era5_sdor_to_t42.py) wrote its latitude coordinate descending, on the mistaken assumption that matched Isca's convention. It doesn't — Isca's own input files (e.g. era-spectral7_T42_64x128.out.nc) are ascending, and mg_drag_init's read_data call trusts the file's row order rather than its coordinate values (same underlying mechanism as (a) above). This silently mirrored the whole topography field north-south. Fixed alongside (a), verified the model's own diagnostic output now shows the topography peak at the correct hemisphere/location.

4. New test cases

Two short (6-month) on/off comparison experiments, isca-experiments/gwd_merge_comparison/:

• cg_drag_qbo_test.py — QBO-configuration non-orographic case (T42L50, RRTM radiation, Rayleigh damping off so cg_drag is the only on/off difference).
• mg_drag_socrates_topo_test.py — Socrates radiation + real topography orographic case (T42L40).

Both confirm the expected physics: all GWD tendency/stress diagnostics (udt_cgwd; udt_gwd, vdt_gwd, taubx, tauby, taus, tdt_diss_gwd, sgsmtn) are populated with physically sane, non-zero magnitudes when their scheme is on, and completely absent (not merely zero) when off — confirming the namelist toggles correctly gate both the physics and diagnostics.

Full verification details, zonal-mean wind and tendency comparison plots (interactive, final 2 months of each run with spin-up excluded), and narrative writeup: see attached gwd_comparison_report.html.

Testing

• Both IscaCodeBase (RRTM) and SocratesCodeBase compile cleanly.
• Both on/off test cases run to completion with no errors.
• mg_drag topography verified against ERA5 ground truth (correct hemisphere/location, zero over open ocean) after both fixes above.
• cg_drag verified bit-identical between NCORES=1 and NCORES=4.
• Branch merged forward from master (111 commits, unrelated work) with zero conflicts outside a single non-overlapping region of idealized_moist_phys.F90; both codebases recompiled cleanly after the merge.

[sit23]

As part of this work I had Claude generate some simple comparison plots of cg_drag on vs off and mg_drag on vs off. This can be found in the attached. This verified in a limited way that both schems are behaving sensibly and in the way that might be expected, but the test cases should still be treated with caution.
gwd_comparison_report.html

[sit23]

Here are the trip test results for this branch:

Trip test results — gwd_merge_orographic_nonorographic

Comparing origin/master (f48bfa1b) against this branch's tip (773a532e) using Isca's own trip test harness (exp/test_cases/trip_test/), covering all 21 registered test cases, including the two new cg_drag_qbo and mg_drag_socrates_topo cases added by this PR.

Summary

Result	Count	Test cases
✅ Pass — bit-identical	15	see below
✅ Fail as expected — proves the fix	2	cg_drag_qbo, mg_drag_socrates_topo
⚪ Fail — pre-existing, unrelated to this PR	4	axisymmetric, giant_planet, grey_mars, radiative_eq_mars

No regressions. Every test case that doesn't touch GWD code is bit-identical between master and this branch. The only two test cases that differ are the two this PR is about, and they differ in exactly the direction expected: broken on master, clean on this branch.

✅ Clean passes (15)

Bit-identical output between master and this branch:

bucket_model, frierson, held_suarez, MiMA, realistic_continents_fixed_sst, realistic_continents_variable_qflux, socrates_aquaplanet, socrates_aquaplanet_cloud, top_down_test, variable_co2_grey, variable_co2_rrtm, ape_aquaplanet, barotropic_vort_eq_stirring, shallow_water_stirring, column_test

(column_test required a fix to the trip test's own comparison logic to complete — see below. Once fixed, it too compares cleanly.)

✅ Expected failures — these are the point of the PR (2)

cg_drag_qbo

• master: fails to even start — FATAL: Unknown namelist, or mistyped namelist variable in namelist cg_drag_nml. Master's cg_drag_nml doesn't recognize kelvin_kludge, weighttop, weightminus1, weightminus2 — direct proof these are genuinely new parameters (Chaim Garfinkel's additions) this branch brings in, not something upstream already has.
• this branch: completely clean run, zero errors, restart archive created.

mg_drag_socrates_topo

• master: crashes with FATAL: lookup_es_1d: saturation vapor pressure table overflow. Master's old, unfixed mg_drag is numerically unstable under the tuned namelist this test case uses.
• this branch: completely clean run, zero errors, restart archive created.

⚪ Pre-existing failures, unrelated to this PR (4)

Each confirmed by checking that both master and this branch crash identically — same error, same location, same values — meaning the failure predates and is independent of this branch:

• axisymmetric — FATAL: interpolator_2D: You must call interpolator_init before calling interpolator, identical on both commits, same PE. Likely specific to running at NCORES=16 (used here for speed; the tool's own default is 4).
• giant_planet — FATAL: temperatures out of valid range, identical on both commits, same grid point, values matching to 3 decimal places. Same likely cause as above.
• grey_mars — FATAL: run length must be multiple of atmosphere time step, identical on both commits. This test case was added to the trip test suite via the 111 unrelated commits merged in from master before this PR's own work began; its dt_atmos doesn't evenly divide the trip test harness's hardcoded 3-day comparison window — a pre-existing harness/test-case mismatch.
• radiative_eq_mars — same run length error as grey_mars, identical on both commits, same root cause.

None of these four touch cg_drag, mg_drag, damping_driver, or idealized_moist_phys — the files this PR actually changes.

Fixes to the trip test harness itself, included in this PR

Two bugs in the trip test infrastructure were found and fixed while running this comparison (both pre-existing, unrelated to GWD, but blocking a clean run):

1. Column model build (cfb5983c) — mg_drag.f90 needed to be mg_drag.F90 for its #ifdef COLUMN_MODEL guard to actually work (mkmf only passes build macros to uppercase .F90 files). Fixed by renaming and updating the six path_names files that reference it.
2. Comparison crash on the column model's output (773a532e) — conduct_comparison_on_test_case subtracted two xarray Datasets directly, which aligns them by coordinate value first; the column model's lonb coordinate has two identical placeholder values (a 1-column grid has no real longitude extent), which pandas can't align on uniquely. Switched to a direct positional (.values) comparison, which is what's actually wanted here anyway since both datasets share the same grid by construction.

How to reproduce

cd exp/test_cases/trip_test
./trip_test_command_line f48bfa1b 773a532e -n 16

(Add -e cg_drag_qbo mg_drag_socrates_topo to check just the two new cases, or -e column_test for that one alone.)