From d99b2b618522d994d01b85f343e97f37fc04142c Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 15 Apr 2026 12:48:44 +0200 Subject: [PATCH 1/2] Misc cleanups and refactors - Bathymetry: minor fix in regrid_bathymetry - DataWrangling: inpainting and restoring improvements - Diagnostics: MLD minor update - EarthSystemModels: simplify time stepping, minor earth_system_model cleanup - InterfaceComputations: minor refactors in sea_ice_ocean fluxes and heat flux formulations Note: Ocean simulation refactoring (flux_and_restoring, Oceans.jl, ocean_simulation.jl) and Diagnostics (interface_fluxes) are handled in PR #155 and #158. --- src/Bathymetry/regrid_bathymetry.jl | 3 ++- src/DataWrangling/inpainting.jl | 13 +++++++++---- src/DataWrangling/restoring.jl | 2 ++ src/Diagnostics/mixed_layer_depth.jl | 1 - .../InterfaceComputations/sea_ice_ocean_fluxes.jl | 12 ++++++------ .../sea_ice_ocean_heat_flux_formulations.jl | 6 +++--- src/EarthSystemModels/earth_system_model.jl | 3 +-- .../time_step_earth_system_model.jl | 9 ++------- 8 files changed, 25 insertions(+), 24 deletions(-) diff --git a/src/Bathymetry/regrid_bathymetry.jl b/src/Bathymetry/regrid_bathymetry.jl index 430cad426..564d531f3 100644 --- a/src/Bathymetry/regrid_bathymetry.jl +++ b/src/Bathymetry/regrid_bathymetry.jl @@ -341,7 +341,8 @@ end # Fix active cells to be at least `-minimum_depth`. active = z < 0 # it's a wet cell - z = ifelse(active, min(z, -minimum_depth), z) + above_minimum_depth = z > -minimum_depth + z = ifelse(active, ifelse(above_minimum_depth, zero(z), z), z) @inbounds target_z[i, j, 1] = z end diff --git a/src/DataWrangling/inpainting.jl b/src/DataWrangling/inpainting.jl index dde0f1f91..a081fa1d7 100644 --- a/src/DataWrangling/inpainting.jl +++ b/src/DataWrangling/inpainting.jl @@ -51,7 +51,12 @@ function propagate_horizontally!(inpainting::NearestNeighborInpainting, field, m iter += 1 end - launch!(arch, grid, size(field), _fill_nans!, field) + # Fill any remaining NaN values with the mean of valid data. + # Using 0 would be catastrophic for fields like salinity (~34 psu). + valid_sum = sum(x -> ifelse(isnan(x), zero(x), x), field; condition=interior(mask)) + valid_count = sum(x -> !isnan(x), field; condition=interior(mask)) + fill_value = convert(eltype(field), valid_sum / valid_count) + launch!(arch, grid, size(field), _fill_nans!, field, fill_value) fill_halo_regions!(field) return field @@ -80,7 +85,7 @@ end end FT_NaN = convert(FT, NaN) - @inbounds substituting_field[i, j, k] = ifelse(value == 0, FT_NaN, value / donors) + @inbounds substituting_field[i, j, k] = ifelse(donors == 0, FT_NaN, value / donors) end @kernel function _substitute_values!(field, substituting_field) @@ -97,9 +102,9 @@ end @inbounds field[i, j, k] = ifelse(mask[i, j, k], FT_NaN, field[i, j, k]) end -@kernel function _fill_nans!(field) +@kernel function _fill_nans!(field, fill_value) i, j, k = @index(Global, NTuple) - @inbounds field[i, j, k] *= !isnan(field[i, j, k]) + @inbounds field[i, j, k] = ifelse(isnan(field[i, j, k]), fill_value, field[i, j, k]) end """ diff --git a/src/DataWrangling/restoring.jl b/src/DataWrangling/restoring.jl index 918edfae0..d26b7c1f4 100644 --- a/src/DataWrangling/restoring.jl +++ b/src/DataWrangling/restoring.jl @@ -12,6 +12,7 @@ using Dates: Second import NumericalEarth: stateindex import Oceananigans.Forcings: materialize_forcing +import Oceananigans.OutputReaders: extract_field_time_series # Variable names for restorable data struct Temperature end @@ -234,6 +235,7 @@ function Base.show(io::IO, dsr::DatasetRestoring) end materialize_forcing(forcing::DatasetRestoring, field, field_name, model_field_names) = forcing +extract_field_time_series(forcing::DatasetRestoring) = forcing.field_time_series ##### ##### Masks for restoring diff --git a/src/Diagnostics/mixed_layer_depth.jl b/src/Diagnostics/mixed_layer_depth.jl index 0b986312b..3faeb750a 100644 --- a/src/Diagnostics/mixed_layer_depth.jl +++ b/src/Diagnostics/mixed_layer_depth.jl @@ -32,7 +32,6 @@ end function compute!(mld::MixedLayerDepthField, time=nothing) compute_mixed_layer_depth!(mld) - #@apply_regionally compute_mixed_layer_depth!(mld) fill_halo_regions!(mld) return mld end diff --git a/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl index a715493d3..5cbffb333 100644 --- a/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl @@ -175,12 +175,12 @@ end qᶠ = δ𝒬ᶠʳᶻ / ℰ @inbounds begin - Tᴺ = Tᵒᶜ[i, j, Nz] - Sᴺ = Sᵒᶜ[i, j, Nz] + Tᴺ = Tᵒᶜ[i, j, Nz] + Sᴺ = Sᵒᶜ[i, j, Nz] Sˢⁱ = ice_salinity[i, j, 1] hˢⁱ = ice_thickness[i, j, 1] - ℵᵢ = ice_concentration[i, j, 1] - hc = ice_consolidation_thickness[i, j, 1] + ℵᵢ = ice_concentration[i, j, 1] + hc = ice_consolidation_thickness[i, j, 1] end # Extract internal temperature (for ConductiveFluxTEF, zero otherwise) @@ -198,8 +198,8 @@ end # ============================================= # Returns interfacial heat flux, melt rate qᵐ, and interface T, S 𝒬ⁱᵒ, qᵐ, Tᵦ, Sᵦ = compute_interface_heat_flux(flux_formulation, - ocean_surface_state, ice_state, - liquidus, ocean_properties, ℰ, u★) + ocean_surface_state, ice_state, + liquidus, ocean_properties, ℰ, u★) # Store interface values and heat flux @inbounds 𝒬ⁱⁿᵗ[i, j, 1] = 𝒬ⁱᵒ diff --git a/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_heat_flux_formulations.jl b/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_heat_flux_formulations.jl index ebf4f4ce8..915fa64b8 100644 --- a/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_heat_flux_formulations.jl +++ b/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_heat_flux_formulations.jl @@ -113,8 +113,8 @@ References - [holland1999modeling](@citet): Holland, D. M., & Jenkins, A. (1999). Modeling thermodynamic ice–ocean interactions at the base of an ice shelf. *Journal of Physical Oceanography*, 29(8), 1787-1800. -- [hieronymus2021comparison](@citet): Hieronymus, M., et al. (2021). A comparison of ocean-ice flux parametrizations. - *Geosci. Model Dev.*, 14, 4891-4908. +- [shi2021sensitivity](@citet): Shi, X., Notz, D., Liu, J., Yang, H., & Lohmann, G. (2021). Sensitivity of Northern + Hemisphere climate to ice-ocean interface heat flux parameterizations. *Geosci. Model Dev.*, 14, 4891-4908. """ struct ThreeEquationHeatFlux{F, T, FT, U} conductive_flux :: F @@ -139,7 +139,7 @@ Adapt.adapt_structure(to, f::ThreeEquationHeatFlux) = Construct a `ThreeEquationHeatFlux` with the specified parameters. -Default values follow [hieronymus2021comparison](@citet) with ``R = \\alpha_h / \\alpha_s = 35``. +Default values follow [shi2021sensitivity](@citet) with ``R = \\alpha_h / \\alpha_s = 35``. Keyword Arguments ================= diff --git a/src/EarthSystemModels/earth_system_model.jl b/src/EarthSystemModels/earth_system_model.jl index 90469dc04..db5ff93dd 100644 --- a/src/EarthSystemModels/earth_system_model.jl +++ b/src/EarthSystemModels/earth_system_model.jl @@ -44,13 +44,12 @@ function Base.show(io::IO, cm::ESM) return nothing end -# Assumption: We have an ocean! architecture(model::ESM) = model.architecture Base.eltype(model::ESM) = Base.eltype(model.interfaces.exchanger.grid) prettytime(model::ESM) = prettytime(model.clock.time) iteration(model::ESM) = model.clock.iteration timestepper(::ESM) = nothing -default_included_properties(::ESM) = tuple() +default_included_properties(::ESM) = [] prognostic_fields(cm::ESM) = nothing fields(::ESM) = NamedTuple() default_clock(TT) = Oceananigans.TimeSteppers.Clock{TT}(0, 0, 1) diff --git a/src/EarthSystemModels/time_step_earth_system_model.jl b/src/EarthSystemModels/time_step_earth_system_model.jl index 0ec04e6e2..c486171b3 100644 --- a/src/EarthSystemModels/time_step_earth_system_model.jl +++ b/src/EarthSystemModels/time_step_earth_system_model.jl @@ -15,13 +15,8 @@ function time_step!(coupled_model::EarthSystemModel, Δt; callbacks=[]) atmosphere = coupled_model.atmosphere # Eventually, split out into OceanOnlyModel - !isnothing(sea_ice) && time_step!(sea_ice, Δt) - - # TODO after ice time-step: - # - Adjust ocean heat flux if the ice completely melts? - !isnothing(ocean) && time_step!(ocean, Δt) - - # Time step the atmosphere + !isnothing(sea_ice) && time_step!(sea_ice, Δt) + !isnothing(ocean) && time_step!(ocean, Δt) !isnothing(atmosphere) && time_step!(atmosphere, Δt) # TODO: From b2052f32adf257db2fa0d87566953f075a312365 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 15 Apr 2026 14:56:30 +0200 Subject: [PATCH 2/2] back to previous filling --- src/DataWrangling/inpainting.jl | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/src/DataWrangling/inpainting.jl b/src/DataWrangling/inpainting.jl index a081fa1d7..dde0f1f91 100644 --- a/src/DataWrangling/inpainting.jl +++ b/src/DataWrangling/inpainting.jl @@ -51,12 +51,7 @@ function propagate_horizontally!(inpainting::NearestNeighborInpainting, field, m iter += 1 end - # Fill any remaining NaN values with the mean of valid data. - # Using 0 would be catastrophic for fields like salinity (~34 psu). - valid_sum = sum(x -> ifelse(isnan(x), zero(x), x), field; condition=interior(mask)) - valid_count = sum(x -> !isnan(x), field; condition=interior(mask)) - fill_value = convert(eltype(field), valid_sum / valid_count) - launch!(arch, grid, size(field), _fill_nans!, field, fill_value) + launch!(arch, grid, size(field), _fill_nans!, field) fill_halo_regions!(field) return field @@ -85,7 +80,7 @@ end end FT_NaN = convert(FT, NaN) - @inbounds substituting_field[i, j, k] = ifelse(donors == 0, FT_NaN, value / donors) + @inbounds substituting_field[i, j, k] = ifelse(value == 0, FT_NaN, value / donors) end @kernel function _substitute_values!(field, substituting_field) @@ -102,9 +97,9 @@ end @inbounds field[i, j, k] = ifelse(mask[i, j, k], FT_NaN, field[i, j, k]) end -@kernel function _fill_nans!(field, fill_value) +@kernel function _fill_nans!(field) i, j, k = @index(Global, NTuple) - @inbounds field[i, j, k] = ifelse(isnan(field[i, j, k]), fill_value, field[i, j, k]) + @inbounds field[i, j, k] *= !isnan(field[i, j, k]) end """