Porting SST and ice physics improvements from uob_fftw_sit23 branch

PHYSICS_PORT_VALIDATION.md

@@ -0,0 +1,253 @@
+# Physics port validation: `uob_fftw_sit23_physics_2026`
+
+This documents validation of the sea-ice/land-SST physics port from
+`uob_fftw_sit23_new_maths2intel` onto current `master`, done on branch
+`uob_fftw_sit23_physics_2026`. Two checks are covered:
+
+1. Does an existing production experiment that depends on this physics
+   (the Maunder Minimum q-flux run) still compile and run on the new branch?
+2. Does it produce the same output as the original, un-ported branch?
+
+## What was ported
+
+Branch `uob_fftw_sit23_physics_2026` (commit `75044a0a`) carries 6 commits
+cherry-picked from `uob_fftw_sit23_new_maths2intel` (commit `a9dc5390`) onto
+current `master`, with original authorship preserved:
+
+- 4 commits adding sea-ice/land-SST patching options to
+  `src/atmos_spectral/driver/solo/mixed_layer.F90` (`specify_sst_over_sea_ice`,
+  `linearly_interpolate_sea_ice_temp_and_sst`, `ice_sst_file`,
+  `specify_sst_over_land_from_separate_file_to_ocean_sst`, `land_sst_file`,
+  `allow_qflux_over_land`, `update_land_mask_from_ice`), hand-reconciled
+  against master's independent ice-albedo work added since the branch forked.
+- 1 compatibility commit, preserving existing bit-reproducibility for
+  `realistic_continents`/`socrates_simcloud` (see below).
+- 1 commit porting `calculate_qflux.py`'s Socrates support and
+  `ignore_ice_for_calculation` option from the original branch.
+
+One deliberate change during the port: the original branch's `binary_ice_albedo`
+switch was dropped, because master independently added an equivalent (and more
+general) `ice_albedo_method` option (`'step_function'` / `'ramp_function'`) in
+the years since the branch forked. The two are algebraically identical for the
+non-default case (`binary_ice_albedo=False` ≡ `ice_albedo_method='ramp_function'`),
+so this is a rename, not a behaviour change.
+
+## Bit-reproducibility of existing test cases
+
+Isca's `trip_test` suite couldn't be used directly to validate the
+`update_land_mask_from_ice` compatibility fix, because it feeds one shared
+namelist file to two different compiled binaries, and the fix specifically
+requires a namelist key the old binary doesn't know about (see the
+"how this was tested" note in the appendix). A standalone comparison giving
+each commit its own correct namelist confirmed **bit-identical** output
+(all fields, 0.0 max abs diff) against `master` for `realistic_continents_fixed_sst`
+and `realistic_continents_variable_qflux` with `update_land_mask_from_ice`
+set explicitly, as well as `frierson`, `held_suarez`, and `bucket_model`
+passing `trip_test` unmodified.
+
+## Does the Maunder Minimum experiment run?
+
+**Yes.** `exp/socrates_fast_land_ice_cmpi5_ozone/maunder_min/socrates_fast_land_fft_new_ice_sbm_do_simple_false_bucket_qflux_mm.py`
+compiles (`SocratesCodeBase`) and runs to completion on
+`uob_fftw_sit23_physics_2026`, continuing from the experiment's existing
+restart file (`input/res0240.tar.gz`, run 240 → 241). One namelist
+adaptation was needed, matching the `binary_ice_albedo` retirement above:
+
+```diff
+- 'binary_ice_albedo': False,
++ 'ice_albedo_method': 'ramp_function',
+```
+
+`cb.enable_fftw3()` remains commented out in the script, as it is on the
+original branch checkout used here — the FFTW/numerics port is a separate,
+not-yet-started piece of work (`uob_fftw_sit23_numerics_2026`), so this
+validation is physics-only and doesn't exercise FFTW either way.
+
+## Does it produce the same output as the original branch?
+
+A direct comparison was run: both `uob_fftw_sit23_new_maths2intel` (the
+original, un-ported branch) and `uob_fftw_sit23_physics_2026` compiled and ran
+**the same experiment, from the same restart file, with the same namelist**
+(modulo the `binary_ice_albedo`→`ice_albedo_method` substitution above)
+for one month (run 241, 30 days).
+
+**Short answer: not bit-identical, and that's expected, not a bug.** Isca is a
+chaotic dynamical system — even a rounding-level difference anywhere in the
+computation (here, unavoidable given the two branches were compiled against
+today's toolchain/library versions rather than whatever was current when the
+original branch was last touched in 2023) amplifies over a month of
+integration into large-looking differences in the instantaneous weather
+state, with no error in the physics involved. Two pieces of evidence support
+that this is exactly what's happening here, not a real regression:
+
+**1. Fields that don't evolve chaotically match exactly.** `albedo`,
+`ice_conc`, `flux_oceanq`, `bk`, `pk`, and `zsurf` — the fields most directly
+exercised by the ported physics (ice-concentration interpolation, the q-flux
+input file) or otherwise fixed/diagnostic — have **0.0 max absolute
+difference** after 30 days. If the port had gotten any of this physics wrong,
+these are exactly the fields that would show it, and they don't.
+
+**2. The dynamical fields diverge with a textbook chaotic-error-growth
+signature, not a step change.** A separate 5-day run with daily output shows
+the difference between the two branches starting essentially at the noise
+floor and growing roughly exponentially:
+
+| Day | max|Δps| (Pa) | max|ΔT| (K) | max|Δu| (m/s) | max|Δv| (m/s) |
+|---|---|---|---|---|
+| 1 | 0.61  | 0.050 | 0.061 | 0.056 |
+| 2 | 7.6   | 0.454 | 0.663 | 0.403 |
+| 3 | 3.5   | 0.443 | 1.007 | 1.021 |
+| 4 | 6.3   | 0.370 | 1.513 | 1.212 |
+| 5 | 12.0  | 0.465 | 3.129 | 2.885 |
+
+Day 1's differences (order 10<sup>-2</sup>&ndash;10<sup>-1</sup> in physical
+units, i.e. parts-per-million relative to typical field magnitudes) are
+consistent with ordinary floating-point rounding differences between two
+different compiles, not a physics discrepancy. By day 30 this has amplified,
+via normal atmospheric chaos, into the full-field differences summarised
+below — large in places, but the same kind of divergence you'd see comparing
+any two bit-non-identical runs of the same physical configuration.
+
+### Full 30-day field comparison (`new_2026` &minus; `old_branch`)
+
+| Field | max&#124;diff&#124; | relative | | Field | max&#124;diff&#124; | relative |
+|---|---|---|---|---|---|---|
+| ps | 916.9 Pa | 0.9% | | flux_lhe | 52.9 W/m&sup2; | 7.4% |
+| temp | 5.03 K | 1.6% | | flux_t | 45.6 W/m&sup2; | 15.8% |
+| ucomp | 14.26 m/s | 18.9% | | soc_surf_flux_lw | 13.80 W/m&sup2; | 8.0% |
+| vcomp | 11.54 m/s | 24.7% | | soc_surf_flux_sw | 2.94 W/m&sup2; | 0.9% |
+| t_surf | 1.74 K | 0.6% | | soc_olr | 4.88 W/m&sup2; | 1.5% |
+| temp_2m | 1.85 K | 0.6% | | soc_toa_sw | 2.49 W/m&sup2; | 0.6% |
+| rh | 17.2 % | 15.8% | | soc_ozone | 6.3&times;10<sup>-8</sup> | 0.4% |
+| precipitation | 4.8&times;10<sup>-5</sup> | 19.7% | | vor | 2.1&times;10<sup>-5</sup> | 33.0% |
+| **albedo** | **0.0** | **0.0%** | | div | 1.5&times;10<sup>-5</sup> | 49.8% |
+| **ice_conc** | **0.0** | **0.0%** | | **bk / pk / zsurf** | **0.0** | **0.0%** |
+| **flux_oceanq** | **0.0** | **0.0%** | | | | |
+
+(Bolded rows are the ones most directly tied to the ported physics; all
+match exactly.)
+
+## Where the divergence actually comes from
+
+The chaotic-growth explanation above is correct as far as it goes, but it
+doesn't answer the sharper question: *is the divergence a generic property of
+comparing any two Isca compiles, or is it caused by something specific?* This
+was checked directly, using a plain `socrates_aquaplanet` control case (no
+ice/land-SST physics involved at all) run the same way (1 month +
+a 5-day daily-output growth check) across several pairings:
+
+| Comparison | Result |
+|---|---|
+| `master` vs. `uob_fftw_sit23_physics_2026` | **Bit-identical.** 0.0 difference in every field, every day, for the full month. |
+| Merge-base (2023) vs. current `master` | **Bit-identical**, bar one ~10<sup>-7</sup>-relative rounding blip in a diagnostic-only flux field that never reaches any dynamical field. |
+| Original `uob_fftw_sit23_new_maths2intel` vs. `physics_2026` | **Diverges** — same signature as the Maunder Minimum comparison above. |
+
+So `master` has been essentially perfectly reproducible for three years and
+372 commits, and the physics port doesn't disturb that at all. The divergence
+is real, and it's specific to the original branch.
+
+**Root cause, isolated by bisection:** every one of the ~35 source files the
+original branch modifies (relative to its merge-base) was reverted to the
+merge-base version in a scratch copy of the branch, one at a time / in
+batches, re-testing against `physics_2026` after each step, until the
+divergence disappeared and reappeared cleanly. Two environment-only patches
+(the `affinity.c` glibc fix and the `spectral_dynamics.F90` format-string fix,
+both already on `master`) were re-applied throughout purely so the branch
+would compile on this machine's current toolchain - they're unrelated to the
+result. This narrowed the entire divergence to a single file:
+**`src/atmos_param/socrates/interface/socrates_interface.F90`** — the
+branch's own (never-merged, unrelated to this port's `mixed_layer.F90`
+changes) modifications to the Socrates radiation coupling code. Reverting
+just that one file, with everything else at the branch's tip, restores
+bit-identical output; restoring just that one file, with everything else at
+merge-base, reproduces the full divergence on its own.
+
+**Why that file specifically:** it isn't one subtle line. Diffing the
+original branch's `socrates_interface.F90` directly against `master`'s
+current version (not against the shared merge-base) turns up **~2400 changed
+lines** out of a ~1650-line file — `master` has independently grown this file
+far beyond where the branch left it, including a full reimplementation of the
+same Mars-orbital diagnostics (`mars_solar_long`, `true_anomaly`, `rrsun`,
+`time_since_ae`, `dec`, `ang`) the branch added back in 2019, plus clear-sky
+diagnostics and SPOOKIE cloud support the branch never had at all. Notably,
+`master`'s version of the shared `astronomy.f90` machinery behind those
+diagnostics explicitly initializes `use_old_r_inv_squared = .TRUE.`, fixing
+an uninitialized-variable bug present in the branch's original version of the
+same logic. In short: this file's divergence is two genuinely different,
+independently-matured implementations of overlapping functionality, not a
+bug to hunt down to one line — and `master`'s is the more complete, already-
+fixed one. This doesn't affect this PR either way:
+**`uob_fftw_sit23_physics_2026` never touches `socrates_interface.F90` at
+all** — confirmed byte-identical to `master`'s current version.
+
+## Confirming the physics port directly, with that confound removed
+
+The comparison above answers "does the whole branch reproduce the whole
+branch's own experiment" — useful, but conflates the `socrates_interface.F90`
+divergence (irrelevant to this port) with the actual physics being merged. A
+cleaner test isolates just the latter: a scratch copy of the original branch
+had every file reverted to match current `master` **except**
+`mixed_layer.F90`, which was left at the branch's own original (un-ported)
+ice/land-SST physics. That "patched" branch was then run through the same
+Maunder Minimum experiment and restart file as above, and compared against
+`physics_2026`.
+
+**Result: bit-identical.** Every field, every day, for the full month —
+`ps`, `temp`, `ucomp`, `vcomp`, `t_surf`, `albedo`, `ice_conc`,
+`flux_oceanq`, all diagnostics, all exactly 0.0 difference, except
+`soc_surf_flux_lw`/`soc_surf_flux_lw_down` at ~10<sup>-7</sup> relative — the
+same pure floating-point noise level seen in the `master`-self-consistency
+check above, not chaotic amplification. This is the direct confirmation that
+the physics actually ported to `uob_fftw_sit23_physics_2026` reproduces
+exactly what the original branch's `mixed_layer.F90` would have produced
+running on top of current `master`, with no gap in what was carried over.
+
+## Bottom line
+
+- The physics port compiles and runs the real Maunder Minimum experiment
+  successfully on top of current `master`.
+- With the original branch patched to match current `master` everywhere
+  except `mixed_layer.F90`, the Maunder Minimum experiment is
+  **bit-identical** between that patched branch and `physics_2026` for the
+  full month — direct confirmation that the ported physics is complete and
+  correct, with no gap between what the branch had and what got carried over.
+- `master` and `physics_2026` are bit-identical to each other on cases the
+  port doesn't touch, and master has been bit-identical to its own 2023
+  merge-base for three years — so the port introduces no detectable
+  regression of its own.
+- The *unpatched* branch's divergence from `physics_2026` is real, but
+  traced to that branch's own `socrates_interface.F90` — which has
+  independently diverged from `master` by ~2400 lines (`master`'s own,
+  more complete reimplementation of the same Mars-diagnostics work, plus
+  clear-sky/SPOOKIE support the branch never had) — a file this port never
+  touches and that isn't part of what's being merged here.
+
+## Appendix: how this was tested
+
+Both branches were compiled and run directly (not via `trip_test`, whose
+shared-namelist-across-commits design can't represent a fix that adds a new
+namelist key only one side understands) via a small standalone script that:
+imports the experiment script's namelist/diagnostics/input-file configuration
+once, builds a `SocratesCodeBase` against each branch's worktree, applies the
+`ice_albedo_method` substitution only for the new branch, compiles, and runs
+`exp.run(241, use_restart=True, ...)` from the existing `res0240.tar.gz`
+restart file, then diffs the resulting `atmos_monthly.nc`/`atmos_daily.nc`
+with `xarray`.
+
+Getting `uob_fftw_sit23_new_maths2intel` to compile at all on this machine's
+current toolchain (gfortran 12.4.0, three years newer than the branch's last
+commit) required two environment-compatibility patches, applied locally to
+the test worktree only, not committed to either branch, and unrelated to
+physics:
+
+- `src/shared/mpp/affinity.c`: the branch's own `static gettid()` declaration
+  conflicts with glibc &ge;2.30's own `gettid()` — master already carries a
+  version-guarded fix for this.
+- `src/atmos_spectral/model/spectral_dynamics.F90`: two `FORMAT` statements
+  (json-logging labels 300/400) are missing a comma between edit descriptors,
+  tolerated by older compilers but a hard runtime error under this gfortran —
+  master already carries the fix (commit `36e7398b`).
+
+Both fixes are already present on `master`/`uob_fftw_sit23_physics_2026`
+through independent, unrelated work; they only needed re-applying to the old
+branch to get a fair comparison running on this specific machine.

exp/test_cases/realistic_continents/namelist_basefile.nml

@@ -58,6 +58,7 @@
     albedo_value = 0.25,                !Ocean albedo value
     land_albedo_prefactor = 1.3,        !What factor to multiply ocean albedo by over land    
     update_albedo_from_ice = .true.,    !Use the simple ice model to update surface albedo
+    update_land_mask_from_ice = .true., !Keep pre-existing behaviour: land_ice_mask includes ice-covered cells
     ice_albedo_value = 0.7,             !What value of albedo to use in regions of ice
     ice_concentration_threshold = 0.5/ !ice concentration threshold above which to make albedo equal to ice_albedo_value        
 

exp/test_cases/simple_clouds/socrates_simcloud.py

@@ -247,6 +247,7 @@
         'specify_sst_over_ocean_only': True,    # Make sure sst only specified in regions of ocean.
         # Copy from realistic_continents namelist
         'update_albedo_from_ice': True,         # Use the simple ice model to update surface albedo
+        'update_land_mask_from_ice': True,      # Keep pre-existing behaviour: land_ice_mask includes ice-covered cells
         'ice_albedo_value': 0.7,                # What value of albedo to use in regions of ice
         #'ice_concentration_threshold': 0.5,    # ice concentration threshold above which to make albedo equal to ice_albedo_value       
         'ice_albedo_method': 'ramp_function',