From 1c09a4d962b77cb00c13d54a1713bf74f23604c4 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 2 Mar 2026 10:54:54 +0100 Subject: [PATCH 01/63] pass grid to sea ice SplitExplicitSolver --- src/SeaIces/sea_ice_simulation.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/SeaIces/sea_ice_simulation.jl b/src/SeaIces/sea_ice_simulation.jl index 90c29e761..0afda344d 100644 --- a/src/SeaIces/sea_ice_simulation.jl +++ b/src/SeaIces/sea_ice_simulation.jl @@ -76,7 +76,7 @@ function sea_ice_dynamics(grid, ocean=nothing; rheology = ElastoViscoPlasticRheology(), coriolis = HydrostaticSphericalCoriolis(; rotation_rate=default_rotation_rate), free_drift = nothing, - solver = SplitExplicitSolver(120)) + solver = SplitExplicitSolver(grid; substeps=120)) SSU, SSV = ocean_surface_velocities(ocean) sea_ice_ocean_drag_coefficient = convert(eltype(grid), sea_ice_ocean_drag_coefficient) @@ -108,7 +108,7 @@ sea_ice_concentration(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.ice_co heat_capacity(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity reference_density(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density -function net_fluxes(sea_ice::Simulation{<:SeaIceModel}) +function net_fluxes(sea_ice::Simulation{<:SeaIceModel}) net_momentum_fluxes = if isnothing(sea_ice.model.dynamics) u = Field{Face, Center, Nothing}(sea_ice.model.grid) v = Field{Center, Face, Nothing}(sea_ice.model.grid) @@ -131,17 +131,17 @@ function default_ai_temperature(sea_ice::Simulation{<:SeaIceModel}) end # Constructor that accepts the sea-ice model -function ThreeEquationHeatFlux(sea_ice::Simulation{<:SeaIceModel}, FT::DataType = Oceananigans.defaults.FloatType; +function ThreeEquationHeatFlux(sea_ice::Simulation{<:SeaIceModel}, FT::DataType = Oceananigans.defaults.FloatType; heat_transfer_coefficient = 0.0095, salt_transfer_coefficient = heat_transfer_coefficient / 35, friction_velocity = convert(FT, 0.002)) conductive_flux = sea_ice.model.ice_thermodynamics.internal_heat_flux.parameters.flux ice_temperature = sea_ice.model.ice_thermodynamics.top_surface_temperature - + return ThreeEquationHeatFlux(conductive_flux, ice_temperature, convert(FT, heat_transfer_coefficient), convert(FT, salt_transfer_coefficient), - friction_velocity) + friction_velocity) end From 2cc2ca91b399f7c769c14907df48270ad633f493 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 2 Mar 2026 10:55:08 +0100 Subject: [PATCH 02/63] bump ClimaSeaIce compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 38b81b4f9..a9558ea54 100644 --- a/Project.toml +++ b/Project.toml @@ -58,7 +58,7 @@ CDSAPI = "2" CFTime = "0.1, 0.2" CUDA = "5.9.5" CUDA_Compiler_jll = "v0.3.1" -ClimaSeaIce = "0.4.2, 0.5" +ClimaSeaIce = "0.4.4, 0.5" CondaPkg = "0.2.33" CopernicusMarine = "0.1.1" CubicSplines = "0.2" From ff5c5a7b1e7601d84ab263307e8de69b6b3b3607 Mon Sep 17 00:00:00 2001 From: Taimoor Sohail Date: Tue, 3 Mar 2026 03:43:44 +1100 Subject: [PATCH 03/63] Added Meridional Heat Transport calculation; subject to changes in the heat flux code --- src/Diagnostics/Diagnostics.jl | 3 ++- src/Diagnostics/MeridionalHeatTransport.jl | 21 +++++++++++++++++++++ 2 files changed, 23 insertions(+), 1 deletion(-) create mode 100644 src/Diagnostics/MeridionalHeatTransport.jl diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl index 76e21b8c0..bde515726 100644 --- a/src/Diagnostics/Diagnostics.jl +++ b/src/Diagnostics/Diagnostics.jl @@ -1,6 +1,6 @@ module Diagnostics -export MixedLayerDepthField, MixedLayerDepthOperand +export MixedLayerDepthField, MixedLayerDepthOperand, Meridional_Heat_Transport using Oceananigans using Oceananigans.Architectures: architecture @@ -14,5 +14,6 @@ using KernelAbstractions: @index, @kernel import Oceananigans.Fields: compute! include("mixed_layer_depth.jl") +include("MeridionalHeatTransport.jl") end # module diff --git a/src/Diagnostics/MeridionalHeatTransport.jl b/src/Diagnostics/MeridionalHeatTransport.jl new file mode 100644 index 000000000..7a9e92163 --- /dev/null +++ b/src/Diagnostics/MeridionalHeatTransport.jl @@ -0,0 +1,21 @@ +""" + Meridional_Heat_Transport(coupled_model) + +Return the cumulative meridional integral of the zonal-depth integrated +temperature minus the zonally integrated interface heat-flux tendency. +""" +function Meridional_Heat_Transport(coupled_model) + ocean = coupled_model.ocean + + T_int = Integral(ocean.model.tracers.T, dims=(1, 3)) + + flux_outputs = InterfaceFluxOutputs(coupled_model; + isolate_sea_ice=false, + units=:physical, + reference_salinity=35) + + heat_flux = haskey(flux_outputs, "heat_flux") ? flux_outputs["heat_flux"] : flux_outputs[:heat_flux] + flux_int = Integral(heat_flux * ocean.model.clock.Δt, dims=(1)) + + return CumulativeIntegral(T_int - flux_int, dims=(2)) +end From ed0fccd73592dfd64f6895d6f6d04ea548563200 Mon Sep 17 00:00:00 2001 From: Taimoor Sohail Date: Tue, 3 Mar 2026 04:41:48 +1100 Subject: [PATCH 04/63] Added pickup utility to MHT calculation --- src/Diagnostics/Diagnostics.jl | 2 +- src/Diagnostics/MeridionalHeatTransport.jl | 107 ++++++++++++++++++-- src/EarthSystemModels/earth_system_model.jl | 8 +- 3 files changed, 107 insertions(+), 10 deletions(-) diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl index bde515726..ebde20608 100644 --- a/src/Diagnostics/Diagnostics.jl +++ b/src/Diagnostics/Diagnostics.jl @@ -1,6 +1,6 @@ module Diagnostics -export MixedLayerDepthField, MixedLayerDepthOperand, Meridional_Heat_Transport +export MixedLayerDepthField, MixedLayerDepthOperand, Meridional_Heat_Transport, reset_meridional_heat_transport_state! using Oceananigans using Oceananigans.Architectures: architecture diff --git a/src/Diagnostics/MeridionalHeatTransport.jl b/src/Diagnostics/MeridionalHeatTransport.jl index 7a9e92163..823749787 100644 --- a/src/Diagnostics/MeridionalHeatTransport.jl +++ b/src/Diagnostics/MeridionalHeatTransport.jl @@ -1,21 +1,112 @@ +import ..EarthSystemModels: EarthSystemModel, checkpoint_auxiliary_state, restore_auxiliary_state! + """ Meridional_Heat_Transport(coupled_model) -Return the cumulative meridional integral of the zonal-depth integrated -temperature minus the zonally integrated interface heat-flux tendency. +Return meridional heat transport diagnosed from the OHC anomaly budget: + +`MHT = CumulativeIntegral((OHC - OHC₀) - ∫ₓ(∫Q dt), dims=(2))`, + +where `OHC = ρₒ cₚ ∫ₓ∫z T`, `OHC₀` is the first sampled OHC profile, and +`Q` is `"heat_flux"` from `InterfaceFluxOutputs`. """ -function Meridional_Heat_Transport(coupled_model) - ocean = coupled_model.ocean +mutable struct MeridionalHeatTransportState + initial_ohc::Any + cumulative_heat_flux::Any + last_time::Float64 +end + +const meridional_heat_transport_states = IdDict{Any, MeridionalHeatTransportState}() - T_int = Integral(ocean.model.tracers.T, dims=(1, 3)) +allocate_storage_like(field) = Field(instantiated_location(field), field.grid; indices=indices(field)) + +""" + reset_meridional_heat_transport_state!(coupled_model) + +Clear cached MHT state (`OHC₀`, cumulative heat flux, last time) for `coupled_model`. +Call this before restarting from a checkpoint or when reusing a model object. +""" +function reset_meridional_heat_transport_state!(coupled_model) + pop!(meridional_heat_transport_states, coupled_model, nothing) + return nothing +end + +function initialize_mht_state!(coupled_model, heat_flux_field, ohc_field, time) + initial_ohc = allocate_storage_like(ohc_field) + set!(initial_ohc, ohc_field) + + cumulative_heat_flux = allocate_storage_like(heat_flux_field) + set!(cumulative_heat_flux, 0) + + state = MeridionalHeatTransportState(initial_ohc, cumulative_heat_flux, time) + meridional_heat_transport_states[coupled_model] = state + return state +end +function current_heat_flux_field(coupled_model) flux_outputs = InterfaceFluxOutputs(coupled_model; isolate_sea_ice=false, units=:physical, reference_salinity=35) - heat_flux = haskey(flux_outputs, "heat_flux") ? flux_outputs["heat_flux"] : flux_outputs[:heat_flux] - flux_int = Integral(heat_flux * ocean.model.clock.Δt, dims=(1)) + heat_flux_raw = haskey(flux_outputs, "heat_flux") ? flux_outputs["heat_flux"] : flux_outputs[:heat_flux] + heat_flux_field = Field(heat_flux_raw) + compute!(heat_flux_field) + return heat_flux_field +end + +function Meridional_Heat_Transport(coupled_model; ρₒ=1035.0, cₚ=3991.86795711963) + ocean = coupled_model.ocean + heat_flux_field = current_heat_flux_field(coupled_model) + ohc_field = Field(ρₒ * cₚ * Integral(ocean.model.tracers.T, dims=(1, 3))) + compute!(ohc_field) + + model_time = Float64(ocean.model.clock.time) + state = get(meridional_heat_transport_states, coupled_model, nothing) + state === nothing && (state = initialize_mht_state!(coupled_model, heat_flux_field, ohc_field, model_time)) + + Δt = max(0.0, model_time - state.last_time) + if Δt == 0.0 && ocean.model.clock.iteration > 0 + Δt = Float64(ocean.model.clock.Δt) + end + state.last_time = model_time + + set!(state.cumulative_heat_flux, state.cumulative_heat_flux + Δt * heat_flux_field) + + Δohc = ohc_field - state.initial_ohc + flux_int = Integral(state.cumulative_heat_flux, dims=(1)) + + return CumulativeIntegral(Δohc - flux_int, dims=(2)) +end + +function checkpoint_auxiliary_state(coupled_model::EarthSystemModel) + state = get(meridional_heat_transport_states, coupled_model, nothing) + state === nothing && return nothing + + return ( + meridional_heat_transport = ( + initial_ohc = Array(interior(state.initial_ohc)), + cumulative_heat_flux = Array(interior(state.cumulative_heat_flux)), + last_time = state.last_time + ), + ) +end + +function restore_auxiliary_state!(coupled_model::EarthSystemModel, auxiliary_state) + auxiliary_state === nothing && return nothing + hasproperty(auxiliary_state, :meridional_heat_transport) || return nothing + + mht_state = auxiliary_state.meridional_heat_transport + mht_state === nothing && return nothing + + heat_flux_field = current_heat_flux_field(coupled_model) + ohc_template = Field(Integral(coupled_model.ocean.model.tracers.T, dims=(1, 3))) + compute!(ohc_template) - return CumulativeIntegral(T_int - flux_int, dims=(2)) + reset_meridional_heat_transport_state!(coupled_model) + state = initialize_mht_state!(coupled_model, heat_flux_field, ohc_template, Float64(mht_state.last_time)) + set!(state.initial_ohc, mht_state.initial_ohc) + set!(state.cumulative_heat_flux, mht_state.cumulative_heat_flux) + state.last_time = Float64(mht_state.last_time) + return nothing end diff --git a/src/EarthSystemModels/earth_system_model.jl b/src/EarthSystemModels/earth_system_model.jl index 931070638..f00f0651c 100644 --- a/src/EarthSystemModels/earth_system_model.jl +++ b/src/EarthSystemModels/earth_system_model.jl @@ -322,12 +322,16 @@ above_freezing_ocean_temperature!(ocean, grid, ::Nothing) = nothing ##### Checkpointing ##### +checkpoint_auxiliary_state(model) = nothing +restore_auxiliary_state!(model, state) = nothing + function prognostic_state(osm::EarthSystemModel) return (clock = prognostic_state(osm.clock), ocean = prognostic_state(osm.ocean), atmosphere = prognostic_state(osm.atmosphere), sea_ice = prognostic_state(osm.sea_ice), - interfaces = prognostic_state(osm.interfaces)) + interfaces = prognostic_state(osm.interfaces), + auxiliary = prognostic_state(checkpoint_auxiliary_state(osm))) end function restore_prognostic_state!(osm::EarthSystemModel, state) @@ -336,6 +340,8 @@ function restore_prognostic_state!(osm::EarthSystemModel, state) restore_prognostic_state!(osm.atmosphere, state.atmosphere) restore_prognostic_state!(osm.sea_ice, state.sea_ice) restore_prognostic_state!(osm.interfaces, state.interfaces) + auxiliary_state = hasproperty(state, :auxiliary) ? state.auxiliary : nothing + restore_auxiliary_state!(osm, auxiliary_state) # Note: we do NOT call update_state! here because: # 1. The checkpoint was saved AFTER update_state! was called at the end of that time step # 2. Calling update_state! would recompute interface fluxes and overwrite restored state From 9fa7f18bda98078add4b1a145f58ae25f71756d6 Mon Sep 17 00:00:00 2001 From: Taimoor Sohail Date: Tue, 3 Mar 2026 19:47:06 +1100 Subject: [PATCH 05/63] Added ability to use the temperature_tendency for OHU --- src/Diagnostics/MeridionalHeatTransport.jl | 101 ++++++++++++++------- 1 file changed, 70 insertions(+), 31 deletions(-) diff --git a/src/Diagnostics/MeridionalHeatTransport.jl b/src/Diagnostics/MeridionalHeatTransport.jl index 823749787..a1da9654f 100644 --- a/src/Diagnostics/MeridionalHeatTransport.jl +++ b/src/Diagnostics/MeridionalHeatTransport.jl @@ -1,19 +1,28 @@ -import ..EarthSystemModels: EarthSystemModel, checkpoint_auxiliary_state, restore_auxiliary_state! +import ..EarthSystemModels: EarthSystemModel, checkpoint_auxiliary_state, restore_auxiliary_state!, + reference_density, heat_capacity """ - Meridional_Heat_Transport(coupled_model) + Meridional_Heat_Transport(coupled_model; method=:ohc_tendency, T_ref=0.0) -Return meridional heat transport diagnosed from the OHC anomaly budget: +Compute meridional heat transport with one of two methods: -`MHT = CumulativeIntegral((OHC - OHC₀) - ∫ₓ(∫Q dt), dims=(2))`, +1. `method = :ohc_tendency` (default): -where `OHC = ρₒ cₚ ∫ₓ∫z T`, `OHC₀` is the first sampled OHC profile, and -`Q` is `"heat_flux"` from `InterfaceFluxOutputs`. + `MHT = CumulativeIntegral((∫ ρₒ cₚ ∫ₓ∫z ∂ₜT dt) - ∫ₓ(∫Q dt), dims=(2))`, + where `Q` is `"heat_flux"` from `InterfaceFluxOutputs`. + +2. `method = :vt_instantaneous`: + + `MHT = ρₒ cₚ ∫ₓ∫z v * (T - T_ref)` + +`ρₒ` and `cₚ` are inferred from `coupled_model.ocean` via +`reference_density` and `heat_capacity`. """ mutable struct MeridionalHeatTransportState - initial_ohc::Any + cumulative_ohc_tendency::Any cumulative_heat_flux::Any last_time::Float64 + last_iteration::Int end const meridional_heat_transport_states = IdDict{Any, MeridionalHeatTransportState}() @@ -23,7 +32,8 @@ allocate_storage_like(field) = Field(instantiated_location(field), field.grid; i """ reset_meridional_heat_transport_state!(coupled_model) -Clear cached MHT state (`OHC₀`, cumulative heat flux, last time) for `coupled_model`. +Clear cached MHT state (cumulative OHC tendency, cumulative heat flux, and clock metadata) +for `coupled_model`. Call this before restarting from a checkpoint or when reusing a model object. """ function reset_meridional_heat_transport_state!(coupled_model) @@ -31,14 +41,14 @@ function reset_meridional_heat_transport_state!(coupled_model) return nothing end -function initialize_mht_state!(coupled_model, heat_flux_field, ohc_field, time) - initial_ohc = allocate_storage_like(ohc_field) - set!(initial_ohc, ohc_field) +function initialize_mht_state!(coupled_model, heat_flux_field, ohc_tendency_field, time, iteration) + cumulative_ohc_tendency = allocate_storage_like(ohc_tendency_field) + set!(cumulative_ohc_tendency, 0) cumulative_heat_flux = allocate_storage_like(heat_flux_field) set!(cumulative_heat_flux, 0) - state = MeridionalHeatTransportState(initial_ohc, cumulative_heat_flux, time) + state = MeridionalHeatTransportState(cumulative_ohc_tendency, cumulative_heat_flux, time, iteration) meridional_heat_transport_states[coupled_model] = state return state end @@ -55,28 +65,51 @@ function current_heat_flux_field(coupled_model) return heat_flux_field end -function Meridional_Heat_Transport(coupled_model; ρₒ=1035.0, cₚ=3991.86795711963) +@inline mht_constants(coupled_model) = reference_density(coupled_model.ocean), heat_capacity(coupled_model.ocean) + +function ohc_tendency_mht(coupled_model, ρₒ, cₚ) ocean = coupled_model.ocean heat_flux_field = current_heat_flux_field(coupled_model) - ohc_field = Field(ρₒ * cₚ * Integral(ocean.model.tracers.T, dims=(1, 3))) - compute!(ohc_field) + ohc_tendency_field = Field(ρₒ * cₚ * Integral(ocean.model.timestepper.Gⁿ.T, dims=(1, 3))) + compute!(ohc_tendency_field) model_time = Float64(ocean.model.clock.time) + model_iteration = Int(ocean.model.clock.iteration) state = get(meridional_heat_transport_states, coupled_model, nothing) - state === nothing && (state = initialize_mht_state!(coupled_model, heat_flux_field, ohc_field, model_time)) - - Δt = max(0.0, model_time - state.last_time) - if Δt == 0.0 && ocean.model.clock.iteration > 0 - Δt = Float64(ocean.model.clock.Δt) + state === nothing && (state = initialize_mht_state!(coupled_model, heat_flux_field, ohc_tendency_field, model_time, model_iteration)) + + if model_iteration != state.last_iteration + Δt = max(0.0, model_time - state.last_time) + if Δt == 0.0 && model_iteration > 0 + Δt = Float64(ocean.model.clock.Δt) + end + + set!(state.cumulative_ohc_tendency, state.cumulative_ohc_tendency + Δt * ohc_tendency_field) + set!(state.cumulative_heat_flux, state.cumulative_heat_flux + Δt * heat_flux_field) + state.last_time = model_time + state.last_iteration = model_iteration end - state.last_time = model_time - - set!(state.cumulative_heat_flux, state.cumulative_heat_flux + Δt * heat_flux_field) - Δohc = ohc_field - state.initial_ohc flux_int = Integral(state.cumulative_heat_flux, dims=(1)) - return CumulativeIntegral(Δohc - flux_int, dims=(2)) + return CumulativeIntegral(state.cumulative_ohc_tendency - flux_int, dims=(2)) +end + +function instantaneous_vt_mht(coupled_model, ρₒ, cₚ, T_ref) + ocean_model = coupled_model.ocean.model + return ρₒ * cₚ * Integral(ocean_model.velocities.v * (ocean_model.tracers.T - T_ref), dims=(1, 3)) +end + +function Meridional_Heat_Transport(coupled_model; method=:ohc_tendency, T_ref=0.0) + ρₒ, cₚ = mht_constants(coupled_model) + + if method === :ohc_tendency + return ohc_tendency_mht(coupled_model, ρₒ, cₚ) + elseif method === :vt_instantaneous + return instantaneous_vt_mht(coupled_model, ρₒ, cₚ, T_ref) + else + throw(ArgumentError("Unknown MHT method=$(repr(method)). Supported methods are :ohc_tendency and :vt_instantaneous.")) + end end function checkpoint_auxiliary_state(coupled_model::EarthSystemModel) @@ -85,9 +118,10 @@ function checkpoint_auxiliary_state(coupled_model::EarthSystemModel) return ( meridional_heat_transport = ( - initial_ohc = Array(interior(state.initial_ohc)), + cumulative_ohc_tendency = Array(interior(state.cumulative_ohc_tendency)), cumulative_heat_flux = Array(interior(state.cumulative_heat_flux)), - last_time = state.last_time + last_time = state.last_time, + last_iteration = state.last_iteration ), ) end @@ -100,13 +134,18 @@ function restore_auxiliary_state!(coupled_model::EarthSystemModel, auxiliary_sta mht_state === nothing && return nothing heat_flux_field = current_heat_flux_field(coupled_model) - ohc_template = Field(Integral(coupled_model.ocean.model.tracers.T, dims=(1, 3))) - compute!(ohc_template) + ohc_tendency_template = Field(Integral(coupled_model.ocean.model.timestepper.Gⁿ.T, dims=(1, 3))) + compute!(ohc_tendency_template) reset_meridional_heat_transport_state!(coupled_model) - state = initialize_mht_state!(coupled_model, heat_flux_field, ohc_template, Float64(mht_state.last_time)) - set!(state.initial_ohc, mht_state.initial_ohc) + state = initialize_mht_state!(coupled_model, + heat_flux_field, + ohc_tendency_template, + Float64(mht_state.last_time), + Int(mht_state.last_iteration)) + set!(state.cumulative_ohc_tendency, mht_state.cumulative_ohc_tendency) set!(state.cumulative_heat_flux, mht_state.cumulative_heat_flux) state.last_time = Float64(mht_state.last_time) + state.last_iteration = Int(mht_state.last_iteration) return nothing end From 0bb0c13bfb7d3c7b9fcbb81d8d61ff84833ca08e Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 3 Mar 2026 11:33:46 +0100 Subject: [PATCH 06/63] rename --- src/Diagnostics/Diagnostics.jl | 2 +- ...{MeridionalHeatTransport.jl => meridional_heat_transport.jl} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/Diagnostics/{MeridionalHeatTransport.jl => meridional_heat_transport.jl} (100%) diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl index ebde20608..3e28356dd 100644 --- a/src/Diagnostics/Diagnostics.jl +++ b/src/Diagnostics/Diagnostics.jl @@ -14,6 +14,6 @@ using KernelAbstractions: @index, @kernel import Oceananigans.Fields: compute! include("mixed_layer_depth.jl") -include("MeridionalHeatTransport.jl") +include("meridional_heat_transport.jl") end # module diff --git a/src/Diagnostics/MeridionalHeatTransport.jl b/src/Diagnostics/meridional_heat_transport.jl similarity index 100% rename from src/Diagnostics/MeridionalHeatTransport.jl rename to src/Diagnostics/meridional_heat_transport.jl From 7108fa57908d150727996d17e6d14868d91b2280 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 3 Mar 2026 11:36:06 +0100 Subject: [PATCH 07/63] docstring on constructor + drop mixed Camel_Case --- src/Diagnostics/meridional_heat_transport.jl | 42 ++++++++++---------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index a1da9654f..7d8d21dad 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -1,8 +1,17 @@ import ..EarthSystemModels: EarthSystemModel, checkpoint_auxiliary_state, restore_auxiliary_state!, - reference_density, heat_capacity + reference_density, heat_capacity + +mutable struct MeridionalHeatTransportState + cumulative_ohc_tendency::Any + cumulative_heat_flux::Any + last_time::Float64 + last_iteration::Int +end + +const meridional_heat_transport_states = IdDict{Any, MeridionalHeatTransportState}() """ - Meridional_Heat_Transport(coupled_model; method=:ohc_tendency, T_ref=0.0) + meridional_heat_transport(esm::EarthSystemModel; method=:ohc_tendency, T_ref=0.0) Compute meridional heat transport with one of two methods: @@ -18,14 +27,17 @@ Compute meridional heat transport with one of two methods: `ρₒ` and `cₚ` are inferred from `coupled_model.ocean` via `reference_density` and `heat_capacity`. """ -mutable struct MeridionalHeatTransportState - cumulative_ohc_tendency::Any - cumulative_heat_flux::Any - last_time::Float64 - last_iteration::Int -end +function meridional_heat_transport(esm::EarthSystemModel; method=:ohc_tendency, T_ref=0.0) + ρₒ, cₚ = mht_constants(esm) -const meridional_heat_transport_states = IdDict{Any, MeridionalHeatTransportState}() + if method === :ohc_tendency + return ohc_tendency_mht(esm, ρₒ, cₚ) + elseif method === :vt_instantaneous + return instantaneous_vt_mht(esm, ρₒ, cₚ, T_ref) + else + throw(ArgumentError("Unknown MHT method=$(repr(method)). Supported methods are :ohc_tendency and :vt_instantaneous.")) + end +end allocate_storage_like(field) = Field(instantiated_location(field), field.grid; indices=indices(field)) @@ -100,18 +112,6 @@ function instantaneous_vt_mht(coupled_model, ρₒ, cₚ, T_ref) return ρₒ * cₚ * Integral(ocean_model.velocities.v * (ocean_model.tracers.T - T_ref), dims=(1, 3)) end -function Meridional_Heat_Transport(coupled_model; method=:ohc_tendency, T_ref=0.0) - ρₒ, cₚ = mht_constants(coupled_model) - - if method === :ohc_tendency - return ohc_tendency_mht(coupled_model, ρₒ, cₚ) - elseif method === :vt_instantaneous - return instantaneous_vt_mht(coupled_model, ρₒ, cₚ, T_ref) - else - throw(ArgumentError("Unknown MHT method=$(repr(method)). Supported methods are :ohc_tendency and :vt_instantaneous.")) - end -end - function checkpoint_auxiliary_state(coupled_model::EarthSystemModel) state = get(meridional_heat_transport_states, coupled_model, nothing) state === nothing && return nothing From 8976ca4029769bab7f1ef1f1941eae4cc64372ce Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 3 Mar 2026 11:46:27 +0100 Subject: [PATCH 08/63] Apply suggestions from code review --- src/Diagnostics/meridional_heat_transport.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 7d8d21dad..a07972aed 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -27,7 +27,7 @@ Compute meridional heat transport with one of two methods: `ρₒ` and `cₚ` are inferred from `coupled_model.ocean` via `reference_density` and `heat_capacity`. """ -function meridional_heat_transport(esm::EarthSystemModel; method=:ohc_tendency, T_ref=0.0) +function meridional_heat_transport(esm::EarthSystemModel; method=:ohc_tendency, reference_temperature=0.0) ρₒ, cₚ = mht_constants(esm) if method === :ohc_tendency @@ -102,7 +102,7 @@ function ohc_tendency_mht(coupled_model, ρₒ, cₚ) state.last_iteration = model_iteration end - flux_int = Integral(state.cumulative_heat_flux, dims=(1)) + flux_int = Integral(state.cumulative_heat_flux, dims=1) return CumulativeIntegral(state.cumulative_ohc_tendency - flux_int, dims=(2)) end From dadec35383fa9906cdae0326a9ed11c32ebdd5cc Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 3 Mar 2026 11:57:26 +0100 Subject: [PATCH 09/63] clearer hopefully --- src/Diagnostics/meridional_heat_transport.jl | 87 ++++++++++---------- 1 file changed, 45 insertions(+), 42 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index a07972aed..a9dda6b66 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -11,29 +11,27 @@ end const meridional_heat_transport_states = IdDict{Any, MeridionalHeatTransportState}() """ - meridional_heat_transport(esm::EarthSystemModel; method=:ohc_tendency, T_ref=0.0) + meridional_heat_transport(esm::EarthSystemModel; method=:ohc_tendency, reference_temperature=0.0) Compute meridional heat transport with one of two methods: 1. `method = :ohc_tendency` (default): - `MHT = CumulativeIntegral((∫ ρₒ cₚ ∫ₓ∫z ∂ₜT dt) - ∫ₓ(∫Q dt), dims=(2))`, + `MHT = CumulativeIntegral((∫ ρᵒᶜ cᵒᶜ ∫ₓ∫z ∂ₜT dt) - ∫ₓ(∫Q dt), dims=(2))`, where `Q` is `"heat_flux"` from `InterfaceFluxOutputs`. 2. `method = :vt_instantaneous`: - `MHT = ρₒ cₚ ∫ₓ∫z v * (T - T_ref)` + `MHT = ρᵒᶜ cᵒᶜ ∫ₓ∫z v * (T - T_ref)` -`ρₒ` and `cₚ` are inferred from `coupled_model.ocean` via -`reference_density` and `heat_capacity`. +`ρᵒᶜ` and `cᵒᶜ` are the ocean reference density and heat capacity respectively +and they are inferred from `esm.ocean`. """ function meridional_heat_transport(esm::EarthSystemModel; method=:ohc_tendency, reference_temperature=0.0) - ρₒ, cₚ = mht_constants(esm) - if method === :ohc_tendency - return ohc_tendency_mht(esm, ρₒ, cₚ) + return ohc_tendency_mht(esm) elseif method === :vt_instantaneous - return instantaneous_vt_mht(esm, ρₒ, cₚ, T_ref) + return instantaneous_vt_mht(esm, reference_temperature) else throw(ArgumentError("Unknown MHT method=$(repr(method)). Supported methods are :ohc_tendency and :vt_instantaneous.")) end @@ -42,18 +40,18 @@ end allocate_storage_like(field) = Field(instantiated_location(field), field.grid; indices=indices(field)) """ - reset_meridional_heat_transport_state!(coupled_model) + reset_meridional_heat_transport_state!(esm) -Clear cached MHT state (cumulative OHC tendency, cumulative heat flux, and clock metadata) -for `coupled_model`. -Call this before restarting from a checkpoint or when reusing a model object. +Clear cached MHT state (cumulative OHC tendency, cumulative heat flux, +and clock metadata) for a coupled `esm`. This method should be called before + restarting from a checkpoint or when reusing a model object. """ -function reset_meridional_heat_transport_state!(coupled_model) - pop!(meridional_heat_transport_states, coupled_model, nothing) +function reset_meridional_heat_transport_state!(esm) + pop!(meridional_heat_transport_states, esm, nothing) return nothing end -function initialize_mht_state!(coupled_model, heat_flux_field, ohc_tendency_field, time, iteration) +function initialize_mht_state!(esm, heat_flux_field, ohc_tendency_field, time, iteration) cumulative_ohc_tendency = allocate_storage_like(ohc_tendency_field) set!(cumulative_ohc_tendency, 0) @@ -61,12 +59,12 @@ function initialize_mht_state!(coupled_model, heat_flux_field, ohc_tendency_fiel set!(cumulative_heat_flux, 0) state = MeridionalHeatTransportState(cumulative_ohc_tendency, cumulative_heat_flux, time, iteration) - meridional_heat_transport_states[coupled_model] = state + meridional_heat_transport_states[esm] = state return state end -function current_heat_flux_field(coupled_model) - flux_outputs = InterfaceFluxOutputs(coupled_model; +function current_heat_flux_field(esm) + flux_outputs = InterfaceFluxOutputs(esm; isolate_sea_ice=false, units=:physical, reference_salinity=35) @@ -77,23 +75,25 @@ function current_heat_flux_field(coupled_model) return heat_flux_field end -@inline mht_constants(coupled_model) = reference_density(coupled_model.ocean), heat_capacity(coupled_model.ocean) +function ohc_tendency_mht(esm) + ρᵒᶜ = reference_density(esm.ocean) + cᵒᶜ = heat_capacity(esm.ocean) + heat_flux = current_heat_flux_field(esm) + T_tendency = esm.ocean.model.timestepper.Gⁿ.T -function ohc_tendency_mht(coupled_model, ρₒ, cₚ) - ocean = coupled_model.ocean - heat_flux_field = current_heat_flux_field(coupled_model) - ohc_tendency_field = Field(ρₒ * cₚ * Integral(ocean.model.timestepper.Gⁿ.T, dims=(1, 3))) - compute!(ohc_tendency_field) + ohc_tendency = Field(ρᵒᶜ * cᵒᶜ * Integral(T_tendency, dims=(1, 3))) + compute!(ohc_tendency) - model_time = Float64(ocean.model.clock.time) - model_iteration = Int(ocean.model.clock.iteration) - state = get(meridional_heat_transport_states, coupled_model, nothing) - state === nothing && (state = initialize_mht_state!(coupled_model, heat_flux_field, ohc_tendency_field, model_time, model_iteration)) + model_time = esm.ocean.model.clock.time + model_iteration = esm.ocean.model.clock.iteration + state = get(meridional_heat_transport_states, esm, nothing) + state === nothing && + (state = initialize_mht_state!(esm, heat_flux, ohc_tendency, model_time, model_iteration)) if model_iteration != state.last_iteration Δt = max(0.0, model_time - state.last_time) if Δt == 0.0 && model_iteration > 0 - Δt = Float64(ocean.model.clock.Δt) + Δt = ocean.model.clock.Δt end set!(state.cumulative_ohc_tendency, state.cumulative_ohc_tendency + Δt * ohc_tendency_field) @@ -102,18 +102,21 @@ function ohc_tendency_mht(coupled_model, ρₒ, cₚ) state.last_iteration = model_iteration end - flux_int = Integral(state.cumulative_heat_flux, dims=1) + ∫heat_flux = Integral(state.cumulative_heat_flux, dims=1) - return CumulativeIntegral(state.cumulative_ohc_tendency - flux_int, dims=(2)) + return CumulativeIntegral(state.cumulative_ohc_tendency - ∫heat_flux, dims=(2)) end -function instantaneous_vt_mht(coupled_model, ρₒ, cₚ, T_ref) - ocean_model = coupled_model.ocean.model - return ρₒ * cₚ * Integral(ocean_model.velocities.v * (ocean_model.tracers.T - T_ref), dims=(1, 3)) +function instantaneous_vt_mht(esm, reference_temperature) + ρᵒᶜ = reference_density(esm.ocean) + cᵒᶜ = heat_capacity(esm.ocean) + T = esm.ocean.model.tracers.T + v = esm.ocean.model.velocities.v + return ρᵒᶜ * cᵒᶜ * Integral(v * (T - reference_temperature), dims=(1, 3)) end -function checkpoint_auxiliary_state(coupled_model::EarthSystemModel) - state = get(meridional_heat_transport_states, coupled_model, nothing) +function checkpoint_auxiliary_state(esm::EarthSystemModel) + state = get(meridional_heat_transport_states, esm, nothing) state === nothing && return nothing return ( @@ -126,19 +129,19 @@ function checkpoint_auxiliary_state(coupled_model::EarthSystemModel) ) end -function restore_auxiliary_state!(coupled_model::EarthSystemModel, auxiliary_state) +function restore_auxiliary_state!(esm::EarthSystemModel, auxiliary_state) auxiliary_state === nothing && return nothing hasproperty(auxiliary_state, :meridional_heat_transport) || return nothing mht_state = auxiliary_state.meridional_heat_transport mht_state === nothing && return nothing - heat_flux_field = current_heat_flux_field(coupled_model) - ohc_tendency_template = Field(Integral(coupled_model.ocean.model.timestepper.Gⁿ.T, dims=(1, 3))) + heat_flux_field = current_heat_flux_field(esm) + ohc_tendency_template = Field(Integral(esm.ocean.model.timestepper.Gⁿ.T, dims=(1, 3))) compute!(ohc_tendency_template) - reset_meridional_heat_transport_state!(coupled_model) - state = initialize_mht_state!(coupled_model, + reset_meridional_heat_transport_state!(esm) + state = initialize_mht_state!(esm, heat_flux_field, ohc_tendency_template, Float64(mht_state.last_time), From c65cde625a11bfe72286ec9a1d4a6ddc15da1553 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 3 Mar 2026 12:12:12 +0100 Subject: [PATCH 10/63] fix docstring --- src/Diagnostics/meridional_heat_transport.jl | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index a9dda6b66..0593ed0eb 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -13,19 +13,24 @@ const meridional_heat_transport_states = IdDict{Any, MeridionalHeatTransportStat """ meridional_heat_transport(esm::EarthSystemModel; method=:ohc_tendency, reference_temperature=0.0) -Compute meridional heat transport with one of two methods: +Return the meridional heat transport for the coupled `esm` using either of two methods: 1. `method = :ohc_tendency` (default): - `MHT = CumulativeIntegral((∫ ρᵒᶜ cᵒᶜ ∫ₓ∫z ∂ₜT dt) - ∫ₓ(∫Q dt), dims=(2))`, - where `Q` is `"heat_flux"` from `InterfaceFluxOutputs`. + ```math + ∫_{y_S}^y dy [(∫dt ∫dx ∫dz ρᵒᶜ cᵒᶜ ∂ₜT) - (∫dt ∫dx Q)] + ``` + + where `Q` is the net heat flux into the ocean. 2. `method = :vt_instantaneous`: - `MHT = ρᵒᶜ cᵒᶜ ∫ₓ∫z v * (T - T_ref)` + ```math + ρᵒᶜ cᵒᶜ ∫dx ∫dz v (T - T_{\\rm reference}) + ``` -`ρᵒᶜ` and `cᵒᶜ` are the ocean reference density and heat capacity respectively -and they are inferred from `esm.ocean`. +Above, ``ρᵒᶜ`` and ``cᵒᶜ`` are the ocean reference density and heat capacity respectively +and they are inferred from the ocean component, `esm.ocean`. """ function meridional_heat_transport(esm::EarthSystemModel; method=:ohc_tendency, reference_temperature=0.0) if method === :ohc_tendency @@ -33,7 +38,7 @@ function meridional_heat_transport(esm::EarthSystemModel; method=:ohc_tendency, elseif method === :vt_instantaneous return instantaneous_vt_mht(esm, reference_temperature) else - throw(ArgumentError("Unknown MHT method=$(repr(method)). Supported methods are :ohc_tendency and :vt_instantaneous.")) + throw(ArgumentError("Unknown method=$(repr(method)). Supported methods are :ohc_tendency and :vt_instantaneous.")) end end From 762fded9457c63c046e7b74862300e75df3b0177 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 3 Mar 2026 12:13:02 +0100 Subject: [PATCH 11/63] fix docstring --- src/Diagnostics/meridional_heat_transport.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 0593ed0eb..73dc7e0a1 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -47,9 +47,10 @@ allocate_storage_like(field) = Field(instantiated_location(field), field.grid; i """ reset_meridional_heat_transport_state!(esm) -Clear cached MHT state (cumulative OHC tendency, cumulative heat flux, -and clock metadata) for a coupled `esm`. This method should be called before - restarting from a checkpoint or when reusing a model object. +Clear cached for the meridional heat transport state +(cumulative OHC tendency, cumulative heat flux, and clock metadata) +for a coupled `esm`. This method should be called before +restarting from a checkpoint or when reusing a model object. """ function reset_meridional_heat_transport_state!(esm) pop!(meridional_heat_transport_states, esm, nothing) From 4b7a2b6f56921f5276cd3e5b10ad66d2c8f614f9 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 3 Mar 2026 12:14:28 +0100 Subject: [PATCH 12/63] fix docstring --- src/Diagnostics/meridional_heat_transport.jl | 28 ++++++++++++-------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 73dc7e0a1..82f5fb92b 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -15,22 +15,28 @@ const meridional_heat_transport_states = IdDict{Any, MeridionalHeatTransportStat Return the meridional heat transport for the coupled `esm` using either of two methods: -1. `method = :ohc_tendency` (default): +Keyword Arguments +================= - ```math - ∫_{y_S}^y dy [(∫dt ∫dx ∫dz ρᵒᶜ cᵒᶜ ∂ₜT) - (∫dt ∫dx Q)] - ``` +* `method`: Determines the method the computation is done. Options are: + 1. `method = :ohc_tendency` (default): - where `Q` is the net heat flux into the ocean. + ```math + ∫_{y_S}^y dy [(∫dt ∫dx ∫dz ρᵒᶜ cᵒᶜ ∂ₜT) - (∫dt ∫dx Q)] + ``` -2. `method = :vt_instantaneous`: + where `Q` is the net heat flux into the ocean. - ```math - ρᵒᶜ cᵒᶜ ∫dx ∫dz v (T - T_{\\rm reference}) - ``` + 2. `method = :vt_instantaneous`: -Above, ``ρᵒᶜ`` and ``cᵒᶜ`` are the ocean reference density and heat capacity respectively -and they are inferred from the ocean component, `esm.ocean`. + ```math + ρᵒᶜ cᵒᶜ ∫dx ∫dz v (T - T_{\\rm reference}) + ``` + + Above, ``ρᵒᶜ`` and ``cᵒᶜ`` are the ocean reference density and heat capacity respectively + and they are inferred from the ocean component, `esm.ocean`. + +* `reference_temperature`: The reference temperature used for `:vt_instantaneous` method. """ function meridional_heat_transport(esm::EarthSystemModel; method=:ohc_tendency, reference_temperature=0.0) if method === :ohc_tendency From bbb39bf020211e7a4cf84a4fcc26d1523f833867 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 3 Mar 2026 14:34:51 +0100 Subject: [PATCH 13/63] simpler --- src/Diagnostics/meridional_heat_transport.jl | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 82f5fb92b..c472bb8d2 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -75,22 +75,10 @@ function initialize_mht_state!(esm, heat_flux_field, ohc_tendency_field, time, i return state end -function current_heat_flux_field(esm) - flux_outputs = InterfaceFluxOutputs(esm; - isolate_sea_ice=false, - units=:physical, - reference_salinity=35) - - heat_flux_raw = haskey(flux_outputs, "heat_flux") ? flux_outputs["heat_flux"] : flux_outputs[:heat_flux] - heat_flux_field = Field(heat_flux_raw) - compute!(heat_flux_field) - return heat_flux_field -end - function ohc_tendency_mht(esm) ρᵒᶜ = reference_density(esm.ocean) cᵒᶜ = heat_capacity(esm.ocean) - heat_flux = current_heat_flux_field(esm) + heat_flux = net_ocean_heat_flux(esm) T_tendency = esm.ocean.model.timestepper.Gⁿ.T ohc_tendency = Field(ρᵒᶜ * cᵒᶜ * Integral(T_tendency, dims=(1, 3))) @@ -148,7 +136,7 @@ function restore_auxiliary_state!(esm::EarthSystemModel, auxiliary_state) mht_state = auxiliary_state.meridional_heat_transport mht_state === nothing && return nothing - heat_flux_field = current_heat_flux_field(esm) + heat_flux_field = net_ocean_heat_flux(esm) ohc_tendency_template = Field(Integral(esm.ocean.model.timestepper.Gⁿ.T, dims=(1, 3))) compute!(ohc_tendency_template) From 73ad7d5cdde9ad88d287716636449415514d20b7 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 3 Mar 2026 15:18:41 +0100 Subject: [PATCH 14/63] changes --- src/Diagnostics/Diagnostics.jl | 9 +- src/Diagnostics/meridional_heat_transport.jl | 121 +++++++++++-------- src/NumericalEarth.jl | 3 +- 3 files changed, 78 insertions(+), 55 deletions(-) diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl index a4f3ee59f..d349b862d 100644 --- a/src/Diagnostics/Diagnostics.jl +++ b/src/Diagnostics/Diagnostics.jl @@ -1,21 +1,18 @@ module Diagnostics -<<<<<<< HEAD -export MixedLayerDepthField, MixedLayerDepthOperand, Meridional_Heat_Transport, reset_meridional_heat_transport_state! -======= export MixedLayerDepthField, MixedLayerDepthOperand +export meridional_heat_transport, OceanHeatContentTendencyMethod, MeridionalHeatFluxMethod 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 ->>>>>>> main using Oceananigans using Oceananigans.Architectures: architecture -using Oceananigans.Models: buoyancy_operation -using Oceananigans.Grids: new_data, inactive_cell, znode using Oceananigans.BoundaryConditions: FieldBoundaryConditions, fill_halo_regions! using Oceananigans.Fields: FieldStatus +using Oceananigans.Grids: new_data, inactive_cell, znode +using Oceananigans.Models: buoyancy_operation using Oceananigans.Utils: launch! using KernelAbstractions: @index, @kernel using NumericalEarth.EarthSystemModels: EarthSystemModel diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index c472bb8d2..00fb165a6 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -1,87 +1,111 @@ -import ..EarthSystemModels: EarthSystemModel, checkpoint_auxiliary_state, restore_auxiliary_state!, - reference_density, heat_capacity - -mutable struct MeridionalHeatTransportState - cumulative_ohc_tendency::Any - cumulative_heat_flux::Any - last_time::Float64 - last_iteration::Int +using ..EarthSystemModels: reference_density, heat_capacity + +import ..EarthSystemModels: EarthSystemModel, checkpoint_auxiliary_state, restore_auxiliary_state! + +struct OceanHeatContentTendencyMethod end +struct MeridionalHeatFluxMethod end + +mutable struct MeridionalHeatTransportState{FT, FL} + cumulative_ohc_tendency :: FL + cumulative_heat_flux :: FL + last_time :: FT + last_iteration :: Int end const meridional_heat_transport_states = IdDict{Any, MeridionalHeatTransportState}() """ - meridional_heat_transport(esm::EarthSystemModel; method=:ohc_tendency, reference_temperature=0.0) + meridional_heat_transport(esm::EarthSystemModel, OceanHeatContentTendencyMethod(); + reference_temperature=0.0) Return the meridional heat transport for the coupled `esm` using either of two methods: -Keyword Arguments -================= +Arguments +========= -* `method`: Determines the method the computation is done. Options are: - 1. `method = :ohc_tendency` (default): +* `esm`: An EarthSystemModel. + +* `OceanHeatContentTendencyMethod()` or `MeridionalHeatFluxMethod()` denoting the method + that the meridional heat transport is computed. + + 1. For `OceanHeatContentTendencyMethod()` (default), the meridional heat transport + is computed via: ```math - ∫_{y_S}^y dy [(∫dt ∫dx ∫dz ρᵒᶜ cᵒᶜ ∂ₜT) - (∫dt ∫dx Q)] + ∫_{y_S}^y dy [(∫dt ∫dx ∫dz ρᵒᶜ cᵒᶜ ∂ₜT) - (∫dt ∫dx 𝒬)] ``` - where `Q` is the net heat flux into the ocean. + where ``y_S`` is the Southern-most latitude of the domain and ``𝒬`` is the + net heat flux into the ocean. - 2. `method = :vt_instantaneous`: + 2. For `MeridionalHeatFluxMethod()`, the meridional heat transport is computed via: ```math - ρᵒᶜ cᵒᶜ ∫dx ∫dz v (T - T_{\\rm reference}) + ρᵒᶜ cᵒᶜ ∫dx ∫dz v (T - T_{\\rm ref}) ``` Above, ``ρᵒᶜ`` and ``cᵒᶜ`` are the ocean reference density and heat capacity respectively - and they are inferred from the ocean component, `esm.ocean`. + and they are inferred from the ocean component, `esm.ocean`, and ``T_{\\rm ref}`` is a + reference temperature. + + +Keyword Arguments +================= * `reference_temperature`: The reference temperature used for `:vt_instantaneous` method. """ -function meridional_heat_transport(esm::EarthSystemModel; method=:ohc_tendency, reference_temperature=0.0) - if method === :ohc_tendency - return ohc_tendency_mht(esm) - elseif method === :vt_instantaneous - return instantaneous_vt_mht(esm, reference_temperature) - else - throw(ArgumentError("Unknown method=$(repr(method)). Supported methods are :ohc_tendency and :vt_instantaneous.")) - end +function meridional_heat_transport(esm::EarthSystemModel, ::OceanHeatContentTendencyMethod; + reference_temperature=0.0) + return meridional_heat_transport_via_ocean_heat_content_tendency(esm) +end + +function meridional_heat_transport(esm::EarthSystemModel, ::MeridionalHeatFluxMethod; + reference_temperature=0.0) + return meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature) end -allocate_storage_like(field) = Field(instantiated_location(field), field.grid; indices=indices(field)) +meridional_heat_transport(esm::EarthSystemModel) = meridional_heat_transport(esm, OceanHeatContentTendencyMethod()) + +meridional_heat_transport(esm::EarthSystemModel, method; reference_temperature=0.0) = + throw(ArgumentError("Unknown method $(method); choose one of OceanHeatContentTendencyMethod() or MeridionalHeatFluxMethod().")) """ reset_meridional_heat_transport_state!(esm) -Clear cached for the meridional heat transport state -(cumulative OHC tendency, cumulative heat flux, and clock metadata) -for a coupled `esm`. This method should be called before -restarting from a checkpoint or when reusing a model object. +Clear cached for the meridional heat transport state (cumulative ocean heat +content tendency, cumulative heat flux, and clock metadata) for a coupled +`esm`. This method should be called before restarting from a checkpoint or +when reusing a model object. """ function reset_meridional_heat_transport_state!(esm) pop!(meridional_heat_transport_states, esm, nothing) return nothing end -function initialize_mht_state!(esm, heat_flux_field, ohc_tendency_field, time, iteration) - cumulative_ohc_tendency = allocate_storage_like(ohc_tendency_field) +function initialize_mht_state!(esm, heat_flux, ohc_tendency, time, iteration) + cumulative_ohc_tendency = deepcopy(ohc_tendency) set!(cumulative_ohc_tendency, 0) - cumulative_heat_flux = allocate_storage_like(heat_flux_field) + cumulative_heat_flux = deepcopy(heat_flux) set!(cumulative_heat_flux, 0) - state = MeridionalHeatTransportState(cumulative_ohc_tendency, cumulative_heat_flux, time, iteration) + FT = eltype(esm) + FL = typeof(cumulative_ohc_tendency) + state = MeridionalHeatTransportState{FT, FL}(cumulative_ohc_tendency, + cumulative_heat_flux, + time, + iteration) meridional_heat_transport_states[esm] = state return state end -function ohc_tendency_mht(esm) +function meridional_heat_transport_via_ocean_heat_content_tendency(esm) ρᵒᶜ = reference_density(esm.ocean) cᵒᶜ = heat_capacity(esm.ocean) heat_flux = net_ocean_heat_flux(esm) - T_tendency = esm.ocean.model.timestepper.Gⁿ.T + ∂T_∂t = esm.ocean.model.timestepper.Gⁿ.T - ohc_tendency = Field(ρᵒᶜ * cᵒᶜ * Integral(T_tendency, dims=(1, 3))) + ohc_tendency = Field(Integral(ρᵒᶜ * cᵒᶜ * ∂T_∂t, dims=(1, 3))) compute!(ohc_tendency) model_time = esm.ocean.model.clock.time @@ -107,7 +131,7 @@ function ohc_tendency_mht(esm) return CumulativeIntegral(state.cumulative_ohc_tendency - ∫heat_flux, dims=(2)) end -function instantaneous_vt_mht(esm, reference_temperature) +function meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature) ρᵒᶜ = reference_density(esm.ocean) cᵒᶜ = heat_capacity(esm.ocean) T = esm.ocean.model.tracers.T @@ -136,19 +160,20 @@ function restore_auxiliary_state!(esm::EarthSystemModel, auxiliary_state) mht_state = auxiliary_state.meridional_heat_transport mht_state === nothing && return nothing - heat_flux_field = net_ocean_heat_flux(esm) - ohc_tendency_template = Field(Integral(esm.ocean.model.timestepper.Gⁿ.T, dims=(1, 3))) - compute!(ohc_tendency_template) + heat_flux = net_ocean_heat_flux(esm) + ∂T_∂t = esm.ocean.model.timestepper.Gⁿ.T + ohc_tendency = Field(Integral(∂T_∂t, dims=(1, 3))) + compute!(ohc_tendency) reset_meridional_heat_transport_state!(esm) state = initialize_mht_state!(esm, - heat_flux_field, - ohc_tendency_template, - Float64(mht_state.last_time), - Int(mht_state.last_iteration)) + heat_flux, + ohc_tendency, + mht_state.last_time, + mht_state.last_iteration) set!(state.cumulative_ohc_tendency, mht_state.cumulative_ohc_tendency) set!(state.cumulative_heat_flux, mht_state.cumulative_heat_flux) - state.last_time = Float64(mht_state.last_time) - state.last_iteration = Int(mht_state.last_iteration) + state.last_time = mht_state.last_time + state.last_iteration = mht_state.last_iteration return nothing end diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl index b6b2850c7..8495aec20 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, OceanHeatContentTendencyMethod, MeridionalHeatFluxMethod using Oceananigans using Oceananigans.Operators: ℑxyᶠᶜᵃ, ℑxyᶜᶠᵃ From c5a256038a7af21f16aab4934066f478a31bea8b Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 3 Mar 2026 15:55:06 +0100 Subject: [PATCH 15/63] use Bounded y in toy example --- test/test_diagnostics.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/test_diagnostics.jl b/test/test_diagnostics.jl index 7d55c7d37..b56087690 100644 --- a/test/test_diagnostics.jl +++ b/test/test_diagnostics.jl @@ -5,6 +5,7 @@ using Oceananigans.Models: buoyancy_operation using NumericalEarth.Diagnostics: MixedLayerDepthField, MixedLayerDepthOperand using SeawaterPolynomials: TEOS10EquationOfState +#= for arch in test_architectures, dataset in (ECCO4Monthly(),) A = typeof(arch) @info "Testing MixedLayerDepthField with $(typeof(dataset)) on $A" @@ -56,16 +57,16 @@ for arch in test_architectures, dataset in (ECCO4Monthly(),) end end end - +=# for arch in test_architectures A = typeof(arch) @info "Testing InterfaceFluxOutputs on $A" @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, From a09ee3a2f15c67b3c748e2e3e18ca057e89e9ceb Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 3 Mar 2026 15:55:06 +0100 Subject: [PATCH 16/63] refactors --- src/Diagnostics/meridional_heat_transport.jl | 70 +++++++++++--------- 1 file changed, 40 insertions(+), 30 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 00fb165a6..3b050a512 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -5,9 +5,9 @@ import ..EarthSystemModels: EarthSystemModel, checkpoint_auxiliary_state, restor struct OceanHeatContentTendencyMethod end struct MeridionalHeatFluxMethod end -mutable struct MeridionalHeatTransportState{FT, FL} - cumulative_ohc_tendency :: FL - cumulative_heat_flux :: FL +mutable struct MeridionalHeatTransportState{FT, OHC, HF} + cumulative_∫ohc_tendency :: OHC + cumulative_∫heat_flux :: HF last_time :: FT last_iteration :: Int end @@ -82,19 +82,20 @@ function reset_meridional_heat_transport_state!(esm) return nothing end -function initialize_mht_state!(esm, heat_flux, ohc_tendency, time, iteration) - cumulative_ohc_tendency = deepcopy(ohc_tendency) - set!(cumulative_ohc_tendency, 0) +function initialize_mht_state!(esm, ∫heat_flux, ∫ohc_tendency, time, iteration) + cumulative_∫ohc_tendency = deepcopy(∫ohc_tendency) + set!(cumulative_∫ohc_tendency, 0) - cumulative_heat_flux = deepcopy(heat_flux) - set!(cumulative_heat_flux, 0) + cumulative_∫heat_flux = deepcopy(∫heat_flux) + set!(cumulative_∫heat_flux, 0) FT = eltype(esm) - FL = typeof(cumulative_ohc_tendency) - state = MeridionalHeatTransportState{FT, FL}(cumulative_ohc_tendency, - cumulative_heat_flux, - time, - iteration) + OHC = typeof(cumulative_∫ohc_tendency) + HF = typeof(cumulative_∫heat_flux) + state = MeridionalHeatTransportState{FT, OHC, HF}(cumulative_∫ohc_tendency, + cumulative_∫heat_flux, + time, + iteration) meridional_heat_transport_states[esm] = state return state end @@ -105,14 +106,17 @@ function meridional_heat_transport_via_ocean_heat_content_tendency(esm) heat_flux = net_ocean_heat_flux(esm) ∂T_∂t = esm.ocean.model.timestepper.Gⁿ.T - ohc_tendency = Field(Integral(ρᵒᶜ * cᵒᶜ * ∂T_∂t, dims=(1, 3))) - compute!(ohc_tendency) + ∫ohc_tendency = Field(Integral(ρᵒᶜ * cᵒᶜ * ∂T_∂t, dims=(1, 3))) + compute!(∫ohc_tendency) + + ∫heat_flux = Field(Integral(heat_flux, dims=1)) + compute!(∫heat_flux) model_time = esm.ocean.model.clock.time model_iteration = esm.ocean.model.clock.iteration state = get(meridional_heat_transport_states, esm, nothing) state === nothing && - (state = initialize_mht_state!(esm, heat_flux, ohc_tendency, model_time, model_iteration)) + (state = initialize_mht_state!(esm, ∫heat_flux, ∫ohc_tendency, model_time, model_iteration)) if model_iteration != state.last_iteration Δt = max(0.0, model_time - state.last_time) @@ -120,15 +124,14 @@ function meridional_heat_transport_via_ocean_heat_content_tendency(esm) Δt = ocean.model.clock.Δt end - set!(state.cumulative_ohc_tendency, state.cumulative_ohc_tendency + Δt * ohc_tendency_field) - set!(state.cumulative_heat_flux, state.cumulative_heat_flux + Δt * heat_flux_field) + set!(state.cumulative_∫ohc_tendency, state.cumulative_ohc_∫tendency + Δt * ∫ohc_tendency) + set!(state.cumulative_∫heat_flux, state.cumulative_∫heat_flux + Δt * ∫heat_flux) + state.last_time = model_time state.last_iteration = model_iteration end - ∫heat_flux = Integral(state.cumulative_heat_flux, dims=1) - - return CumulativeIntegral(state.cumulative_ohc_tendency - ∫heat_flux, dims=(2)) + return Field(CumulativeIntegral(state.cumulative_∫ohc_tendency - ∫heat_flux, dims=2)) end function meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature) @@ -136,7 +139,8 @@ function meridional_heat_transport_via_meridional_heat_flux(esm; reference_tempe cᵒᶜ = heat_capacity(esm.ocean) T = esm.ocean.model.tracers.T v = esm.ocean.model.velocities.v - return ρᵒᶜ * cᵒᶜ * Integral(v * (T - reference_temperature), dims=(1, 3)) + mht = Field(Integral(ρᵒᶜ * cᵒᶜ * v * (T - reference_temperature), dims=(1, 3))) + return mht end function checkpoint_auxiliary_state(esm::EarthSystemModel) @@ -145,8 +149,8 @@ function checkpoint_auxiliary_state(esm::EarthSystemModel) return ( meridional_heat_transport = ( - cumulative_ohc_tendency = Array(interior(state.cumulative_ohc_tendency)), - cumulative_heat_flux = Array(interior(state.cumulative_heat_flux)), + cumulative_ohc_tendency = Array(interior(state.cumulative_∫ohc_tendency)), + cumulative_heat_flux = Array(interior(state.cumulative_∫heat_flux)), last_time = state.last_time, last_iteration = state.last_iteration ), @@ -160,19 +164,25 @@ function restore_auxiliary_state!(esm::EarthSystemModel, auxiliary_state) mht_state = auxiliary_state.meridional_heat_transport mht_state === nothing && return nothing + ρᵒᶜ = reference_density(esm.ocean) + cᵒᶜ = heat_capacity(esm.ocean) heat_flux = net_ocean_heat_flux(esm) ∂T_∂t = esm.ocean.model.timestepper.Gⁿ.T - ohc_tendency = Field(Integral(∂T_∂t, dims=(1, 3))) - compute!(ohc_tendency) + + ∫ohc_tendency = Field(Integral(ρᵒᶜ * cᵒᶜ * ∂T_∂t, dims=(1, 3))) + compute!(∫ohc_tendency) + + ∫heat_flux = Field(Integral(heat_flux, dims=1)) + compute!(∫heat_flux) reset_meridional_heat_transport_state!(esm) state = initialize_mht_state!(esm, - heat_flux, - ohc_tendency, + ∫heat_flux, + ∫ohc_tendency, mht_state.last_time, mht_state.last_iteration) - set!(state.cumulative_ohc_tendency, mht_state.cumulative_ohc_tendency) - set!(state.cumulative_heat_flux, mht_state.cumulative_heat_flux) + set!(state.cumulative_∫ohc_tendency, mht_state.∫cumulative_ohc_tendency) + set!(state.cumulative_∫heat_flux, mht_state.∫cumulative_heat_flux) state.last_time = mht_state.last_time state.last_iteration = mht_state.last_iteration return nothing From 56322f4f3b9578f62499c78ff190425f803f53d5 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 3 Mar 2026 15:59:36 +0100 Subject: [PATCH 17/63] split using and imports --- src/Diagnostics/meridional_heat_transport.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 3b050a512..acde11e18 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -1,6 +1,6 @@ -using ..EarthSystemModels: reference_density, heat_capacity +using ..EarthSystemModels: EarthSystemModel, reference_density, heat_capacity -import ..EarthSystemModels: EarthSystemModel, checkpoint_auxiliary_state, restore_auxiliary_state! +import ..EarthSystemModels: checkpoint_auxiliary_state, restore_auxiliary_state! struct OceanHeatContentTendencyMethod end struct MeridionalHeatFluxMethod end From c5da1907f93a65a7242cdab4117535c3c04ab630 Mon Sep 17 00:00:00 2001 From: Taimoor Sohail Date: Tue, 3 Mar 2026 23:17:43 +0100 Subject: [PATCH 18/63] Update src/Diagnostics/meridional_heat_transport.jl Co-authored-by: Gregory L. Wagner --- src/Diagnostics/meridional_heat_transport.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index acde11e18..457c27cf3 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -119,7 +119,7 @@ function meridional_heat_transport_via_ocean_heat_content_tendency(esm) (state = initialize_mht_state!(esm, ∫heat_flux, ∫ohc_tendency, model_time, model_iteration)) if model_iteration != state.last_iteration - Δt = max(0.0, model_time - state.last_time) + Δt = max(0, model_time - state.last_time) if Δt == 0.0 && model_iteration > 0 Δt = ocean.model.clock.Δt end From 2a23b79fdf54f6e34499e1daa4111c6b152324e6 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Wed, 4 Mar 2026 14:56:31 +0100 Subject: [PATCH 19/63] fix OHC diag; return AbstractOperations --- src/Diagnostics/meridional_heat_transport.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index acde11e18..630195ba8 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -32,7 +32,7 @@ Arguments is computed via: ```math - ∫_{y_S}^y dy [(∫dt ∫dx ∫dz ρᵒᶜ cᵒᶜ ∂ₜT) - (∫dt ∫dx 𝒬)] + ∫_{y_S}^y \\mathrm{d}y \\left( ∫\\mathrm{d}t ∫\\mathrm{d}x ∫\\mathrm{d}z \\, ρᵒᶜ cᵒᶜ ∂ₜT - ∫\\mathrm{d}t ∫\\mathrm{d}x \\, 𝒬 \\right) ``` where ``y_S`` is the Southern-most latitude of the domain and ``𝒬`` is the @@ -41,7 +41,7 @@ Arguments 2. For `MeridionalHeatFluxMethod()`, the meridional heat transport is computed via: ```math - ρᵒᶜ cᵒᶜ ∫dx ∫dz v (T - T_{\\rm ref}) + ρᵒᶜ cᵒᶜ ∫\\mathrm{d}x ∫\\mathrm{d}z \\, v (T - T_{\\rm ref}) ``` Above, ``ρᵒᶜ`` and ``cᵒᶜ`` are the ocean reference density and heat capacity respectively @@ -52,7 +52,7 @@ Arguments Keyword Arguments ================= -* `reference_temperature`: The reference temperature used for `:vt_instantaneous` method. +* `reference_temperature`: The reference temperature used for `MeridionalHeatFluxMethod()`. """ function meridional_heat_transport(esm::EarthSystemModel, ::OceanHeatContentTendencyMethod; reference_temperature=0.0) @@ -106,7 +106,7 @@ function meridional_heat_transport_via_ocean_heat_content_tendency(esm) heat_flux = net_ocean_heat_flux(esm) ∂T_∂t = esm.ocean.model.timestepper.Gⁿ.T - ∫ohc_tendency = Field(Integral(ρᵒᶜ * cᵒᶜ * ∂T_∂t, dims=(1, 3))) + ∫ohc_tendency = Field(sum(Field(Integral(ρᵒᶜ * cᵒᶜ * ∂T_∂t, dims=1)), dims=3)) compute!(∫ohc_tendency) ∫heat_flux = Field(Integral(heat_flux, dims=1)) @@ -131,7 +131,7 @@ function meridional_heat_transport_via_ocean_heat_content_tendency(esm) state.last_iteration = model_iteration end - return Field(CumulativeIntegral(state.cumulative_∫ohc_tendency - ∫heat_flux, dims=2)) + return CumulativeIntegral(state.cumulative_∫ohc_tendency - state.cumulative_∫heat_flux, dims=2) end function meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature) @@ -139,7 +139,7 @@ function meridional_heat_transport_via_meridional_heat_flux(esm; reference_tempe cᵒᶜ = heat_capacity(esm.ocean) T = esm.ocean.model.tracers.T v = esm.ocean.model.velocities.v - mht = Field(Integral(ρᵒᶜ * cᵒᶜ * v * (T - reference_temperature), dims=(1, 3))) + mht = Integral(ρᵒᶜ * cᵒᶜ * v * (T - reference_temperature), dims=(1, 3)) return mht end From 841607a9e533bd54394d92f320ed17a179c98b0b Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Wed, 4 Mar 2026 14:56:46 +0100 Subject: [PATCH 20/63] add example to compute MHT --- examples/meridional_heat_transport_ecco.jl | 106 +++++++++++++++++++++ 1 file changed, 106 insertions(+) create mode 100644 examples/meridional_heat_transport_ecco.jl diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl new file mode 100644 index 000000000..8d23dcf78 --- /dev/null +++ b/examples/meridional_heat_transport_ecco.jl @@ -0,0 +1,106 @@ +using NumericalEarth +using Oceananigans +using Oceananigans.Units +using Dates +using Statistics +using Printf + +arch = CPU() +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=90days) + + +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(", extrema(To): (%.1f, %.1f) ᵒC, mean(To(z=0)): %.1f ᵒC", Tmin, Tmax, Tavg) + msg4 = @sprintf(", max(e): %.2f m² s⁻²", emax) + msg5 = @sprintf(", wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg3 * msg4 * msg5 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(simulation, progress, TimeInterval(1days)) + +mht_OHC = meridional_heat_transport(esm, OceanHeatContentTendencyMethod()) |> Field +mht_vT = meridional_heat_transport(esm, MeridionalHeatFluxMethod()) |> Field + +mht_outputs = (; mht_OHC, mht_vT) + +ocean.output_writers[:mth] = JLD2Writer(ocean.model, mht_outputs; + schedule = TimeInterval(1hours), + filename = "ocean_one_degree_mht", + overwrite_existing = true) + + +run!(simulation) + +## + +using GLMakie + +fig = Figure() +ax1 = Axis(fig[1, 1]) +ax2 = Axis(fig[2, 1]) + +φ = φnodes(grid, Face()) +lines!(ax1, φ, mht_method1 / 1e15, linewidth=4) +lines!(ax2, φ, mht_method2 / 1e15, linewidth=4) +current_figure() + +display(fig) From c815495f113a3a2aefcab777e407ea78c4bb8463 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 5 Mar 2026 04:57:58 -0500 Subject: [PATCH 21/63] updates --- examples/meridional_heat_transport_ecco.jl | 37 +++++++++++++++------- 1 file changed, 26 insertions(+), 11 deletions(-) diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl index 8d23dcf78..9a7fdcf3a 100644 --- a/examples/meridional_heat_transport_ecco.jl +++ b/examples/meridional_heat_transport_ecco.jl @@ -5,7 +5,7 @@ using Dates using Statistics using Printf -arch = CPU() +arch = GPU() Nx = 360 Ny = 180 Nz = 50 @@ -45,7 +45,7 @@ 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=90days) +simulation = Simulation(esm; Δt=20minutes, stop_time=2*365days) wall_time = Ref(time_ns()) @@ -83,24 +83,39 @@ mht_vT = meridional_heat_transport(esm, MeridionalHeatFluxMethod()) |> Field mht_outputs = (; mht_OHC, mht_vT) ocean.output_writers[:mth] = JLD2Writer(ocean.model, mht_outputs; - schedule = TimeInterval(1hours), + schedule = TimeInterval(3hours), filename = "ocean_one_degree_mht", overwrite_existing = true) run!(simulation) -## -using GLMakie +mht_OCH = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht_OHC"; backend = OnDisk()) +mht_vT = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht_vT"; backend = OnDisk()) + +times = mht_OCH.times +Nt = length(times) + +mht_OCH_mean = deepcopy(mht_OCH[1][1, :, 1]) +mht_vT_mean = deepcopy(mht_vT[1][1, :, 1]) + +for j in 1:Nt + mht_OCH_mean += mht_OCH[j][1, :, 1] + mht_vT_mean += mht_vT[j][1, :, 1] +end +@. mht_OCH_mean = mht_OCH_mean / Nt +@. mht_vT_mean = mht_vT_mean / Nt + + +using CairoMakie fig = Figure() -ax1 = Axis(fig[1, 1]) -ax2 = Axis(fig[2, 1]) +ax1 = Axis(fig[1, 1], xlabel="latitude (deg)", ylabel="MHT (PW)") +ax2 = Axis(fig[2, 1], xlabel="latitude (deg)", ylabel="MHT (PW)") φ = φnodes(grid, Face()) -lines!(ax1, φ, mht_method1 / 1e15, linewidth=4) -lines!(ax2, φ, mht_method2 / 1e15, linewidth=4) -current_figure() +lines!(ax1, φnodes(grid, Center()), diff(mht_OCH_mean[1:Ny+1]) / 1e15, linewidth=4) +lines!(ax2, φ, mht_vT_mean[1:Ny+1] / 1e15, linewidth=4) -display(fig) +save("mht.png", fig) From 77d518a40dc70f1e9ca5ff268947e2b52b9665b5 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 6 Mar 2026 01:25:45 -0500 Subject: [PATCH 22/63] simpler OHC-based MHT diag --- src/Diagnostics/meridional_heat_transport.jl | 224 ++++++++----------- 1 file changed, 90 insertions(+), 134 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 9b6d86c6c..56e80a832 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -1,3 +1,4 @@ +using Oceananigans.TimeSteppers: SplitRungeKuttaTimeStepper, QuasiAdamsBashforth2TimeStepper using ..EarthSystemModels: EarthSystemModel, reference_density, heat_capacity import ..EarthSystemModels: checkpoint_auxiliary_state, restore_auxiliary_state! @@ -5,133 +6,133 @@ import ..EarthSystemModels: checkpoint_auxiliary_state, restore_auxiliary_state! struct OceanHeatContentTendencyMethod end struct MeridionalHeatFluxMethod end -mutable struct MeridionalHeatTransportState{FT, OHC, HF} - cumulative_∫ohc_tendency :: OHC - cumulative_∫heat_flux :: HF - last_time :: FT - last_iteration :: Int -end - -const meridional_heat_transport_states = IdDict{Any, MeridionalHeatTransportState}() - """ - meridional_heat_transport(esm::EarthSystemModel, OceanHeatContentTendencyMethod(); - reference_temperature=0.0) + meridional_heat_transport(esm::EarthSystemModel, MeridionalHeatFluxMethod(); + reference_temperature = 0) -Return the meridional heat transport for the coupled `esm` using either of two methods: +Return the meridional heat transport for the coupled `esm::EarthSystemModel` using either +two methods: either by directly computing the meridional heat flux or indirectly using +the total ocean heat content and the ocean heat uptake. Arguments ========= * `esm`: An EarthSystemModel. -* `OceanHeatContentTendencyMethod()` or `MeridionalHeatFluxMethod()` denoting the method - that the meridional heat transport is computed. +* The method for the computation. Available options are: `MeridionalHeatFluxMethod()` (default) + and `OceanHeatContentTendencyMethod()`. + + `MeridionalHeatFluxMethod()` computes the meridional heat transport directly by summing + the meridional heat flux; `OceanHeatContentTendencyMethod()` computes the meridional heat + transport indirectly using the ocean heat content. + + 1. For `MeridionalHeatFluxMethod()`, 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 heat capacity respectively. + + 2. For `OceanHeatContentTendencyMethod()` we have: + + Let ``T` be three-dimensional (potential) temperature, `ρᵒᶜ` a reference density, + `cᵒᶜ` a heat capacity, `H` the resting depth, and `η` the free-surface elevation. + + The column heat content per unit horizontal area (units of J m⁻²) is: - 1. For `OceanHeatContentTendencyMethod()` (default), the meridional heat transport - is computed via: + ```math + ℋ = ρᵒᶜ cᵒᶜ ∫_{-H}^{η} T \\, \\mathrm{d}z + ``` + + The vertically-integrated heat budget is + + ```math + \\frac{∂ℋ}{∂t} = - \\boldsymbol{∇}_h \\cdot \\boldsymbol{F}_h - 𝒬_{\\rm net} + ℛ + ``` + + where + + * ``\\boldsymbol{F}_h`` is the depth-integrated horizontal heat flux vector (units W m⁻¹), + that includes advection and any parameterized lateral/eddy fluxes, + * ``- 𝒬_{\\rm net}`` is the net surface heat flux into the ocean at the surface + (units W m⁻²), and + * ``ℛ`` is the residual sources/sinks and non-closed terms (e.g. numerics, unaccounted + physics, mass/volume effects). + + The total ocean heat content (OHC) South of latitude ``φ`` is: ```math - ∫_{y_S}^y \\mathrm{d}y \\left( ∫\\mathrm{d}t ∫\\mathrm{d}x ∫\\mathrm{d}z \\, ρᵒᶜ cᵒᶜ ∂ₜT - ∫\\mathrm{d}t ∫\\mathrm{d}x \\, 𝒬 \\right) + \\mathrm{OHC}_S(φ, t) ≡ ∫_{A(φ)} ℋ \\, \\mathrm{d}A ``` - where ``y_S`` is the Southern-most latitude of the domain and ``𝒬`` is the - net heat flux into the ocean. + where `A(φ)` is the ocean area South of latitude ``φ``. - 2. For `MeridionalHeatFluxMethod()`, the meridional heat transport is computed via: + Integrating the vertically-integrated budget over `A(φ)` and using the divergence + theorem we get ```math - ρᵒᶜ cᵒᶜ ∫\\mathrm{d}x ∫\\mathrm{d}z \\, v (T - T_{\\rm ref}) + \\frac{d}{dt} \\, \\mathrm{OHC}_S(φ, t) = + - ∮_{∂A(φ)} \\boldsymbol{F}_h \\cdot \\hat{\\boldsymbol{n}} \\, \\mathrm{d}ℓ + + ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A + + ∫_{A(φ)} ℛ \\, \\mathrm{d}A ``` - Above, ``ρᵒᶜ`` and ``cᵒᶜ`` are the ocean reference density and heat capacity respectively - and they are inferred from the ocean component, `esm.ocean`, and ``T_{\\rm ref}`` is a - reference temperature. + The northward meridional heat transport across latitude ``φ`` is + + ```math + \\mathrm{MHT}(φ, t) ≡ - ∮_{\\mathrm{lat}=φ} \\boldsymbol{F}_h \\cdot \\hat{\boldsymbol{n}} \\, \\mathrm{d}ℓ + ``` + with the sign convention that ``\\mathrm{MHT} > 0`` is northward. + + Ignoring the residual ``ℛ``, the OHC-based diagnostic relation is + + ```math + \\mathrm{MHT} = ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A - \\frac{d}{dt} \\, \\mathrm{OHC}_S + ``` Keyword Arguments ================= -* `reference_temperature`: The reference temperature used for `MeridionalHeatFluxMethod()`. +* `reference_temperature`: The reference temperature (in ᵒC) used for `MeridionalHeatFluxMethod()`; default: 0 ᵒC. """ function meridional_heat_transport(esm::EarthSystemModel, ::OceanHeatContentTendencyMethod; - reference_temperature=0.0) - return meridional_heat_transport_via_ocean_heat_content_tendency(esm) + reference_temperature=0) + return meridional_heat_transport_via_ocean_heat_content(esm) end function meridional_heat_transport(esm::EarthSystemModel, ::MeridionalHeatFluxMethod; - reference_temperature=0.0) - return meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature) + reference_temperature=0) + return meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature = eltype(esm)(reference_temperature)) end meridional_heat_transport(esm::EarthSystemModel) = meridional_heat_transport(esm, OceanHeatContentTendencyMethod()) meridional_heat_transport(esm::EarthSystemModel, method; reference_temperature=0.0) = - throw(ArgumentError("Unknown method $(method); choose one of OceanHeatContentTendencyMethod() or MeridionalHeatFluxMethod().")) - -""" - reset_meridional_heat_transport_state!(esm) + throw(ArgumentError("Unknown method $(method); choose either MeridionalHeatFluxMethod() or OceanHeatContentTendencyMethod().")) -Clear cached for the meridional heat transport state (cumulative ocean heat -content tendency, cumulative heat flux, and clock metadata) for a coupled -`esm`. This method should be called before restarting from a checkpoint or -when reusing a model object. -""" -function reset_meridional_heat_transport_state!(esm) - pop!(meridional_heat_transport_states, esm, nothing) - return nothing -end - -function initialize_mht_state!(esm, ∫heat_flux, ∫ohc_tendency, time, iteration) - cumulative_∫ohc_tendency = deepcopy(∫ohc_tendency) - set!(cumulative_∫ohc_tendency, 0) - - cumulative_∫heat_flux = deepcopy(∫heat_flux) - set!(cumulative_∫heat_flux, 0) - - FT = eltype(esm) - OHC = typeof(cumulative_∫ohc_tendency) - HF = typeof(cumulative_∫heat_flux) - state = MeridionalHeatTransportState{FT, OHC, HF}(cumulative_∫ohc_tendency, - cumulative_∫heat_flux, - time, - iteration) - meridional_heat_transport_states[esm] = state - return state -end - -function meridional_heat_transport_via_ocean_heat_content_tendency(esm) +function meridional_heat_transport_via_ocean_heat_content(esm) ρᵒᶜ = reference_density(esm.ocean) cᵒᶜ = heat_capacity(esm.ocean) - heat_flux = net_ocean_heat_flux(esm) - ∂T_∂t = esm.ocean.model.timestepper.Gⁿ.T - - ∫ohc_tendency = Field(sum(Field(Integral(ρᵒᶜ * cᵒᶜ * ∂T_∂t, dims=1)), dims=3)) - compute!(∫ohc_tendency) + ∂T_∂t = temperature_evolution_tendency(esm.ocean.model.timestepper) + heat_flux = net_ocean_heat_flux(esm) |> Field - ∫heat_flux = Field(Integral(heat_flux, dims=1)) - compute!(∫heat_flux) + ∂ℋ_∂t = Integral(ρᵒᶜ * cᵒᶜ * ∂T_∂t, dims=3) |> Field + ∫sum_dx = Integral(heat_flux + ∂ℋ_∂t, dims=1) |> Field + MHT = CumulativeIntegral(- ∫sum_dx, dims=2) - model_time = esm.ocean.model.clock.time - model_iteration = esm.ocean.model.clock.iteration - state = get(meridional_heat_transport_states, esm, nothing) - state === nothing && - (state = initialize_mht_state!(esm, ∫heat_flux, ∫ohc_tendency, model_time, model_iteration)) - - if model_iteration != state.last_iteration - Δt = max(0, model_time - state.last_time) - if Δt == 0.0 && model_iteration > 0 - Δt = ocean.model.clock.Δt - end - - set!(state.cumulative_∫ohc_tendency, state.cumulative_ohc_∫tendency + Δt * ∫ohc_tendency) - set!(state.cumulative_∫heat_flux, state.cumulative_∫heat_flux + Δt * ∫heat_flux) + return MHT +end - state.last_time = model_time - state.last_iteration = model_iteration - end +temperature_evolution_tendency(timestepper::SplitRungeKuttaTimeStepper) = timestepper.Gⁿ.T - return CumulativeIntegral(state.cumulative_∫ohc_tendency - state.cumulative_∫heat_flux, dims=2) +function temperature_evolution_tendency(timestepper::QuasiAdamsBashforth2TimeStepper) + Gⁿ = timestepper.Gⁿ.T + G⁻ = timestepper.G⁻.T + χ = timestepper.χ + return (3/2 + χ) * Gⁿ - (1/2 + χ) * G⁻ end function meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature) @@ -139,51 +140,6 @@ function meridional_heat_transport_via_meridional_heat_flux(esm; reference_tempe 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 - -function checkpoint_auxiliary_state(esm::EarthSystemModel) - state = get(meridional_heat_transport_states, esm, nothing) - state === nothing && return nothing - - return ( - meridional_heat_transport = ( - cumulative_ohc_tendency = Array(interior(state.cumulative_∫ohc_tendency)), - cumulative_heat_flux = Array(interior(state.cumulative_∫heat_flux)), - last_time = state.last_time, - last_iteration = state.last_iteration - ), - ) -end - -function restore_auxiliary_state!(esm::EarthSystemModel, auxiliary_state) - auxiliary_state === nothing && return nothing - hasproperty(auxiliary_state, :meridional_heat_transport) || return nothing - - mht_state = auxiliary_state.meridional_heat_transport - mht_state === nothing && return nothing - - ρᵒᶜ = reference_density(esm.ocean) - cᵒᶜ = heat_capacity(esm.ocean) - heat_flux = net_ocean_heat_flux(esm) - ∂T_∂t = esm.ocean.model.timestepper.Gⁿ.T - - ∫ohc_tendency = Field(Integral(ρᵒᶜ * cᵒᶜ * ∂T_∂t, dims=(1, 3))) - compute!(∫ohc_tendency) - - ∫heat_flux = Field(Integral(heat_flux, dims=1)) - compute!(∫heat_flux) - - reset_meridional_heat_transport_state!(esm) - state = initialize_mht_state!(esm, - ∫heat_flux, - ∫ohc_tendency, - mht_state.last_time, - mht_state.last_iteration) - set!(state.cumulative_∫ohc_tendency, mht_state.∫cumulative_ohc_tendency) - set!(state.cumulative_∫heat_flux, mht_state.∫cumulative_heat_flux) - state.last_time = mht_state.last_time - state.last_iteration = mht_state.last_iteration - return nothing + MHT = Integral(ρᵒᶜ * cᵒᶜ * v * (T - reference_temperature), dims=(1, 3)) + return MHT end From 59073a2aa2a14556a335cdd480f3ed8916e53269 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 6 Mar 2026 01:31:47 -0500 Subject: [PATCH 23/63] minor --- src/Diagnostics/meridional_heat_transport.jl | 22 +++++++++++--------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 56e80a832..29dc4f533 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -105,7 +105,9 @@ end function meridional_heat_transport(esm::EarthSystemModel, ::MeridionalHeatFluxMethod; reference_temperature=0) - return meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature = eltype(esm)(reference_temperature)) + FT = eltype(esm) + reference_temperature = convert(reference_temperature, FT) + return meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature) end meridional_heat_transport(esm::EarthSystemModel) = meridional_heat_transport(esm, OceanHeatContentTendencyMethod()) @@ -113,6 +115,15 @@ meridional_heat_transport(esm::EarthSystemModel) = meridional_heat_transport(esm meridional_heat_transport(esm::EarthSystemModel, method; reference_temperature=0.0) = throw(ArgumentError("Unknown method $(method); choose either MeridionalHeatFluxMethod() or OceanHeatContentTendencyMethod().")) +function meridional_heat_transport_via_meridional_heat_flux(esm; 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 + function meridional_heat_transport_via_ocean_heat_content(esm) ρᵒᶜ = reference_density(esm.ocean) cᵒᶜ = heat_capacity(esm.ocean) @@ -134,12 +145,3 @@ function temperature_evolution_tendency(timestepper::QuasiAdamsBashforth2TimeSte χ = timestepper.χ return (3/2 + χ) * Gⁿ - (1/2 + χ) * G⁻ end - -function meridional_heat_transport_via_meridional_heat_flux(esm; 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 From 81d1cffce71e4c1ec2ccc6b9ac77600c9b2cc92b Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 6 Mar 2026 03:14:33 -0500 Subject: [PATCH 24/63] reorder --- src/Diagnostics/Diagnostics.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl index d349b862d..7ec08091c 100644 --- a/src/Diagnostics/Diagnostics.jl +++ b/src/Diagnostics/Diagnostics.jl @@ -20,7 +20,7 @@ using NumericalEarth.EarthSystemModels: EarthSystemModel import Oceananigans.Fields: compute! include("mixed_layer_depth.jl") -include("interface_fluxes.jl") include("meridional_heat_transport.jl") +include("interface_fluxes.jl") end # module From cafa9c6185a169a7940f1ba5057948c2362b487c Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 6 Mar 2026 03:14:52 -0500 Subject: [PATCH 25/63] tweaks --- examples/meridional_heat_transport_ecco.jl | 41 +++++++++++++--------- 1 file changed, 24 insertions(+), 17 deletions(-) diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl index 9a7fdcf3a..d5a39fe96 100644 --- a/examples/meridional_heat_transport_ecco.jl +++ b/examples/meridional_heat_transport_ecco.jl @@ -5,6 +5,8 @@ using Dates using Statistics using Printf +using CUDA; CUDA.device!(3) + arch = GPU() Nx = 360 Ny = 180 @@ -45,8 +47,7 @@ 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=2*365days) - +simulation = Simulation(esm; Δt=20minutes, stop_time=5*365days) wall_time = Ref(time_ns()) @@ -75,47 +76,53 @@ function progress(sim) end # And add it as a callback to the simulation. -add_callback!(simulation, progress, TimeInterval(1days)) - -mht_OHC = meridional_heat_transport(esm, OceanHeatContentTendencyMethod()) |> Field -mht_vT = meridional_heat_transport(esm, MeridionalHeatFluxMethod()) |> Field +add_callback!(simulation, progress, IterationInterval(200)) -mht_outputs = (; mht_OHC, mht_vT) +mht_OHC= meridional_heat_transport(esm, OceanHeatContentTendencyMethod()) +mht_OHC = Field(mht_OHC) +mht_vT = Field(meridional_heat_transport(esm, MeridionalHeatFluxMethod())) -ocean.output_writers[:mth] = JLD2Writer(ocean.model, mht_outputs; +ocean.output_writers[:mth] = JLD2Writer(ocean.model, (; mht_OHC, mht_vT); schedule = TimeInterval(3hours), filename = "ocean_one_degree_mht", overwrite_existing = true) - run!(simulation) +## -mht_OCH = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht_OHC"; backend = OnDisk()) +using Oceananigans + +mht_OHC = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht_OHC"; backend = OnDisk()) mht_vT = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht_vT"; backend = OnDisk()) times = mht_OCH.times Nt = length(times) +grid = mht_OHC.grid +Ny = size(mht_OHC.grid, 2) + mht_OCH_mean = deepcopy(mht_OCH[1][1, :, 1]) -mht_vT_mean = deepcopy(mht_vT[1][1, :, 1]) +mht_vT_mean = deepcopy(mht_vT[1][1, :, 1]) for j in 1:Nt - mht_OCH_mean += mht_OCH[j][1, :, 1] + @info "iteration $j out of $Nt" + mht_OHC_mean += mht_OHC[j][1, :, 1] mht_vT_mean += mht_vT[j][1, :, 1] end -@. mht_OCH_mean = mht_OCH_mean / Nt + +@. mht_OHC_mean = mht_OHC_mean / Nt @. mht_vT_mean = mht_vT_mean / Nt using CairoMakie fig = Figure() -ax1 = Axis(fig[1, 1], xlabel="latitude (deg)", ylabel="MHT (PW)") -ax2 = Axis(fig[2, 1], xlabel="latitude (deg)", ylabel="MHT (PW)") +ax = Axis(fig[1, 1], xlabel="latitude (deg)", ylabel="MHT (PW)") φ = φnodes(grid, Face()) -lines!(ax1, φnodes(grid, Center()), diff(mht_OCH_mean[1:Ny+1]) / 1e15, linewidth=4) -lines!(ax2, φ, mht_vT_mean[1:Ny+1] / 1e15, linewidth=4) + +lines!(ax, φ, mht_OHC_mean[1:Ny+1] / 1e15, linewidth=4) +lines!(ax, φ, mht_vT_mean[1:Ny+1] / 1e15, linewidth=4) save("mht.png", fig) From eb920b8ac68cf78347630461d06774ec0e811a37 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 6 Mar 2026 03:18:56 -0500 Subject: [PATCH 26/63] uncomment tests --- test/test_diagnostics.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/test_diagnostics.jl b/test/test_diagnostics.jl index b56087690..ebf523e55 100644 --- a/test/test_diagnostics.jl +++ b/test/test_diagnostics.jl @@ -5,7 +5,6 @@ using Oceananigans.Models: buoyancy_operation using NumericalEarth.Diagnostics: MixedLayerDepthField, MixedLayerDepthOperand using SeawaterPolynomials: TEOS10EquationOfState -#= for arch in test_architectures, dataset in (ECCO4Monthly(),) A = typeof(arch) @info "Testing MixedLayerDepthField with $(typeof(dataset)) on $A" @@ -57,7 +56,7 @@ for arch in test_architectures, dataset in (ECCO4Monthly(),) end end end -=# + for arch in test_architectures A = typeof(arch) @info "Testing InterfaceFluxOutputs on $A" From 84ab3368c563d4f276a33d341a7a03023eff6e69 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 6 Mar 2026 03:37:32 -0500 Subject: [PATCH 27/63] cleanup --- src/Diagnostics/meridional_heat_transport.jl | 29 ++++++++++---------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 29dc4f533..4b1996573 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -98,23 +98,22 @@ Keyword Arguments * `reference_temperature`: The reference temperature (in ᵒC) used for `MeridionalHeatFluxMethod()`; default: 0 ᵒC. """ -function meridional_heat_transport(esm::EarthSystemModel, ::OceanHeatContentTendencyMethod; - reference_temperature=0) - return meridional_heat_transport_via_ocean_heat_content(esm) +function meridional_heat_transport(esm::EarthSystemModel, method=MeridionalHeatFluxMethod(); reference_temperature=0) + + esm.ocean.model.grid.underlying_grid isa OrthogonalSphericalShellGrid && + throw(ArgumentError("meridional_heat_transport diagnostic does not work on OrthogonalSphericalShellGrid at the moment; use LatitudeLongitudeGrid.")) + + if method isa MeridionalHeatFluxMethod + FT = eltype(esm) + reference_temperature = convert(reference_temperature, FT) + return meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature) + elseif method isa OceanHeatContentTendencyMethod + return meridional_heat_transport_via_ocean_heat_content(esm) + else + throw(ArgumentError("Unknown method $(method); choose either MeridionalHeatFluxMethod() or OceanHeatContentTendencyMethod().")) + end end -function meridional_heat_transport(esm::EarthSystemModel, ::MeridionalHeatFluxMethod; - reference_temperature=0) - FT = eltype(esm) - reference_temperature = convert(reference_temperature, FT) - return meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature) -end - -meridional_heat_transport(esm::EarthSystemModel) = meridional_heat_transport(esm, OceanHeatContentTendencyMethod()) - -meridional_heat_transport(esm::EarthSystemModel, method; reference_temperature=0.0) = - throw(ArgumentError("Unknown method $(method); choose either MeridionalHeatFluxMethod() or OceanHeatContentTendencyMethod().")) - function meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature) ρᵒᶜ = reference_density(esm.ocean) cᵒᶜ = heat_capacity(esm.ocean) From 5260063dad2ee3da3597c68e3772b6cc3c5f2c21 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 6 Mar 2026 03:49:42 -0500 Subject: [PATCH 28/63] fix convert --- src/Diagnostics/meridional_heat_transport.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 4b1996573..5f1596c32 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -105,7 +105,7 @@ function meridional_heat_transport(esm::EarthSystemModel, method=MeridionalHeatF if method isa MeridionalHeatFluxMethod FT = eltype(esm) - reference_temperature = convert(reference_temperature, FT) + reference_temperature = convert(FT, reference_temperature) return meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature) elseif method isa OceanHeatContentTendencyMethod return meridional_heat_transport_via_ocean_heat_content(esm) From 411e5b33844b3810b95caa431ae9b4eadc2bd101 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 6 Mar 2026 03:51:57 -0500 Subject: [PATCH 29/63] simplify --- examples/meridional_heat_transport_ecco.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl index d5a39fe96..6eb12f0fa 100644 --- a/examples/meridional_heat_transport_ecco.jl +++ b/examples/meridional_heat_transport_ecco.jl @@ -78,8 +78,7 @@ end # And add it as a callback to the simulation. add_callback!(simulation, progress, IterationInterval(200)) -mht_OHC= meridional_heat_transport(esm, OceanHeatContentTendencyMethod()) -mht_OHC = Field(mht_OHC) +mht_OHC= Field(meridional_heat_transport(esm, OceanHeatContentTendencyMethod())) mht_vT = Field(meridional_heat_transport(esm, MeridionalHeatFluxMethod())) ocean.output_writers[:mth] = JLD2Writer(ocean.model, (; mht_OHC, mht_vT); @@ -96,13 +95,13 @@ using Oceananigans mht_OHC = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht_OHC"; backend = OnDisk()) mht_vT = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht_vT"; backend = OnDisk()) -times = mht_OCH.times +times = mht_OHC.times Nt = length(times) grid = mht_OHC.grid Ny = size(mht_OHC.grid, 2) -mht_OCH_mean = deepcopy(mht_OCH[1][1, :, 1]) +mht_OHC_mean = deepcopy(mht_OHC[1][1, :, 1]) mht_vT_mean = deepcopy(mht_vT[1][1, :, 1]) for j in 1:Nt From c646272632ee478841833c094fe136f9130327b3 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 6 Mar 2026 11:30:42 +0100 Subject: [PATCH 30/63] Update earth_system_model.jl --- src/EarthSystemModels/earth_system_model.jl | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/EarthSystemModels/earth_system_model.jl b/src/EarthSystemModels/earth_system_model.jl index f00f0651c..ab0ae2309 100644 --- a/src/EarthSystemModels/earth_system_model.jl +++ b/src/EarthSystemModels/earth_system_model.jl @@ -322,16 +322,12 @@ above_freezing_ocean_temperature!(ocean, grid, ::Nothing) = nothing ##### Checkpointing ##### -checkpoint_auxiliary_state(model) = nothing -restore_auxiliary_state!(model, state) = nothing - function prognostic_state(osm::EarthSystemModel) return (clock = prognostic_state(osm.clock), ocean = prognostic_state(osm.ocean), atmosphere = prognostic_state(osm.atmosphere), sea_ice = prognostic_state(osm.sea_ice), - interfaces = prognostic_state(osm.interfaces), - auxiliary = prognostic_state(checkpoint_auxiliary_state(osm))) + interfaces = prognostic_state(osm.interfaces)) end function restore_prognostic_state!(osm::EarthSystemModel, state) From 8f1340ff9d32bf3c54df525e9439d8688d075057 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 6 Mar 2026 05:33:31 -0500 Subject: [PATCH 31/63] add grid compat warning --- src/Diagnostics/meridional_heat_transport.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 5f1596c32..dd0395835 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -14,6 +14,11 @@ Return the meridional heat transport for the coupled `esm::EarthSystemModel` usi two methods: either by directly computing the meridional heat flux or indirectly using the total ocean heat content and the ocean heat uptake. +!!! warning "Only works on LatitudeLongitudeGrid" + + The `meridional_heat_transport` diagnostic currently is only supported only on + `LongitudeLatitudeGrid`s. + Arguments ========= From 84e1a99fa034db876ad14a9b5f2d8a1e531bf401 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 9 Mar 2026 10:24:59 +0200 Subject: [PATCH 32/63] simpler mean over fts --- examples/meridional_heat_transport_ecco.jl | 14 ++------------ 1 file changed, 2 insertions(+), 12 deletions(-) diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl index 6eb12f0fa..3e2d02420 100644 --- a/examples/meridional_heat_transport_ecco.jl +++ b/examples/meridional_heat_transport_ecco.jl @@ -101,18 +101,8 @@ Nt = length(times) grid = mht_OHC.grid Ny = size(mht_OHC.grid, 2) -mht_OHC_mean = deepcopy(mht_OHC[1][1, :, 1]) -mht_vT_mean = deepcopy(mht_vT[1][1, :, 1]) - -for j in 1:Nt - @info "iteration $j out of $Nt" - mht_OHC_mean += mht_OHC[j][1, :, 1] - mht_vT_mean += mht_vT[j][1, :, 1] -end - -@. mht_OHC_mean = mht_OHC_mean / Nt -@. mht_vT_mean = mht_vT_mean / Nt - +mht_OHC_mean = first(mean(mht_OHC, dims=4)) +mht_vT_mean = first(mean(mht_vT, dims=4)) using CairoMakie From 17cdf75ca4c253442db5611c6d03ceacd7475a4b Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 9 Mar 2026 10:28:14 +0200 Subject: [PATCH 33/63] simpler progress --- examples/meridional_heat_transport_ecco.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl index 3e2d02420..b1ba85a7f 100644 --- a/examples/meridional_heat_transport_ecco.jl +++ b/examples/meridional_heat_transport_ecco.jl @@ -64,11 +64,10 @@ function progress(sim) msg1 = @sprintf("time: %s, iter: %d", prettytime(sim), iteration(sim)) msg2 = @sprintf(", max|uo|: (%.1e, %.1e, %.1e) m s⁻¹", umax...) - msg3 = @sprintf(", extrema(To): (%.1f, %.1f) ᵒC, mean(To(z=0)): %.1f ᵒC", Tmin, Tmax, Tavg) - msg4 = @sprintf(", max(e): %.2f m² s⁻²", emax) - msg5 = @sprintf(", wall time: %s \n", prettytime(step_time)) + msg3 = @sprintf(", max(e): %.2f m² s⁻²", emax) + msg4 = @sprintf(", wall time: %s \n", prettytime(step_time)) - @info msg1 * msg2 * msg3 * msg4 * msg5 + @info msg1 * msg2 * msg3 * msg4 wall_time[] = time_ns() From 856301959195f9961daa996a8dea2f9bad1844db Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 9 Mar 2026 10:32:22 +0200 Subject: [PATCH 34/63] fix latex render --- src/Diagnostics/meridional_heat_transport.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index dd0395835..3f3a5902b 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -42,8 +42,8 @@ Arguments 2. For `OceanHeatContentTendencyMethod()` we have: - Let ``T` be three-dimensional (potential) temperature, `ρᵒᶜ` a reference density, - `cᵒᶜ` a heat capacity, `H` the resting depth, and `η` the free-surface elevation. + Let ``T`` be three-dimensional (potential) temperature, ``ρᵒᶜ`` a reference density, + ``cᵒᶜ`` the heat capacity, ``H`` the resting depth, and ``η`` the free-surface elevation. The column heat content per unit horizontal area (units of J m⁻²) is: @@ -87,7 +87,7 @@ Arguments The northward meridional heat transport across latitude ``φ`` is ```math - \\mathrm{MHT}(φ, t) ≡ - ∮_{\\mathrm{lat}=φ} \\boldsymbol{F}_h \\cdot \\hat{\boldsymbol{n}} \\, \\mathrm{d}ℓ + \\mathrm{MHT}(φ, t) ≡ - ∮_{\\mathrm{lat}=φ} \\boldsymbol{F}_h \\cdot \\hat{\\boldsymbol{n}} \\, \\mathrm{d}ℓ ``` with the sign convention that ``\\mathrm{MHT} > 0`` is northward. @@ -95,7 +95,7 @@ Arguments Ignoring the residual ``ℛ``, the OHC-based diagnostic relation is ```math - \\mathrm{MHT} = ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A - \\frac{d}{dt} \\, \\mathrm{OHC}_S + \\mathrm{MHT} = ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A - \\frac{\\mathrm{d}}{\\mathrm{d}t} \\, \\mathrm{OHC}_S ``` Keyword Arguments From 59f1bdbe961ed8c2c980fb3e1496997618619bce Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 9 Mar 2026 10:39:39 +0200 Subject: [PATCH 35/63] fix signs --- examples/one_degree_simulation.jl | 0 src/Diagnostics/meridional_heat_transport.jl | 12 ++++++------ 2 files changed, 6 insertions(+), 6 deletions(-) mode change 100644 => 100755 examples/one_degree_simulation.jl diff --git a/examples/one_degree_simulation.jl b/examples/one_degree_simulation.jl old mode 100644 new mode 100755 diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 3f3a5902b..d2733c4b1 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -61,7 +61,7 @@ Arguments * ``\\boldsymbol{F}_h`` is the depth-integrated horizontal heat flux vector (units W m⁻¹), that includes advection and any parameterized lateral/eddy fluxes, - * ``- 𝒬_{\\rm net}`` is the net surface heat flux into the ocean at the surface + * ``𝒬_{\\rm net}`` is the [net ocean surface heat flux](@ref NumericalEarth.net_ocean_heat_flux) (units W m⁻²), and * ``ℛ`` is the residual sources/sinks and non-closed terms (e.g. numerics, unaccounted physics, mass/volume effects). @@ -72,22 +72,22 @@ Arguments \\mathrm{OHC}_S(φ, t) ≡ ∫_{A(φ)} ℋ \\, \\mathrm{d}A ``` - where `A(φ)` is the ocean area South of latitude ``φ``. + where ``A(φ)`` is the ocean area South of latitude ``φ``. - Integrating the vertically-integrated budget over `A(φ)` and using the divergence + Integrating the vertically-integrated budget over ``A(φ)`` and using the divergence theorem we get ```math \\frac{d}{dt} \\, \\mathrm{OHC}_S(φ, t) = - ∮_{∂A(φ)} \\boldsymbol{F}_h \\cdot \\hat{\\boldsymbol{n}} \\, \\mathrm{d}ℓ - + ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A + - ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A + ∫_{A(φ)} ℛ \\, \\mathrm{d}A ``` The northward meridional heat transport across latitude ``φ`` is ```math - \\mathrm{MHT}(φ, t) ≡ - ∮_{\\mathrm{lat}=φ} \\boldsymbol{F}_h \\cdot \\hat{\\boldsymbol{n}} \\, \\mathrm{d}ℓ + \\mathrm{MHT}(φ, t) ≡ ∮_{\\mathrm{lat}=φ} \\boldsymbol{F}_h \\cdot \\hat{\\boldsymbol{n}} \\, \\mathrm{d}ℓ ``` with the sign convention that ``\\mathrm{MHT} > 0`` is northward. @@ -95,7 +95,7 @@ Arguments Ignoring the residual ``ℛ``, the OHC-based diagnostic relation is ```math - \\mathrm{MHT} = ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A - \\frac{\\mathrm{d}}{\\mathrm{d}t} \\, \\mathrm{OHC}_S + \\mathrm{MHT} = - ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A - \\frac{\\mathrm{d}}{\\mathrm{d}t} \\, \\mathrm{OHC}_S ``` Keyword Arguments From 576d0bc5236e8dc8870a451cb208e8231864df28 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 9 Mar 2026 04:58:40 -0400 Subject: [PATCH 36/63] mean(fts, dims=4) does not work with immersed boundaries --- examples/meridional_heat_transport_ecco.jl | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl index b1ba85a7f..260b66e7b 100644 --- a/examples/meridional_heat_transport_ecco.jl +++ b/examples/meridional_heat_transport_ecco.jl @@ -100,8 +100,17 @@ Nt = length(times) grid = mht_OHC.grid Ny = size(mht_OHC.grid, 2) -mht_OHC_mean = first(mean(mht_OHC, dims=4)) -mht_vT_mean = first(mean(mht_vT, dims=4)) +mht_OHC_mean = deepcopy(mht_OHC[1][1, :, 1]) +mht_vT_mean = deepcopy(mht_vT[1][1, :, 1]) + +for iter in 1:Nt + @info "iteration $j out of $Nt" + mht_OHC_mean += mht_OHC[iter][1, :, 1] + mht_vT_mean += mht_vT[iter][1, :, 1] +end + +@. mht_OHC_mean = mht_OHC_mean / Nt +@. mht_vT_mean = mht_vT_mean / Nt using CairoMakie From 9b18f261622177a4e03b63d23e3d559148c0bd1e Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 9 Mar 2026 05:02:08 -0400 Subject: [PATCH 37/63] fix info msg --- examples/meridional_heat_transport_ecco.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl index 260b66e7b..333536078 100644 --- a/examples/meridional_heat_transport_ecco.jl +++ b/examples/meridional_heat_transport_ecco.jl @@ -104,7 +104,7 @@ mht_OHC_mean = deepcopy(mht_OHC[1][1, :, 1]) mht_vT_mean = deepcopy(mht_vT[1][1, :, 1]) for iter in 1:Nt - @info "iteration $j out of $Nt" + @info "iteration $iter out of $Nt" mht_OHC_mean += mht_OHC[iter][1, :, 1] mht_vT_mean += mht_vT[iter][1, :, 1] end From 34332ed42a0733b483d75eb7486eedd9d33e9acb Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 10 Mar 2026 06:41:05 +0200 Subject: [PATCH 38/63] fix latex rendering --- src/Diagnostics/meridional_heat_transport.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index d2733c4b1..b846623b3 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -78,7 +78,7 @@ Arguments theorem we get ```math - \\frac{d}{dt} \\, \\mathrm{OHC}_S(φ, t) = + \\frac{\\mathrm{d}}{\\mathrm{d}t} \\, \\mathrm{OHC}_S(φ, t) = - ∮_{∂A(φ)} \\boldsymbol{F}_h \\cdot \\hat{\\boldsymbol{n}} \\, \\mathrm{d}ℓ - ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A + ∫_{A(φ)} ℛ \\, \\mathrm{d}A From d1b8174b3b73fe0f28077a52943e524c5ff3033f Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 10 Mar 2026 06:42:37 +0200 Subject: [PATCH 39/63] clarify G^T --- src/Diagnostics/meridional_heat_transport.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index b846623b3..3235c6ad1 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -144,8 +144,8 @@ end temperature_evolution_tendency(timestepper::SplitRungeKuttaTimeStepper) = timestepper.Gⁿ.T function temperature_evolution_tendency(timestepper::QuasiAdamsBashforth2TimeStepper) - Gⁿ = timestepper.Gⁿ.T - G⁻ = timestepper.G⁻.T - χ = timestepper.χ - return (3/2 + χ) * Gⁿ - (1/2 + χ) * G⁻ + Gᵀⁿ = timestepper.Gⁿ.T + Gᵀ⁻ = timestepper.G⁻.T + χ = timestepper.χ + return (3/2 + χ) * Gᵀⁿ - (1/2 + χ) * Gᵀ⁻ end From 2d741856f6353f521fbab82ea9d8b49d8c41b9c2 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 10 Mar 2026 06:43:03 +0200 Subject: [PATCH 40/63] remove stray time --- src/Diagnostics/meridional_heat_transport.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 3235c6ad1..edd3ffeae 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -137,7 +137,6 @@ function meridional_heat_transport_via_ocean_heat_content(esm) ∂ℋ_∂t = Integral(ρᵒᶜ * cᵒᶜ * ∂T_∂t, dims=3) |> Field ∫sum_dx = Integral(heat_flux + ∂ℋ_∂t, dims=1) |> Field MHT = CumulativeIntegral(- ∫sum_dx, dims=2) - return MHT end From dda8eb22830f5ac9b1d41994542bce3219775cb6 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 10 Mar 2026 06:45:54 +0200 Subject: [PATCH 41/63] fix docs ref --- src/Diagnostics/meridional_heat_transport.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index edd3ffeae..b7d24aeed 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -61,7 +61,7 @@ Arguments * ``\\boldsymbol{F}_h`` is the depth-integrated horizontal heat flux vector (units W m⁻¹), that includes advection and any parameterized lateral/eddy fluxes, - * ``𝒬_{\\rm net}`` is the [net ocean surface heat flux](@ref NumericalEarth.net_ocean_heat_flux) + * ``𝒬_{\\rm net}`` is the [net ocean surface heat flux](@ref NumericalEarth.Diagnostics.net_ocean_heat_flux) (units W m⁻²), and * ``ℛ`` is the residual sources/sinks and non-closed terms (e.g. numerics, unaccounted physics, mass/volume effects). From 6e4794bcc32cd5456d4503af683fc80c0357be16 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 10 Mar 2026 06:46:37 +0200 Subject: [PATCH 42/63] align code --- src/Diagnostics/meridional_heat_transport.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index b7d24aeed..d2bf372c7 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -78,10 +78,10 @@ Arguments theorem we get ```math - \\frac{\\mathrm{d}}{\\mathrm{d}t} \\, \\mathrm{OHC}_S(φ, t) = - - ∮_{∂A(φ)} \\boldsymbol{F}_h \\cdot \\hat{\\boldsymbol{n}} \\, \\mathrm{d}ℓ - - ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A - + ∫_{A(φ)} ℛ \\, \\mathrm{d}A + \\frac{\\mathrm{d}}{\\mathrm{d}t} \\, \\mathrm{OHC}_S(φ, t) = + - ∮_{∂A(φ)} \\boldsymbol{F}_h \\cdot \\hat{\\boldsymbol{n}} \\, \\mathrm{d}ℓ + - ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A + + ∫_{A(φ)} ℛ \\, \\mathrm{d}A ``` The northward meridional heat transport across latitude ``φ`` is From b7b924a6c80c54df225b57e1d6dea0e1db3fadcf Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 10 Mar 2026 09:07:17 +0200 Subject: [PATCH 43/63] temperature_evolution_tendency -> temperature_tendency --- src/Diagnostics/meridional_heat_transport.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index d2bf372c7..77940ceba 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -131,7 +131,7 @@ end function meridional_heat_transport_via_ocean_heat_content(esm) ρᵒᶜ = reference_density(esm.ocean) cᵒᶜ = heat_capacity(esm.ocean) - ∂T_∂t = temperature_evolution_tendency(esm.ocean.model.timestepper) + ∂T_∂t = temperature_tendency(esm.ocean.model.timestepper) heat_flux = net_ocean_heat_flux(esm) |> Field ∂ℋ_∂t = Integral(ρᵒᶜ * cᵒᶜ * ∂T_∂t, dims=3) |> Field @@ -140,9 +140,9 @@ function meridional_heat_transport_via_ocean_heat_content(esm) return MHT end -temperature_evolution_tendency(timestepper::SplitRungeKuttaTimeStepper) = timestepper.Gⁿ.T +temperature_tendency(timestepper::SplitRungeKuttaTimeStepper) = timestepper.Gⁿ.T -function temperature_evolution_tendency(timestepper::QuasiAdamsBashforth2TimeStepper) +function temperature_tendency(timestepper::QuasiAdamsBashforth2TimeStepper) Gᵀⁿ = timestepper.Gⁿ.T Gᵀ⁻ = timestepper.G⁻.T χ = timestepper.χ From d505092ac7017e5b9b60a0de04e5721157ae2137 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 10 Mar 2026 09:12:38 +0200 Subject: [PATCH 44/63] clarify docstring --- src/Diagnostics/meridional_heat_transport.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 77940ceba..a7494ad0b 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -84,13 +84,15 @@ Arguments + ∫_{A(φ)} ℛ \\, \\mathrm{d}A ``` + where ``∂A(φ)`` is the boundary of ``A(φ)``, i.e., the circle at latitude ``φ``. + The northward meridional heat transport across latitude ``φ`` is ```math \\mathrm{MHT}(φ, t) ≡ ∮_{\\mathrm{lat}=φ} \\boldsymbol{F}_h \\cdot \\hat{\\boldsymbol{n}} \\, \\mathrm{d}ℓ ``` - with the sign convention that ``\\mathrm{MHT} > 0`` is northward. + with the understanding that ``\\mathrm{MHT} > 0`` means northward heat transport. Ignoring the residual ``ℛ``, the OHC-based diagnostic relation is From ad87f027bba2015a63668c6727a8483646b6fed5 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 10 Mar 2026 15:34:54 +0200 Subject: [PATCH 45/63] cleanup + example --- src/Diagnostics/meridional_heat_transport.jl | 51 +++++++++++++++++--- 1 file changed, 43 insertions(+), 8 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index a7494ad0b..d2bf87756 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -1,8 +1,6 @@ using Oceananigans.TimeSteppers: SplitRungeKuttaTimeStepper, QuasiAdamsBashforth2TimeStepper using ..EarthSystemModels: EarthSystemModel, reference_density, heat_capacity -import ..EarthSystemModels: checkpoint_auxiliary_state, restore_auxiliary_state! - struct OceanHeatContentTendencyMethod end struct MeridionalHeatFluxMethod end @@ -38,12 +36,13 @@ Arguments ``` Above, ``T_{\\rm ref}`` is a reference temperature and ``ρᵒᶜ`` and ``cᵒᶜ`` are the - ocean reference density and heat capacity respectively. + ocean reference density and specific heat capacity respectively. 2. For `OceanHeatContentTendencyMethod()` we have: - Let ``T`` be three-dimensional (potential) temperature, ``ρᵒᶜ`` a reference density, - ``cᵒᶜ`` the heat capacity, ``H`` the resting depth, and ``η`` the free-surface elevation. + Let ``T`` be three-dimensional (potential) temperature, ``ρᵒᶜ`` the ocean reference + density, ``cᵒᶜ`` the specific heat capacity, ``H`` the resting depth, and ``η`` the + free-surface elevation. The column heat content per unit horizontal area (units of J m⁻²) is: @@ -92,22 +91,58 @@ Arguments \\mathrm{MHT}(φ, t) ≡ ∮_{\\mathrm{lat}=φ} \\boldsymbol{F}_h \\cdot \\hat{\\boldsymbol{n}} \\, \\mathrm{d}ℓ ``` - with the understanding that ``\\mathrm{MHT} > 0`` means northward heat transport. + with the understanding that ``\\mathrm{MHT} > 0`` implies northward heat transport. Ignoring the residual ``ℛ``, the OHC-based diagnostic relation is ```math - \\mathrm{MHT} = - ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A - \\frac{\\mathrm{d}}{\\mathrm{d}t} \\, \\mathrm{OHC}_S + \\mathrm{MHT} = - ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A + - \\frac{\\mathrm{d}}{\\mathrm{d}t} \\, \\mathrm{OHC}_S ``` Keyword Arguments ================= * `reference_temperature`: The reference temperature (in ᵒC) used for `MeridionalHeatFluxMethod()`; default: 0 ᵒC. + +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, method=MeridionalHeatFluxMethod(); reference_temperature=0) - esm.ocean.model.grid.underlying_grid isa OrthogonalSphericalShellGrid && + 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.")) if method isa MeridionalHeatFluxMethod From 20dea9a3980001350a20d0db35ff4d4cdffd0ae1 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 10 Mar 2026 15:56:23 +0200 Subject: [PATCH 46/63] breathing space --- examples/meridional_heat_transport_ecco.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl index 333536078..48ef04221 100644 --- a/examples/meridional_heat_transport_ecco.jl +++ b/examples/meridional_heat_transport_ecco.jl @@ -77,10 +77,10 @@ end # And add it as a callback to the simulation. add_callback!(simulation, progress, IterationInterval(200)) -mht_OHC= Field(meridional_heat_transport(esm, OceanHeatContentTendencyMethod())) +mht_OHC = Field(meridional_heat_transport(esm, OceanHeatContentTendencyMethod())) mht_vT = Field(meridional_heat_transport(esm, MeridionalHeatFluxMethod())) -ocean.output_writers[:mth] = JLD2Writer(ocean.model, (; mht_OHC, mht_vT); +ocean.output_writers[:mth] = JLD2Writer(ocean.model, (; mht_vT, mht_OHC); schedule = TimeInterval(3hours), filename = "ocean_one_degree_mht", overwrite_existing = true) From 61bded3c498aca07ed20d816d65e2fa070bbee73 Mon Sep 17 00:00:00 2001 From: Taimoor Sohail Date: Wed, 11 Mar 2026 15:31:37 +1100 Subject: [PATCH 47/63] Update documentation for reference_temperature argument Clarify the relevance of reference temperature in meridional heat transport calculations. --- src/Diagnostics/meridional_heat_transport.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index d2bf87756..d8b762112 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -103,7 +103,9 @@ Arguments Keyword Arguments ================= -* `reference_temperature`: The reference temperature (in ᵒC) used for `MeridionalHeatFluxMethod()`; default: 0 ᵒC. +* `reference_temperature`: The reference temperature (in ᵒC) used for `MeridionalHeatFluxMethod()`; default: 0 ᵒC. +Note that reference temperature is only relevant where there is a net volume transport through the section, which for the case of global latitude circles there is not. +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 ======= From 032c9c9d5d6525a459a6a183d27b43ad6ce5d425 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Wed, 11 Mar 2026 07:26:51 +0200 Subject: [PATCH 48/63] add admonition --- src/Diagnostics/meridional_heat_transport.jl | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index d8b762112..068ba9040 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -103,9 +103,14 @@ Arguments Keyword Arguments ================= -* `reference_temperature`: The reference temperature (in ᵒC) used for `MeridionalHeatFluxMethod()`; default: 0 ᵒC. -Note that reference temperature is only relevant where there is a net volume transport through the section, which for the case of global latitude circles there is not. -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. +* `reference_temperature`: The reference temperature (in ᵒC) used for `MeridionalHeatFluxMethod()`; default: 0 ᵒC. + + !!! info "Repeat-year forcing" + + The reference temperature is only relevant where there is a net volume transport through the section, + which for the case of global latitude circles there is not. 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 ======= From 88635816bc6ae72f478ec120e6b3308c110bafa8 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 13 Mar 2026 11:03:02 +0200 Subject: [PATCH 49/63] notation --- src/Diagnostics/meridional_heat_transport.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 068ba9040..19da5804b 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -175,11 +175,11 @@ end function meridional_heat_transport_via_ocean_heat_content(esm) ρᵒᶜ = reference_density(esm.ocean) cᵒᶜ = heat_capacity(esm.ocean) - ∂T_∂t = temperature_tendency(esm.ocean.model.timestepper) + ∂t_T = temperature_tendency(esm.ocean.model.timestepper) heat_flux = net_ocean_heat_flux(esm) |> Field - ∂ℋ_∂t = Integral(ρᵒᶜ * cᵒᶜ * ∂T_∂t, dims=3) |> Field - ∫sum_dx = Integral(heat_flux + ∂ℋ_∂t, dims=1) |> Field + ∂t_ℋ = Integral(ρᵒᶜ * cᵒᶜ * ∂t_T, dims=3) |> Field + ∫sum_dx = Integral(heat_flux + ∂t_ℋ, dims=1) |> Field MHT = CumulativeIntegral(- ∫sum_dx, dims=2) return MHT end From ba466fd2366f3bfbb5e1f94563da44b5cff857a8 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 13 Mar 2026 11:28:33 +0200 Subject: [PATCH 50/63] notation --- src/Diagnostics/meridional_heat_transport.jl | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 19da5804b..6af7dcfb4 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -50,17 +50,19 @@ Arguments ℋ = ρᵒᶜ cᵒᶜ ∫_{-H}^{η} T \\, \\mathrm{d}z ``` + and its time-derivative is a heat flux, ``𝒬 = ∂ℋ/∂t``. + The vertically-integrated heat budget is ```math - \\frac{∂ℋ}{∂t} = - \\boldsymbol{∇}_h \\cdot \\boldsymbol{F}_h - 𝒬_{\\rm net} + ℛ + 𝒬 = \\frac{∂ℋ}{∂t} = - \\boldsymbol{∇}_h \\cdot \\boldsymbol{F}_h - 𝒬ᵃᵒ_{\\rm net} + ℛ ``` where * ``\\boldsymbol{F}_h`` is the depth-integrated horizontal heat flux vector (units W m⁻¹), that includes advection and any parameterized lateral/eddy fluxes, - * ``𝒬_{\\rm net}`` is the [net ocean surface heat flux](@ref NumericalEarth.Diagnostics.net_ocean_heat_flux) + * ``𝒬ᵃᵒ_{\\rm net}`` is the [net ocean surface heat flux](@ref NumericalEarth.Diagnostics.net_ocean_heat_flux) (units W m⁻²), and * ``ℛ`` is the residual sources/sinks and non-closed terms (e.g. numerics, unaccounted physics, mass/volume effects). @@ -68,7 +70,7 @@ Arguments The total ocean heat content (OHC) South of latitude ``φ`` is: ```math - \\mathrm{OHC}_S(φ, t) ≡ ∫_{A(φ)} ℋ \\, \\mathrm{d}A + \\mathrm{OHC}_S(φ, t) ≡ ∫_{A(φ)} 𝒬 \\, \\mathrm{d}A ``` where ``A(φ)`` is the ocean area South of latitude ``φ``. @@ -79,7 +81,7 @@ Arguments ```math \\frac{\\mathrm{d}}{\\mathrm{d}t} \\, \\mathrm{OHC}_S(φ, t) = - ∮_{∂A(φ)} \\boldsymbol{F}_h \\cdot \\hat{\\boldsymbol{n}} \\, \\mathrm{d}ℓ - - ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A + - ∫_{A(φ)} 𝒬ᵃᵒ_{\\rm net} \\, \\mathrm{d}A + ∫_{A(φ)} ℛ \\, \\mathrm{d}A ``` @@ -96,7 +98,7 @@ Arguments Ignoring the residual ``ℛ``, the OHC-based diagnostic relation is ```math - \\mathrm{MHT} = - ∫_{A(φ)} 𝒬_{\\rm net} \\, \\mathrm{d}A + \\mathrm{MHT} = - ∫_{A(φ)} 𝒬ᵃᵒ_{\\rm net} \\, \\mathrm{d}A - \\frac{\\mathrm{d}}{\\mathrm{d}t} \\, \\mathrm{OHC}_S ``` @@ -176,11 +178,11 @@ function meridional_heat_transport_via_ocean_heat_content(esm) ρᵒᶜ = reference_density(esm.ocean) cᵒᶜ = heat_capacity(esm.ocean) ∂t_T = temperature_tendency(esm.ocean.model.timestepper) - heat_flux = net_ocean_heat_flux(esm) |> Field + 𝒬ᵃᵒₙₑₜ = net_ocean_heat_flux(esm) |> Field - ∂t_ℋ = Integral(ρᵒᶜ * cᵒᶜ * ∂t_T, dims=3) |> Field - ∫sum_dx = Integral(heat_flux + ∂t_ℋ, dims=1) |> Field - MHT = CumulativeIntegral(- ∫sum_dx, dims=2) + 𝒬 = Integral(ρᵒᶜ * cᵒᶜ * ∂t_T, dims=3) |> Field + ∫Σ𝒬_dx = Integral(𝒬ᵃᵒₙₑₜ + 𝒬, dims=1) |> Field + MHT = CumulativeIntegral(- ∫Σ𝒬_dx, dims=2) return MHT end From 545be9b29ab3364b24927f77f34d725fa400f69d Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 13 Mar 2026 14:07:11 +0200 Subject: [PATCH 51/63] enhance reference temperature admonition --- src/Diagnostics/meridional_heat_transport.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 6af7dcfb4..a2cc19c19 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -107,12 +107,14 @@ Keyword Arguments * `reference_temperature`: The reference temperature (in ᵒC) used for `MeridionalHeatFluxMethod()`; default: 0 ᵒC. - !!! info "Repeat-year forcing" - - The reference temperature is only relevant where there is a net volume transport through the section, - which for the case of global latitude circles there is not. 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. + !!! 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 ======= From d4e1d1a3e1cb714ba4216dc94630df1f61867e9d Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 13 Mar 2026 14:08:40 +0200 Subject: [PATCH 52/63] remove stray time --- src/Diagnostics/meridional_heat_transport.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index a2cc19c19..32afebceb 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -75,7 +75,7 @@ Arguments where ``A(φ)`` is the ocean area South of latitude ``φ``. - Integrating the vertically-integrated budget over ``A(φ)`` and using the divergence + Integrating the vertically-integrated heat budget over ``A(φ)`` and using the divergence theorem we get ```math From cadcdd8df8e4c18667373e0aa3646718cfc6860a Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 13 Mar 2026 14:14:13 +0200 Subject: [PATCH 53/63] connect MHT to scrQ --- src/Diagnostics/meridional_heat_transport.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 32afebceb..e685c52d9 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -70,7 +70,7 @@ Arguments The total ocean heat content (OHC) South of latitude ``φ`` is: ```math - \\mathrm{OHC}_S(φ, t) ≡ ∫_{A(φ)} 𝒬 \\, \\mathrm{d}A + \\mathrm{OHC}_S(φ, t) ≡ ∫_{A(φ)} ℋ \\, \\mathrm{d}A ``` where ``A(φ)`` is the ocean area South of latitude ``φ``. @@ -98,8 +98,11 @@ Arguments Ignoring the residual ``ℛ``, the OHC-based diagnostic relation is ```math - \\mathrm{MHT} = - ∫_{A(φ)} 𝒬ᵃᵒ_{\\rm net} \\, \\mathrm{d}A - - \\frac{\\mathrm{d}}{\\mathrm{d}t} \\, \\mathrm{OHC}_S + \\begin{align*} + \\mathrm{MHT} & = - ∫_{A(φ)} 𝒬ᵃᵒ_{\\rm net} \\, \\mathrm{d}A + - \\frac{\\mathrm{d}}{\\mathrm{d}t} \\, \\mathrm{OHC}_S \\\\ + & = - ∫_{A(φ)} \\left( 𝒬ᵃᵒ_{\\rm net} + 𝒬 \\right) \\, \\mathrm{d}A + \\end{align*} ``` Keyword Arguments From bde7d0804c2fb97155ed995403a5b63e8c2721a5 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 13 Mar 2026 14:15:21 +0200 Subject: [PATCH 54/63] define --- src/Diagnostics/meridional_heat_transport.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index e685c52d9..1eade1243 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -50,7 +50,7 @@ Arguments ℋ = ρᵒᶜ cᵒᶜ ∫_{-H}^{η} T \\, \\mathrm{d}z ``` - and its time-derivative is a heat flux, ``𝒬 = ∂ℋ/∂t``. + and its time-derivative is a heat flux, ``𝒬 ≡ ∂ℋ/∂t``. The vertically-integrated heat budget is From 46efae2a8eb68cf7c517e7ceaa8ff6158db020fb Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 13 Mar 2026 13:29:25 -0400 Subject: [PATCH 55/63] enhance mht fig --- examples/meridional_heat_transport_ecco.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl index 48ef04221..7c49725d8 100644 --- a/examples/meridional_heat_transport_ecco.jl +++ b/examples/meridional_heat_transport_ecco.jl @@ -119,7 +119,9 @@ ax = Axis(fig[1, 1], xlabel="latitude (deg)", ylabel="MHT (PW)") φ = φnodes(grid, Face()) -lines!(ax, φ, mht_OHC_mean[1:Ny+1] / 1e15, linewidth=4) -lines!(ax, φ, mht_vT_mean[1:Ny+1] / 1e15, linewidth=4) +lines!(ax, φ, mht_OHC_mean[1:Ny+1] / 1e15, linewidth=4, label="via OHC") +lines!(ax, φ, mht_vT_mean[1:Ny+1] / 1e15, linewidth=4, label="via vT") +Legend(fig[2, :], ax, orientation=:horizontal) +Label(fig[0, :], "Meridional heat transport", fontsize=16, tellwidth=false) save("mht.png", fig) From 45eb5c61af434752a9aefc9a590b98cfc8f139f9 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 30 Mar 2026 15:27:02 +1100 Subject: [PATCH 56/63] fix filename --- test/{tet_diagnostics_2.jl => test_diagnostics_2.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename test/{tet_diagnostics_2.jl => test_diagnostics_2.jl} (100%) diff --git a/test/tet_diagnostics_2.jl b/test/test_diagnostics_2.jl similarity index 100% rename from test/tet_diagnostics_2.jl rename to test/test_diagnostics_2.jl From debc4040f364fe83b7d1bf9fd9c7777922eec68f Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 30 Mar 2026 15:37:14 +1100 Subject: [PATCH 57/63] add budget warning for OceanHeatContentTendencyMethod --- src/Diagnostics/meridional_heat_transport.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 1eade1243..f773ae9b7 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -164,6 +164,11 @@ function meridional_heat_transport(esm::EarthSystemModel, method=MeridionalHeatF reference_temperature = convert(FT, reference_temperature) return meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature) elseif method isa OceanHeatContentTendencyMethod + @warn """ + Computing the meridional heat transport via the OceanHeatContentTendencyMethod + does not ensure the heat budget is closed! + If you require a closed heat budget, then use the MeridionalHeatFluxMethod. + """ return meridional_heat_transport_via_ocean_heat_content(esm) else throw(ArgumentError("Unknown method $(method); choose either MeridionalHeatFluxMethod() or OceanHeatContentTendencyMethod().")) From 368d6a19070e6e5eca1bc9952e62e3a12fa35b6f Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 30 Mar 2026 16:41:26 +1100 Subject: [PATCH 58/63] Update meridional_heat_transport_ecco.jl --- examples/meridional_heat_transport_ecco.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl index 7c49725d8..a896f7e63 100644 --- a/examples/meridional_heat_transport_ecco.jl +++ b/examples/meridional_heat_transport_ecco.jl @@ -15,7 +15,7 @@ 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 = 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, From 898b8115b5064b38e5fbaa5961e75c9d9de67469 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 11 Apr 2026 20:58:22 +1000 Subject: [PATCH 59/63] remove unwanted imports --- src/Diagnostics/Diagnostics.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl index 7ec08091c..87965b309 100644 --- a/src/Diagnostics/Diagnostics.jl +++ b/src/Diagnostics/Diagnostics.jl @@ -11,8 +11,6 @@ using Oceananigans using Oceananigans.Architectures: architecture using Oceananigans.BoundaryConditions: FieldBoundaryConditions, fill_halo_regions! using Oceananigans.Fields: FieldStatus -using Oceananigans.Grids: new_data, inactive_cell, znode -using Oceananigans.Models: buoyancy_operation using Oceananigans.Utils: launch! using KernelAbstractions: @index, @kernel using NumericalEarth.EarthSystemModels: EarthSystemModel From c369138e62897ca9f004f3351b2cfb41dd6bfa1b Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 11 Apr 2026 20:58:59 +1000 Subject: [PATCH 60/63] undo remove unwanted imports --- src/Diagnostics/Diagnostics.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl index 87965b309..7ec08091c 100644 --- a/src/Diagnostics/Diagnostics.jl +++ b/src/Diagnostics/Diagnostics.jl @@ -11,6 +11,8 @@ using Oceananigans using Oceananigans.Architectures: architecture using Oceananigans.BoundaryConditions: FieldBoundaryConditions, fill_halo_regions! using Oceananigans.Fields: FieldStatus +using Oceananigans.Grids: new_data, inactive_cell, znode +using Oceananigans.Models: buoyancy_operation using Oceananigans.Utils: launch! using KernelAbstractions: @index, @kernel using NumericalEarth.EarthSystemModels: EarthSystemModel From 9c608e27e1c41859d6b20d163c39043d87c855ab Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 11 Apr 2026 21:00:10 +1000 Subject: [PATCH 61/63] remove unecessary --- src/EarthSystemModels/earth_system_model.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/EarthSystemModels/earth_system_model.jl b/src/EarthSystemModels/earth_system_model.jl index 448412daf..90469dc04 100644 --- a/src/EarthSystemModels/earth_system_model.jl +++ b/src/EarthSystemModels/earth_system_model.jl @@ -335,8 +335,6 @@ function restore_prognostic_state!(osm::EarthSystemModel, state) restore_prognostic_state!(osm.atmosphere, state.atmosphere) restore_prognostic_state!(osm.sea_ice, state.sea_ice) restore_prognostic_state!(osm.interfaces, state.interfaces) - auxiliary_state = hasproperty(state, :auxiliary) ? state.auxiliary : nothing - restore_auxiliary_state!(osm, auxiliary_state) # Note: we do NOT call update_state! here because: # 1. The checkpoint was saved AFTER update_state! was called at the end of that time step # 2. Calling update_state! would recompute interface fluxes and overwrite restored state From 658923ded4d06d9178c07fbce68baece3fe73bcd Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 12 Apr 2026 10:19:46 +1000 Subject: [PATCH 62/63] cleaner method names --- examples/meridional_heat_transport_ecco.jl | 4 +-- src/Diagnostics/Diagnostics.jl | 2 +- src/Diagnostics/meridional_heat_transport.jl | 32 ++++++++++---------- src/NumericalEarth.jl | 2 +- 4 files changed, 20 insertions(+), 20 deletions(-) diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl index a896f7e63..ee4e60fb8 100644 --- a/examples/meridional_heat_transport_ecco.jl +++ b/examples/meridional_heat_transport_ecco.jl @@ -77,8 +77,8 @@ end # And add it as a callback to the simulation. add_callback!(simulation, progress, IterationInterval(200)) -mht_OHC = Field(meridional_heat_transport(esm, OceanHeatContentTendencyMethod())) -mht_vT = Field(meridional_heat_transport(esm, MeridionalHeatFluxMethod())) +mht_OHC = Field(meridional_heat_transport(esm, TendencyMethod())) +mht_vT = Field(meridional_heat_transport(esm, MeridionalFluxMethod())) ocean.output_writers[:mth] = JLD2Writer(ocean.model, (; mht_vT, mht_OHC); schedule = TimeInterval(3hours), diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl index 7ec08091c..8db2d120b 100644 --- a/src/Diagnostics/Diagnostics.jl +++ b/src/Diagnostics/Diagnostics.jl @@ -1,7 +1,7 @@ module Diagnostics export MixedLayerDepthField, MixedLayerDepthOperand -export meridional_heat_transport, OceanHeatContentTendencyMethod, MeridionalHeatFluxMethod +export meridional_heat_transport, TendencyMethod, MeridionalFluxMethod 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, diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 75ade860b..a4c583fd9 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -1,11 +1,11 @@ using Oceananigans.TimeSteppers: SplitRungeKuttaTimeStepper, QuasiAdamsBashforth2TimeStepper using ..EarthSystemModels: EarthSystemModel, reference_density, heat_capacity -struct OceanHeatContentTendencyMethod end -struct MeridionalHeatFluxMethod end +struct MeridionalFluxMethod end +struct TendencyMethod end """ - meridional_heat_transport(esm::EarthSystemModel, MeridionalHeatFluxMethod(); + meridional_heat_transport(esm::EarthSystemModel, MeridionalFluxMethod(); reference_temperature = 0) Return the meridional heat transport for the coupled `esm::EarthSystemModel` using either @@ -22,14 +22,14 @@ Arguments * `esm`: An EarthSystemModel. -* The method for the computation. Available options are: `MeridionalHeatFluxMethod()` (default) - and `OceanHeatContentTendencyMethod()`. +* The method for the computation. Available options are: `MeridionalFluxMethod()` (default) + and `TendencyMethod()`. - `MeridionalHeatFluxMethod()` computes the meridional heat transport directly by summing - the meridional heat flux; `OceanHeatContentTendencyMethod()` computes the meridional heat + `MeridionalFluxMethod()` computes the meridional heat transport directly by summing + the meridional heat flux; `TendencyMethod()` computes the meridional heat transport indirectly using the ocean heat content. - 1. For `MeridionalHeatFluxMethod()`, the meridional heat transport is computed via: + 1. For `MeridionalFluxMethod()`, the meridional heat transport is computed via: ```math \\mathrm{MHT} ≡ ρᵒᶜ cᵒᶜ ∫ v (T - T_{\\rm ref}) \\, \\mathrm{d}x \\, \\mathrm{d}z @@ -38,7 +38,7 @@ Arguments Above, ``T_{\\rm ref}`` is a reference temperature and ``ρᵒᶜ`` and ``cᵒᶜ`` are the ocean reference density and specific heat capacity respectively. - 2. For `OceanHeatContentTendencyMethod()` we have: + 2. For `TendencyMethod()` we have: Let ``T`` be three-dimensional (potential) temperature, ``ρᵒᶜ`` the ocean reference density, ``cᵒᶜ`` the specific heat capacity, ``H`` the resting depth, and ``η`` the @@ -108,7 +108,7 @@ Arguments Keyword Arguments ================= -* `reference_temperature`: The reference temperature (in ᵒC) used for `MeridionalHeatFluxMethod()`; default: 0 ᵒC. +* `reference_temperature`: The reference temperature (in ᵒC) used for `MeridionalFluxMethod()`; default: 0 ᵒC. !!! info "Reference temperature" @@ -150,7 +150,7 @@ Integral of BinaryOperation at (Center, Face, Center) over dims (1, 3) └── grid: 4×5×2 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×2 halo ``` """ -function meridional_heat_transport(esm::EarthSystemModel, method=MeridionalHeatFluxMethod(); +function meridional_heat_transport(esm::EarthSystemModel, method=MeridionalFluxMethod(); reference_temperature=0) grid = esm.ocean.model.grid @@ -160,19 +160,19 @@ function meridional_heat_transport(esm::EarthSystemModel, method=MeridionalHeatF grid isa OrthogonalSphericalShellGrid && throw(ArgumentError("meridional_heat_transport diagnostic does not work on OrthogonalSphericalShellGrid at the moment; use LatitudeLongitudeGrid.")) - if method isa MeridionalHeatFluxMethod + if method isa MeridionalFluxMethod FT = eltype(esm) reference_temperature = convert(FT, reference_temperature) return meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature) - elseif method isa OceanHeatContentTendencyMethod + elseif method isa TendencyMethod @warn """ - Computing the meridional heat transport via the OceanHeatContentTendencyMethod + Computing the meridional heat transport via the TendencyMethod does not ensure the heat budget is closed! - If you require a closed heat budget, then use the MeridionalHeatFluxMethod. + If you require a closed heat budget, then use the MeridionalFluxMethod. """ return meridional_heat_transport_via_ocean_heat_content(esm) else - throw(ArgumentError("Unknown method $(method); choose either MeridionalHeatFluxMethod() or OceanHeatContentTendencyMethod().")) + throw(ArgumentError("Unknown method $(method); choose either MeridionalFluxMethod() or TendencyMethod().")) end end diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl index 301f990f0..b7e9886cb 100644 --- a/src/NumericalEarth.jl +++ b/src/NumericalEarth.jl @@ -56,7 +56,7 @@ export 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, - meridional_heat_transport, OceanHeatContentTendencyMethod, MeridionalHeatFluxMethod + meridional_heat_transport, TendencyMethod, MeridionalFluxMethod using Oceananigans using Oceananigans.Operators: ℑxyᶠᶜᵃ, ℑxyᶜᶠᵃ From 8db34f8cf63db4cb043cd341ca2a75e0ad2c0974 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 12 Apr 2026 10:21:34 +1000 Subject: [PATCH 63/63] reorder --- examples/meridional_heat_transport_ecco.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl index ee4e60fb8..6570a7652 100644 --- a/examples/meridional_heat_transport_ecco.jl +++ b/examples/meridional_heat_transport_ecco.jl @@ -77,8 +77,8 @@ end # And add it as a callback to the simulation. add_callback!(simulation, progress, IterationInterval(200)) -mht_OHC = Field(meridional_heat_transport(esm, TendencyMethod())) mht_vT = Field(meridional_heat_transport(esm, MeridionalFluxMethod())) +mht_OHC = Field(meridional_heat_transport(esm, TendencyMethod())) ocean.output_writers[:mth] = JLD2Writer(ocean.model, (; mht_vT, mht_OHC); schedule = TimeInterval(3hours), @@ -91,22 +91,22 @@ run!(simulation) using Oceananigans -mht_OHC = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht_OHC"; backend = OnDisk()) mht_vT = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht_vT"; backend = OnDisk()) +mht_OHC = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht_OHC"; backend = OnDisk()) -times = mht_OHC.times +times = mht_vT.times Nt = length(times) -grid = mht_OHC.grid -Ny = size(mht_OHC.grid, 2) +grid = mht_vT.grid +Ny = size(mht_vT.grid, 2) -mht_OHC_mean = deepcopy(mht_OHC[1][1, :, 1]) mht_vT_mean = deepcopy(mht_vT[1][1, :, 1]) +mht_OHC_mean = deepcopy(mht_OHC[1][1, :, 1]) for iter in 1:Nt @info "iteration $iter out of $Nt" - mht_OHC_mean += mht_OHC[iter][1, :, 1] mht_vT_mean += mht_vT[iter][1, :, 1] + mht_OHC_mean += mht_OHC[iter][1, :, 1] end @. mht_OHC_mean = mht_OHC_mean / Nt @@ -119,8 +119,8 @@ ax = Axis(fig[1, 1], xlabel="latitude (deg)", ylabel="MHT (PW)") φ = φnodes(grid, Face()) -lines!(ax, φ, mht_OHC_mean[1:Ny+1] / 1e15, linewidth=4, label="via OHC") lines!(ax, φ, mht_vT_mean[1:Ny+1] / 1e15, linewidth=4, label="via vT") +lines!(ax, φ, mht_OHC_mean[1:Ny+1] / 1e15, linewidth=4, label="via OHC") Legend(fig[2, :], ax, orientation=:horizontal) Label(fig[0, :], "Meridional heat transport", fontsize=16, tellwidth=false)