From d99b2b618522d994d01b85f343e97f37fc04142c Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 15 Apr 2026 12:48:44 +0200 Subject: [PATCH 1/9] 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 97392107afcdc6b1c93bb231226382881f236301 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 15 Apr 2026 12:48:44 +0200 Subject: [PATCH 2/9] Misc cleanups and refactors - Bathymetry: minor fix in regrid_bathymetry - DataWrangling: inpainting and restoring improvements - Diagnostics: interface flux diagnostics and MLD updates - Oceans: refactor ocean simulation, extract flux_and_restoring - EarthSystemModels: simplify time stepping, minor earth_system_model cleanup - InterfaceComputations: minor refactors in sea_ice_ocean fluxes and heat flux formulations --- src/Diagnostics/Diagnostics.jl | 2 + src/Diagnostics/interface_fluxes.jl | 9 ++-- src/Oceans/Oceans.jl | 19 ++++--- src/Oceans/flux_and_restoring.jl | 46 +++++++++++++++++ src/Oceans/ocean_simulation.jl | 77 ++++++++++++++++++----------- 5 files changed, 115 insertions(+), 38 deletions(-) create mode 100644 src/Oceans/flux_and_restoring.jl diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl index 86cb7446a..4f698a4b1 100644 --- a/src/Diagnostics/Diagnostics.jl +++ b/src/Diagnostics/Diagnostics.jl @@ -14,7 +14,9 @@ using Oceananigans.BoundaryConditions: FieldBoundaryConditions, fill_halo_region using Oceananigans.Fields: FieldStatus using Oceananigans.Utils: launch! using KernelAbstractions: @index, @kernel +using Oceananigans.BoundaryConditions: DiscreteBoundaryFunction using NumericalEarth.EarthSystemModels: EarthSystemModel +using NumericalEarth.Oceans: FluxAndRestoring import Oceananigans.Fields: compute! diff --git a/src/Diagnostics/interface_fluxes.jl b/src/Diagnostics/interface_fluxes.jl index 4d144d84c..2d5b80758 100644 --- a/src/Diagnostics/interface_fluxes.jl +++ b/src/Diagnostics/interface_fluxes.jl @@ -1,4 +1,8 @@ +@inline flux_field(condition) = condition +@inline flux_field(bc::FluxAndRestoring) = bc.flux_field +@inline flux_field(bc::DiscreteBoundaryFunction) = flux_field(bc.func) + ########################### ### Temperature fluxes ########################### @@ -21,12 +25,11 @@ end Return the net temperature flux (K m s⁻¹) at the ocean's surface in a coupled `esm`. """ function net_ocean_temperature_flux(esm::EarthSystemModel) - Jᵀ = esm.ocean.model.tracers.T.boundary_conditions.top.condition + Jᵀ = flux_field(esm.ocean.model.tracers.T.boundary_conditions.top.condition) net_ocean_temperature_flux = Jᵀ + frazil_temperature_flux(esm) return Field(net_ocean_temperature_flux) end - """ sea_ice_ocean_temperature_flux(esm::EarthSystemModel) @@ -116,7 +119,7 @@ end Return the net salinity flux (g/kg m s⁻¹) at the ocean's surface in a coupled `esm`. """ function net_ocean_salinity_flux(esm::EarthSystemModel) - Jˢ = esm.ocean.model.tracers.S.boundary_conditions.top.condition + Jˢ = flux_field(esm.ocean.model.tracers.S.boundary_conditions.top.condition) return Jˢ end diff --git a/src/Oceans/Oceans.jl b/src/Oceans/Oceans.jl index 93b1d2727..87093e37f 100644 --- a/src/Oceans/Oceans.jl +++ b/src/Oceans/Oceans.jl @@ -1,13 +1,13 @@ module Oceans -export ocean_simulation, SlabOcean +export ocean_simulation, SlabOcean, FluxAndRestoring using Oceananigans using Oceananigans.Units using Oceananigans.Utils using Oceananigans.Utils: with_tracers using Oceananigans.Advection: FluxFormAdvection -using Oceananigans.BoundaryConditions: DefaultBoundaryCondition +using Oceananigans.BoundaryConditions: DefaultBoundaryCondition, DiscreteBoundaryFunction using Oceananigans.ImmersedBoundaries: immersed_peripheral_node, inactive_node, MutableGridOfSomeKind using Oceananigans.OrthogonalSphericalShellGrids using Oceananigans.Operators @@ -62,6 +62,7 @@ default_or_override(override, alternative_default=nothing) = override include("slab_ocean.jl") include("barotropic_potential_forcing.jl") include("radiative_forcing.jl") +include("flux_and_restoring.jl") include("ocean_simulation.jl") include("assemble_net_ocean_fluxes.jl") @@ -74,12 +75,12 @@ ocean_temperature(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = ocean.mode function ocean_surface_salinity(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) kᴺ = size(ocean.model.grid, 3) - return interior(ocean.model.tracers.S, :, :, kᴺ:kᴺ) + return view(ocean.model.tracers.S.data, :, :, kᴺ:kᴺ) end function ocean_surface_temperature(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) kᴺ = size(ocean.model.grid, 3) - return interior(ocean.model.tracers.T, :, :, kᴺ:kᴺ) + return view(ocean.model.tracers.T.data, :, :, kᴺ:kᴺ) end function ocean_surface_velocities(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) @@ -109,14 +110,18 @@ function ComponentExchanger(ocean::Simulation{<:HydrostaticFreeSurfaceModel}, gr return ComponentExchanger((; u, v, T, S), nothing) end +@inline net_flux(condition) = condition +@inline net_flux(bc::FluxAndRestoring) = bc.flux_field +@inline net_flux(bc::DiscreteBoundaryFunction) = net_flux(bc.func) + function net_fluxes(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) # TODO: Generalize this to work with any ocean model - τˣ = ocean.model.velocities.u.boundary_conditions.top.condition - τʸ = ocean.model.velocities.v.boundary_conditions.top.condition + τˣ = net_flux(ocean.model.velocities.u.boundary_conditions.top.condition) + τʸ = net_flux(ocean.model.velocities.v.boundary_conditions.top.condition) net_ocean_surface_fluxes = (; u=τˣ, v=τʸ) tracers = ocean.model.tracers - ocean_surface_tracer_fluxes = NamedTuple(name => tracers[name].boundary_conditions.top.condition for name in keys(tracers)) + ocean_surface_tracer_fluxes = NamedTuple(name => net_flux(tracers[name].boundary_conditions.top.condition) for name in keys(tracers)) return merge(ocean_surface_tracer_fluxes, net_ocean_surface_fluxes) end diff --git a/src/Oceans/flux_and_restoring.jl b/src/Oceans/flux_and_restoring.jl new file mode 100644 index 000000000..8da53a7f0 --- /dev/null +++ b/src/Oceans/flux_and_restoring.jl @@ -0,0 +1,46 @@ +using Oceananigans.Operators: Δzᶜᶜᶜ + +using Adapt + +""" + FluxAndRestoring(flux_field, restoring) + +A boundary-condition condition (intended to be wrapped in a discrete-form +`FluxBoundaryCondition`) that combines two contributions at a tracer's top +boundary: + +1. `flux_field`: a 2D `Field{Center, Center, Nothing}` that an external flux + solver (e.g. the OMIP coupled atmosphere/sea-ice solver) writes into each + step. This is read at `(i, j, 1)`. + +2. `restoring`: a callable with signature `(i, j, k, grid, clock, fields)` that + returns a tendency in the top cell — typically a `DatasetRestoring`, + evaluating to `r * μ * (ψ_dataset - ψ)`. The tendency is converted to a + surface flux by multiplying by `-Δz` at the top cell, consistent with the + Oceananigans top-flux sign convention (top cell tendency contribution is + `-J / Δz`). + +This lets the coupled flux solver and a dataset restoring share the same top +boundary condition without one clobbering the other. +""" +struct FluxAndRestoring{F, R} <: Function + flux_field :: F + restoring :: R +end + +Adapt.adapt_structure(to, fr::FluxAndRestoring) = + FluxAndRestoring(adapt(to, fr.flux_field), + adapt(to, fr.restoring)) + +@inline function (fr::FluxAndRestoring)(i, j, grid, clock, fields) + Nz = grid.Nz + @inbounds J = fr.flux_field[i, j, 1] + + # Restoring accessed as a tendency forcing (compatible with DatasetRestoring) + G = fr.restoring(i, j, Nz, grid, clock, fields) + + # Top BC convention: tendency contribution = -J / Δz, so to inject + # `G` in the top cell the flux is `-G * Δz`. + Δz = Δzᶜᶜᶜ(i, j, Nz, grid) + return J - G * Δz +end diff --git a/src/Oceans/ocean_simulation.jl b/src/Oceans/ocean_simulation.jl index 88627bb1b..71a24f3a3 100644 --- a/src/Oceans/ocean_simulation.jl +++ b/src/Oceans/ocean_simulation.jl @@ -19,15 +19,18 @@ using Statistics: mean ##### @inline ϕ²(i, j, k, grid, ϕ) = @inbounds ϕ[i, j, k]^2 -@inline spᶠᶜᶜ(i, j, k, grid, Φ) = @inbounds sqrt(Φ.u[i, j, k]^2 + ℑxyᶠᶜᵃ(i, j, k, grid, ϕ², Φ.v)) -@inline spᶜᶠᶜ(i, j, k, grid, Φ) = @inbounds sqrt(Φ.v[i, j, k]^2 + ℑxyᶜᶠᵃ(i, j, k, grid, ϕ², Φ.u)) +@inline spᶠᶜᶜ(i, j, k, grid, Φ, Uᴮ) = @inbounds sqrt(Φ.u[i, j, k]^2 + ℑxyᶠᶜᵃ(i, j, k, grid, ϕ², Φ.v) + Uᴮ^2) +@inline spᶜᶠᶜ(i, j, k, grid, Φ, Uᴮ) = @inbounds sqrt(Φ.v[i, j, k]^2 + ℑxyᶜᶠᵃ(i, j, k, grid, ϕ², Φ.u) + Uᴮ^2) -@inline u_quadratic_bottom_drag(i, j, grid, c, Φ, μ) = @inbounds - μ * Φ.u[i, j, 1] * spᶠᶜᶜ(i, j, 1, grid, Φ) -@inline v_quadratic_bottom_drag(i, j, grid, c, Φ, μ) = @inbounds - μ * Φ.v[i, j, 1] * spᶜᶠᶜ(i, j, 1, grid, Φ) +@inline u_quadratic_bottom_drag(i, j, grid, c, Φ, p) = @inbounds - p.μ * Φ.u[i, j, 1] * spᶠᶜᶜ(i, j, 1, grid, Φ, p.Uᴮ) +@inline v_quadratic_bottom_drag(i, j, grid, c, Φ, p) = @inbounds - p.μ * Φ.v[i, j, 1] * spᶜᶠᶜ(i, j, 1, grid, Φ, p.Uᴮ) # Keep a constant linear drag parameter independent on vertical level -@inline u_immersed_bottom_drag(i, j, k, grid, clock, Φ, μ) = @inbounds - μ * Φ.u[i, j, k] * spᶠᶜᶜ(i, j, k, grid, Φ) -@inline v_immersed_bottom_drag(i, j, k, grid, clock, Φ, μ) = @inbounds - μ * Φ.v[i, j, k] * spᶜᶠᶜ(i, j, k, grid, Φ) +@inline u_immersed_bottom_drag(i, j, k, grid, clock, Φ, p) = @inbounds - p.μ * Φ.u[i, j, k] * spᶠᶜᶜ(i, j, k, grid, Φ, p.Uᴮ) +@inline v_immersed_bottom_drag(i, j, k, grid, clock, Φ, p) = @inbounds - p.μ * Φ.v[i, j, k] * spᶜᶠᶜ(i, j, k, grid, Φ, p.Uᴮ) + +@inline build_top_tracer_bc(flux_field, ::Nothing) = FluxBoundaryCondition(flux_field) +@inline build_top_tracer_bc(flux_field, restoring) = FluxBoundaryCondition(FluxAndRestoring(flux_field, restoring); discrete_form=true) ##### ##### Defaults @@ -100,6 +103,7 @@ end """ ocean_simulation(grid; Δt = estimate_maximum_Δt(grid), + clock = Clock(grid), closure = default_ocean_closure(), tracers = (:T, :S), free_surface = default_free_surface(grid), @@ -107,7 +111,9 @@ end rotation_rate = default_planet_rotation_rate, gravitational_acceleration = default_gravitational_acceleration, bottom_drag_coefficient = Default(0.003), + drag_bulk_velocity = Default(0.1), forcing = NamedTuple(), + surface_restoring = NamedTuple(), biogeochemistry = nothing, timestepper = :SplitRungeKutta3, coriolis = Default(HydrostaticSphericalCoriolis(; rotation_rate)), @@ -126,7 +132,6 @@ consistent defaults for advection, closures, the equation of state, surface flux barotropic pressure–gradient forcing, boundary conditions, and optional biogeochemistry. It then wraps the model into an Oceananigans's `Simulation` with the specified timestepping options. - ## Behaviour and automatic configuration ### Coriolis @@ -161,6 +166,7 @@ defaults on a per-field basis. ## Keyword Arguments - `Δt`: Timestep used by the `Simulation`. Defaults to the maximum stable timestep estimated from the `grid`. +- `clock`: Clock object. Defaults to `Clock(grid)`. - `closure`: A turbulence or mixing closure. Defaults to `default_ocean_closure()`. - `tracers`: Tuple of tracer names. Defaults to `(:T, :S)`. - `free_surface`: Free–surface solver. Defaults to `default_free_surface(grid)`. @@ -168,7 +174,9 @@ defaults on a per-field basis. - `rotation_rate`: Planetary rotation rate used for Coriolis forcing. - `gravitational_acceleration`: Gravitational acceleration, passed to buoyancy. - `bottom_drag_coefficient`: Bottom drag coefficient. May be a `Default` wrapper. +- `drag_bulk_velocity`: a minimum velocity for the bottom drag. - `forcing`: Named tuple of additional forcing(s) for individual fields. +- `surface_restoring`: Named tuple of dataset restorings to apply as part of the tracer top boundary condition. - `biogeochemistry`: A biogeochemical model or `nothing`. - `timestepper`: Time-stepping scheme; options are `:SplitRungeKutta3` (default), or `:QuasiAdamsBashforth2`. - `coriolis`: Coriolis object or `Default(...)` wrapper. @@ -182,6 +190,7 @@ defaults on a per-field basis. """ function ocean_simulation(grid; Δt = estimate_maximum_Δt(grid), + clock = Clock(grid), closure = default_ocean_closure(), tracers = (:T, :S), free_surface = default_free_surface(grid), @@ -189,7 +198,10 @@ function ocean_simulation(grid; rotation_rate = default_planet_rotation_rate, gravitational_acceleration = default_gravitational_acceleration, bottom_drag_coefficient = Default(0.003), + drag_bulk_velocity = Default(0.05), + use_barotropic_potential = true, forcing = NamedTuple(), + surface_restoring = NamedTuple(), biogeochemistry = nothing, timestepper = :SplitRungeKutta3, coriolis = Default(HydrostaticSphericalCoriolis(; rotation_rate)), @@ -198,11 +210,12 @@ function ocean_simulation(grid; equation_of_state = TEOS10EquationOfState(; reference_density), boundary_conditions::NamedTuple = NamedTuple(), radiative_forcing = default_radiative_forcing(grid), + materialize_buoyancy_gradients = true, warn = true, verbose = false) FT = eltype(grid) - + if grid isa RectilinearGrid # turn off Coriolis unless user-supplied coriolis = default_or_override(coriolis, nothing) else @@ -216,6 +229,8 @@ function ocean_simulation(grid; if single_column_simulation # Let users put a bottom drag if they want bottom_drag_coefficient = default_or_override(bottom_drag_coefficient, zero(grid)) + drag_bulk_velocity = default_or_override(drag_bulk_velocity, zero(grid)) + drag_parameters = (μ = bottom_drag_coefficient, Uᴮ = drag_bulk_velocity) # Don't let users use advection in a single column model tracer_advection = nothing @@ -236,21 +251,27 @@ function ocean_simulation(grid; end bottom_drag_coefficient = default_or_override(bottom_drag_coefficient) + drag_bulk_velocity = default_or_override(drag_bulk_velocity) + bottom_drag_coefficient = convert(FT, bottom_drag_coefficient) + drag_bulk_velocity = convert(FT, drag_bulk_velocity) + drag_parameters = (μ = bottom_drag_coefficient, Uᴮ = drag_bulk_velocity) - u_immersed_drag = FluxBoundaryCondition(u_immersed_bottom_drag, discrete_form=true, parameters=bottom_drag_coefficient) - v_immersed_drag = FluxBoundaryCondition(v_immersed_bottom_drag, discrete_form=true, parameters=bottom_drag_coefficient) + u_immersed_drag = FluxBoundaryCondition(u_immersed_bottom_drag, discrete_form=true, parameters=drag_parameters) + v_immersed_drag = FluxBoundaryCondition(v_immersed_bottom_drag, discrete_form=true, parameters=drag_parameters) u_immersed_bc = ImmersedBoundaryCondition(bottom=u_immersed_drag) v_immersed_bc = ImmersedBoundaryCondition(bottom=v_immersed_drag) - # Forcing for u, v - barotropic_potential = Field{Center, Center, Nothing}(grid) - u_forcing = BarotropicPotentialForcing(XDirection(), barotropic_potential) - v_forcing = BarotropicPotentialForcing(YDirection(), barotropic_potential) + if use_barotropic_potential + # Forcing for u, v + barotropic_potential = Field{Center, Center, Nothing}(grid) + u_forcing = BarotropicPotentialForcing(XDirection(), barotropic_potential) + v_forcing = BarotropicPotentialForcing(YDirection(), barotropic_potential) - :u ∈ keys(forcing) && (u_forcing = (u_forcing, forcing[:u])) - :v ∈ keys(forcing) && (v_forcing = (v_forcing, forcing[:v])) - forcing = merge(forcing, (u=u_forcing, v=v_forcing)) + :u ∈ keys(forcing) && (u_forcing = (u_forcing, forcing[:u])) + :v ∈ keys(forcing) && (v_forcing = (v_forcing, forcing[:v])) + forcing = merge(forcing, (u=u_forcing, v=v_forcing)) + end end if !isnothing(radiative_forcing) @@ -262,22 +283,20 @@ function ocean_simulation(grid; forcing = merge(forcing, (; T=T_forcing)) end - bottom_drag_coefficient = convert(FT, bottom_drag_coefficient) - # Set up boundary conditions using Field - top_zonal_momentum_flux = τˣ = Field{Face, Center, Nothing}(grid) - top_meridional_momentum_flux = τʸ = Field{Center, Face, Nothing}(grid) + top_zonal_momentum_flux = τˣ = Field{Face, Center, Nothing}(grid) + top_meridional_momentum_flux = τʸ = Field{Center, Face, Nothing}(grid) top_ocean_heat_flux = Jᵀ = Field{Center, Center, Nothing}(grid) top_salt_flux = Jˢ = Field{Center, Center, Nothing}(grid) # Construct ocean boundary conditions including surface forcing and bottom drag u_top_bc = FluxBoundaryCondition(τˣ) v_top_bc = FluxBoundaryCondition(τʸ) - T_top_bc = FluxBoundaryCondition(Jᵀ) - S_top_bc = FluxBoundaryCondition(Jˢ) + T_top_bc = build_top_tracer_bc(Jᵀ, get(surface_restoring, :T, nothing)) + S_top_bc = build_top_tracer_bc(Jˢ, get(surface_restoring, :S, nothing)) - u_bot_bc = FluxBoundaryCondition(u_quadratic_bottom_drag, discrete_form=true, parameters=bottom_drag_coefficient) - v_bot_bc = FluxBoundaryCondition(v_quadratic_bottom_drag, discrete_form=true, parameters=bottom_drag_coefficient) + u_bot_bc = FluxBoundaryCondition(u_quadratic_bottom_drag, discrete_form=true, parameters=drag_parameters) + v_bot_bc = FluxBoundaryCondition(v_quadratic_bottom_drag, discrete_form=true, parameters=drag_parameters) default_boundary_conditions = (u = FieldBoundaryConditions(top=u_top_bc, bottom=u_bot_bc, immersed=u_immersed_bc), v = FieldBoundaryConditions(top=v_top_bc, bottom=v_bot_bc, immersed=v_immersed_bc), @@ -289,7 +308,8 @@ function ocean_simulation(grid; # conditions even when a user-bc is supplied). boundary_conditions = merge(default_boundary_conditions, boundary_conditions) buoyancy = SeawaterBuoyancy(; gravitational_acceleration, equation_of_state) - + buoyancy = Oceananigans.BuoyancyFormulations.BuoyancyForce(grid, buoyancy; materialize_gradients=materialize_buoyancy_gradients) + if tracer_advection isa NamedTuple tracer_advection = with_tracers(tracers, tracer_advection, default_tracer_advection()) else @@ -297,12 +317,13 @@ function ocean_simulation(grid; end if hasclosure(closure, CATKEVerticalDiffusivity) - # Turn off CATKE tracer advection - tke_advection = (; e=nothing) + # Use the same advection as for temperature + tke_advection = (; e=tracer_advection[1]) tracer_advection = merge(tracer_advection, tke_advection) end ocean_model = HydrostaticFreeSurfaceModel(grid; + clock, buoyancy, closure, biogeochemistry, From c1345d7dd31de06a707fece2288ac9d4dbd00d08 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 15 Apr 2026 12:50:01 +0200 Subject: [PATCH 3/9] Add temperature/snow-dependent sea ice albedo (CCSM3 scheme) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New type SeaIceAlbedo implementing the CCSM3 albedo parameterization (Briegleb et al. 2004): - Broadband albedos: bare ice 0.54, snow-covered 0.83 - Temperature-dependent reduction near melting (implicit melt ponds) - Thin-ice transition to ocean albedo below 0.5 m - Snow cover blending (full snow albedo at hs > 0.02 m) - Handles nothing snow_thickness gracefully Also evaluates state-dependent albedo to scalar via stateindex at the top of the atmosphere-sea ice flux kernel, fixing a potential struct-in-arithmetic bug when the albedo is not a plain number. Changes to atmosphere_sea_ice_fluxes.jl: - Evaluate radiation properties (α, ϵ) via stateindex before iteration - Pass local_interface_properties with scalar values to compute_interface_state - Add hs (snow thickness) to local_interior_state - Rename h → hi for ice thickness in interior state References: - Briegleb et al. (2004): NCAR Tech Note - Briegleb & Light (2007): NCAR/TN-472+STR --- .../InterfaceComputations.jl | 10 ++ .../atmosphere_sea_ice_fluxes.jl | 24 +++- .../InterfaceComputations/sea_ice_albedo.jl | 133 ++++++++++++++++++ 3 files changed, 161 insertions(+), 6 deletions(-) create mode 100644 src/EarthSystemModels/InterfaceComputations/sea_ice_albedo.jl diff --git a/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl b/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl index ed93f8c9c..0f543cf39 100644 --- a/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl +++ b/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl @@ -9,14 +9,21 @@ export Radiation, ComponentInterfaces, LatitudeDependentAlbedo, + SeaIceAlbedo, SimilarityTheoryFluxes, MomentumRoughnessLength, ScalarRoughnessLength, + NCARMomentumRoughnessLength, + NCARScalarRoughnessLength, + NCARBulkFluxes, CoefficientBasedFluxes, SkinTemperature, BulkTemperature, + LinearStableStabilityFunction, + COARELogarithmicSimilarityProfile, atmosphere_ocean_stability_functions, atmosphere_sea_ice_stability_functions, + ncar_stability_functions, compute_atmosphere_ocean_fluxes!, compute_atmosphere_sea_ice_fluxes!, compute_sea_ice_ocean_fluxes!, @@ -67,13 +74,16 @@ end include("radiation.jl") include("latitude_dependent_albedo.jl") include("tabulated_albedo.jl") +include("sea_ice_albedo.jl") # Turbulent fluxes include("roughness_lengths.jl") +include("ncar_roughness_lengths.jl") include("interface_states.jl") include("compute_interface_state.jl") include("similarity_theory_turbulent_fluxes.jl") include("coefficient_based_turbulent_fluxes.jl") +include("ncar_bulk_fluxes.jl") # State exchanger and interfaces include("state_exchanger.jl") diff --git a/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl index 1292d6311..667a92804 100644 --- a/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl @@ -84,18 +84,30 @@ end # Sea ice properties uˢⁱ = zero(FT) # ℑxᶜᵃᵃ(i, j, 1, grid, interior_state.u) vˢⁱ = zero(FT) # ℑyᵃᶜᵃ(i, j, 1, grid, interior_state.v) - hˢⁱ = interior_state.h[i, j, 1] + hˢⁱ = interior_state.hi[i, j, 1] + hˢⁿ = interior_state.hs[i, j, 1] hc = interior_state.hc[i, j, 1] ℵᵢ = interior_state.ℵ[i, j, 1] Tₛ = interface_temperature[i, j, 1] Tₛ = convert_to_kelvin(sea_ice_properties.temperature_units, Tₛ) end + # Evaluate state-dependent radiation properties at this grid point. + # The albedo may be a struct (e.g., CCSM3SeaIceAlbedo) that reads model fields; + # we evaluate it here so the iteration uses a scalar. + time = Time(clock.time) + σ = interface_properties.radiation.σ + α = stateindex(interface_properties.radiation.α, i, j, kᴺ, grid, time, CCC) + ϵ = stateindex(interface_properties.radiation.ϵ, i, j, kᴺ, grid, time, CCC) + local_radiation = (; σ, α, ϵ) + local_interface_properties = InterfaceProperties(local_radiation, + interface_properties.specific_humidity_formulation, + interface_properties.temperature_formulation, + interface_properties.velocity_formulation) + # Build thermodynamic and dynamic states in the atmosphere and interface. - # Notation: - # ⋅ 𝒰 ≡ "dynamic" state vector (thermodynamics + reference height + velocity) ℂᵃᵗ = atmosphere_properties.thermodynamics_parameters - zᵃᵗ = atmosphere_properties.surface_layer_height # elevation of atmos variables relative to interface + zᵃᵗ = atmosphere_properties.surface_layer_height local_atmosphere_state = (z = zᵃᵗ, u = uᵃᵗ, @@ -106,7 +118,7 @@ end h_bℓ = atmosphere_state.h_bℓ) downwelling_radiation = (; ℐꜜˢʷ, ℐꜜˡʷ) - local_interior_state = (u=uˢⁱ, v=vˢⁱ, T=Tᵒᶜ, S=Sᵒᶜ, h=hˢⁱ, hc=hc) + local_interior_state = (u=uˢⁱ, v=vˢⁱ, T=Tᵒᶜ, S=Sᵒᶜ, hi=hˢⁱ, hs=hˢⁿ, hc=hc) # Estimate initial interface state (FP32 compatible) u★ = convert(FT, 1f-4) @@ -132,7 +144,7 @@ end local_atmosphere_state, local_interior_state, downwelling_radiation, - interface_properties, + local_interface_properties, atmosphere_properties, sea_ice_properties) end diff --git a/src/EarthSystemModels/InterfaceComputations/sea_ice_albedo.jl b/src/EarthSystemModels/InterfaceComputations/sea_ice_albedo.jl new file mode 100644 index 000000000..ab2fc8442 --- /dev/null +++ b/src/EarthSystemModels/InterfaceComputations/sea_ice_albedo.jl @@ -0,0 +1,133 @@ +""" + SeaIceAlbedo{FT, HI, HS, TS} + +Sea ice albedo parameterization following the CCSM3 scheme (Briegleb et al. 2004). + +Computes broadband albedo as a function of ice thickness, snow depth, and surface +temperature. The scheme blends between bare ice and snow-covered albedos, with +a temperature-dependent reduction near the melting point to implicitly represent +melt pond formation. + +Algorithm: +1. Base cold albedos: bare ice (0.53) and snow-covered (0.82) +2. Temperature reduction within 1C of melting: Δα_ice = 0.075, Δα_snow = 0.10 +3. Thin-ice transition to ocean albedo below h_amin = 0.5 m +4. Snow cover interpolation: full snow albedo at h_snow > h_smin = 0.02 m + +References: +- Briegleb, B.P., C.M. Bitz, E.C. Hunke, W.H. Lipscomb, and M.M. Schramm (2004): + Scientific description of the sea ice component in CCSM3. NCAR Tech Note. +- Briegleb, B.P. and B. Light (2007): NCAR/TN-472+STR. +""" +struct SeaIceAlbedo{FT, HI, HS, TS} + # Cold base albedos (broadband, approx 0.52 * vis + 0.48 * nir) + ice_albedo :: FT # 0.52*0.73 + 0.48*0.33 = 0.538 ≈ 0.54 + snow_albedo :: FT # 0.52*0.96 + 0.48*0.68 = 0.825 ≈ 0.83 + # Melt reduction + ice_melt_reduction :: FT # 0.075 + snow_melt_reduction :: FT # 0.10 + melting_temperature :: FT # 0 C + temperature_range :: FT # 1 C + # Thickness scales + ocean_albedo :: FT # 0.06 + minimum_ice_thickness :: FT # 0.5 m + minimum_snow_depth :: FT # 0.02 m + # References to model fields + ice_thickness :: HI + snow_thickness :: HS + surface_temperature :: TS +end + +Adapt.adapt_structure(to, α::SeaIceAlbedo) = + SeaIceAlbedo(α.ice_albedo, + α.snow_albedo, + α.ice_melt_reduction, + α.snow_melt_reduction, + α.melting_temperature, + α.temperature_range, + α.ocean_albedo, + α.minimum_ice_thickness, + α.minimum_snow_depth, + Adapt.adapt(to, α.ice_thickness), + Adapt.adapt(to, α.snow_thickness), + Adapt.adapt(to, α.surface_temperature)) + +""" + SeaIceAlbedo(ice_thickness, snow_thickness, surface_temperature; + ice_albedo = 0.54, + snow_albedo = 0.83, + ice_melt_reduction = 0.075, + snow_melt_reduction = 0.10, + melting_temperature = 0.0, + temperature_range = 1.0, + ocean_albedo = 0.06, + minimum_ice_thickness = 0.5, + minimum_snow_depth = 0.02) + +Construct a CCSM3 sea ice albedo parameterization. Requires references to the sea ice +model's ice thickness, snow thickness, and surface temperature fields. + +Broadband albedos are approximate averages of the visible and near-IR bands +weighted by solar spectrum (52% visible, 48% near-IR): +- ice_albedo ≈ 0.52 x 0.73 + 0.48 x 0.33 ≈ 0.54 +- snow_albedo ≈ 0.52 x 0.96 + 0.48 x 0.68 ≈ 0.83 +""" +function SeaIceAlbedo(ice_thickness, snow_thickness, surface_temperature; + FT = Float64, + ice_albedo = 0.54, + snow_albedo = 0.83, + ice_melt_reduction = 0.075, + snow_melt_reduction = 0.10, + melting_temperature = 0.0, + temperature_range = 1.0, + ocean_albedo = 0.06, + minimum_ice_thickness = 0.5, + minimum_snow_depth = 0.02) + + return SeaIceAlbedo(convert(FT, ice_albedo), + convert(FT, snow_albedo), + convert(FT, ice_melt_reduction), + convert(FT, snow_melt_reduction), + convert(FT, melting_temperature), + convert(FT, temperature_range), + convert(FT, ocean_albedo), + convert(FT, minimum_ice_thickness), + convert(FT, minimum_snow_depth), + ice_thickness, + snow_thickness, + surface_temperature) +end + +Base.summary(::SeaIceAlbedo{FT}) where FT = "SeaIceAlbedo{$FT}" +Base.show(io::IO, α::SeaIceAlbedo{FT}) where FT = + print(io, "SeaIceAlbedo{$FT}(ice=", α.ice_albedo, + ", snow=", α.snow_albedo, ")") + +@inline function stateindex(α::SeaIceAlbedo, i, j, k, grid, time, loc, args...) + @inbounds hi = α.ice_thickness[i, j, 1] + @inbounds Ts = α.surface_temperature[i, j, 1] + + # Snow thickness: may be nothing (no snow model) + hs = get_snow_thickness(α.snow_thickness, i, j) + + # Temperature-dependent reduction (implicit melt ponds) + Tm = α.melting_temperature + ΔT = α.temperature_range + fT = clamp((Ts - Tm + ΔT) / ΔT, zero(Ts), one(Ts)) + + αi = α.ice_albedo - α.ice_melt_reduction * fT + αs = α.snow_albedo - α.snow_melt_reduction * fT + + # Thin ice → transition to ocean albedo + αo = α.ocean_albedo + fh = clamp(hi / α.minimum_ice_thickness, zero(hi), one(hi)) + αi = αo + (αi - αo) * fh + + # Snow cover blending + fs = clamp(hs / α.minimum_snow_depth, zero(hs), one(hs)) + return fs * αs + (1 - fs) * αi +end + +# Helper to handle nothing snow thickness (no snow model) +@inline get_snow_thickness(hs::Nothing, i, j, grid) = zero(grid) +@inline get_snow_thickness(hs, i, j, grid) = @inbounds hs[i, j, 1] From ac1b6a20fd635328336cbca8910f9f74374549a7 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 15 Apr 2026 14:31:17 +0200 Subject: [PATCH 4/9] Add Adapt for InterfaceProperties, fix get_snow_thickness call - Add Adapt.adapt_structure for InterfaceProperties so that SeaIceAlbedo Fields are properly adapted for GPU kernels - Fix get_snow_thickness call to pass grid argument Co-Authored-By: Claude Opus 4.6 (1M context) --- .../InterfaceComputations/interface_states.jl | 6 ++++++ .../InterfaceComputations/sea_ice_albedo.jl | 2 +- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/EarthSystemModels/InterfaceComputations/interface_states.jl b/src/EarthSystemModels/InterfaceComputations/interface_states.jl index 4f1ae6c2d..22b8cb1d2 100644 --- a/src/EarthSystemModels/InterfaceComputations/interface_states.jl +++ b/src/EarthSystemModels/InterfaceComputations/interface_states.jl @@ -16,6 +16,12 @@ struct InterfaceProperties{R, Q, T, V} velocity_formulation :: V end +Adapt.adapt_structure(to, p::InterfaceProperties) = + InterfaceProperties(Adapt.adapt(to, p.radiation), + Adapt.adapt(to, p.specific_humidity_formulation), + Adapt.adapt(to, p.temperature_formulation), + Adapt.adapt(to, p.velocity_formulation)) + ##### ##### Interface specific humidity formulations ##### diff --git a/src/EarthSystemModels/InterfaceComputations/sea_ice_albedo.jl b/src/EarthSystemModels/InterfaceComputations/sea_ice_albedo.jl index ab2fc8442..9ffc88417 100644 --- a/src/EarthSystemModels/InterfaceComputations/sea_ice_albedo.jl +++ b/src/EarthSystemModels/InterfaceComputations/sea_ice_albedo.jl @@ -108,7 +108,7 @@ Base.show(io::IO, α::SeaIceAlbedo{FT}) where FT = @inbounds Ts = α.surface_temperature[i, j, 1] # Snow thickness: may be nothing (no snow model) - hs = get_snow_thickness(α.snow_thickness, i, j) + hs = get_snow_thickness(α.snow_thickness, i, j, grid) # Temperature-dependent reduction (implicit melt ponds) Tm = α.melting_temperature From b2052f32adf257db2fa0d87566953f075a312365 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 15 Apr 2026 14:56:30 +0200 Subject: [PATCH 5/9] 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 """ From 5027bb7ed347c094f1924d4f126a139c1363f301 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 15 Apr 2026 15:12:49 +0200 Subject: [PATCH 6/9] Remove unrelated changes, add SeaIceAlbedo tests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Strip out changes that belong to other PRs: FluxAndRestoring, bottom drag rework, barotropic potential toggle, meridional heat transport deletion, bathymetry/inpainting/restoring changes, NCAR bulk flux includes/exports, snow-specific h→hi rename and hs additions. Keep only sea ice albedo: SeaIceAlbedo struct, stateindex evaluation in atmosphere-sea-ice kernel, Adapt for InterfaceProperties, export/include. Add test_sea_ice_albedo.jl covering: constructor, cold/melting ice, thin ice→ocean transition, snow cover blending, nothing snow thickness, Adapt.adapt_structure for both SeaIceAlbedo and InterfaceProperties. Co-Authored-By: Claude Opus 4.6 (1M context) --- examples/meridional_heat_transport_ecco.jl | 119 +++++++++++++++ src/Bathymetry/regrid_bathymetry.jl | 3 +- src/DataWrangling/inpainting.jl | 13 +- src/DataWrangling/restoring.jl | 2 - src/Diagnostics/Diagnostics.jl | 4 +- src/Diagnostics/interface_fluxes.jl | 9 +- src/Diagnostics/meridional_heat_transport.jl | 95 ++++++++++++ src/Diagnostics/mixed_layer_depth.jl | 1 + .../InterfaceComputations.jl | 8 - .../atmosphere_sea_ice_fluxes.jl | 5 +- .../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 +- src/NumericalEarth.jl | 3 +- src/Oceans/Oceans.jl | 19 +-- src/Oceans/flux_and_restoring.jl | 46 ------ src/Oceans/ocean_simulation.jl | 77 ++++------ ...diagnostics_2.jl => test_diagnostics_2.jl} | 4 +- test/test_sea_ice_albedo.jl | 141 ++++++++++++++++++ 20 files changed, 425 insertions(+), 154 deletions(-) create mode 100755 examples/meridional_heat_transport_ecco.jl create mode 100644 src/Diagnostics/meridional_heat_transport.jl delete mode 100644 src/Oceans/flux_and_restoring.jl rename test/{tet_diagnostics_2.jl => test_diagnostics_2.jl} (97%) create mode 100644 test/test_sea_ice_albedo.jl diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl new file mode 100755 index 000000000..58b3e7849 --- /dev/null +++ b/examples/meridional_heat_transport_ecco.jl @@ -0,0 +1,119 @@ +using NumericalEarth +using Oceananigans +using Oceananigans.Units +using Dates +using Statistics +using Printf + +using CUDA; CUDA.device!(3) + +arch = GPU() +Nx = 360 +Ny = 180 +Nz = 50 + +depth = 5000meters +z = ExponentialDiscretization(Nz, -depth, 0; scale = depth/4) + +underlying_grid = TripolarGrid(arch; size = (Nx, Ny, Nz), halo = (5, 5, 4), z) +underlying_grid = LatitudeLongitudeGrid(arch; size = (Nx, Ny, Nz), halo = (5, 5, 4), z, longitude = (0, 360), latitude = (-80, 80)) +bottom_height = regrid_bathymetry(underlying_grid; + minimum_depth = 10, + interpolation_passes = 10, + major_basins = 2) +grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height); + active_cells_map=true) + +free_surface = SplitExplicitFreeSurface(grid; substeps=70) +momentum_advection = WENOVectorInvariant(order=5) +tracer_advection = WENO(order=5) +vertical_mixing = NumericalEarth.Oceans.default_ocean_closure() +ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surface, + closure=(vertical_mixing,)) +sea_ice = sea_ice_simulation(grid, ocean; advection=tracer_advection) + +date = DateTime(1993, 1, 1) +dataset = ECCO4Monthly() +ecco_temperature = Metadatum(:temperature; date, dataset) +ecco_salinity = Metadatum(:salinity; date, dataset) +ecco_sea_ice_thickness = Metadatum(:sea_ice_thickness; date, dataset) +ecco_sea_ice_concentration = Metadatum(:sea_ice_concentration; date, dataset) + +set!(ocean.model, T=ecco_temperature, S=ecco_salinity) +set!(sea_ice.model, h=ecco_sea_ice_thickness, ℵ=ecco_sea_ice_concentration) + +radiation = Radiation(arch) +atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(80), + include_rivers_and_icebergs = false) +esm = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + +simulation = Simulation(esm; Δt=20minutes, stop_time=5*365days) + +wall_time = Ref(time_ns()) + +function progress(sim) + ocean = sim.model.ocean + u, v, w = ocean.model.velocities + T = ocean.model.tracers.T + e = ocean.model.tracers.e + Tmin, Tmax, Tavg = minimum(T), maximum(T), mean(view(T, :, :, ocean.model.grid.Nz)) + emax = maximum(e) + umax = (maximum(abs, u), maximum(abs, v), maximum(abs, w)) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iter: %d", prettytime(sim), iteration(sim)) + msg2 = @sprintf(", max|uo|: (%.1e, %.1e, %.1e) m s⁻¹", umax...) + msg3 = @sprintf(", max(e): %.2f m² s⁻²", emax) + msg4 = @sprintf(", wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg3 * msg4 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(simulation, progress, IterationInterval(200)) + +mht = Field(meridional_heat_transport(esm)) + +ocean.output_writers[:mth] = JLD2Writer(ocean.model, (; mht); + schedule = TimeInterval(3hours), + filename = "ocean_one_degree_mht", + overwrite_existing = true) + +run!(simulation) + +## + +using Oceananigans + +mht = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht"; backend = OnDisk()) + +times = mht.times +Nt = length(times) + +grid = mht.grid +Ny = size(mht.grid, 2) + +mht_mean = deepcopy(mht[1][1, :, 1]) + +for iter in 1:Nt + @info "iteration $iter out of $Nt" + mht_mean += mht[iter][1, :, 1] +end + +@. mht_mean = mht_mean / Nt + +using CairoMakie + +fig = Figure() +ax = Axis(fig[1, 1], xlabel="latitude (deg)", ylabel="MHT (PW)") + +φ = φnodes(grid, Face()) + +lines!(ax, φ, mht_mean[1:Ny+1] / 1e15, linewidth=4) + +save("mht.png", fig) diff --git a/src/Bathymetry/regrid_bathymetry.jl b/src/Bathymetry/regrid_bathymetry.jl index 564d531f3..430cad426 100644 --- a/src/Bathymetry/regrid_bathymetry.jl +++ b/src/Bathymetry/regrid_bathymetry.jl @@ -341,8 +341,7 @@ end # Fix active cells to be at least `-minimum_depth`. active = z < 0 # it's a wet cell - above_minimum_depth = z > -minimum_depth - z = ifelse(active, ifelse(above_minimum_depth, zero(z), z), z) + z = ifelse(active, min(z, -minimum_depth), z) @inbounds target_z[i, j, 1] = z end 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 """ diff --git a/src/DataWrangling/restoring.jl b/src/DataWrangling/restoring.jl index d26b7c1f4..918edfae0 100644 --- a/src/DataWrangling/restoring.jl +++ b/src/DataWrangling/restoring.jl @@ -12,7 +12,6 @@ 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 @@ -235,7 +234,6 @@ 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/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl index 4f698a4b1..3c7783ae0 100644 --- a/src/Diagnostics/Diagnostics.jl +++ b/src/Diagnostics/Diagnostics.jl @@ -1,6 +1,7 @@ module Diagnostics export MixedLayerDepthField, MixedLayerDepthOperand +export meridional_heat_transport export frazil_temperature_flux, net_ocean_temperature_flux, sea_ice_ocean_temperature_flux, atmosphere_ocean_temperature_flux, frazil_heat_flux, net_ocean_heat_flux, sea_ice_ocean_heat_flux, atmosphere_ocean_heat_flux, net_ocean_salinity_flux, sea_ice_ocean_salinity_flux, atmosphere_ocean_salinity_flux, @@ -14,13 +15,12 @@ using Oceananigans.BoundaryConditions: FieldBoundaryConditions, fill_halo_region using Oceananigans.Fields: FieldStatus using Oceananigans.Utils: launch! using KernelAbstractions: @index, @kernel -using Oceananigans.BoundaryConditions: DiscreteBoundaryFunction using NumericalEarth.EarthSystemModels: EarthSystemModel -using NumericalEarth.Oceans: FluxAndRestoring import Oceananigans.Fields: compute! include("mixed_layer_depth.jl") +include("meridional_heat_transport.jl") include("interface_fluxes.jl") end # module diff --git a/src/Diagnostics/interface_fluxes.jl b/src/Diagnostics/interface_fluxes.jl index 2d5b80758..4d144d84c 100644 --- a/src/Diagnostics/interface_fluxes.jl +++ b/src/Diagnostics/interface_fluxes.jl @@ -1,8 +1,4 @@ -@inline flux_field(condition) = condition -@inline flux_field(bc::FluxAndRestoring) = bc.flux_field -@inline flux_field(bc::DiscreteBoundaryFunction) = flux_field(bc.func) - ########################### ### Temperature fluxes ########################### @@ -25,11 +21,12 @@ end Return the net temperature flux (K m s⁻¹) at the ocean's surface in a coupled `esm`. """ function net_ocean_temperature_flux(esm::EarthSystemModel) - Jᵀ = flux_field(esm.ocean.model.tracers.T.boundary_conditions.top.condition) + Jᵀ = esm.ocean.model.tracers.T.boundary_conditions.top.condition net_ocean_temperature_flux = Jᵀ + frazil_temperature_flux(esm) return Field(net_ocean_temperature_flux) end + """ sea_ice_ocean_temperature_flux(esm::EarthSystemModel) @@ -119,7 +116,7 @@ end Return the net salinity flux (g/kg m s⁻¹) at the ocean's surface in a coupled `esm`. """ function net_ocean_salinity_flux(esm::EarthSystemModel) - Jˢ = flux_field(esm.ocean.model.tracers.S.boundary_conditions.top.condition) + Jˢ = esm.ocean.model.tracers.S.boundary_conditions.top.condition return Jˢ end diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl new file mode 100644 index 000000000..6bc6338db --- /dev/null +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -0,0 +1,95 @@ +using ..EarthSystemModels: EarthSystemModel, reference_density, heat_capacity + +""" + meridional_heat_transport(esm::EarthSystemModel; + reference_temperature = 0) + +Return the meridional heat transport for the coupled `esm::EarthSystemModel` by computing +the meridional heat flux. + +The meridional heat transport is computed via: + +```math +\\mathrm{MHT} ≡ ρᵒᶜ cᵒᶜ ∫ v (T - T_{\\rm ref}) \\, \\mathrm{d}x \\, \\mathrm{d}z +``` + +Above, ``T_{\\rm ref}`` is a reference temperature and ``ρᵒᶜ`` and ``cᵒᶜ`` are the +ocean reference density and specific heat capacity respectively. + +!!! warning "Only works on LatitudeLongitudeGrid" + + The `meridional_heat_transport` diagnostic currently is only supported only on + `LongitudeLatitudeGrid`s. + +Arguments +========= + +* `esm`: An EarthSystemModel. + + +Keyword Arguments +================= + +* `reference_temperature`: The reference temperature (in ᵒC) used for the calculation; default: 0 ᵒC. + + !!! info "Reference temperature" + + The reference temperature is only relevant when we compute the meridional heat transport over a section + where there is a net volume transport. If we are computing the diagnostic globally, i.e., around a whole + latitude circle, then by necessity there is no net volume transport and thus the reference temperature + value is irrelevant. Section-averaged transport could also be considered as a reference temperature to + remove residual barotropic volume fluxes in basin-scale/regional analyses where a net volume transport + is present. + +Example +======= + +```jldoctest +using NumericalEarth +using Oceananigans + +grid = RectilinearGrid(size = (4, 5, 2), extent = (1, 1, 1), + topology = (Periodic, Bounded, Bounded)) + +ocean = ocean_simulation(grid; + momentum_advection = nothing, + tracer_advection = nothing, + closure = nothing, + coriolis = nothing) + +sea_ice = sea_ice_simulation(grid, ocean) + +atmosphere = PrescribedAtmosphere(grid, [0.0]) + +esm = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation = Radiation()) + +mht = meridional_heat_transport(esm) + +# output + +Integral of BinaryOperation at (Center, Face, Center) over dims (1, 3) +└── operand: BinaryOperation at (Center, Face, Center) + └── grid: 4×5×2 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×2 halo +``` +""" +function meridional_heat_transport(esm::EarthSystemModel; reference_temperature=0) + + grid = esm.ocean.model.grid + + validation_grid = grid isa ImmersedBoundaryGrid ? grid.underlying_grid : grid + + grid isa OrthogonalSphericalShellGrid && + throw(ArgumentError("meridional_heat_transport diagnostic does not work on OrthogonalSphericalShellGrid at the moment; use LatitudeLongitudeGrid.")) + + FT = eltype(esm) + reference_temperature = convert(FT, reference_temperature) + + ρᵒᶜ = reference_density(esm.ocean) + cᵒᶜ = heat_capacity(esm.ocean) + + T = esm.ocean.model.tracers.T + v = esm.ocean.model.velocities.v + + MHT = Integral(ρᵒᶜ * cᵒᶜ * v * (T - reference_temperature), dims=(1, 3)) + return MHT +end diff --git a/src/Diagnostics/mixed_layer_depth.jl b/src/Diagnostics/mixed_layer_depth.jl index 3faeb750a..0b986312b 100644 --- a/src/Diagnostics/mixed_layer_depth.jl +++ b/src/Diagnostics/mixed_layer_depth.jl @@ -32,6 +32,7 @@ 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/InterfaceComputations.jl b/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl index 0f543cf39..0306aca4b 100644 --- a/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl +++ b/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl @@ -13,17 +13,11 @@ export SimilarityTheoryFluxes, MomentumRoughnessLength, ScalarRoughnessLength, - NCARMomentumRoughnessLength, - NCARScalarRoughnessLength, - NCARBulkFluxes, CoefficientBasedFluxes, SkinTemperature, BulkTemperature, - LinearStableStabilityFunction, - COARELogarithmicSimilarityProfile, atmosphere_ocean_stability_functions, atmosphere_sea_ice_stability_functions, - ncar_stability_functions, compute_atmosphere_ocean_fluxes!, compute_atmosphere_sea_ice_fluxes!, compute_sea_ice_ocean_fluxes!, @@ -78,12 +72,10 @@ include("sea_ice_albedo.jl") # Turbulent fluxes include("roughness_lengths.jl") -include("ncar_roughness_lengths.jl") include("interface_states.jl") include("compute_interface_state.jl") include("similarity_theory_turbulent_fluxes.jl") include("coefficient_based_turbulent_fluxes.jl") -include("ncar_bulk_fluxes.jl") # State exchanger and interfaces include("state_exchanger.jl") diff --git a/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl index 667a92804..f18078a2d 100644 --- a/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl @@ -84,8 +84,7 @@ end # Sea ice properties uˢⁱ = zero(FT) # ℑxᶜᵃᵃ(i, j, 1, grid, interior_state.u) vˢⁱ = zero(FT) # ℑyᵃᶜᵃ(i, j, 1, grid, interior_state.v) - hˢⁱ = interior_state.hi[i, j, 1] - hˢⁿ = interior_state.hs[i, j, 1] + hˢⁱ = interior_state.h[i, j, 1] hc = interior_state.hc[i, j, 1] ℵᵢ = interior_state.ℵ[i, j, 1] Tₛ = interface_temperature[i, j, 1] @@ -118,7 +117,7 @@ end h_bℓ = atmosphere_state.h_bℓ) downwelling_radiation = (; ℐꜜˢʷ, ℐꜜˡʷ) - local_interior_state = (u=uˢⁱ, v=vˢⁱ, T=Tᵒᶜ, S=Sᵒᶜ, hi=hˢⁱ, hs=hˢⁿ, hc=hc) + local_interior_state = (u=uˢⁱ, v=vˢⁱ, T=Tᵒᶜ, S=Sᵒᶜ, h=hˢⁱ, hc=hc) # Estimate initial interface state (FP32 compatible) u★ = convert(FT, 1f-4) diff --git a/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl index 5cbffb333..a715493d3 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 915fa64b8..ebf4f4ce8 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. -- [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. +- [hieronymus2021comparison](@citet): Hieronymus, M., et al. (2021). A comparison of ocean-ice flux parametrizations. + *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 [shi2021sensitivity](@citet) with ``R = \\alpha_h / \\alpha_s = 35``. +Default values follow [hieronymus2021comparison](@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 db5ff93dd..90469dc04 100644 --- a/src/EarthSystemModels/earth_system_model.jl +++ b/src/EarthSystemModels/earth_system_model.jl @@ -44,12 +44,13 @@ 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) = [] +default_included_properties(::ESM) = tuple() 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 c486171b3..0ec04e6e2 100644 --- a/src/EarthSystemModels/time_step_earth_system_model.jl +++ b/src/EarthSystemModels/time_step_earth_system_model.jl @@ -15,8 +15,13 @@ 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) - !isnothing(ocean) && time_step!(ocean, Δt) + !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(atmosphere) && time_step!(atmosphere, Δt) # TODO: diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl index 8b23f354a..1421eadc5 100644 --- a/src/NumericalEarth.jl +++ b/src/NumericalEarth.jl @@ -55,7 +55,8 @@ export frazil_temperature_flux, net_ocean_temperature_flux, sea_ice_ocean_temperature_flux, atmosphere_ocean_temperature_flux, frazil_heat_flux, net_ocean_heat_flux, sea_ice_ocean_heat_flux, atmosphere_ocean_heat_flux, net_ocean_salinity_flux, sea_ice_ocean_salinity_flux, atmosphere_ocean_salinity_flux, - net_ocean_freshwater_flux, sea_ice_ocean_freshwater_flux, atmosphere_ocean_freshwater_flux + net_ocean_freshwater_flux, sea_ice_ocean_freshwater_flux, atmosphere_ocean_freshwater_flux, + meridional_heat_transport using Oceananigans using Oceananigans.Operators: ℑxyᶠᶜᵃ, ℑxyᶜᶠᵃ diff --git a/src/Oceans/Oceans.jl b/src/Oceans/Oceans.jl index 87093e37f..93b1d2727 100644 --- a/src/Oceans/Oceans.jl +++ b/src/Oceans/Oceans.jl @@ -1,13 +1,13 @@ module Oceans -export ocean_simulation, SlabOcean, FluxAndRestoring +export ocean_simulation, SlabOcean using Oceananigans using Oceananigans.Units using Oceananigans.Utils using Oceananigans.Utils: with_tracers using Oceananigans.Advection: FluxFormAdvection -using Oceananigans.BoundaryConditions: DefaultBoundaryCondition, DiscreteBoundaryFunction +using Oceananigans.BoundaryConditions: DefaultBoundaryCondition using Oceananigans.ImmersedBoundaries: immersed_peripheral_node, inactive_node, MutableGridOfSomeKind using Oceananigans.OrthogonalSphericalShellGrids using Oceananigans.Operators @@ -62,7 +62,6 @@ default_or_override(override, alternative_default=nothing) = override include("slab_ocean.jl") include("barotropic_potential_forcing.jl") include("radiative_forcing.jl") -include("flux_and_restoring.jl") include("ocean_simulation.jl") include("assemble_net_ocean_fluxes.jl") @@ -75,12 +74,12 @@ ocean_temperature(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = ocean.mode function ocean_surface_salinity(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) kᴺ = size(ocean.model.grid, 3) - return view(ocean.model.tracers.S.data, :, :, kᴺ:kᴺ) + return interior(ocean.model.tracers.S, :, :, kᴺ:kᴺ) end function ocean_surface_temperature(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) kᴺ = size(ocean.model.grid, 3) - return view(ocean.model.tracers.T.data, :, :, kᴺ:kᴺ) + return interior(ocean.model.tracers.T, :, :, kᴺ:kᴺ) end function ocean_surface_velocities(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) @@ -110,18 +109,14 @@ function ComponentExchanger(ocean::Simulation{<:HydrostaticFreeSurfaceModel}, gr return ComponentExchanger((; u, v, T, S), nothing) end -@inline net_flux(condition) = condition -@inline net_flux(bc::FluxAndRestoring) = bc.flux_field -@inline net_flux(bc::DiscreteBoundaryFunction) = net_flux(bc.func) - function net_fluxes(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) # TODO: Generalize this to work with any ocean model - τˣ = net_flux(ocean.model.velocities.u.boundary_conditions.top.condition) - τʸ = net_flux(ocean.model.velocities.v.boundary_conditions.top.condition) + τˣ = ocean.model.velocities.u.boundary_conditions.top.condition + τʸ = ocean.model.velocities.v.boundary_conditions.top.condition net_ocean_surface_fluxes = (; u=τˣ, v=τʸ) tracers = ocean.model.tracers - ocean_surface_tracer_fluxes = NamedTuple(name => net_flux(tracers[name].boundary_conditions.top.condition) for name in keys(tracers)) + ocean_surface_tracer_fluxes = NamedTuple(name => tracers[name].boundary_conditions.top.condition for name in keys(tracers)) return merge(ocean_surface_tracer_fluxes, net_ocean_surface_fluxes) end diff --git a/src/Oceans/flux_and_restoring.jl b/src/Oceans/flux_and_restoring.jl deleted file mode 100644 index 8da53a7f0..000000000 --- a/src/Oceans/flux_and_restoring.jl +++ /dev/null @@ -1,46 +0,0 @@ -using Oceananigans.Operators: Δzᶜᶜᶜ - -using Adapt - -""" - FluxAndRestoring(flux_field, restoring) - -A boundary-condition condition (intended to be wrapped in a discrete-form -`FluxBoundaryCondition`) that combines two contributions at a tracer's top -boundary: - -1. `flux_field`: a 2D `Field{Center, Center, Nothing}` that an external flux - solver (e.g. the OMIP coupled atmosphere/sea-ice solver) writes into each - step. This is read at `(i, j, 1)`. - -2. `restoring`: a callable with signature `(i, j, k, grid, clock, fields)` that - returns a tendency in the top cell — typically a `DatasetRestoring`, - evaluating to `r * μ * (ψ_dataset - ψ)`. The tendency is converted to a - surface flux by multiplying by `-Δz` at the top cell, consistent with the - Oceananigans top-flux sign convention (top cell tendency contribution is - `-J / Δz`). - -This lets the coupled flux solver and a dataset restoring share the same top -boundary condition without one clobbering the other. -""" -struct FluxAndRestoring{F, R} <: Function - flux_field :: F - restoring :: R -end - -Adapt.adapt_structure(to, fr::FluxAndRestoring) = - FluxAndRestoring(adapt(to, fr.flux_field), - adapt(to, fr.restoring)) - -@inline function (fr::FluxAndRestoring)(i, j, grid, clock, fields) - Nz = grid.Nz - @inbounds J = fr.flux_field[i, j, 1] - - # Restoring accessed as a tendency forcing (compatible with DatasetRestoring) - G = fr.restoring(i, j, Nz, grid, clock, fields) - - # Top BC convention: tendency contribution = -J / Δz, so to inject - # `G` in the top cell the flux is `-G * Δz`. - Δz = Δzᶜᶜᶜ(i, j, Nz, grid) - return J - G * Δz -end diff --git a/src/Oceans/ocean_simulation.jl b/src/Oceans/ocean_simulation.jl index 71a24f3a3..88627bb1b 100644 --- a/src/Oceans/ocean_simulation.jl +++ b/src/Oceans/ocean_simulation.jl @@ -19,18 +19,15 @@ using Statistics: mean ##### @inline ϕ²(i, j, k, grid, ϕ) = @inbounds ϕ[i, j, k]^2 -@inline spᶠᶜᶜ(i, j, k, grid, Φ, Uᴮ) = @inbounds sqrt(Φ.u[i, j, k]^2 + ℑxyᶠᶜᵃ(i, j, k, grid, ϕ², Φ.v) + Uᴮ^2) -@inline spᶜᶠᶜ(i, j, k, grid, Φ, Uᴮ) = @inbounds sqrt(Φ.v[i, j, k]^2 + ℑxyᶜᶠᵃ(i, j, k, grid, ϕ², Φ.u) + Uᴮ^2) +@inline spᶠᶜᶜ(i, j, k, grid, Φ) = @inbounds sqrt(Φ.u[i, j, k]^2 + ℑxyᶠᶜᵃ(i, j, k, grid, ϕ², Φ.v)) +@inline spᶜᶠᶜ(i, j, k, grid, Φ) = @inbounds sqrt(Φ.v[i, j, k]^2 + ℑxyᶜᶠᵃ(i, j, k, grid, ϕ², Φ.u)) -@inline u_quadratic_bottom_drag(i, j, grid, c, Φ, p) = @inbounds - p.μ * Φ.u[i, j, 1] * spᶠᶜᶜ(i, j, 1, grid, Φ, p.Uᴮ) -@inline v_quadratic_bottom_drag(i, j, grid, c, Φ, p) = @inbounds - p.μ * Φ.v[i, j, 1] * spᶜᶠᶜ(i, j, 1, grid, Φ, p.Uᴮ) +@inline u_quadratic_bottom_drag(i, j, grid, c, Φ, μ) = @inbounds - μ * Φ.u[i, j, 1] * spᶠᶜᶜ(i, j, 1, grid, Φ) +@inline v_quadratic_bottom_drag(i, j, grid, c, Φ, μ) = @inbounds - μ * Φ.v[i, j, 1] * spᶜᶠᶜ(i, j, 1, grid, Φ) # Keep a constant linear drag parameter independent on vertical level -@inline u_immersed_bottom_drag(i, j, k, grid, clock, Φ, p) = @inbounds - p.μ * Φ.u[i, j, k] * spᶠᶜᶜ(i, j, k, grid, Φ, p.Uᴮ) -@inline v_immersed_bottom_drag(i, j, k, grid, clock, Φ, p) = @inbounds - p.μ * Φ.v[i, j, k] * spᶜᶠᶜ(i, j, k, grid, Φ, p.Uᴮ) - -@inline build_top_tracer_bc(flux_field, ::Nothing) = FluxBoundaryCondition(flux_field) -@inline build_top_tracer_bc(flux_field, restoring) = FluxBoundaryCondition(FluxAndRestoring(flux_field, restoring); discrete_form=true) +@inline u_immersed_bottom_drag(i, j, k, grid, clock, Φ, μ) = @inbounds - μ * Φ.u[i, j, k] * spᶠᶜᶜ(i, j, k, grid, Φ) +@inline v_immersed_bottom_drag(i, j, k, grid, clock, Φ, μ) = @inbounds - μ * Φ.v[i, j, k] * spᶜᶠᶜ(i, j, k, grid, Φ) ##### ##### Defaults @@ -103,7 +100,6 @@ end """ ocean_simulation(grid; Δt = estimate_maximum_Δt(grid), - clock = Clock(grid), closure = default_ocean_closure(), tracers = (:T, :S), free_surface = default_free_surface(grid), @@ -111,9 +107,7 @@ end rotation_rate = default_planet_rotation_rate, gravitational_acceleration = default_gravitational_acceleration, bottom_drag_coefficient = Default(0.003), - drag_bulk_velocity = Default(0.1), forcing = NamedTuple(), - surface_restoring = NamedTuple(), biogeochemistry = nothing, timestepper = :SplitRungeKutta3, coriolis = Default(HydrostaticSphericalCoriolis(; rotation_rate)), @@ -132,6 +126,7 @@ consistent defaults for advection, closures, the equation of state, surface flux barotropic pressure–gradient forcing, boundary conditions, and optional biogeochemistry. It then wraps the model into an Oceananigans's `Simulation` with the specified timestepping options. + ## Behaviour and automatic configuration ### Coriolis @@ -166,7 +161,6 @@ defaults on a per-field basis. ## Keyword Arguments - `Δt`: Timestep used by the `Simulation`. Defaults to the maximum stable timestep estimated from the `grid`. -- `clock`: Clock object. Defaults to `Clock(grid)`. - `closure`: A turbulence or mixing closure. Defaults to `default_ocean_closure()`. - `tracers`: Tuple of tracer names. Defaults to `(:T, :S)`. - `free_surface`: Free–surface solver. Defaults to `default_free_surface(grid)`. @@ -174,9 +168,7 @@ defaults on a per-field basis. - `rotation_rate`: Planetary rotation rate used for Coriolis forcing. - `gravitational_acceleration`: Gravitational acceleration, passed to buoyancy. - `bottom_drag_coefficient`: Bottom drag coefficient. May be a `Default` wrapper. -- `drag_bulk_velocity`: a minimum velocity for the bottom drag. - `forcing`: Named tuple of additional forcing(s) for individual fields. -- `surface_restoring`: Named tuple of dataset restorings to apply as part of the tracer top boundary condition. - `biogeochemistry`: A biogeochemical model or `nothing`. - `timestepper`: Time-stepping scheme; options are `:SplitRungeKutta3` (default), or `:QuasiAdamsBashforth2`. - `coriolis`: Coriolis object or `Default(...)` wrapper. @@ -190,7 +182,6 @@ defaults on a per-field basis. """ function ocean_simulation(grid; Δt = estimate_maximum_Δt(grid), - clock = Clock(grid), closure = default_ocean_closure(), tracers = (:T, :S), free_surface = default_free_surface(grid), @@ -198,10 +189,7 @@ function ocean_simulation(grid; rotation_rate = default_planet_rotation_rate, gravitational_acceleration = default_gravitational_acceleration, bottom_drag_coefficient = Default(0.003), - drag_bulk_velocity = Default(0.05), - use_barotropic_potential = true, forcing = NamedTuple(), - surface_restoring = NamedTuple(), biogeochemistry = nothing, timestepper = :SplitRungeKutta3, coriolis = Default(HydrostaticSphericalCoriolis(; rotation_rate)), @@ -210,12 +198,11 @@ function ocean_simulation(grid; equation_of_state = TEOS10EquationOfState(; reference_density), boundary_conditions::NamedTuple = NamedTuple(), radiative_forcing = default_radiative_forcing(grid), - materialize_buoyancy_gradients = true, warn = true, verbose = false) FT = eltype(grid) - + if grid isa RectilinearGrid # turn off Coriolis unless user-supplied coriolis = default_or_override(coriolis, nothing) else @@ -229,8 +216,6 @@ function ocean_simulation(grid; if single_column_simulation # Let users put a bottom drag if they want bottom_drag_coefficient = default_or_override(bottom_drag_coefficient, zero(grid)) - drag_bulk_velocity = default_or_override(drag_bulk_velocity, zero(grid)) - drag_parameters = (μ = bottom_drag_coefficient, Uᴮ = drag_bulk_velocity) # Don't let users use advection in a single column model tracer_advection = nothing @@ -251,27 +236,21 @@ function ocean_simulation(grid; end bottom_drag_coefficient = default_or_override(bottom_drag_coefficient) - drag_bulk_velocity = default_or_override(drag_bulk_velocity) - bottom_drag_coefficient = convert(FT, bottom_drag_coefficient) - drag_bulk_velocity = convert(FT, drag_bulk_velocity) - drag_parameters = (μ = bottom_drag_coefficient, Uᴮ = drag_bulk_velocity) - u_immersed_drag = FluxBoundaryCondition(u_immersed_bottom_drag, discrete_form=true, parameters=drag_parameters) - v_immersed_drag = FluxBoundaryCondition(v_immersed_bottom_drag, discrete_form=true, parameters=drag_parameters) + u_immersed_drag = FluxBoundaryCondition(u_immersed_bottom_drag, discrete_form=true, parameters=bottom_drag_coefficient) + v_immersed_drag = FluxBoundaryCondition(v_immersed_bottom_drag, discrete_form=true, parameters=bottom_drag_coefficient) u_immersed_bc = ImmersedBoundaryCondition(bottom=u_immersed_drag) v_immersed_bc = ImmersedBoundaryCondition(bottom=v_immersed_drag) - if use_barotropic_potential - # Forcing for u, v - barotropic_potential = Field{Center, Center, Nothing}(grid) - u_forcing = BarotropicPotentialForcing(XDirection(), barotropic_potential) - v_forcing = BarotropicPotentialForcing(YDirection(), barotropic_potential) + # Forcing for u, v + barotropic_potential = Field{Center, Center, Nothing}(grid) + u_forcing = BarotropicPotentialForcing(XDirection(), barotropic_potential) + v_forcing = BarotropicPotentialForcing(YDirection(), barotropic_potential) - :u ∈ keys(forcing) && (u_forcing = (u_forcing, forcing[:u])) - :v ∈ keys(forcing) && (v_forcing = (v_forcing, forcing[:v])) - forcing = merge(forcing, (u=u_forcing, v=v_forcing)) - end + :u ∈ keys(forcing) && (u_forcing = (u_forcing, forcing[:u])) + :v ∈ keys(forcing) && (v_forcing = (v_forcing, forcing[:v])) + forcing = merge(forcing, (u=u_forcing, v=v_forcing)) end if !isnothing(radiative_forcing) @@ -283,20 +262,22 @@ function ocean_simulation(grid; forcing = merge(forcing, (; T=T_forcing)) end + bottom_drag_coefficient = convert(FT, bottom_drag_coefficient) + # Set up boundary conditions using Field - top_zonal_momentum_flux = τˣ = Field{Face, Center, Nothing}(grid) - top_meridional_momentum_flux = τʸ = Field{Center, Face, Nothing}(grid) + top_zonal_momentum_flux = τˣ = Field{Face, Center, Nothing}(grid) + top_meridional_momentum_flux = τʸ = Field{Center, Face, Nothing}(grid) top_ocean_heat_flux = Jᵀ = Field{Center, Center, Nothing}(grid) top_salt_flux = Jˢ = Field{Center, Center, Nothing}(grid) # Construct ocean boundary conditions including surface forcing and bottom drag u_top_bc = FluxBoundaryCondition(τˣ) v_top_bc = FluxBoundaryCondition(τʸ) - T_top_bc = build_top_tracer_bc(Jᵀ, get(surface_restoring, :T, nothing)) - S_top_bc = build_top_tracer_bc(Jˢ, get(surface_restoring, :S, nothing)) + T_top_bc = FluxBoundaryCondition(Jᵀ) + S_top_bc = FluxBoundaryCondition(Jˢ) - u_bot_bc = FluxBoundaryCondition(u_quadratic_bottom_drag, discrete_form=true, parameters=drag_parameters) - v_bot_bc = FluxBoundaryCondition(v_quadratic_bottom_drag, discrete_form=true, parameters=drag_parameters) + u_bot_bc = FluxBoundaryCondition(u_quadratic_bottom_drag, discrete_form=true, parameters=bottom_drag_coefficient) + v_bot_bc = FluxBoundaryCondition(v_quadratic_bottom_drag, discrete_form=true, parameters=bottom_drag_coefficient) default_boundary_conditions = (u = FieldBoundaryConditions(top=u_top_bc, bottom=u_bot_bc, immersed=u_immersed_bc), v = FieldBoundaryConditions(top=v_top_bc, bottom=v_bot_bc, immersed=v_immersed_bc), @@ -308,8 +289,7 @@ function ocean_simulation(grid; # conditions even when a user-bc is supplied). boundary_conditions = merge(default_boundary_conditions, boundary_conditions) buoyancy = SeawaterBuoyancy(; gravitational_acceleration, equation_of_state) - buoyancy = Oceananigans.BuoyancyFormulations.BuoyancyForce(grid, buoyancy; materialize_gradients=materialize_buoyancy_gradients) - + if tracer_advection isa NamedTuple tracer_advection = with_tracers(tracers, tracer_advection, default_tracer_advection()) else @@ -317,13 +297,12 @@ function ocean_simulation(grid; end if hasclosure(closure, CATKEVerticalDiffusivity) - # Use the same advection as for temperature - tke_advection = (; e=tracer_advection[1]) + # Turn off CATKE tracer advection + tke_advection = (; e=nothing) tracer_advection = merge(tracer_advection, tke_advection) end ocean_model = HydrostaticFreeSurfaceModel(grid; - clock, buoyancy, closure, biogeochemistry, diff --git a/test/tet_diagnostics_2.jl b/test/test_diagnostics_2.jl similarity index 97% rename from test/tet_diagnostics_2.jl rename to test/test_diagnostics_2.jl index d9912c756..59e6d1c68 100644 --- a/test/tet_diagnostics_2.jl +++ b/test/test_diagnostics_2.jl @@ -11,9 +11,9 @@ for arch in test_architectures @testset "InterfaceFluxOutputs on $A" begin grid = RectilinearGrid(arch; - size = (4, 4, 2), + size = (4, 5, 2), extent = (1, 1, 1), - topology = (Periodic, Periodic, Bounded)) + topology = (Periodic, Bounded, Bounded)) ocean = ocean_simulation(grid; momentum_advection = nothing, diff --git a/test/test_sea_ice_albedo.jl b/test/test_sea_ice_albedo.jl new file mode 100644 index 000000000..e80b0bbcb --- /dev/null +++ b/test/test_sea_ice_albedo.jl @@ -0,0 +1,141 @@ +include("runtests_setup.jl") + +using NumericalEarth: stateindex +using NumericalEarth.EarthSystemModels.InterfaceComputations: SeaIceAlbedo, + InterfaceProperties, + Radiation + +using Oceananigans.Units: Time + +@testset "SeaIceAlbedo" begin + for arch in test_architectures + A = typeof(arch) + + grid = RectilinearGrid(arch; + size = (4, 4, 1), + extent = (1, 1, 1), + topology = (Periodic, Periodic, Bounded)) + + hi = Field{Center, Center, Nothing}(grid) + hs = Field{Center, Center, Nothing}(grid) + Ts = Field{Center, Center, Nothing}(grid) + + @testset "Constructor [$A]" begin + α = SeaIceAlbedo(hi, hs, Ts) + @test α isa SeaIceAlbedo{Float64} + @test α.ice_albedo == 0.54 + @test α.snow_albedo == 0.83 + @test α.ocean_albedo == 0.06 + + α32 = SeaIceAlbedo(hi, hs, Ts; FT = Float32) + @test α32 isa SeaIceAlbedo{Float32} + end + + @testset "Nothing snow thickness [$A]" begin + α = SeaIceAlbedo(hi, nothing, Ts) + @test α.snow_thickness === nothing + end + + @testset "Cold thick ice, no snow [$A]" begin + # Thick ice, well below melting: should give cold bare-ice albedo + set!(hi, 2.0) # 2 m (>> minimum_ice_thickness = 0.5) + set!(hs, 0.0) + set!(Ts, -20.0) # well below melting + + α = SeaIceAlbedo(hi, hs, Ts) + time = Time(0.0) + loc = (Center, Center, Center) + + @allowscalar begin + val = stateindex(α, 1, 1, 1, grid, time, loc) + # With Ts = -20, fT = clamp((-20 - 0 + 1)/1, 0, 1) = 0 + # αi = 0.54 - 0 = 0.54, no thin-ice effect (fh=1), no snow (fs=0) + @test val ≈ 0.54 + end + end + + @testset "Cold thick ice, deep snow [$A]" begin + set!(hi, 2.0) + set!(hs, 0.1) # >> minimum_snow_depth = 0.02 + set!(Ts, -20.0) + + α = SeaIceAlbedo(hi, hs, Ts) + time = Time(0.0) + loc = (Center, Center, Center) + + @allowscalar begin + val = stateindex(α, 1, 1, 1, grid, time, loc) + # fs = clamp(0.1/0.02, 0, 1) = 1 → full snow albedo + @test val ≈ 0.83 + end + end + + @testset "Melting ice, no snow [$A]" begin + set!(hi, 2.0) + set!(hs, 0.0) + set!(Ts, 0.0) # at melting point + + α = SeaIceAlbedo(hi, hs, Ts) + time = Time(0.0) + loc = (Center, Center, Center) + + @allowscalar begin + val = stateindex(α, 1, 1, 1, grid, time, loc) + # fT = clamp((0 - 0 + 1)/1, 0, 1) = 1 + # αi = 0.54 - 0.075 = 0.465 + @test val ≈ 0.54 - 0.075 + end + end + + @testset "Thin ice transitions to ocean albedo [$A]" begin + set!(hi, 0.0) # no ice + set!(hs, 0.0) + set!(Ts, -10.0) + + α = SeaIceAlbedo(hi, hs, Ts) + time = Time(0.0) + loc = (Center, Center, Center) + + @allowscalar begin + val = stateindex(α, 1, 1, 1, grid, time, loc) + # fh = clamp(0/0.5, 0, 1) = 0 → pure ocean albedo + @test val ≈ 0.06 + end + end + + @testset "Partial snow cover [$A]" begin + set!(hi, 2.0) + set!(hs, 0.01) # half of minimum_snow_depth + set!(Ts, -20.0) + + α = SeaIceAlbedo(hi, hs, Ts) + time = Time(0.0) + loc = (Center, Center, Center) + + @allowscalar begin + val = stateindex(α, 1, 1, 1, grid, time, loc) + # fs = clamp(0.01/0.02, 0, 1) = 0.5 + # αi = 0.54, αs = 0.83 + @test val ≈ 0.5 * 0.83 + 0.5 * 0.54 + end + end + + @testset "Adapt.adapt_structure [$A]" begin + α = SeaIceAlbedo(hi, hs, Ts) + # Adapt should not error + adapted = Adapt.adapt_structure(Array, α) + @test adapted isa SeaIceAlbedo + @test adapted.ice_albedo == α.ice_albedo + end + + @testset "InterfaceProperties Adapt [$A]" begin + σ = 5.67e-8 + α = SeaIceAlbedo(hi, hs, Ts) + ϵ = 1.0 + radiation = (; σ, α, ϵ) + props = InterfaceProperties(radiation, nothing, nothing, nothing) + adapted = Adapt.adapt_structure(Array, props) + @test adapted.radiation.α isa SeaIceAlbedo + end + end +end From 467c0d5892769507008812478d281616eed5a923 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 20 Apr 2026 14:48:03 +0200 Subject: [PATCH 7/9] fix tests --- test/test_sea_ice_albedo.jl | 18 ------------------ 1 file changed, 18 deletions(-) diff --git a/test/test_sea_ice_albedo.jl b/test/test_sea_ice_albedo.jl index e80b0bbcb..6969cc556 100644 --- a/test/test_sea_ice_albedo.jl +++ b/test/test_sea_ice_albedo.jl @@ -119,23 +119,5 @@ using Oceananigans.Units: Time @test val ≈ 0.5 * 0.83 + 0.5 * 0.54 end end - - @testset "Adapt.adapt_structure [$A]" begin - α = SeaIceAlbedo(hi, hs, Ts) - # Adapt should not error - adapted = Adapt.adapt_structure(Array, α) - @test adapted isa SeaIceAlbedo - @test adapted.ice_albedo == α.ice_albedo - end - - @testset "InterfaceProperties Adapt [$A]" begin - σ = 5.67e-8 - α = SeaIceAlbedo(hi, hs, Ts) - ϵ = 1.0 - radiation = (; σ, α, ϵ) - props = InterfaceProperties(radiation, nothing, nothing, nothing) - adapted = Adapt.adapt_structure(Array, props) - @test adapted.radiation.α isa SeaIceAlbedo - end end end From 3ec70d374d3922405b6a1bd656bd77effc3112e8 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 30 Apr 2026 09:07:18 +0200 Subject: [PATCH 8/9] Remove comments on albedo struct usage Removed comments about albedo struct evaluation. --- .../InterfaceComputations/atmosphere_sea_ice_fluxes.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl index f18078a2d..a5228f879 100644 --- a/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl @@ -92,8 +92,6 @@ end end # Evaluate state-dependent radiation properties at this grid point. - # The albedo may be a struct (e.g., CCSM3SeaIceAlbedo) that reads model fields; - # we evaluate it here so the iteration uses a scalar. time = Time(clock.time) σ = interface_properties.radiation.σ α = stateindex(interface_properties.radiation.α, i, j, kᴺ, grid, time, CCC) From d128172376342446edd847dd7b88df84c2be2ee6 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 30 Apr 2026 09:08:23 +0200 Subject: [PATCH 9/9] Clarify elevation of atmosphere variables MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Added comment to clarify the purpose of the variable zᵃᵗ. --- .../InterfaceComputations/atmosphere_sea_ice_fluxes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl index a5228f879..5bd715e63 100644 --- a/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl @@ -104,7 +104,7 @@ end # Build thermodynamic and dynamic states in the atmosphere and interface. ℂᵃᵗ = atmosphere_properties.thermodynamics_parameters - zᵃᵗ = atmosphere_properties.surface_layer_height + zᵃᵗ = atmosphere_properties.surface_layer_height # elevation of atmos variables relative to interface local_atmosphere_state = (z = zᵃᵗ, u = uᵃᵗ,