diff --git a/PHYSICS_PORT_VALIDATION.md b/PHYSICS_PORT_VALIDATION.md
new file mode 100644
index 000000000..6a5a924e8
--- /dev/null
+++ b/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-2–10-1 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` − `old_branch`)
+
+| Field | max|diff| | relative | | Field | max|diff| | relative |
+|---|---|---|---|---|---|---|
+| ps | 916.9 Pa | 0.9% | | flux_lhe | 52.9 W/m² | 7.4% |
+| temp | 5.03 K | 1.6% | | flux_t | 45.6 W/m² | 15.8% |
+| ucomp | 14.26 m/s | 18.9% | | soc_surf_flux_lw | 13.80 W/m² | 8.0% |
+| vcomp | 11.54 m/s | 24.7% | | soc_surf_flux_sw | 2.94 W/m² | 0.9% |
+| t_surf | 1.74 K | 0.6% | | soc_olr | 4.88 W/m² | 1.5% |
+| temp_2m | 1.85 K | 0.6% | | soc_toa_sw | 2.49 W/m² | 0.6% |
+| rh | 17.2 % | 15.8% | | soc_ozone | 6.3×10-8 | 0.4% |
+| precipitation | 4.8×10-5 | 19.7% | | vor | 2.1×10-5 | 33.0% |
+| **albedo** | **0.0** | **0.0%** | | div | 1.5×10-5 | 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-7-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-7 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 ≥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.
diff --git a/exp/test_cases/realistic_continents/namelist_basefile.nml b/exp/test_cases/realistic_continents/namelist_basefile.nml
index b4f5de343..5d10ca74a 100644
--- a/exp/test_cases/realistic_continents/namelist_basefile.nml
+++ b/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
diff --git a/exp/test_cases/simple_clouds/socrates_simcloud.py b/exp/test_cases/simple_clouds/socrates_simcloud.py
index 7fdaba0fa..9a1162b7a 100644
--- a/exp/test_cases/simple_clouds/socrates_simcloud.py
+++ b/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',
diff --git a/src/atmos_spectral/driver/solo/mixed_layer.F90 b/src/atmos_spectral/driver/solo/mixed_layer.F90
index ce7e1541d..6ecfb6150 100644
--- a/src/atmos_spectral/driver/solo/mixed_layer.F90
+++ b/src/atmos_spectral/driver/solo/mixed_layer.F90
@@ -133,7 +133,17 @@ module mixed_layer_mod
real :: ice_albedo_value = 0.7
real :: ice_concentration_threshold = 0.5
logical :: update_albedo_from_ice = .false.
+logical :: update_land_mask_from_ice = .false.
character(len=256) :: ice_albedo_method = 'step_function'
+logical :: allow_qflux_over_land = .false.
+
+logical :: specify_sst_over_sea_ice = .false. !Problem with specifying SSTs when Isca has no ice model is poles are too warm. This adds the option to have a separate SST specification over SEA ice only.
+logical :: linearly_interpolate_sea_ice_temp_and_sst = .true. !Mix SST input and ice sst linearly based on ice concentration
+character(len=256) :: ice_sst_file = 'temp_2m_input'
+
+logical :: specify_sst_over_land_from_separate_file_to_ocean_sst = .false. !Isca's lack of land model introduces temperature biases. Try forcibly correcting them by specifying ssts there.
+character(len=256) :: land_sst_file = 'temp_2m_input'
+
logical :: add_latent_heat_flux_anom = .false.
character(len=256) :: flux_lhe_anom_file_name = 'INPUT/flux_lhe_anom.nc'
@@ -153,11 +163,17 @@ module mixed_layer_mod
land_albedo_prefactor, & !s
load_qflux,qflux_file_name,time_varying_qflux, &
update_albedo_from_ice, ice_file_name, &
+ update_land_mask_from_ice, &
ice_albedo_value, specify_sst_over_ocean_only, &
ice_concentration_threshold, ice_albedo_method,&
add_latent_heat_flux_anom,flux_lhe_anom_file_name,&
flux_lhe_anom_field_name, specify_constant_sst,&
sst_prescribed_constant, &
+ 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, &
do_ape_sst, qflux_field_name
!=================================================================================================================================
@@ -206,6 +222,9 @@ module mixed_layer_mod
zsurf, & ! mj know about topography
land_sea_heat_capacity,&
sst_new, & ! mj input SST
+ sst_input, & ! mj input SST
+ ice_sst_new, &
+ land_sst_new, &
albedo_initial
logical, allocatable, dimension(:,:) :: land_mask
@@ -214,6 +233,8 @@ module mixed_layer_mod
type(interpolate_type),save :: sst_interp
type(interpolate_type),save :: qflux_interp
type(interpolate_type),save :: ice_interp
+ type(interpolate_type),save :: ice_sst_interp
+ type(interpolate_type),save :: land_sst_interp
type(interpolate_type),save :: flux_lhe_anom_interp
real inv_cp_air
@@ -294,7 +315,10 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax
allocate(land_sea_heat_capacity (is:ie, js:je))
allocate(zsurf (is:ie, js:je))
allocate(sst_new (is:ie, js:je))
-allocate(land_mask (is:ie, js:je)); land_mask=land
+allocate(sst_input (is:ie, js:je))
+allocate(ice_sst_new (is:ie, js:je))
+allocate(land_sst_new (is:ie, js:je))
+allocate(land_mask (is:ie, js:je)); land_mask=land
!
!see if restart file exists for the surface temperature
!
@@ -323,8 +347,17 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax
call interpolator_init( sst_interp, trim(sst_file)//'.nc', rad_lonb_2d, rad_latb_2d, data_out_of_bounds=(/CONSTANT/) )
endif
+ if (specify_sst_over_sea_ice) then
+ call interpolator_init( ice_sst_interp, trim(ice_sst_file)//'.nc', rad_lonb_2d, rad_latb_2d, data_out_of_bounds=(/CONSTANT/) )
+ endif
+ if (specify_sst_over_land_from_separate_file_to_ocean_sst) then
+ call interpolator_init( land_sst_interp, trim(land_sst_file)//'.nc', rad_lonb_2d, rad_latb_2d, data_out_of_bounds=(/CONSTANT/) )
+ endif
+ if (update_albedo_from_ice) then
+ call interpolator_init( ice_interp, trim(ice_file_name)//'.nc', rad_lonb_2d, rad_latb_2d, data_out_of_bounds=(/CONSTANT/) )
+ endif
if (file_exist('INPUT/mixed_layer.res.nc')) then
@@ -346,6 +379,23 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax
call interpolator( sst_interp, Time, t_surf, trim(sst_file) )
+ if(specify_sst_over_sea_ice) then
+ call interpolator( ice_sst_interp, Time, ice_sst_new, trim(ice_sst_file) )
+
+ if(linearly_interpolate_sea_ice_temp_and_sst) then
+ call read_ice_conc(Time)
+ where(.not. land_mask) t_surf = t_surf + (ice_concentration) *(ice_sst_new - t_surf) !Only do this where sea ice.
+ endif
+ endif
+
+ if(specify_sst_over_land_from_separate_file_to_ocean_sst) then
+ call interpolator( land_sst_interp, Time, land_sst_new, trim(land_sst_file) )
+ where (land_mask) t_surf = land_sst_new
+ endif
+
+
+
+
elseif (prescribe_initial_dist) then
! call error_mesg('mixed_layer','mixed_layer restart file not found - initializing from prescribed distribution', WARNING)
@@ -488,7 +538,6 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax
albedo_initial=albedo
if (update_albedo_from_ice) then
- call interpolator_init( ice_interp, trim(ice_file_name)//'.nc', rad_lonb_2d, rad_latb_2d, data_out_of_bounds=(/CONSTANT/) )
call read_ice_conc(Time)
call albedo_calc(albedo,Time)
else
@@ -611,14 +660,14 @@ subroutine mixed_layer ( &
call error_mesg('mixed_layer','mixed_layer module is not initialized',FATAL)
endif
-if(update_albedo_from_ice) then
- call read_ice_conc(Time_next)
- land_ice_mask=.false.
- where(land_mask.or.(ice_concentration.gt.ice_concentration_threshold))
- land_ice_mask=.true.
- end where
+if(update_land_mask_from_ice) then
+ call read_ice_conc(Time_next)
+ land_ice_mask=.false.
+ where(land_mask.or.(ice_concentration.gt.ice_concentration_threshold))
+ land_ice_mask=.true.
+ end where
else
- land_ice_mask=land_mask
+ land_ice_mask=land_mask
endif
call albedo_calc(albedo_out,Time_next)
@@ -664,9 +713,11 @@ subroutine mixed_layer ( &
if(load_qflux.and.time_varying_qflux) then
call interpolator( qflux_interp, Time, ocean_qflux, trim(qflux_file_name) )
- if(update_albedo_from_ice) then
- where (land_ice_mask) ocean_qflux=0.
- endif
+ if(.not.allow_qflux_over_land) then
+ if(update_albedo_from_ice) then
+ where (land_ice_mask) ocean_qflux=0.
+ endif
+ endif
endif
@@ -686,18 +737,36 @@ subroutine mixed_layer ( &
if(do_sc_sst) then !mj sst read from input file
! read at the new time, as that is what we are stepping to
if (specify_constant_sst) then
- sst_new = sst_prescribed_constant
+ sst_input = sst_prescribed_constant
else
- call interpolator( sst_interp, Time_next, sst_new, trim(sst_file) )
+ call interpolator( sst_interp, Time_next, sst_input, trim(sst_file) )
endif
- if(specify_sst_over_ocean_only) then
+ sst_new = sst_input
+
+ if(specify_sst_over_sea_ice) then
+ call interpolator( ice_sst_interp, Time_next, ice_sst_new, trim(ice_sst_file) )
+
+ if(linearly_interpolate_sea_ice_temp_and_sst) then
+ call read_ice_conc(Time_next)
+ where (.not.land_mask) sst_new = sst_input + (ice_concentration) *(ice_sst_new - sst_input) !Only do this where sea ice.
+ where (land_mask) sst_new = sst_input
+ endif
+ endif
+
+ if(specify_sst_over_land_from_separate_file_to_ocean_sst) then
+ call interpolator( land_sst_interp, Time_next, land_sst_new, trim(land_sst_file) )
+ where (land_mask) sst_new = land_sst_new
+ endif
+
+
+ if(specify_sst_over_ocean_only) then
where (.not.land_ice_mask) delta_t_surf = sst_new - t_surf
- where (.not.land_ice_mask) t_surf = t_surf + delta_t_surf
- else
+ where (.not.land_ice_mask) t_surf = t_surf + delta_t_surf
+ else
delta_t_surf = sst_new - t_surf
t_surf = t_surf + delta_t_surf
- endif
+ endif
end if
diff --git a/src/extra/python/scripts/calculate_qflux/calculate_qflux.py b/src/extra/python/scripts/calculate_qflux/calculate_qflux.py
index ae2df9bbe..6214c02da 100644
--- a/src/extra/python/scripts/calculate_qflux/calculate_qflux.py
+++ b/src/extra/python/scripts/calculate_qflux/calculate_qflux.py
@@ -15,10 +15,10 @@
__author__='Stephen Thomson'
-def qflux_calc(dataset, model_params, output_file_name, ice_file_name=None, groupby_name='months'):
+def qflux_calc(dataset, model_params, output_file_name, ice_file_name=None, groupby_name='months', ignore_ice_for_calculation=False):
if groupby_name=='months':
- time_varying_ice = ice_mask_calculation(dataset, dataset.land, ice_file_name)
+ time_varying_ice = ice_mask_calculation(dataset, dataset.land, ice_file_name, ignore_ice_for_calculation=ignore_ice_for_calculation)
upper_ocean_heat_content(dataset, model_params, time_varying_ice)
net_surf_energy_flux(dataset, model_params)
deep_ocean_heat_content(dataset, model_params)
@@ -27,7 +27,7 @@ def qflux_calc(dataset, model_params, output_file_name, ice_file_name=None, grou
output_dict={'manual_grid_option':False, 'is_thd':False, 'num_years':1., 'time_spacing_days':12, 'file_name':output_file_name+'.nc', 'var_name':output_file_name} #Have specified that var name is the same as file name as this is what the fortran assumes.
elif groupby_name=='dayofyear':
- time_varying_ice = ice_mask_calculation(dataset, dataset.land, ice_file_name)
+ time_varying_ice = ice_mask_calculation(dataset, dataset.land, ice_file_name, ignore_ice_for_calculation=ignore_ice_for_calculation)
upper_ocean_heat_content(dataset, model_params, time_varying_ice, dayofyear_or_months='dayofyear')
net_surf_energy_flux(dataset, model_params)
deep_ocean_heat_content(dataset, model_params)
@@ -36,7 +36,7 @@ def qflux_calc(dataset, model_params, output_file_name, ice_file_name=None, grou
output_dict={'manual_grid_option':False, 'is_thd':False, 'num_years':1., 'time_spacing_days':12, 'file_name':output_file_name+'.nc', 'var_name':output_file_name}
elif groupby_name=='all_time':
- time_varying_ice = ice_mask_calculation(dataset, dataset.land, ice_file_name, dayofyear_or_months=groupby_name)
+ time_varying_ice = ice_mask_calculation(dataset, dataset.land, ice_file_name, dayofyear_or_months=groupby_name, ignore_ice_for_calculation=ignore_ice_for_calculation)
upper_ocean_heat_content(dataset, model_params, time_varying_ice, dayofyear_or_months=groupby_name)
net_surf_energy_flux(dataset, model_params, dayofyear_or_months=groupby_name)
deep_ocean_heat_content(dataset, model_params, dayofyear_or_months=groupby_name)
@@ -53,30 +53,35 @@ def time_gradient(data_in, delta_t):
return data_out
-def ice_mask_calculation(dataset, land_array, ice_file_name, dayofyear_or_months='months'):
-
- try:
- ice_climatology=dataset['ice_conc'].groupby(dayofyear_or_months).mean('time').load()
- ice_array=ice_climatology.values
- ice_idx=ice_array !=0.
- ice_array[ice_idx]=1.0
- time_varying_ice=True
- print('have gotten ice concentration from climatology')
- except KeyError:
+def ice_mask_calculation(dataset, land_array, ice_file_name, dayofyear_or_months='months', ignore_ice_for_calculation=False):
+ if ignore_ice_for_calculation:
+ ice_array=np.zeros_like(land_array)
+ time_varying_ice = False
+ print('Ignoring any ice climatology')
+ else:
try:
- ice_data_temp = xarray.open_dataset(ice_file_name, decode_times=False)
- albedo_constant_value = ice_data_temp['albedo']
- albedo_array = albedo_constant_value.values.squeeze()
- ice_idx = np.round(albedo_array, decimals=2) == 0.7
- ice_array=np.zeros_like(albedo_array)
- ice_array[ice_idx] = 1.0
- time_varying_ice = False
- print('have gotten ice concentration from one month albedo')
- except TypeError:
- ice_array=np.zeros_like(land_array)
- time_varying_ice = False
- print('no ice climatology')
+ ice_climatology=dataset['ice_conc'].groupby(dayofyear_or_months).mean('time').load()
+ ice_array=ice_climatology.values
+ ice_idx=ice_array !=0.
+ ice_array[ice_idx]=1.0
+ time_varying_ice=True
+ print('have gotten ice concentration from climatology')
+ except KeyError:
+
+ try:
+ ice_data_temp = xarray.open_dataset(ice_file_name, decode_times=False)
+ albedo_constant_value = ice_data_temp['albedo']
+ albedo_array = albedo_constant_value.values.squeeze()
+ ice_idx = np.round(albedo_array, decimals=2) == 0.7
+ ice_array=np.zeros_like(albedo_array)
+ ice_array[ice_idx] = 1.0
+ time_varying_ice = False
+ print('have gotten ice concentration from one month albedo')
+ except TypeError:
+ ice_array=np.zeros_like(land_array)
+ time_varying_ice = False
+ print('no ice climatology')
land_ice_mask=np.zeros_like(ice_array)
@@ -191,7 +196,12 @@ def net_surf_energy_flux(dataset, model_params, dayofyear_or_months='months'):
# scaling_factor=1.0
- print('using scale factor for SW of '+ str(scaling_factor))
+ try:
+ scaling_factor_string = str(scaling_factor.values)
+ except:
+ scaling_factor_string = str(scaling_factor)
+
+ print('using scale factor for SW of '+ scaling_factor_string)
net_surf_energy_fl=(scaling_factor*dataset['flux_sw_clim']+dataset['flux_lw_clim']-(model_params['sigma_sb']*dataset['sst_clim']**4.0)-dataset['flux_t_clim']-dataset['flux_lhe_clim'])*(1.0-dataset['land_ice_mask'])
@@ -233,38 +243,53 @@ def regrid_in_time(dataset, groupby_name):
dataset['masked_ocean_transport'] = (('months_ax','lat','lon'), dataset_to_output)
-def check_surface_flux_dims(dataset):
+def check_surface_flux_dims(dataset, rad_scheme_type):
''' This surface flux checker is designed to decide if we're using grey rad or not. If we're using grey rad then the definition
of flux_sw and flux_lw are different to RRTM. The script was written to use RRTM output, so it changes variable names etc to be
equivalent to RRTM definitions.
'''
- flux_dims = dataset['flux_sw'].dims
+ lower_case_rad_scheme = rad_scheme_type.lower()
- if 'phalf' in flux_dims:
- dataset.rename({'flux_sw':'flux_sw'+'_3d'}, inplace=True)
- max_pressure = dataset.phalf.max()
- flux_at_bottom_phalf_level = dataset['flux_sw_3d'].sel(phalf=max_pressure)
- new_dims = ('time','lat','lon')
- dataset['flux_sw'] = (new_dims, flux_at_bottom_phalf_level)
+ if lower_case_rad_scheme not in ['socrates', 'rrtm', 'grey']:
+ raise NotImplementedError(f"Radiation scheme should be one of 'socrates', 'rrtm' or 'grey'. The value entered was {lower_case_rad_scheme}")
- flux_dims_lw = dataset['flux_lw'].dims
+ if lower_case_rad_scheme=='rrtm':
+ print('No renaming of variables required when using RRTM')
+ elif lower_case_rad_scheme=='grey':
+ print('Renaming variables to enable use of grey radiation')
- if 'phalf' in flux_dims_lw:
- dataset.rename({'flux_lw':'flux_lw'+'_3d'}, inplace=True)
- try:
- # Script assumes flux_lw is the surface lw down (i.e. not a net flux). This is the case with RRTM, but with
- # grey radiation 'flux_lw' is the net lw flux in 3D. So we take the lwdn_sfc output from grey rad and rename it
- # flux_lw.
- dataset['lwdn_sfc']
- dataset.rename({'lwdn_sfc':'flux_lw'}, inplace=True)
- except:
- #If lwdn_sfc is not available, then we re-calculate it from flux_lw by adding back sigma*t_surf**4, then call it flux_lw
- print('lwdn_sfc not present when using grey radiation, so re-calculating it from flux_lw.')
+ flux_dims = dataset['flux_sw'].dims
+
+ if 'phalf' in flux_dims:
+ dataset.rename({'flux_sw':'flux_sw'+'_3d'}, inplace=True)
max_pressure = dataset.phalf.max()
- lwdn_sfc = dataset.flux_lw_3d.sel(phalf=max_pressure) + sigma_sb*dataset.t_surf**4.
+ flux_at_bottom_phalf_level = dataset['flux_sw_3d'].sel(phalf=max_pressure)
new_dims = ('time','lat','lon')
- dataset['flux_lw'] = (new_dims, lwdn_sfc)
+ dataset['flux_sw'] = (new_dims, flux_at_bottom_phalf_level)
+
+ flux_dims_lw = dataset['flux_lw'].dims
+
+ if 'phalf' in flux_dims_lw:
+ dataset.rename({'flux_lw':'flux_lw'+'_3d'}, inplace=True)
+ try:
+ # Script assumes flux_lw is the surface lw down (i.e. not a net flux). This is the case with RRTM, but with
+ # grey radiation 'flux_lw' is the net lw flux in 3D. So we take the lwdn_sfc output from grey rad and rename it
+ # flux_lw.
+ dataset['lwdn_sfc']
+ dataset.rename({'lwdn_sfc':'flux_lw'}, inplace=True)
+ except:
+ #If lwdn_sfc is not available, then we re-calculate it from flux_lw by adding back sigma*t_surf**4, then call it flux_lw
+ print('lwdn_sfc not present when using grey radiation, so re-calculating it from flux_lw.')
+ max_pressure = dataset.phalf.max()
+ lwdn_sfc = dataset.flux_lw_3d.sel(phalf=max_pressure) + sigma_sb*dataset.t_surf**4.
+ new_dims = ('time','lat','lon')
+ dataset['flux_lw'] = (new_dims, lwdn_sfc)
+
+ elif lower_case_rad_scheme=='socrates':
+ print('Renaming variables to enable use of socrates radiation')
+ dataset.rename({'soc_surf_flux_sw':'flux_sw'}, inplace=True)
+ dataset.rename({'soc_surf_flux_lw_down':'flux_lw'}, inplace=True)
if __name__ == "__main__":
@@ -281,19 +306,22 @@ def check_surface_flux_dims(dataset):
input_dir=GFDL_BASE
- base_dir=GFDL_DATA
- land_file='input/land.nc'
- base_exp_name='annual_mean_ice_post_princeton_fixed_sst/' #Folder containing the python script and input files that ran the experiment
- exp_name='annual_mean_ice_post_princeton_fixed_sst_1' #Folder within the data directory where the files can be found
+ base_dir='/disca/share/sit204/data_from_isca_cpu/cssp_perturb_exps/control/'
+ land_file='input/era-spectral7_T42_64x128.out.nc'
+ base_exp_name='socrates_fast_land/' #Folder containing the python script and input files that ran the experiment
+ exp_name='soc_ga3_do_simple_false_cmip_o3_bucket/grouped_nc_files/' #Folder within the data directory where the files can be found
# ice_file_name=base_dir+'annual_mean_ice_albedo_change_test_mk2_4320_dt_rad_4/'+'run360/'+'atmos_monthly.nc'
- ice_file_name = '/scratch/sit204/data_isca/realistic_continents_fixed_sst_test_experiment_albedo/run0001/atmos_daily.nc'
- output_file_name='ami_test_interp' #Proposed name of your output qflux file. Will also be qflux field name in q-flux netcdf file as the fortran assumes file-name = field name. No need to add '.nc' or any file paths in this variable as otherwise they will end up in the field name too. Output file will be stored in the same directory as this script.
+ ice_file_name = None
+ ignore_ice=True
+ output_file_name='soc_ga3_bucket_jra_55_ice_temps' #Proposed name of your output qflux file. Will also be qflux field name in q-flux netcdf file as the fortran assumes file-name = field name. No need to add '.nc' or any file paths in this variable as otherwise they will end up in the field name too. Output file will be stored in the same directory as this script.
- start_file=240
- end_file=360
+ start_file=121
+ end_file=3000
land_present=True
use_interpolated_pressure_level_data = False #Conditions the script on whether to expect data on sigma levels (if False) or pressure levels (if True). Script should be insensitive to this choice if both sets of files exist.
+ use_clim = True
+
#Set time increments of input files (e.g. `monthly` for `atmos_monthly` files.
avg_or_daily='monthly'
@@ -302,13 +330,21 @@ def check_surface_flux_dims(dataset):
model_params = sagp.model_params_set(input_dir, delta_t=720., ml_depth=20., res=42)
- dataset, time_arr, size_list = io.read_data( base_dir,exp_name,start_file,end_file,avg_or_daily,use_interpolated_pressure_level_data)
+ if use_clim:
+ clim_str = '_clim'
+ else:
+ clim_str = ''
+ file_name = 'atmos_'+avg_or_daily+'_together_interp_new_height_temp_not_below_ground'+str(start_file)+'_'+str(end_file)+clim_str
+
+ dataset, time_arr, size_list = io.read_data( base_dir,exp_name,start_file,end_file,avg_or_daily,use_interpolated_pressure_level_data, file_name=file_name)
land_array, topo_array = io.read_land(input_dir,base_exp_name,land_present,use_interpolated_pressure_level_data,size_list,land_file)
dataset['land'] = (('lat','lon'),land_array)
- check_surface_flux_dims(dataset)
+ which_radiation_scheme_are_you_using = 'socrates'
+
+ check_surface_flux_dims(dataset, which_radiation_scheme_are_you_using)
- qflux_calc(dataset, model_params, output_file_name, ice_file_name, groupby_name=time_divisions_of_qflux_to_be_calculated)
+ qflux_calc(dataset, model_params, output_file_name, ice_file_name, groupby_name=time_divisions_of_qflux_to_be_calculated, ignore_ice_for_calculation=ignore_ice)
diff --git a/src/extra/python/scripts/calculate_qflux/nc_file_io_xarray.py b/src/extra/python/scripts/calculate_qflux/nc_file_io_xarray.py
index d7416d747..d0daa4a88 100644
--- a/src/extra/python/scripts/calculate_qflux/nc_file_io_xarray.py
+++ b/src/extra/python/scripts/calculate_qflux/nc_file_io_xarray.py
@@ -13,35 +13,40 @@
def read_data( base_dir, exp_name, start_file, end_file, avg_or_daily, use_interpolated_pressure_level_data, model='fms13', file_name=None):
if model=='fms13':
+ if file_name is None:
- possible_format_strs = [[base_dir+'/'+exp_name+'/run%03d' % m for m in range(start_file, end_file+1)],
- [base_dir+'/'+exp_name+'/run%04d' % m for m in range(start_file, end_file+1)],
- [base_dir+'/'+exp_name+'/run%d' % m for m in range(start_file, end_file+1)]]
+ possible_format_strs = [[base_dir+'/'+exp_name+'/run%03d' % m for m in range(start_file, end_file+1)],
+ [base_dir+'/'+exp_name+'/run%04d' % m for m in range(start_file, end_file+1)],
+ [base_dir+'/'+exp_name+'/run%d' % m for m in range(start_file, end_file+1)]]
- if(use_interpolated_pressure_level_data):
- if avg_or_daily == 'monthly':
-# extra='_interp.nc'
- extra='_interp_new_height.nc'
+ if(use_interpolated_pressure_level_data):
+ if avg_or_daily == 'monthly':
+ # extra='_interp.nc'
+ extra='_interp_new_height.nc'
+ else:
+ extra='_interp_new_height_temp.nc'
else:
- extra='_interp_new_height_temp.nc'
- else:
- extra='.nc'
+ extra='.nc'
- thd_string = '/atmos_'+avg_or_daily+extra
+ thd_string = '/atmos_'+avg_or_daily+extra
- for format_str_files in possible_format_strs:
- files_temp = format_str_files
+ for format_str_files in possible_format_strs:
+ files_temp = format_str_files
- thd_files = [s + thd_string for s in files_temp]
+ thd_files = [s + thd_string for s in files_temp]
- thd_files_exist=[os.path.isfile(s) for s in thd_files]
+ thd_files_exist=[os.path.isfile(s) for s in thd_files]
+
+ if thd_files_exist[0]:
+ break
- if thd_files_exist[0]:
- break
+ if not thd_files_exist[0] and possible_format_strs.index(format_str_files)==(len(possible_format_strs)-1):
+ raise EOFError('EXITING BECAUSE NO APPROPRIATE FORMAT STR', [thd_files[elem] for elem in [0] if not thd_files_exist[elem]])
- if not thd_files_exist[0] and possible_format_strs.index(format_str_files)==(len(possible_format_strs)-1):
- raise EOFError('EXITING BECAUSE NO APPROPRIATE FORMAT STR', [thd_files[elem] for elem in [0] if not thd_files_exist[elem]])
+ else:
+ thd_files=[base_dir+'/'+exp_name+'/'+file_name+'.nc']
+ thd_files_exist=[os.path.isfile(s) for s in thd_files]
print(thd_files[0])
diff --git a/src/extra/python/scripts/cell_area.py b/src/extra/python/scripts/cell_area.py
index e8f596a72..6a9fb0c04 100644
--- a/src/extra/python/scripts/cell_area.py
+++ b/src/extra/python/scripts/cell_area.py
@@ -1,5 +1,7 @@
import numpy as np
from netCDF4 import Dataset
+import gauss_grid as gg
+import pdb
def cell_area_all(t_res,base_dir, radius=6376.0e3):
"""read in grid from approriate file, and return 2D array of grid cell areas in metres**2."""
@@ -11,7 +13,7 @@ def cell_area_all(t_res,base_dir, radius=6376.0e3):
lonb = resolution_file.variables['lonb'][:]
latb = resolution_file.variables['latb'][:]
- area_array,xsize_array,ysize_array = cell_area_calculate(lons, lats, lonb, latb, radius)
+ area_array,xsize_array,ysize_array = cell_area_calculate(lons, lats, lonb, latb)
return area_array,xsize_array,ysize_array
@@ -49,23 +51,31 @@ def cell_area_from_xar(dataset, lat_name='lat', lon_name = 'lon', latb_name='lat
latb = dataset[latb_name].values
lonb = dataset[lonb_name].values
except KeyError:
- delta_lat=(lats[1]-lats[0])
- if np.all((lats[1:10]-lats[0:9]) == delta_lat):
+ delta_lat=np.round(lats[1]-lats[0],4)
+ if np.all(np.round(lats[1:10]-lats[0:9],4) == delta_lat):
latb = np.zeros((len(lats)+1))
for latb_idx in range(len(lats)):
- latb[latb_idx] = lats[latb_idx]-delta_lat / 2.
- latb[-1] = lats[-1] + delta_lat / 2.
-
- delta_lon=(lons[1]-lons[0])
- if np.all((lons[1:10]-lons[0:9]) == delta_lon):
+ latb[latb_idx] = np.round(lats[latb_idx],4)-delta_lat / 2.
+ latb[-1] = np.round(lats[-1],4) + delta_lat / 2.
+ else:
+ n_lat_model = lats.shape[0]
+ model_grid_lats = gg.gaussian_latitudes(int(n_lat_model/2.))[0]
+ if np.all(np.around(model_grid_lats,2)==np.around(lats,2)):
+ model_grid_latbs_bounds = gg.gaussian_latitudes(int(n_lat_model/2.))[1]
+ latb = [model_grid_latbs_bounds[i][0] for i in range(n_lat_model)]
+ latb.append(90.)
+ dataset['latb'] = (('latb'), latb)
+
+ delta_lon=np.round(lons[1]-lons[0],4)
+ if np.all(np.round(lons[1:10]-lons[0:9],4) == delta_lon):
lonb = np.zeros((len(lons)+1))
for lonb_idx in range(len(lons)):
- lonb[lonb_idx] = lons[lonb_idx]-delta_lon / 2.
- lonb[-1] = lons[-1] + delta_lon / 2.
-
-
+ lonb[lonb_idx] = np.round(lons[lonb_idx],4)-delta_lon / 2.
+ lonb[-1] = np.round(lons[-1],4) + delta_lon / 2.
+
+ dataset['lonb'] = (('lonb'), lonb)
area_array,xsize_array,ysize_array = cell_area_calculate(lons, lats, lonb, latb, radius)