Skip to content

Pull request to put the Isca Mars configurations into Isca's master - #299

Merged
sit23 merged 108 commits into
ExeClim:masterfrom
sit23:mars_dev_gfort_local_heating_2026
Aug 13, 2026
Merged

Pull request to put the Isca Mars configurations into Isca's master#299
sit23 merged 108 commits into
ExeClim:masterfrom
sit23:mars_dev_gfort_local_heating_2026

Conversation

@sit23

@sit23 sit23 commented Aug 4, 2026

Copy link
Copy Markdown
Contributor

The Isca Mars configurations haven't ever been merged into the master, and I'd like to make sure this is done ASAP.

In this P/R I have put in my grey-radiation Mars and Socrates Mars configurations without dust as test cases that can be run with the trip tests. These are the ones that are described in this paper:

https://www.mdpi.com/2073-4433/10/12/803

I'd also like to get in Emily Ball's Isca-Mars configurations that built on my Mars stuff, but I wanted to do this first. This P/R is ready for merging (pending passing the trip tests).

sit23 and others added 30 commits May 11, 2018 12:18
…my module in order that we can correctly calculate rrsun using the true anomaly rather than the mean anomaly. Have also added multiplication by rrsun to two-stream-grey-rad, which makes no difference to results in initial tests when ecc = 0.0
…ue anomaly, and therefore be able to do mars time-telling.
… expects orbital_period in seconds, but I had supplied it in days. This is because Alex's code expects it in days. Will fix Alex code to make it consistent.
…e of surface optical depths and albedos, assuming transparancy in the visible. Useful for working out rough values for optical depths and albedos for a simple mars model.
…ay is not the same as one mars day. Have altered rotation rate and orbital rate such that each can be set as an integer, and we end up with a length of sol which is also an integer. This can then be used as the averaging period.
…dependently. Then also only calculating rh when rh is asked for. Model now runs and looks vaguely mars-like, but dates of equinox etc still not right.
…nml, and ability to use specified temperatures for evaporation calculations in surface flux. This means that, when using a dry model where temperatures are outside standard range for sat specific humidity calculations, then false temps can be used to stop it failing.
Cores on a Mac cannot be tied in the same way that they can on linux, so core affinity is not possible.  Here I am just stubbing out the core affinity functions to mimic the linux case where it cannot get the affinity info.

More info in links inline in `affinity.c`.
…xis. For some reason, the interpolator will get rid of such variables when used on its own, so I modified the python to add them back.
…o idealized moist phys back in teh old GFDLMoistmodel repo on the local_heating_dev branch. This seems to work now, and should hopefully be useful. Only tried the Isidoro option, and not the input file option, but this should be fine.
Adding mac functionality to local heating branch
…cal prescribed heating input files. Seems to work alright, but problem is that local heating code always reads zeros no matter what I seem to do. Reading the created files into RRTM as an ozone does work, and reading ozone-1990 into local heating does not. Very odd. Tried everything I can think of, but we'll have to carry on testing.
…not being passed time, and therefore was not being fed to the correct interpolator within the interface type structure.
Conflicts:
	src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90
sit23 added 3 commits August 4, 2026 11:32
…p dry heating test case and associated script for generating heating rate. This now produced a heating array of zeros purely for test purposes.
sit23 and others added 17 commits August 4, 2026 13:38
…ase. Confirmed results bitwise identical, and have made it so that the max number of processors used for parallel compilation is 8, as this gives 32% speedup on serial compilation, but more cores does not give compilation speedup.
…oducible as random seed was not set. Now added as namelist parameter for both test cases, and verified that runs are reproducible.
…so non-Titan cases stay bit-identical to master

Both update_tracers's grid-tracer branch (spectral_dynamics.F90) and the
bucket-water diffusion call (idealized_moist_phys.F90) unconditionally ran a
grid->spherical->grid round-trip that was only meant to support new,
opt-in Titan features (do_spec_tracer_filter, damping_coeff_bucket damping).
Even with those options left at their defaults (off), the round-trip isn't a
no-op - it silently low-pass-filters real grid-scale structure (e.g. via
floating-point reassociation in the tracer tendency, or truncation of
land/ocean boundaries in the bucket case) - which broke bit-reproducibility
against master for frierson, bucket_model, and realistic_continents_variable_qflux.

Both call sites are now branched on their respective flags: the Titan path
runs unchanged when requested, and everything else falls back to the exact
pre-existing formulation. Verified bit-identical to master (eb64615) for
frierson, bucket_model, and realistic_continents_variable_qflux.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
7f4201d ("Important fix for two-stream") corrected frac_of_day to account
for total elapsed seconds rather than just the current day's remainder,
needed for solar day lengths > 86400s (Mars, Titan) where the model's
internal "days" bookkeeping (fixed 86400s calendar units) and the actual
solar day length diverge. Master never needed this fix and still uses the
original formulation.

For day_in_s == 86400 (Earth-standard), the original formulation was
already correct, but the code ran the total-elapsed-seconds path
unconditionally. Since that path needs mod() to strip the accumulated whole-days
component back out, it introduces a small precision loss vs directly using
r_seconds - breaking bit-reproducibility against master for any do_seasonal
Earth-day-length case (e.g. variable_co2_grey).

Branch on day_in_s == 86400 instead: Earth-standard runs get the exact
pre-7f4201dd formulation back, non-Earth day lengths keep 7f4201d's fix
unchanged. Verified bit-identical to master for variable_co2_grey (and
frierson, bucket_model together).

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
@sit23

sit23 commented Aug 12, 2026

Copy link
Copy Markdown
Contributor Author

Having put all of the Mars-related modifications into this P/R, nearly all of the trip tests. However, I worked with Claude to find the causes of these trip test failures, and have identified and fixed all of them. The summary of this investigation (generated by Claude) is given below:

Summary

Three genuine numerical regressions between mars_dev_gfort_local_heating_2026 and master were tracked down and fixed. All three shared the same shape: a Titan-motivated code change was added as the only code path instead of being gated behind its own (already-existing, default-off) namelist option, so it silently ran for every experiment — including ones that never asked for it — and broke bit-reproducibility against master.

  • frierson and realistic_continents_variable_qflux diverged from master (tiny for frierson, substantial for realistic_continents_variable_qflux — ps off by up to 1.7).
  • bucket_model diverged substantially (ps off by 4.5, temp by 0.5).
  • variable_co2_grey diverged moderately (ps off by 0.023).

All three are now bit-identical to master. Socrates-based trip tests (socrates_aquaplanet, socrates_aquaplanet_cloud, ape_aquaplanet), previously untestable locally due to an unrelated environment issue, now run and pass as well.

Fix 1 — spectral_dynamics.F90: grid-tracer accumulation order (commit 3cdb2e93)

update_tracers's grid-tracer branch was restructured to support a new Titan option (do_spec_tracer_filter, spectral damping of the grid tracer) by accumulating the horizontal+vertical advection tendencies into dt_tr and recomputing tr_future from it once at the end, instead of incrementing tr_future directly after each advection step. Mathematically equivalent, but a different floating-point summation order — and that changed bit-level results for every moist run, not just Titan ones, since do_spec_tracer_filter defaults to .false..

Fix: branch on do_spec_tracer_filter. The new accumulate-then-filter path runs unchanged when the option is actually requested; everything else falls back to the original direct-increment formulation, byte-for-byte.

Fix 2 — idealized_moist_phys.F90: unconditional bucket-water spectral diffusion (commit 3cdb2e93)

Same underlying pattern, bigger effect. The bucket-water tendency block gained an unconditional call to a new diffuse_surf_water() subroutine, regardless of damping_coeff_bucket (which defaults to 0., i.e. "no damping requested"). The subroutine round-trips the tendency field through a spherical harmonic transform and back — and even at zero damping, that round-trip isn't a no-op: it truncates the field to the spectral resolution, silently low-pass-filtering real grid-scale structure (e.g. land/ocean boundaries in bucket_depth). That's a genuine physical change, not roundoff, which is why bucket_model's divergence was far larger than the others.

Fix: only call diffuse_surf_water() when damping_coeff_bucket > 0.

Fix 3 — two_stream_gray_rad.F90: seasonal frac-of-day precision loss (commit ea8e155b)

Different in kind from the first two: not an accidental regression. Commit 7f4201dd ("Important fix for two-stream") correctly fixed frac_of_day to account for total elapsed seconds rather than just the current day's remainder — needed for solar day lengths > 86400s (Mars, Titan), where the model's internal "days" bookkeeping (fixed 86400s calendar units) and the actual solar day length diverge. Master never needed this and still uses the original formulation. For the standard Earth day length, the new formulation still runs and still gets the right answer, but only after mod()-ing away an accumulated whole-days component — introducing a small but real precision loss relative to just using the bounded seconds-of-day value directly. That's what broke variable_co2_grey (which sets do_seasonal: True) against master.

Fix: branch on day_in_s == 86400.. Earth-standard day length gets the exact pre-7f4201dd formulation back, bit-for-bit; non-Earth day lengths keep 7f4201dd's correction unchanged, preserving correct Mars/Titan seasonal timing.

Socrates trip tests

Previously blocked locally by an unrelated environment issue: GFDL_SOC points to a nonexistent path, and scratch codebase clones (which don't inherit the untracked trunk symlink our real checkout already has set up correctly) picked up a dangling symlink to that broken location. No code change — just needed trunk repointed at the real local Socrates source per clone, same method already in use in this checkout. All three Socrates test cases (socrates_aquaplanet, socrates_aquaplanet_cloud, ape_aquaplanet) pass against master.

Known remaining issue (not fixed in this PR)

top_down_test shows a small, stable divergence from master (ucomp ~1.2e-7, vor ~1.1e-13) that is not caused by any of the above — confirmed by reverting all candidate files from the relevant merge and finding no change. Bisecting it further is blocked by an unrelated double free or corruption crash in old commits (likely in mixed_layer.F90's heat-capacity spinup code), which would need real memory debugging rather than a mechanical fix to work around. Given this test case sees little use, it's left as the one known open item rather than pursued further here.

Test plan

  • frierson vs master: bit-identical
  • bucket_model vs master: bit-identical
  • realistic_continents_variable_qflux vs master: bit-identical
  • variable_co2_grey vs master: bit-identical
  • socrates_aquaplanet, socrates_aquaplanet_cloud, ape_aquaplanet vs master: bit-identical
  • held_suarez, MiMA, realistic_continents_fixed_sst, variable_co2_rrtm, column_test vs master: unaffected, still bit-identical
  • top_down_test vs master: still diverges (pre-existing, unrelated, documented above — not addressed in this PR)

@sit23

sit23 commented Aug 12, 2026

Copy link
Copy Markdown
Contributor Author

The plan for this branch now is to attempt to track down the cause of the top-down test divergence bug. If this can be found, a fix will be applied, and then this branch can be merged.

…actually set

update_orbit's theta calculation gained an equinox_day offset and a modulo()
wraparound to support a non-default calendar alignment. equinox_day defaults
to 0, and theta is 2*pi-periodic either way (sin(theta) is invariant to
adding whole orbits), so the new formulation is mathematically equivalent
for the default case - but it reorders the multiply/divide and adds a
wrap-around that's never actually needed there, which isn't bit-identical
to the original even when no wrapping occurs. This broke bit-reproducibility
against master for any equilibrium_t_option='top_down' run (e.g.
top_down_test), the only trip test case that exercises update_orbit.

Branch on equinox_day == 0: the default case gets the exact original
formulation back, bit-for-bit; a non-zero equinox_day keeps the new
wraparound behavior unchanged.

Verified bit-identical to master for top_down_test, together with frierson,
bucket_model, variable_co2_grey, and realistic_continents_variable_qflux.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
@sit23

sit23 commented Aug 13, 2026

Copy link
Copy Markdown
Contributor Author

Claude found the fix for top-down test:

Update: found and fixed the last remaining divergence (top_down_test)

Since the description above, the last known open item — top_down_test's small divergence from master (ucomp ~1.2e-7, vor ~1.1e-13) — has also been found and fixed.

Fix 4 — hs_forcing.F90: orbital-phase wraparound (commit a89b5c83)

top_down_test is a dry test case (idealized_moist_model=False) using equilibrium_t_option='top_down', a code path (hs_forcing.F90, update_orbit) none of the other three fixes touch. update_orbit's orbital-phase angle theta gained a new equinox_day offset (default 0.) and a modulo() wraparound:

- theta = 2*pi*current_time/(orbital_period*86400)
+ theta = 2*pi*modulo((current_time/(orbital_period*86400))-equinox_day, 1.0)

Same class of issue as Fix 1: theta is 2π-periodic either way (sin(theta) is invariant to adding whole orbits), so the two formulations are mathematically equivalent at the default equinox_day=0 — but the new one reorders the multiply/divide and adds a wraparound that's never actually needed there, which isn't bit-identical to the original even when no wrapping occurs.

Fix: branch on equinox_day == 0.. The default case gets the exact original formulation back, bit-for-bit; a non-zero equinox_day keeps the new wraparound behavior unchanged.

What this means for the PR

Every genuine numerical divergence found during this investigation is now fixed. The complete local trip test suite was re-run one final time against the finished branch head with all four fixes together: every test that was ever a real divergence now passes, and the only remaining failures are a cluster of pre-existing, unrelated crashes (see updated test plan below) that have nothing to do with this PR — a structural mismatch from comparing master's older Fortran against this branch's current Python test-case scripts, not a code bug.

Updated test plan

  • frierson vs master: bit-identical
  • bucket_model vs master: bit-identical
  • realistic_continents_variable_qflux vs master: bit-identical
  • variable_co2_grey vs master: bit-identical
  • top_down_test vs master: bit-identical (new — previously the one open item)
  • socrates_aquaplanet, socrates_aquaplanet_cloud, ape_aquaplanet vs master: bit-identical
  • held_suarez, MiMA, realistic_continents_fixed_sst, variable_co2_rrtm, column_test vs master: unaffected, still bit-identical
  • Full local trip test suite re-run against the final branch head with all four fixes together: confirms no regressions
  • axisymmetric, giant_planet, barotropic_vort_eq_stirring, shallow_water_stirring, grey_mars, radiative_eq_mars vs master: still crash — pre-existing, unrelated namelist/version-skew mismatch between master and this branch's current test-case scripts, not addressed in this PR (not a numerical divergence)

Final trip test results (all cases, vs master)

Test case Result Notes
frierson ✅ Pass Fixed in this PR (Fix 1)
bucket_model ✅ Pass Fixed in this PR (Fix 2)
realistic_continents_variable_qflux ✅ Pass Fixed in this PR (Fix 1)
variable_co2_grey ✅ Pass Fixed in this PR (Fix 3)
top_down_test ✅ Pass Fixed in this PR (Fix 4)
held_suarez ✅ Pass Unaffected, always passed
MiMA ✅ Pass Unaffected, always passed
realistic_continents_fixed_sst ✅ Pass Unaffected, always passed
variable_co2_rrtm ✅ Pass Unaffected, always passed
column_test ✅ Pass Unaffected, always passed
socrates_aquaplanet ✅ Pass Previously untestable locally (env issue, not a code bug); now runs and passes
socrates_aquaplanet_cloud ✅ Pass Previously untestable locally (env issue, not a code bug); now runs and passes
ape_aquaplanet ✅ Pass Previously untestable locally (env issue, not a code bug); now runs and passes
axisymmetric ❌ Crash Pre-existing, unrelated - master-side namelist/version-skew mismatch
giant_planet ❌ Crash Pre-existing, unrelated - master-side crash
barotropic_vort_eq_stirring ❌ Crash Pre-existing, unrelated - master's Fortran rejects stirring_nml option only the branch's current test scripts write
shallow_water_stirring ❌ Crash Pre-existing, unrelated - same stirring_nml mismatch as above
grey_mars ❌ Crash Pre-existing, unrelated - master's Fortran rejects constants_nml options only the branch's current test scripts write
radiative_eq_mars ❌ Crash Pre-existing, unrelated - same constants_nml mismatch as above

19 test cases total: 13 pass (5 newly fixed by this PR, 3 previously blocked by an unrelated local env issue, 5 always passing), 6 fail with pre-existing crashes unrelated to this PR.

@sit23

sit23 commented Aug 13, 2026

Copy link
Copy Markdown
Contributor Author

Once the CI passes, I'm happy and ready to merge this in. After that, for Mars things, I'm going to move onto merging in Emily Ball's Isca-Mars work so that we have a test case for that too.

@sit23

sit23 commented Aug 13, 2026

Copy link
Copy Markdown
Contributor Author

Here's Claude's summary of the contents of this P/R:

High-level summary

This PR brings several years of accumulated Mars and Titan development work into master: new radiation schemes, a generalized (non-water) condensible framework, orbital mechanics improvements needed for eccentric and non-Earth-length orbits, an alternative surface-flux scheme, local/prescribed heating support, four new trip test cases, and a round of gfortran/compiler-compatibility and build-performance work. It also includes four targeted bit-reproducibility fixes (documented separately) that make every new capability opt-in rather than silently changing results for existing, unrelated experiments.

New physics capabilities

Mars radiation and dynamics

  • Grey-gas Mars radiation scheme, built out from an early prototype into a working, tested configuration (two_stream_gray_rad.F90).
  • Mars-specific rotation/orbital rate handling so the length of a Martian sol comes out as a clean integer, usable directly as an averaging period.
  • MOLA Mars topography import script (mars_topo_maker.py).
  • N₂ added as a selectable background-gas option for Socrates (Mars' atmosphere is mostly CO₂/N₂), plus new Martian diagnostics in the Socrates interface.

Titan radiation and surface

  • Titan-specific grey radiation scheme (from Schneider et al. 2012), including an option to override net SW flux with measured Huygens probe data.
  • Generalized the condensible-species framework so it isn't hard-coded to water: saturation vapor pressure calculation and the vapor-related constants (RVGAS, HLV, HLF, TFREEZE, ES0, DENS_VAPOR, TRIPLE_POINT_PRES) are now namelist-configurable, enabling methane as a working fluid.
  • Spectral (grid-tracer) damping option (do_spec_tracer_filter) and surface-liquid diffusion (diffuse_surf_water, damping_coeff_bucket), both from Tapio Schneider's Titan model, for numerically stabilizing methane transport and surface reservoirs.

Local / prescribed heating

  • Resurrected and integrated a local prescribed-heating capability into idealized_moist_phys.F90 (local_heating_option), supporting both an analytic ("Isidoro") profile and reading a heating rate from file.
  • New script for generating local-heating input files (create_local_heating_input_file.py).

Alternative surface flux scheme

  • Added the Frierson Monin-Obukhov surface-flux formulation from MiMA as a new, selectable module (frierson_monin_obukhov.F90, use_frierson_mo_drag) — a simplified alternative to the full Monin-Obukhov similarity calculation.

Orbital mechanics / astronomy

  • Proper eccentric- and true-anomaly calculation for orbital distance (rrsun), replacing a mean-anomaly approximation — needed for accurate insolation on eccentric orbits like Mars'.
  • equinox_day added to hs_forcing, making it straightforward to align the model's seasonal cycle with an Earth-like calendar.
  • Fixed a genuine bug in the two-stream radiation scheme's day-fraction calculation that broke insolation timing (fixed sub-stellar longitude) for solar days longer than 86400s — relevant to both Mars and Titan.
  • Fixed a diurnal_solar segfault triggered when called without its optional declination/hour-angle arguments.

New test cases

Test case Description
grey_mars Full Mars configuration using the grey-gas Mars radiation scheme
radiative_eq_mars Dry radiative-equilibrium Mars experiment, from the Isca Atmospheres paper
socrates_mars Mars configuration using the Socrates radiation scheme (N₂/CO₂ atmosphere)
frierson_dry_heating Frierson-configuration test exercising the new local/prescribed heating capability

All four are wired into the trip test suite alongside the existing cases.

Build, compiler compatibility, and performance

  • gfortran compatibility fixes across the codebase (this branch's namesake): logical-comparison operators (==.eqv.), glibc/gettid conflicts in affinity.c, and related toolchain issues.
  • Mac (Homebrew) and Docker build environments added.
  • Parallel compilation support, capped at 8 cores by default (~32% faster serial-to-parallel speedup on the barotropic test case), verified to produce bitwise-identical results to serial compilation.

Test infrastructure improvements

  • Test case scripts renamed (*_test.py*_test_case.py) and a pytest.ini added, so pytest no longer auto-discovers and tries to run them as unit tests.
  • column_test moved from the automated compile-test suite into the trip test suite, alongside the Mars test cases.
  • Stirring test cases (barotropic and shallow-water) were non-reproducible because their random seed wasn't set; a namelist parameter now fixes this, and reproducibility has been verified.
  • Various Socrates test-case configuration and compilation bugs fixed (missing codebase selection, missing path-names entries for new modules, socrates namelist initialization).

Bit-reproducibility fixes (this investigation)

Four of the new, default-off features above were being exercised unconditionally rather than only when requested, which silently changed results for unrelated experiments and broke bit-reproducibility against master. All four are now correctly gated behind their own namelist options, with the affected trip test cases (frierson, bucket_model, realistic_continents_variable_qflux, variable_co2_grey, top_down_test) verified bit-identical to master. Full details are in the separate fix write-up and PR update comment.

@sit23

sit23 commented Aug 13, 2026

Copy link
Copy Markdown
Contributor Author

OK - I'm going to merge this in now. The summary being that Isca Mars configurations as in the atmospheres paper that Geoff and I had are now available in the Isca master.

@sit23
sit23 merged commit f48bfa1 into ExeClim:master Aug 13, 2026
1 check passed
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.

2 participants