From 16655a26183cf32258ea342f4c4558e66326a7fd Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 14:04:09 +0100 Subject: [PATCH 01/14] Revert "Revert "add a prescribed ocean"" This reverts commit 71d02b840dedf8378d43429af23a709d7d793e8b. --- .../flux_climatology/flux_climatology.jl | 55 ------------- src/ClimaOcean.jl | 10 +++ src/OceanSeaIceModels/OceanSeaIceModels.jl | 1 - src/OceanSeaIceModels/ocean_sea_ice_model.jl | 19 ----- src/OceanSimulations/OceanSimulations.jl | 12 ++- .../prescribed_ocean_model.jl | 81 +++++++++++++++++++ src/SeaIceSimulations.jl | 5 ++ test/test_ocean_sea_ice_model.jl | 19 +++++ 8 files changed, 126 insertions(+), 76 deletions(-) create mode 100644 src/OceanSimulations/prescribed_ocean_model.jl diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index 94a1a951e..f0f766663 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -115,32 +115,6 @@ end ##### A prescribed ocean... ##### -struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} - architecture :: A - grid :: G - clock :: Clock{C} - velocities :: U - tracers :: T - timeseries :: F -end - -function PrescribedOcean(timeseries; - grid, - clock=Clock{Float64}(time = 0)) - - τx = Field{Face, Center, Nothing}(grid) - τy = Field{Center, Face, Nothing}(grid) - Jᵀ = Field{Center, Center, Nothing}(grid) - Jˢ = Field{Center, Center, Nothing}(grid) - - u = XFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Face, Center, Center), top = FluxBoundaryCondition(τx))) - v = YFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Face, Center), top = FluxBoundaryCondition(τy))) - T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) - S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) - - PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) -end - # ...with prescribed velocity and tracer fields version = ECCO4Monthly() dates = all_ECCO_dates(version)[1:24] @@ -157,35 +131,6 @@ grid = ECCO.ECCO_immersed_grid(arch) ocean_model = PrescribedOcean((; u, v, T, S); grid) ocean = Simulation(ocean_model, Δt=3hour, stop_time=365days) -##### -##### Need to extend a couple of methods -##### - -function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) - tick!(model.clock, Δt) - time = Time(model.clock.time) - - possible_fts = merge(model.velocities, model.tracers) - time_series_tuple = extract_field_time_series(possible_fts) - - for fts in time_series_tuple - update_field_time_series!(fts, time) - end - - parent(model.velocities.u) .= parent(model.timeseries.u[time]) - parent(model.velocities.v) .= parent(model.timeseries.v[time]) - parent(model.tracers.T) .= parent(model.timeseries.T[time]) - parent(model.tracers.S) .= parent(model.timeseries.S[time]) - - return nothing -end - -update_state!(::PrescribedOcean) = nothing -timestepper(::PrescribedOcean) = nothing - -reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6 -heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6 - ##### ##### A prescribed atmosphere... ##### diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl index 2a90d73e8..96a543673 100644 --- a/src/ClimaOcean.jl +++ b/src/ClimaOcean.jl @@ -41,6 +41,16 @@ using DataDeps using Oceananigans.OutputReaders: GPUAdaptedFieldTimeSeries, FieldTimeSeries using Oceananigans.Grids: node +# Does not really matter if there is no model +reference_density(::Nothing) = 0 +heat_capacity(::Nothing) = 0 + +reference_density(unsupported) = + throw(ArgumentError("Cannot extract reference density from $(typeof(unsupported))")) + +heat_capacity(unsupported) = + throw(ArgumentError("Cannot deduce the heat capacity from $(typeof(unsupported))")) + const SomeKindOfFieldTimeSeries = Union{FieldTimeSeries, GPUAdaptedFieldTimeSeries} diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl index 8041f497d..4dee01e38 100644 --- a/src/OceanSeaIceModels/OceanSeaIceModels.jl +++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl @@ -22,7 +22,6 @@ using ClimaSeaIce: SeaIceModel using ClimaSeaIce.SeaIceThermodynamics: melting_temperature using ClimaOcean: stateindex - using KernelAbstractions: @kernel, @index using KernelAbstractions.Extras.LoopInfo: @unroll diff --git a/src/OceanSeaIceModels/ocean_sea_ice_model.jl b/src/OceanSeaIceModels/ocean_sea_ice_model.jl index d2344dabb..9c97a6cb5 100644 --- a/src/OceanSeaIceModels/ocean_sea_ice_model.jl +++ b/src/OceanSeaIceModels/ocean_sea_ice_model.jl @@ -73,25 +73,6 @@ function initialize!(model::OSIM) return nothing end -reference_density(unsupported) = - throw(ArgumentError("Cannot extract reference density from $(typeof(unsupported))")) - -heat_capacity(unsupported) = - throw(ArgumentError("Cannot deduce the heat capacity from $(typeof(unsupported))")) - -reference_density(ocean::Simulation) = reference_density(ocean.model.buoyancy.formulation) -reference_density(buoyancy_formulation::SeawaterBuoyancy) = reference_density(buoyancy_formulation.equation_of_state) -reference_density(eos::TEOS10EquationOfState) = eos.reference_density -reference_density(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density - -heat_capacity(ocean::Simulation) = heat_capacity(ocean.model.buoyancy.formulation) -heat_capacity(buoyancy_formulation::SeawaterBuoyancy) = heat_capacity(buoyancy_formulation.equation_of_state) -heat_capacity(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity - -# Does not really matter if there is no model -reference_density(::Nothing) = 0 -heat_capacity(::Nothing) = 0 - function heat_capacity(::TEOS10EquationOfState{FT}) where FT cₚ⁰ = SeawaterPolynomials.TEOS10.teos10_reference_heat_capacity return convert(FT, cₚ⁰) diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 32423c690..49d48aa6e 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -1,6 +1,6 @@ module OceanSimulations -export ocean_simulation +export ocean_simulation, PrescribedOceanModel using Oceananigans using Oceananigans.Units @@ -21,6 +21,15 @@ using Oceananigans.BuoyancyFormulations: g_Earth using Oceananigans.Coriolis: Ω_Earth using Oceananigans.Operators +import ClimaOcean: reference_density, heat_capacity + +reference_density(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = reference_density(ocean.model.buoyancy.formulation) +reference_density(buoyancy_formulation::SeawaterBuoyancy) = reference_density(buoyancy_formulation.equation_of_state) +reference_density(eos::TEOS10EquationOfState) = eos.reference_density + +heat_capacity(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = heat_capacity(ocean.model.buoyancy.formulation) +heat_capacity(buoyancy_formulation::SeawaterBuoyancy) = heat_capacity(buoyancy_formulation.equation_of_state) + struct Default{V} value :: V end @@ -41,5 +50,6 @@ default_or_override(override, alternative_default=nothing) = override include("barotropic_potential_forcing.jl") include("ocean_simulation.jl") +include("prescribed_ocean_model.jl") end # module diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl new file mode 100644 index 000000000..a5309e4a8 --- /dev/null +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -0,0 +1,81 @@ +using Oceananigans.OutputReaders: extract_field_time_series, update_field_time_series! + +import Oceananigans.TimeSteppers: time_step!, update_state!, reset!, tick! +import Oceananigans.Models: timestepper, update_model_field_time_series! + +import ClimaOcean: reference_density, heat_capacity +import Oceananigans.Architectures: on_architecture + +##### +##### A prescribed ocean... +##### + +struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} + architecture :: A + grid :: G + clock :: Clock{C} + velocities :: U + tracers :: T + timeseries :: F +end + +function PrescribedOcean(timeseries; + grid, + clock=Clock{Float64}(time = 0)) + + # Make sure all elements of the timeseries are on the same grid + for (k, v) in timeseries + if !isa(v, FieldTimeSeries) + throw(ArgumentError("All variables in the `timeseries` argument must be `FieldTimeSeries`")) + end + if v.grid != grid + throw(ArgumentError("All variables in the timeseries reside on the provided grid")) + end + end + + τx = Field{Face, Center, Nothing}(grid) + τy = Field{Center, Face, Nothing}(grid) + Jᵀ = Field{Center, Center, Nothing}(grid) + Jˢ = Field{Center, Center, Nothing}(grid) + + u = XFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Face, Center, Center), top = FluxBoundaryCondition(τx))) + v = YFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Face, Center), top = FluxBoundaryCondition(τy))) + T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) + S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) + + PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) +end + +##### +##### Need to extend a couple of methods +##### + +function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) + tick!(model.clock, Δt) + time = Time(model.clock.time) + + possible_fts = merge(model.velocities, model.tracers) + time_series_tuple = extract_field_time_series(possible_fts) + + for fts in time_series_tuple + update_field_time_series!(fts, time) + end + + update_u_velocity = haskey(model.timeseries, :u) + update_v_velocity = haskey(model.timeseries, :v) + update_temperature = haskey(model.timeseries, :T) + update_salinity = haskey(model.timeseries, :S) + + update_u_velocity && iterpolate!(model.velocities.u, model.timeseries.u[time]) + update_v_velocity && iterpolate!(model.velocities.v, model.timeseries.v[time]) + update_temperature && iterpolate!(model.tracers.T, model.timeseries.T[time]) + update_salinity && iterpolate!(model.tracers.S, model.timeseries.S[time]) + + return nothing +end + +update_state!(::PrescribedOcean) = nothing +timestepper(::PrescribedOcean) = nothing + +reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6 +heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6 diff --git a/src/SeaIceSimulations.jl b/src/SeaIceSimulations.jl index 2b7dc17b8..510680651 100644 --- a/src/SeaIceSimulations.jl +++ b/src/SeaIceSimulations.jl @@ -20,6 +20,11 @@ using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium using ClimaOcean.OceanSimulations: Default +import ClimaOcean: reference_density, heat_capacity + +reference_density(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density +heat_capacity(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity + function sea_ice_simulation(grid; Δt = 5minutes, ice_salinity = 0, # psu diff --git a/test/test_ocean_sea_ice_model.jl b/test/test_ocean_sea_ice_model.jl index 55b6ac3f5..16d55fc9d 100644 --- a/test/test_ocean_sea_ice_model.jl +++ b/test/test_ocean_sea_ice_model.jl @@ -100,5 +100,24 @@ using ClimaSeaIce.Rheologies true end end + + @testset "Test a prescribed ocean model" begin + grid = LatitudeLongitudeGrid(size = (10, 10, 10), latitude = (20, 30), longitude = (40, 50), z = (-100, 0)) + T = ECCOFieldTimeSeries(:temperature, grid; time_indices_in_memory = 2) + S = ECCOFieldTimeSeries(:salinity, grid; time_indices_in_memory = 2) + ocean_model = PrescribedOcean((; T, S); grid) + + time_step!(ocean_model, T.times[2]) + + @test ocean_model.clock.time == T.times[2] + @test ocean_model.tracers.T.data == T[2].data + @test ocean_model.tracers.S.data == S[2].data + + time_step!(ocean_model, 10days) + + @test ocean_model.clock.time == 10days + @test T.backend.start == 2 + @test ocean_model.tracers.T.data != T[3].data + end end From 1a00396ddb397a16107817590427bcc3ad0ad5c2 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 14:09:13 +0100 Subject: [PATCH 02/14] better like this --- .../prescribed_ocean_model.jl | 23 ++++++++++--------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl index a5309e4a8..e2bf596fb 100644 --- a/src/OceanSimulations/prescribed_ocean_model.jl +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -19,17 +19,18 @@ struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} timeseries :: F end -function PrescribedOcean(timeseries; +function PrescribedOcean(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) # Make sure all elements of the timeseries are on the same grid - for (k, v) in timeseries - if !isa(v, FieldTimeSeries) - throw(ArgumentError("All variables in the `timeseries` argument must be `FieldTimeSeries`")) - end - if v.grid != grid - throw(ArgumentError("All variables in the timeseries reside on the provided grid")) + # If we decide to pass a timeseries + if !isempty(timeseries) + for (k, v) in timeseries + isa(v, FieldTimeSeries) || + throw(ArgumentError("All variables in the `timeseries` argument must be `FieldTimeSeries`")) + v.grid == grid || + throw(ArgumentError("All variables in the timeseries reside on the provided grid")) end end @@ -66,10 +67,10 @@ function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) update_temperature = haskey(model.timeseries, :T) update_salinity = haskey(model.timeseries, :S) - update_u_velocity && iterpolate!(model.velocities.u, model.timeseries.u[time]) - update_v_velocity && iterpolate!(model.velocities.v, model.timeseries.v[time]) - update_temperature && iterpolate!(model.tracers.T, model.timeseries.T[time]) - update_salinity && iterpolate!(model.tracers.S, model.timeseries.S[time]) + update_u_velocity && parent(model.velocities.u) .= parent(model.timeseries.u[time]) + update_v_velocity && parent(model.velocities.v) .= parent(model.timeseries.v[time]) + update_temperature && parent(model.tracers.T) .= parent(model.timeseries.T[time]) + update_salinity && parent(model.tracers.S) .= parent(model.timeseries.S[time]) return nothing end From 7e5babd0c494a741babc0b49170b45bd5a27b802 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 14:11:37 +0100 Subject: [PATCH 03/14] add a docstring --- src/OceanSimulations/prescribed_ocean_model.jl | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl index e2bf596fb..5d0c42792 100644 --- a/src/OceanSimulations/prescribed_ocean_model.jl +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -10,8 +10,8 @@ import Oceananigans.Architectures: on_architecture ##### A prescribed ocean... ##### -struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} - architecture :: A +struct PrescribedOcean{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} + architecture :: Arch grid :: G clock :: Clock{C} velocities :: U @@ -19,6 +19,18 @@ struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} timeseries :: F end +""" + PrescribedOcean(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) + +Create a prescribed ocean model to be used in combination with ClimaOcean's `OceanSeaIceModel` +on a `grid` with a `clock`. + +Arguments +========= +- `timeseries::NamedTuple`: A named tuple containing time series data for various fields. The named tuple can be empty or + include any combination of the following fields: `u`, `v`, `T`, `S`. + All elements provided must be of type `FieldTimeSeries` and reside on the provided `grid`. +""" function PrescribedOcean(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) @@ -44,7 +56,7 @@ function PrescribedOcean(timeseries=NamedTuple(); T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) - PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) + return PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) end ##### From 83b1a8ab8df52551d8c5f0821517b84d9b2e3c8c Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 14:12:24 +0100 Subject: [PATCH 04/14] better comment --- src/OceanSimulations/prescribed_ocean_model.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl index 5d0c42792..15a73ca33 100644 --- a/src/OceanSimulations/prescribed_ocean_model.jl +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -27,9 +27,10 @@ on a `grid` with a `clock`. Arguments ========= -- `timeseries::NamedTuple`: A named tuple containing time series data for various fields. The named tuple can be empty or - include any combination of the following fields: `u`, `v`, `T`, `S`. - All elements provided must be of type `FieldTimeSeries` and reside on the provided `grid`. +- `timeseries`: A `NamedTuple` containing time series data for various fields. The named tuple can be empty + (meaning that the ocean does not evolve in time) or include any combination of the + following fields: `u`, `v`, `T`, `S`. All elements provided must be of type `FieldTimeSeries` + and reside on the provided `grid`. """ function PrescribedOcean(timeseries=NamedTuple(); grid, From 6a924458ed8b1bc336768bb1ea03f72735d614df Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 14:22:08 +0100 Subject: [PATCH 05/14] let's go --- src/OceanSimulations/prescribed_ocean_model.jl | 1 + src/SeaIceSimulations.jl | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl index 15a73ca33..8fbac108e 100644 --- a/src/OceanSimulations/prescribed_ocean_model.jl +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -1,3 +1,4 @@ +using Oceananigans.Models: AbstractModel using Oceananigans.OutputReaders: extract_field_time_series, update_field_time_series! import Oceananigans.TimeSteppers: time_step!, update_state!, reset!, tick! diff --git a/src/SeaIceSimulations.jl b/src/SeaIceSimulations.jl index 510680651..59576aa0c 100644 --- a/src/SeaIceSimulations.jl +++ b/src/SeaIceSimulations.jl @@ -22,8 +22,8 @@ using ClimaOcean.OceanSimulations: Default import ClimaOcean: reference_density, heat_capacity -reference_density(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density -heat_capacity(sea_ice::SeaIceSimulation) = 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 +heat_capacity(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity function sea_ice_simulation(grid; Δt = 5minutes, From 60a0c65ae9c58f1d57c88bec74a27a03b052b2d8 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 15:01:47 +0100 Subject: [PATCH 06/14] it's a model --- src/OceanSimulations/prescribed_ocean_model.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl index 8fbac108e..54d79c578 100644 --- a/src/OceanSimulations/prescribed_ocean_model.jl +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -11,7 +11,7 @@ import Oceananigans.Architectures: on_architecture ##### A prescribed ocean... ##### -struct PrescribedOcean{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} +struct PrescribedOceanModel{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} architecture :: Arch grid :: G clock :: Clock{C} @@ -21,7 +21,7 @@ struct PrescribedOcean{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} end """ - PrescribedOcean(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) + PrescribedOceanModel(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) Create a prescribed ocean model to be used in combination with ClimaOcean's `OceanSeaIceModel` on a `grid` with a `clock`. @@ -33,7 +33,7 @@ Arguments following fields: `u`, `v`, `T`, `S`. All elements provided must be of type `FieldTimeSeries` and reside on the provided `grid`. """ -function PrescribedOcean(timeseries=NamedTuple(); +function PrescribedOceanModel(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) @@ -58,14 +58,14 @@ function PrescribedOcean(timeseries=NamedTuple(); T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) - return PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) + return PrescribedOceanModel(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) end ##### ##### Need to extend a couple of methods ##### -function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) +function time_step!(model::PrescribedOceanModel, Δt; callbacks=[], euler=true) tick!(model.clock, Δt) time = Time(model.clock.time) @@ -89,8 +89,8 @@ function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) return nothing end -update_state!(::PrescribedOcean) = nothing -timestepper(::PrescribedOcean) = nothing +update_state!(::PrescribedOceanModel) = nothing +timestepper(::PrescribedOceanModel) = nothing -reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6 -heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6 +reference_density(ocean::Simulation{<:PrescribedOceanModel}) = 1025.6 +heat_capacity(ocean::Simulation{<:PrescribedOceanModel}) = 3995.6 From c1441c2c91c1066bfe3b4ce17eabac636f34e987 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 15:22:03 +0100 Subject: [PATCH 07/14] update --- .../prescribed_ocean_model.jl | 32 ++++++++++++------- 1 file changed, 21 insertions(+), 11 deletions(-) diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl index 54d79c578..e1ff82296 100644 --- a/src/OceanSimulations/prescribed_ocean_model.jl +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -1,4 +1,5 @@ using Oceananigans.Models: AbstractModel +using Oceananigans.Fields: ZeroField using Oceananigans.OutputReaders: extract_field_time_series, update_field_time_series! import Oceananigans.TimeSteppers: time_step!, update_state!, reset!, tick! @@ -40,10 +41,11 @@ function PrescribedOceanModel(timeseries=NamedTuple(); # Make sure all elements of the timeseries are on the same grid # If we decide to pass a timeseries if !isempty(timeseries) - for (k, v) in timeseries - isa(v, FieldTimeSeries) || + for k in keys(timeseries) + f = timeseries[k] + isa(f, FieldTimeSeries) || throw(ArgumentError("All variables in the `timeseries` argument must be `FieldTimeSeries`")) - v.grid == grid || + f.grid == grid || throw(ArgumentError("All variables in the timeseries reside on the provided grid")) end end @@ -76,15 +78,23 @@ function time_step!(model::PrescribedOceanModel, Δt; callbacks=[], euler=true) update_field_time_series!(fts, time) end - update_u_velocity = haskey(model.timeseries, :u) - update_v_velocity = haskey(model.timeseries, :v) - update_temperature = haskey(model.timeseries, :T) - update_salinity = haskey(model.timeseries, :S) + # Time stepping the model! - update_u_velocity && parent(model.velocities.u) .= parent(model.timeseries.u[time]) - update_v_velocity && parent(model.velocities.v) .= parent(model.timeseries.v[time]) - update_temperature && parent(model.tracers.T) .= parent(model.timeseries.T[time]) - update_salinity && parent(model.tracers.S) .= parent(model.timeseries.S[time]) + if haskey(model.timeseries, :u) + arent(model.velocities.u) .= parent(model.timeseries.u[time]) + end + + if haskey(model.timeseries, :v) + parent(model.velocities.v) .= parent(model.timeseries.v[time]) + end + + if haskey(model.timeseries, :T) + parent(model.tracers.T) .= parent(model.timeseries.T[time]) + end + + if haskey(model.timeseries, :S) + parent(model.tracers.S) .= parent(model.timeseries.S[time]) + end return nothing end From fb76671aba1d4dae8a54144fa91aba4828c92670 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 15:26:30 +0100 Subject: [PATCH 08/14] import correct stuff --- src/OceanSeaIceModels/OceanSeaIceModels.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl index 4dee01e38..8316a157f 100644 --- a/src/OceanSeaIceModels/OceanSeaIceModels.jl +++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl @@ -25,10 +25,10 @@ using ClimaOcean: stateindex using KernelAbstractions: @kernel, @index using KernelAbstractions.Extras.LoopInfo: @unroll +import ClimaOcean: reference_density, heat_capacity + function downwelling_radiation end function freshwater_flux end -function reference_density end -function heat_capacity end const default_gravitational_acceleration = 9.80665 const default_freshwater_density = 1000 From 100f96cf5c2ce734d448ba4d76dc569ec089d868 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Sun, 16 Nov 2025 10:20:58 +0100 Subject: [PATCH 09/14] prescribed ocean model to prescribed ocean --- src/OceanSimulations/OceanSimulations.jl | 2 +- ...bed_ocean_model.jl => prescribed_ocean.jl} | 39 ++++++++++++------- 2 files changed, 27 insertions(+), 14 deletions(-) rename src/OceanSimulations/{prescribed_ocean_model.jl => prescribed_ocean.jl} (72%) diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 10f425c2f..402f62034 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -51,6 +51,6 @@ default_or_override(override, alternative_default=nothing) = override include("barotropic_potential_forcing.jl") include("radiative_forcing.jl") include("ocean_simulation.jl") -include("prescribed_ocean_model.jl") +include("prescribed_ocean.jl") end # module diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean.jl similarity index 72% rename from src/OceanSimulations/prescribed_ocean_model.jl rename to src/OceanSimulations/prescribed_ocean.jl index e1ff82296..4f377b95e 100644 --- a/src/OceanSimulations/prescribed_ocean_model.jl +++ b/src/OceanSimulations/prescribed_ocean.jl @@ -12,17 +12,18 @@ import Oceananigans.Architectures: on_architecture ##### A prescribed ocean... ##### -struct PrescribedOceanModel{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} - architecture :: Arch +struct PrescribedOcean{G, C, U, T, F, FT} grid :: G clock :: Clock{C} velocities :: U tracers :: T timeseries :: F + reference_density :: FT + heat_capacity :: FT end """ - PrescribedOceanModel(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) + PrescribedOcean(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) Create a prescribed ocean model to be used in combination with ClimaOcean's `OceanSeaIceModel` on a `grid` with a `clock`. @@ -34,9 +35,11 @@ Arguments following fields: `u`, `v`, `T`, `S`. All elements provided must be of type `FieldTimeSeries` and reside on the provided `grid`. """ -function PrescribedOceanModel(timeseries=NamedTuple(); +function PrescribedOcean(timeseries=NamedTuple(); grid, - clock=Clock{Float64}(time = 0)) + clock = Clock{Float64}(time = 0), + reference_density = 1029, + heat_capacity = 3998) # Make sure all elements of the timeseries are on the same grid # If we decide to pass a timeseries @@ -60,14 +63,23 @@ function PrescribedOceanModel(timeseries=NamedTuple(); T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) - return PrescribedOceanModel(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) + reference_density = convert(eltype(grid), reference_density) + heat_capacity = convert(eltype(grid), heat_capacity) + + return PrescribedOcean(grid, + clock, + (; u, v, w=ZeroField()), + (; T, S), + timeseries, + reference_density, + heat_capacity) end ##### ##### Need to extend a couple of methods ##### -function time_step!(model::PrescribedOceanModel, Δt; callbacks=[], euler=true) +function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) tick!(model.clock, Δt) time = Time(model.clock.time) @@ -79,9 +91,8 @@ function time_step!(model::PrescribedOceanModel, Δt; callbacks=[], euler=true) end # Time stepping the model! - if haskey(model.timeseries, :u) - arent(model.velocities.u) .= parent(model.timeseries.u[time]) + parent(model.velocities.u) .= parent(model.timeseries.u[time]) end if haskey(model.timeseries, :v) @@ -99,8 +110,10 @@ function time_step!(model::PrescribedOceanModel, Δt; callbacks=[], euler=true) return nothing end -update_state!(::PrescribedOceanModel) = nothing -timestepper(::PrescribedOceanModel) = nothing +get_ocean_state() + +update_state!(::PrescribedOcean) = nothing +timestepper(::PrescribedOcean) = nothing -reference_density(ocean::Simulation{<:PrescribedOceanModel}) = 1025.6 -heat_capacity(ocean::Simulation{<:PrescribedOceanModel}) = 3995.6 +reference_density(ocean::PrescribedOcean) = ocean.reference_density +heat_capacity(ocean::PrescribedOcean) = ocean.heat_capacity From 44c76a1512aba3e3f643c236bff30fe6fe891854 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Sun, 16 Nov 2025 11:10:27 +0100 Subject: [PATCH 10/14] simplify --- src/OceanSimulations/prescribed_ocean.jl | 40 +++++------------------- 1 file changed, 8 insertions(+), 32 deletions(-) diff --git a/src/OceanSimulations/prescribed_ocean.jl b/src/OceanSimulations/prescribed_ocean.jl index 4f377b95e..46c580c80 100644 --- a/src/OceanSimulations/prescribed_ocean.jl +++ b/src/OceanSimulations/prescribed_ocean.jl @@ -35,8 +35,9 @@ Arguments following fields: `u`, `v`, `T`, `S`. All elements provided must be of type `FieldTimeSeries` and reside on the provided `grid`. """ -function PrescribedOcean(timeseries=NamedTuple(); - grid, +function PrescribedOcean(grid; + velocities = default_ocean_velocities(grid), + tracers = default_ocean_tracers(grid), clock = Clock{Float64}(time = 0), reference_density = 1029, heat_capacity = 3998) @@ -53,28 +54,20 @@ function PrescribedOcean(timeseries=NamedTuple(); end end - τx = Field{Face, Center, Nothing}(grid) - τy = Field{Center, Face, Nothing}(grid) - Jᵀ = Field{Center, Center, Nothing}(grid) - Jˢ = Field{Center, Center, Nothing}(grid) - - u = XFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Face, Center, Center), top = FluxBoundaryCondition(τx))) - v = YFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Face, Center), top = FluxBoundaryCondition(τy))) - T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) - S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) - reference_density = convert(eltype(grid), reference_density) heat_capacity = convert(eltype(grid), heat_capacity) return PrescribedOcean(grid, clock, - (; u, v, w=ZeroField()), - (; T, S), - timeseries, + velocities, + tracers, reference_density, heat_capacity) end +default_ocean_velocities(grid) = (u = XFaceField(grid), v = YFaceField(grid)) +default_ocean_tracers(grid) = (T = CenterField(grid), S = CenterField(grid)) + ##### ##### Need to extend a couple of methods ##### @@ -90,23 +83,6 @@ function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) update_field_time_series!(fts, time) end - # Time stepping the model! - if haskey(model.timeseries, :u) - parent(model.velocities.u) .= parent(model.timeseries.u[time]) - end - - if haskey(model.timeseries, :v) - parent(model.velocities.v) .= parent(model.timeseries.v[time]) - end - - if haskey(model.timeseries, :T) - parent(model.tracers.T) .= parent(model.timeseries.T[time]) - end - - if haskey(model.timeseries, :S) - parent(model.tracers.S) .= parent(model.timeseries.S[time]) - end - return nothing end From 80fae2ca420a991c375e4d0faf078ebb514e79c7 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Sun, 16 Nov 2025 12:00:16 +0100 Subject: [PATCH 11/14] improve --- src/OceanSimulations/prescribed_ocean.jl | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/OceanSimulations/prescribed_ocean.jl b/src/OceanSimulations/prescribed_ocean.jl index 46c580c80..c335d01a1 100644 --- a/src/OceanSimulations/prescribed_ocean.jl +++ b/src/OceanSimulations/prescribed_ocean.jl @@ -42,18 +42,6 @@ function PrescribedOcean(grid; reference_density = 1029, heat_capacity = 3998) - # Make sure all elements of the timeseries are on the same grid - # If we decide to pass a timeseries - if !isempty(timeseries) - for k in keys(timeseries) - f = timeseries[k] - isa(f, FieldTimeSeries) || - throw(ArgumentError("All variables in the `timeseries` argument must be `FieldTimeSeries`")) - f.grid == grid || - throw(ArgumentError("All variables in the timeseries reside on the provided grid")) - end - end - reference_density = convert(eltype(grid), reference_density) heat_capacity = convert(eltype(grid), heat_capacity) From 546e8448dbdf27eb8b9170808390c01d4ed0e2c2 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Sun, 16 Nov 2025 12:01:37 +0100 Subject: [PATCH 12/14] start with one fix --- src/OceanSimulations/prescribed_ocean.jl | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/OceanSimulations/prescribed_ocean.jl b/src/OceanSimulations/prescribed_ocean.jl index c335d01a1..4223c8b61 100644 --- a/src/OceanSimulations/prescribed_ocean.jl +++ b/src/OceanSimulations/prescribed_ocean.jl @@ -2,8 +2,7 @@ using Oceananigans.Models: AbstractModel using Oceananigans.Fields: ZeroField using Oceananigans.OutputReaders: extract_field_time_series, update_field_time_series! -import Oceananigans.TimeSteppers: time_step!, update_state!, reset!, tick! -import Oceananigans.Models: timestepper, update_model_field_time_series! +import Oceananigans.TimeSteppers: time_step! import ClimaOcean: reference_density, heat_capacity import Oceananigans.Architectures: on_architecture @@ -74,10 +73,8 @@ function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) return nothing end -get_ocean_state() - -update_state!(::PrescribedOcean) = nothing -timestepper(::PrescribedOcean) = nothing +# TODO: fix this one... +get_ocean_state(::PrescribedOcean) = .... reference_density(ocean::PrescribedOcean) = ocean.reference_density heat_capacity(ocean::PrescribedOcean) = ocean.heat_capacity From 3c8b099d2b9ec36404261b5b292e2002bfbc3600 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Sun, 16 Nov 2025 12:07:56 +0100 Subject: [PATCH 13/14] start the restructuring --- .../PrescribedAtmospheres.jl | 3 + src/ClimaOcean.jl | 3 +- src/DataWrangling/restoring.jl | 6 +- .../assemble_net_fluxes.jl | 23 +++-- .../atmosphere_ocean_fluxes.jl | 20 +++-- .../component_interfaces.jl | 87 +++++++++---------- .../interpolate_atmospheric_state.jl | 6 +- .../sea_ice_ocean_fluxes.jl | 16 ++-- src/OceanSeaIceModels/OceanSeaIceModels.jl | 3 + src/OceanSimulations/OceanSimulations.jl | 12 +-- src/OceanSimulations/ocean_simulation.jl | 11 +++ src/OceanSimulations/prescribed_ocean.jl | 7 +- .../SeaIceSimulations.jl | 0 13 files changed, 106 insertions(+), 91 deletions(-) rename src/{OceanSeaIceModels => AtmosphereSimulations}/PrescribedAtmospheres.jl (99%) rename src/{ => SeaIceSimulations}/SeaIceSimulations.jl (100%) diff --git a/src/OceanSeaIceModels/PrescribedAtmospheres.jl b/src/AtmosphereSimulations/PrescribedAtmospheres.jl similarity index 99% rename from src/OceanSeaIceModels/PrescribedAtmospheres.jl rename to src/AtmosphereSimulations/PrescribedAtmospheres.jl index e7a3ca18c..15764b14c 100644 --- a/src/OceanSeaIceModels/PrescribedAtmospheres.jl +++ b/src/AtmosphereSimulations/PrescribedAtmospheres.jl @@ -5,6 +5,7 @@ using Oceananigans.Fields: Center using Oceananigans.Grids: grid_name using Oceananigans.OutputReaders: FieldTimeSeries, update_field_time_series!, extract_field_time_series using Oceananigans.TimeSteppers: Clock, tick! +using Oceananigans.Simulations: TimeStepWizard using Oceananigans.Utils: prettysummary, Time using Adapt @@ -371,6 +372,8 @@ end return nothing end +(wizard::TimeStepWizard)(atmos::PrescribedAtmosphere) = Inf + @inline thermodynamics_parameters(atmos::Nothing) = nothing @inline thermodynamics_parameters(atmos::PrescribedAtmosphere) = atmos.thermodynamics_parameters @inline surface_layer_height(atmos::PrescribedAtmosphere) = atmos.surface_layer_height diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl index 76407c547..e0b7a2af1 100644 --- a/src/ClimaOcean.jl +++ b/src/ClimaOcean.jl @@ -89,7 +89,8 @@ end end include("OceanSimulations/OceanSimulations.jl") -include("SeaIceSimulations.jl") +include("SeaIceSimulations/SeaIceSimulations.jl") +include("AtmosphereSimulations/AtmosphereSimulations.jl") include("OceanSeaIceModels/OceanSeaIceModels.jl") include("InitialConditions/InitialConditions.jl") include("DataWrangling/DataWrangling.jl") diff --git a/src/DataWrangling/restoring.jl b/src/DataWrangling/restoring.jl index 489efb1d0..5355acfbb 100644 --- a/src/DataWrangling/restoring.jl +++ b/src/DataWrangling/restoring.jl @@ -11,7 +11,7 @@ using NCDatasets using Dates: Second import ClimaOcean: stateindex -import Oceananigans.Forcings: materialize_forcing +import Oceananigans.Forcings: regularize_forcing # Variable names for restorable data struct Temperature end @@ -230,7 +230,7 @@ function Base.show(io::IO, dsr::DatasetRestoring) "└── native_grid: ", summary(dsr.native_grid)) end -materialize_forcing(forcing::DatasetRestoring, field, field_name, model_field_names) = forcing +regularize_forcing(forcing::DatasetRestoring, field, field_name, model_field_names) = forcing ##### ##### Masks for restoring @@ -293,4 +293,4 @@ end LX, LY, LZ = loc λ, φ, z = node(i, j, k, grid, LX(), LY(), LZ()) return mask(φ, z) -end +end \ No newline at end of file diff --git a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl index 1532a3be7..1fefa089f 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl @@ -21,18 +21,22 @@ using ClimaOcean.OceanSeaIceModels: sea_ice_concentration return zero(Iˢʷ) end -get_radiative_forcing(FT) = FT -function get_radiative_forcing(FT::MultipleForcings) +@inline get_radiative_forcing(ocean::OceananigansSimulation) = get_radiative_forcing(ocean.model.forcing.T) +@inline get_radiative_forcing(FT) = FT + +@inline function get_radiative_forcing(FT::MultipleForcings) for forcing in FT.forcings forcing isa TwoColorRadiation && return forcing end return nothing end +# No need to do this for an Oceananigans Simulation +fill_net_fluxes!(ocean, net_ocean_fluxes) = nothing + function compute_net_ocean_fluxes!(coupled_model) - ocean = coupled_model.ocean sea_ice = coupled_model.sea_ice - grid = ocean.model.grid + grid = coupled_model.interfaces.exchanger.exchange_grid arch = architecture(grid) clock = coupled_model.clock @@ -56,13 +60,14 @@ function compute_net_ocean_fluxes!(coupled_model) freshwater_flux = atmosphere_fields.Mp.data ice_concentration = sea_ice_concentration(sea_ice) - ocean_salinity = ocean.model.tracers.S + ocean_state = get_ocean_state(coupled_model.ocean, coupled_model.interfaces.exchanger) + ocean_salinity = ocean_state.S atmos_ocean_properties = coupled_model.interfaces.atmosphere_ocean_interface.properties ocean_properties = coupled_model.interfaces.ocean_properties kernel_parameters = interface_kernel_parameters(grid) ocean_surface_temperature = coupled_model.interfaces.atmosphere_ocean_interface.temperature - penetrating_radiation = get_radiative_forcing(ocean.model.forcing.T) + penetrating_radiation = get_radiative_forcing(coupled_model.ocean) launch!(arch, grid, kernel_parameters, _assemble_net_ocean_fluxes!, @@ -80,6 +85,8 @@ function compute_net_ocean_fluxes!(coupled_model) atmos_ocean_properties, ocean_properties) + fill_net_fluxes!(coupled_model.ocean, net_ocean_fluxes) + return nothing end @@ -261,8 +268,8 @@ end @inbounds begin Ts = surface_temperature[i, j, kᴺ] Ts = convert_to_kelvin(sea_ice_properties.temperature_units, Ts) - ℵi = ice_concentration[i, j, 1] - + ℵi = ice_concentration[i, j, kᴺ] + Qs = downwelling_radiation.Qs[i, j, 1] Qℓ = downwelling_radiation.Qℓ[i, j, 1] Qc = atmosphere_sea_ice_fluxes.sensible_heat[i, j, 1] # sensible or "conductive" heat flux diff --git a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl index 319ba5366..041e81a43 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl @@ -4,19 +4,21 @@ using ClimaOcean.OceanSeaIceModels.PrescribedAtmospheres: thermodynamics_paramet surface_layer_height, boundary_layer_height -function compute_atmosphere_ocean_fluxes!(coupled_model) - ocean = coupled_model.ocean - atmosphere = coupled_model.atmosphere - grid = ocean.model.grid - arch = architecture(grid) - clock = coupled_model.clock - - ocean_state = (u = ocean.model.velocities.u, +@inline get_ocean_state(ocean::OceananigansSimulation, coupled_model) = + (u = ocean.model.velocities.u, v = ocean.model.velocities.v, T = ocean.model.tracers.T, S = ocean.model.tracers.S) - atmosphere_fields = coupled_model.interfaces.exchanger.exchange_atmosphere_state +function compute_atmosphere_ocean_fluxes!(coupled_model) + ocean = coupled_model.ocean + atmosphere = coupled_model.atmosphere + exchanger = coupled_model.interfaces.exchanger + grid = exchanger.exchange_grid + arch = architecture(grid) + clock = coupled_model.clock + ocean_state = get_ocean_state(ocean, exchanger) + atmosphere_fields = exchanger.exchange_atmosphere_state # Simplify NamedTuple to reduce parameter space consumption. # See https://github.com/CliMA/ClimaOcean.jl/issues/116. diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 592794f0e..63e42f231 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -8,7 +8,8 @@ using ..OceanSeaIceModels: reference_density, sea_ice_thickness, downwelling_radiation, freshwater_flux, - SeaIceSimulation + SeaIceSimulation, + OceananigansSimulation using ..OceanSeaIceModels.PrescribedAtmospheres: PrescribedAtmosphere, @@ -17,7 +18,7 @@ using ..OceanSeaIceModels.PrescribedAtmospheres: using ClimaSeaIce: SeaIceModel using Oceananigans: HydrostaticFreeSurfaceModel, architecture -using Oceananigans.Grids: inactive_node, node, topology +using Oceananigans.Grids: inactive_node, node, topology, AbstractGrid using Oceananigans.BoundaryConditions: fill_halo_regions! using Oceananigans.Fields: ConstantField, interpolate, FractionalIndices using Oceananigans.Utils: launch!, Time, KernelParameters @@ -86,10 +87,9 @@ ExchangeAtmosphereState(grid) = ExchangeAtmosphereState(Field{Center, Center, No fractional_index_type(FT, Topo) = FT fractional_index_type(FT, ::Flat) = Nothing +state_exchanger(ocean::Simulation, ::Nothing) = nothing -StateExchanger(ocean::Simulation, ::Nothing) = nothing - -function StateExchanger(ocean::Simulation, atmosphere) +function state_exchanger(ocean::Simulation, atmosphere) # TODO: generalize this exchange_grid = ocean.model.grid exchange_atmosphere_state = ExchangeAtmosphereState(exchange_grid) @@ -116,12 +116,11 @@ function atmosphere_exchanger(atmosphere::PrescribedAtmosphere, exchange_grid) end initialize!(exchanger::StateExchanger, ::Nothing) = nothing +initialize!(exchanger::StateExchanger, atmosphere) = initialize!(exchanger.atmosphere_exchanger, exchanger.exchange_grid, atmosphere) -function initialize!(exchanger::StateExchanger, atmosphere) +function initialize!(frac_indices, exchange_grid::AbstractGrid, atmosphere) atmos_grid = atmosphere.grid - exchange_grid = exchanger.exchange_grid arch = architecture(exchange_grid) - frac_indices = exchanger.atmosphere_exchanger kernel_parameters = interface_kernel_parameters(exchange_grid) launch!(arch, exchange_grid, kernel_parameters, _compute_fractional_indices!, frac_indices, exchange_grid, atmos_grid) @@ -165,26 +164,28 @@ Base.summary(crf::ComponentInterfaces) = "ComponentInterfaces" Base.show(io::IO, crf::ComponentInterfaces) = print(io, summary(crf)) atmosphere_ocean_interface(::Nothing, args...) = nothing +atmosphere_ocean_interface(ocean::OceananigansSimulation, args...) = + atmosphere_ocean_interface(ocean.model.grid, args...) -function atmosphere_ocean_interface(atmos, - ocean, +function atmosphere_ocean_interface(grid::AbstractGrid, + atmos, radiation, ao_flux_formulation, temperature_formulation, velocity_formulation, specific_humidity_formulation) - water_vapor = Field{Center, Center, Nothing}(ocean.model.grid) - latent_heat = Field{Center, Center, Nothing}(ocean.model.grid) - sensible_heat = Field{Center, Center, Nothing}(ocean.model.grid) - x_momentum = Field{Center, Center, Nothing}(ocean.model.grid) - y_momentum = Field{Center, Center, Nothing}(ocean.model.grid) - friction_velocity = Field{Center, Center, Nothing}(ocean.model.grid) - temperature_scale = Field{Center, Center, Nothing}(ocean.model.grid) - water_vapor_scale = Field{Center, Center, Nothing}(ocean.model.grid) - upwelling_longwave = Field{Center, Center, Nothing}(ocean.model.grid) - downwelling_longwave = Field{Center, Center, Nothing}(ocean.model.grid) - downwelling_shortwave = Field{Center, Center, Nothing}(ocean.model.grid) + water_vapor = Field{Center, Center, Nothing}(grid) + latent_heat = Field{Center, Center, Nothing}(grid) + sensible_heat = Field{Center, Center, Nothing}(grid) + x_momentum = Field{Center, Center, Nothing}(grid) + y_momentum = Field{Center, Center, Nothing}(grid) + friction_velocity = Field{Center, Center, Nothing}(grid) + temperature_scale = Field{Center, Center, Nothing}(grid) + water_vapor_scale = Field{Center, Center, Nothing}(grid) + upwelling_longwave = Field{Center, Center, Nothing}(grid) + downwelling_longwave = Field{Center, Center, Nothing}(grid) + downwelling_shortwave = Field{Center, Center, Nothing}(grid) ao_fluxes = (; latent_heat, sensible_heat, @@ -208,7 +209,7 @@ function atmosphere_ocean_interface(atmos, temperature_formulation, velocity_formulation) - interface_temperature = Field{Center, Center, Nothing}(ocean.model.grid) + interface_temperature = Field{Center, Center, Nothing}(grid) return AtmosphereInterface(ao_fluxes, ao_flux_formulation, interface_temperature, ao_properties) end @@ -250,16 +251,18 @@ function atmosphere_sea_ice_interface(atmos, return AtmosphereInterface(fluxes, ai_flux_formulation, interface_temperature, properties) end -sea_ice_ocean_interface(sea_ice, ocean) = nothing +sea_ice_ocean_interface(ocean, sea_ice) = nothing +sea_ice_ocean_interface(ocean, sea_ice::SeaIceSimulation; kw...) = + sea_ice_ocean_interface(ocean.grid, sea_ice; kw...) -function sea_ice_ocean_interface(sea_ice::SeaIceSimulation, ocean; +function sea_ice_ocean_interface(grid::AbstractGrid, sea_ice::SeaIceSimulation; characteristic_melting_speed = 1e-5) - io_bottom_heat_flux = Field{Center, Center, Nothing}(ocean.model.grid) - io_frazil_heat_flux = Field{Center, Center, Nothing}(ocean.model.grid) - io_salt_flux = Field{Center, Center, Nothing}(ocean.model.grid) - x_momentum = Field{Face, Center, Nothing}(ocean.model.grid) - y_momentum = Field{Center, Face, Nothing}(ocean.model.grid) + io_bottom_heat_flux = Field{Center, Center, Nothing}(grid) + io_frazil_heat_flux = Field{Center, Center, Nothing}(grid) + io_salt_flux = Field{Center, Center, Nothing}(grid) + x_momentum = Field{Face, Center, Nothing}(grid) + y_momentum = Field{Center, Face, Nothing}(grid) @assert io_frazil_heat_flux isa Field{Center, Center, Nothing} @assert io_bottom_heat_flux isa Field{Center, Center, Nothing} @@ -284,7 +287,7 @@ function default_ai_temperature(sea_ice::SeaIceSimulation) end function default_ao_specific_humidity(ocean) - FT = eltype(ocean.model.grid) + FT = eltype(ocean) phase = AtmosphericThermodynamics.Liquid() x_H₂O = convert(FT, 0.98) return ImpureSaturationSpecificHumidity(phase, x_H₂O) @@ -310,8 +313,8 @@ end function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; radiation = Radiation(), freshwater_density = default_freshwater_density, - atmosphere_ocean_fluxes = SimilarityTheoryFluxes(eltype(ocean.model.grid)), - atmosphere_sea_ice_fluxes = SimilarityTheoryFluxes(eltype(ocean.model.grid)), + atmosphere_ocean_fluxes = SimilarityTheoryFluxes(), + atmosphere_sea_ice_fluxes = SimilarityTheoryFluxes(), atmosphere_ocean_interface_temperature = BulkTemperature(), atmosphere_ocean_velocity_difference = RelativeVelocity(), atmosphere_ocean_interface_specific_humidity = default_ao_specific_humidity(ocean), @@ -325,8 +328,7 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; sea_ice_heat_capacity = heat_capacity(sea_ice), gravitational_acceleration = default_gravitational_acceleration) - ocean_grid = ocean.model.grid - FT = eltype(ocean_grid) + FT = eltype(ocean) ocean_reference_density = convert(FT, ocean_reference_density) ocean_heat_capacity = convert(FT, ocean_heat_capacity) @@ -342,8 +344,8 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; freshwater_density = freshwater_density, temperature_units = ocean_temperature_units) - ao_interface = atmosphere_ocean_interface(atmosphere, - ocean, + ao_interface = atmosphere_ocean_interface(ocean, + atmosphere, radiation, atmosphere_ocean_fluxes, atmosphere_ocean_interface_temperature, @@ -384,23 +386,14 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; net_bottom_sea_ice_fluxes = nothing end - τx = surface_flux(ocean.model.velocities.u) - τy = surface_flux(ocean.model.velocities.v) - tracers = ocean.model.tracers - ρₒ = ocean_reference_density - cₒ = ocean_heat_capacity - Qₒ = ρₒ * cₒ * surface_flux(ocean.model.tracers.T) - net_ocean_surface_fluxes = (u=τx, v=τy, Q=Qₒ) - - ocean_surface_tracer_fluxes = NamedTuple(name => surface_flux(tracers[name]) for name in keys(tracers)) - net_ocean_surface_fluxes = merge(ocean_surface_tracer_fluxes, net_ocean_surface_fluxes) + net_ocean_surface_fluxes = ocean_surface_fluxes(ocean, ocean_reference_density, ocean_heat_capacity) # Total interface fluxes net_fluxes = (ocean_surface = net_ocean_surface_fluxes, sea_ice_top = net_top_sea_ice_fluxes, sea_ice_bottom = net_bottom_sea_ice_fluxes) - exchanger = StateExchanger(ocean, atmosphere) + exchanger = state_exchanger(ocean, atmosphere) properties = (; gravitational_acceleration) diff --git a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl index d79a9caf9..63baf19f2 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl @@ -26,7 +26,7 @@ function interpolate_atmosphere_state!(interfaces, atmosphere::PrescribedAtmosph atmosphere_grid = atmosphere.grid # Basic model properties - grid = ocean.model.grid + grid = coupled_model.interfaces.exchanger.exchange_grid arch = architecture(grid) clock = coupled_model.clock @@ -125,8 +125,8 @@ function interpolate_atmosphere_state!(interfaces, atmosphere::PrescribedAtmosph # # TODO: find a better design for this that doesn't have redundant # arrays for the barotropic potential - u_potential = forcing_barotropic_potential(ocean.model.forcing.u) - v_potential = forcing_barotropic_potential(ocean.model.forcing.v) + u_potential = forcing_barotropic_potential(ocean) + v_potential = forcing_barotropic_potential(ocean) ρₒ = coupled_model.interfaces.ocean_properties.reference_density if !isnothing(u_potential) diff --git a/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl index 1cc1f7cb8..8d5d805a0 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl @@ -23,10 +23,13 @@ function compute_sea_ice_ocean_fluxes!(sea_ice_ocean_fluxes, ocean, sea_ice, mel ℵᵢ = sea_ice.model.ice_concentration hᵢ = sea_ice.model.ice_thickness Gh = sea_ice.model.ice_thermodynamics.thermodynamic_tendency + Δt = sea_ice.Δt + + ocean_state = get_ocean_state(ocean, coupled_model) liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus - grid = ocean.model.grid - clock = ocean.model.clock + grid = sea_ice.model.grid + clock = sea_ice.model.clock arch = architecture(grid) uᵢ, vᵢ = sea_ice.model.velocities @@ -41,7 +44,7 @@ function compute_sea_ice_ocean_fluxes!(sea_ice_ocean_fluxes, ocean, sea_ice, mel # What about the latent heat removed from the ocean when ice forms? # Is it immediately removed from the ocean? Or is it stored in the ice? launch!(arch, grid, :xy, _compute_sea_ice_ocean_fluxes!, - sea_ice_ocean_fluxes, grid, clock, hᵢ, ℵᵢ, Sᵢ, Gh, Tₒ, Sₒ, uᵢ, vᵢ, + sea_ice_ocean_fluxes, grid, clock, hᵢ, ℵᵢ, Sᵢ, Gh, ocean_state, uᵢ, vᵢ, τs, liquidus, ocean_properties, melting_speed, Δt) return nothing @@ -54,8 +57,7 @@ end ice_concentration, ice_salinity, thermodynamic_tendency, - ocean_temperature, - ocean_salinity, + ocean_state, sea_ice_u_velocity, sea_ice_v_velocity, sea_ice_ocean_stresses, @@ -74,8 +76,8 @@ end τy = sea_ice_ocean_fluxes.y_momentum uᵢ = sea_ice_u_velocity vᵢ = sea_ice_v_velocity - Tₒ = ocean_temperature - Sₒ = ocean_salinity + Tₒ = ocean_state.T + Sₒ = ocean_state.S Sᵢ = ice_salinity hᵢ = ice_thickness ℵᵢ = ice_concentration diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl index 87e2129e5..677cfeb7f 100644 --- a/src/OceanSeaIceModels/OceanSeaIceModels.jl +++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl @@ -39,6 +39,9 @@ const default_gravitational_acceleration = Oceananigans.defaults.gravitational_a const default_freshwater_density = 1000 # kg m⁻³ const SeaIceSimulation = Simulation{<:SeaIceModel} +const OceananigansSimulation = Simulation{<:HydrostaticFreeSurfaceModel} + +Base.eltype(ocean::OceananigansSimulation) = eltype(ocean.model.grid) sea_ice_thickness(::Nothing) = ZeroField() sea_ice_thickness(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thickness diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 402f62034..9c3bfabca 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -1,6 +1,6 @@ module OceanSimulations -export ocean_simulation, PrescribedOceanModel +export ocean_simulation using Oceananigans using Oceananigans.Units @@ -21,15 +21,6 @@ using SeawaterPolynomials.TEOS10: TEOS10EquationOfState default_gravitational_acceleration = Oceananigans.defaults.gravitational_acceleration default_planet_rotation_rate = Oceananigans.defaults.planet_rotation_rate -import ClimaOcean: reference_density, heat_capacity - -reference_density(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = reference_density(ocean.model.buoyancy.formulation) -reference_density(buoyancy_formulation::SeawaterBuoyancy) = reference_density(buoyancy_formulation.equation_of_state) -reference_density(eos::TEOS10EquationOfState) = eos.reference_density - -heat_capacity(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = heat_capacity(ocean.model.buoyancy.formulation) -heat_capacity(buoyancy_formulation::SeawaterBuoyancy) = heat_capacity(buoyancy_formulation.equation_of_state) - struct Default{V} value :: V end @@ -51,6 +42,5 @@ default_or_override(override, alternative_default=nothing) = override include("barotropic_potential_forcing.jl") include("radiative_forcing.jl") include("ocean_simulation.jl") -include("prescribed_ocean.jl") end # module diff --git a/src/OceanSimulations/ocean_simulation.jl b/src/OceanSimulations/ocean_simulation.jl index 9a34ccfbf..1de430db3 100644 --- a/src/OceanSimulations/ocean_simulation.jl +++ b/src/OceanSimulations/ocean_simulation.jl @@ -267,3 +267,14 @@ end hasclosure(closure, ClosureType) = closure isa ClosureType hasclosure(closure_tuple::Tuple, ClosureType) = any(hasclosure(c, ClosureType) for c in closure_tuple) + +function ocean_surface_fluxes(ocean::OceananigansSimulation, ρₒ, cₒ) + τx = surface_flux(ocean.model.velocities.u) + τy = surface_flux(ocean.model.velocities.v) + tracers = ocean.model.tracers + Qₒ = ρₒ * cₒ * surface_flux(ocean.model.tracers.T) + net_ocean_surface_fluxes = (u=τx, v=τy, Q=Qₒ) + + ocean_surface_tracer_fluxes = NamedTuple(name => surface_flux(tracers[name]) for name in keys(tracers)) + return merge(ocean_surface_tracer_fluxes, net_ocean_surface_fluxes) +end diff --git a/src/OceanSimulations/prescribed_ocean.jl b/src/OceanSimulations/prescribed_ocean.jl index 4223c8b61..70276cd83 100644 --- a/src/OceanSimulations/prescribed_ocean.jl +++ b/src/OceanSimulations/prescribed_ocean.jl @@ -73,8 +73,11 @@ function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) return nothing end -# TODO: fix this one... -get_ocean_state(::PrescribedOcean) = .... reference_density(ocean::PrescribedOcean) = ocean.reference_density heat_capacity(ocean::PrescribedOcean) = ocean.heat_capacity + +# TODO: fix these two... +get_ocean_state(::PrescribedOcean) = .... +function ocean_surface_fluxes(ocean::PrescribedOcean, ρₒ, cₒ) +end diff --git a/src/SeaIceSimulations.jl b/src/SeaIceSimulations/SeaIceSimulations.jl similarity index 100% rename from src/SeaIceSimulations.jl rename to src/SeaIceSimulations/SeaIceSimulations.jl From 5dc4887aaf9ac915cbbcd0b3112b8f6fe5bad08b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 18 Nov 2025 10:22:58 +0100 Subject: [PATCH 14/14] add a prescribed atmosphere in its own module --- .../AtmosphereSimulations.jl | 20 ++ ...> atmosphere_thermodynamics_parameters.jl} | 186 ------------------ .../prescribed_atmosphere.jl | 163 +++++++++++++++ src/ClimaOcean.jl | 3 + src/OceanSeaIceModels/OceanSeaIceModels.jl | 6 +- 5 files changed, 191 insertions(+), 187 deletions(-) create mode 100644 src/AtmosphereSimulations/AtmosphereSimulations.jl rename src/AtmosphereSimulations/{PrescribedAtmospheres.jl => atmosphere_thermodynamics_parameters.jl} (64%) create mode 100644 src/AtmosphereSimulations/prescribed_atmosphere.jl diff --git a/src/AtmosphereSimulations/AtmosphereSimulations.jl b/src/AtmosphereSimulations/AtmosphereSimulations.jl new file mode 100644 index 000000000..777101d65 --- /dev/null +++ b/src/AtmosphereSimulations/AtmosphereSimulations.jl @@ -0,0 +1,20 @@ +module AtmosphereSimulations + +using Oceananigans +using Oceananigans.Fields: Center +using Oceananigans.Grids: grid_name +using Oceananigans.OutputReaders: FieldTimeSeries, update_field_time_series!, extract_field_time_series +using Oceananigans.TimeSteppers: Clock, tick! +using Oceananigans.Simulations: TimeStepWizard +using Oceananigans.Utils: prettysummary, Time + +using Adapt +using Thermodynamics.Parameters: AbstractThermodynamicsParameters + +import Oceananigans.TimeSteppers: time_step!, update_state! +import ClimaOcean: compute_net_atmosphere_fluxes! + +include("atmospheric_thermodynamics_parameters.jl") +include("prescribed_atmosphere.jl") + +end # module diff --git a/src/AtmosphereSimulations/PrescribedAtmospheres.jl b/src/AtmosphereSimulations/atmosphere_thermodynamics_parameters.jl similarity index 64% rename from src/AtmosphereSimulations/PrescribedAtmospheres.jl rename to src/AtmosphereSimulations/atmosphere_thermodynamics_parameters.jl index 1d6d4d57e..aa1ba4da0 100644 --- a/src/AtmosphereSimulations/PrescribedAtmospheres.jl +++ b/src/AtmosphereSimulations/atmosphere_thermodynamics_parameters.jl @@ -1,18 +1,3 @@ -module PrescribedAtmospheres - -using Oceananigans -using Oceananigans.Fields: Center -using Oceananigans.Grids: grid_name -using Oceananigans.OutputReaders: FieldTimeSeries, update_field_time_series!, extract_field_time_series -using Oceananigans.TimeSteppers: Clock, tick! -using Oceananigans.Simulations: TimeStepWizard -using Oceananigans.Utils: prettysummary, Time - -using Adapt -using Thermodynamics.Parameters: AbstractThermodynamicsParameters - -import Oceananigans.TimeSteppers: time_step!, update_state! - import Thermodynamics.Parameters: gas_constant, # molmass_dryair, # Molar mass of dry air (without moisture) @@ -41,11 +26,6 @@ import Thermodynamics.Parameters: pow_icenuc # "Power parameter" that controls liquid/ice condensate partitioning # during partial ice nucleation -import ..OceanSeaIceModels: - downwelling_radiation, - freshwater_flux, - compute_net_atmosphere_fluxes! - ##### ##### Atmospheric thermodynamics parameters ##### @@ -287,169 +267,3 @@ const ATP = AtmosphereThermodynamicsParameters @inline cp_d(p::ATP) = R_d(p) / kappa_d(p) @inline cv_d(p::ATP) = cp_d(p) - R_d(p) @inline cv_v(p::ATP) = cp_v(p) - R_v(p) - -##### -##### Prescribed atmosphere (as opposed to dynamically evolving / prognostic) -##### - -mutable struct PrescribedAtmosphere{FT, G, T, U, P, C, F, I, R, TP, TI} - grid :: G - clock :: Clock{T} - velocities :: U - pressure :: P - tracers :: C - freshwater_flux :: F - auxiliary_freshwater_flux :: I - downwelling_radiation :: R - thermodynamics_parameters :: TP - times :: TI - surface_layer_height :: FT - boundary_layer_height :: FT -end - -function Base.summary(pa::PrescribedAtmosphere{FT}) where FT - Nx, Ny, Nz = size(pa.grid) - Nt = length(pa.times) - sz_str = string(Nx, "×", Ny, "×", Nz, "×", Nt) - return string(sz_str, " PrescribedAtmosphere{$FT}") -end - -function Base.show(io::IO, pa::PrescribedAtmosphere) - print(io, summary(pa), " on ", grid_name(pa.grid), ":", '\n') - print(io, "├── times: ", prettysummary(pa.times), '\n') - print(io, "├── surface_layer_height: ", prettysummary(pa.surface_layer_height), '\n') - print(io, "└── boundary_layer_height: ", prettysummary(pa.boundary_layer_height)) -end - -function default_atmosphere_velocities(grid, times) - ua = FieldTimeSeries{Center, Center, Nothing}(grid, times) - va = FieldTimeSeries{Center, Center, Nothing}(grid, times) - return (u=ua, v=va) -end - -function default_atmosphere_tracers(grid, times) - Ta = FieldTimeSeries{Center, Center, Nothing}(grid, times) - qa = FieldTimeSeries{Center, Center, Nothing}(grid, times) - parent(Ta) .= 273.15 + 20 - return (T=Ta, q=qa) -end - -function default_downwelling_radiation(grid, times) - Qℓ = FieldTimeSeries{Center, Center, Nothing}(grid, times) - Qs = FieldTimeSeries{Center, Center, Nothing}(grid, times) - return TwoBandDownwellingRadiation(shortwave=Qs, longwave=Qℓ) -end - -function default_freshwater_flux(grid, times) - rain = FieldTimeSeries{Center, Center, Nothing}(grid, times) - snow = FieldTimeSeries{Center, Center, Nothing}(grid, times) - return (; rain, snow) -end - -""" The standard unit of atmospheric pressure; 1 standard atmosphere (atm) = 101,325 Pascals (Pa) -in SI units. This is approximately equal to the mean sea-level atmospheric pressure on Earth. """ -function default_atmosphere_pressure(grid, times) - pa = FieldTimeSeries{Center, Center, Nothing}(grid, times) - parent(pa) .= 101325 - return pa -end - - -@inline function update_state!(atmos::PrescribedAtmosphere) - time = Time(atmos.clock.time) - ftses = extract_field_time_series(atmos) - - for fts in ftses - update_field_time_series!(fts, time) - end - return nothing -end - -@inline function time_step!(atmos::PrescribedAtmosphere, Δt) - tick!(atmos.clock, Δt) - - update_state!(atmos) - - return nothing -end - -(wizard::TimeStepWizard)(atmos::PrescribedAtmosphere) = Inf - -@inline thermodynamics_parameters(atmos::Nothing) = nothing -@inline thermodynamics_parameters(atmos::PrescribedAtmosphere) = atmos.thermodynamics_parameters -@inline surface_layer_height(atmos::PrescribedAtmosphere) = atmos.surface_layer_height -@inline boundary_layer_height(atmos::PrescribedAtmosphere) = atmos.boundary_layer_height - -# No need to compute anything here... -compute_net_atmosphere_fluxes!(coupled_model, ::PrescribedAtmosphere) = nothing - -""" - PrescribedAtmosphere(grid, times=[zero(grid)]; - clock = Clock{Float64}(time = 0), - surface_layer_height = 10, # meters - boundary_layer_height = 512 # meters, - thermodynamics_parameters = AtmosphereThermodynamicsParameters(eltype(grid)), - auxiliary_freshwater_flux = nothing, - velocities = default_atmosphere_velocities(grid, times), - tracers = default_atmosphere_tracers(grid, times), - pressure = default_atmosphere_pressure(grid, times), - freshwater_flux = default_freshwater_flux(grid, times), - downwelling_radiation = default_downwelling_radiation(grid, times)) - -Return a representation of a prescribed time-evolving atmospheric -state with data given at `times`. -""" -function PrescribedAtmosphere(grid, times=[zero(grid)]; - clock = Clock{Float64}(time = 0), - surface_layer_height = 10, - boundary_layer_height = 512, - thermodynamics_parameters = AtmosphereThermodynamicsParameters(eltype(grid)), - auxiliary_freshwater_flux = nothing, - velocities = default_atmosphere_velocities(grid, times), - tracers = default_atmosphere_tracers(grid, times), - pressure = default_atmosphere_pressure(grid, times), - freshwater_flux = default_freshwater_flux(grid, times), - downwelling_radiation = default_downwelling_radiation(grid, times)) - - FT = eltype(grid) - if isnothing(thermodynamics_parameters) - thermodynamics_parameters = AtmosphereThermodynamicsParameters(FT) - end - - atmosphere = PrescribedAtmosphere(grid, - clock, - velocities, - pressure, - tracers, - freshwater_flux, - auxiliary_freshwater_flux, - downwelling_radiation, - thermodynamics_parameters, - times, - convert(FT, surface_layer_height), - convert(FT, boundary_layer_height)) - update_state!(atmosphere) - - return atmosphere -end - -struct TwoBandDownwellingRadiation{SW, LW} - shortwave :: SW - longwave :: LW -end - -""" - TwoBandDownwellingRadiation(shortwave=nothing, longwave=nothing) - -Return a two-band model for downwelling radiation (split into a shortwave band -and a longwave band) that passes through the atmosphere and arrives at the surface of ocean -or sea ice. -""" -TwoBandDownwellingRadiation(; shortwave=nothing, longwave=nothing) = - TwoBandDownwellingRadiation(shortwave, longwave) - -Adapt.adapt_structure(to, tsdr::TwoBandDownwellingRadiation) = - TwoBandDownwellingRadiation(adapt(to, tsdr.shortwave), - adapt(to, tsdr.longwave)) - -end # module diff --git a/src/AtmosphereSimulations/prescribed_atmosphere.jl b/src/AtmosphereSimulations/prescribed_atmosphere.jl new file mode 100644 index 000000000..cd4c385c8 --- /dev/null +++ b/src/AtmosphereSimulations/prescribed_atmosphere.jl @@ -0,0 +1,163 @@ +##### +##### Prescribed atmosphere (as opposed to dynamically evolving / prognostic) +##### + +mutable struct PrescribedAtmosphere{FT, G, T, U, P, C, F, I, R, TP, TI} + grid :: G + clock :: Clock{T} + velocities :: U + pressure :: P + tracers :: C + freshwater_flux :: F + auxiliary_freshwater_flux :: I + downwelling_radiation :: R + thermodynamics_parameters :: TP + times :: TI + surface_layer_height :: FT + boundary_layer_height :: FT +end + +function Base.summary(pa::PrescribedAtmosphere{FT}) where FT + Nx, Ny, Nz = size(pa.grid) + Nt = length(pa.times) + sz_str = string(Nx, "×", Ny, "×", Nz, "×", Nt) + return string(sz_str, " PrescribedAtmosphere{$FT}") +end + +function Base.show(io::IO, pa::PrescribedAtmosphere) + print(io, summary(pa), " on ", grid_name(pa.grid), ":", '\n') + print(io, "├── times: ", prettysummary(pa.times), '\n') + print(io, "├── surface_layer_height: ", prettysummary(pa.surface_layer_height), '\n') + print(io, "└── boundary_layer_height: ", prettysummary(pa.boundary_layer_height)) +end + +function default_atmosphere_velocities(grid, times) + ua = FieldTimeSeries{Center, Center, Nothing}(grid, times) + va = FieldTimeSeries{Center, Center, Nothing}(grid, times) + return (u=ua, v=va) +end + +function default_atmosphere_tracers(grid, times) + Ta = FieldTimeSeries{Center, Center, Nothing}(grid, times) + qa = FieldTimeSeries{Center, Center, Nothing}(grid, times) + parent(Ta) .= 273.15 + 20 + return (T=Ta, q=qa) +end + +function default_downwelling_radiation(grid, times) + Qℓ = FieldTimeSeries{Center, Center, Nothing}(grid, times) + Qs = FieldTimeSeries{Center, Center, Nothing}(grid, times) + return TwoBandDownwellingRadiation(shortwave=Qs, longwave=Qℓ) +end + +function default_freshwater_flux(grid, times) + rain = FieldTimeSeries{Center, Center, Nothing}(grid, times) + snow = FieldTimeSeries{Center, Center, Nothing}(grid, times) + return (; rain, snow) +end + +""" The standard unit of atmospheric pressure; 1 standard atmosphere (atm) = 101,325 Pascals (Pa) +in SI units. This is approximately equal to the mean sea-level atmospheric pressure on Earth. """ +function default_atmosphere_pressure(grid, times) + pa = FieldTimeSeries{Center, Center, Nothing}(grid, times) + parent(pa) .= 101325 + return pa +end + + +@inline function update_state!(atmos::PrescribedAtmosphere) + time = Time(atmos.clock.time) + ftses = extract_field_time_series(atmos) + + for fts in ftses + update_field_time_series!(fts, time) + end + return nothing +end + +@inline function time_step!(atmos::PrescribedAtmosphere, Δt) + tick!(atmos.clock, Δt) + + update_state!(atmos) + + return nothing +end + +(wizard::TimeStepWizard)(atmos::PrescribedAtmosphere) = Inf + +@inline thermodynamics_parameters(atmos::Nothing) = nothing +@inline thermodynamics_parameters(atmos::PrescribedAtmosphere) = atmos.thermodynamics_parameters +@inline surface_layer_height(atmos::PrescribedAtmosphere) = atmos.surface_layer_height +@inline boundary_layer_height(atmos::PrescribedAtmosphere) = atmos.boundary_layer_height + +# No need to compute anything here... +compute_net_atmosphere_fluxes!(coupled_model, ::PrescribedAtmosphere) = nothing + +""" + PrescribedAtmosphere(grid, times=[zero(grid)]; + clock = Clock{Float64}(time = 0), + surface_layer_height = 10, # meters + boundary_layer_height = 512 # meters, + thermodynamics_parameters = AtmosphereThermodynamicsParameters(eltype(grid)), + auxiliary_freshwater_flux = nothing, + velocities = default_atmosphere_velocities(grid, times), + tracers = default_atmosphere_tracers(grid, times), + pressure = default_atmosphere_pressure(grid, times), + freshwater_flux = default_freshwater_flux(grid, times), + downwelling_radiation = default_downwelling_radiation(grid, times)) + +Return a representation of a prescribed time-evolving atmospheric +state with data given at `times`. +""" +function PrescribedAtmosphere(grid, times=[zero(grid)]; + clock = Clock{Float64}(time = 0), + surface_layer_height = 10, + boundary_layer_height = 512, + thermodynamics_parameters = AtmosphereThermodynamicsParameters(eltype(grid)), + auxiliary_freshwater_flux = nothing, + velocities = default_atmosphere_velocities(grid, times), + tracers = default_atmosphere_tracers(grid, times), + pressure = default_atmosphere_pressure(grid, times), + freshwater_flux = default_freshwater_flux(grid, times), + downwelling_radiation = default_downwelling_radiation(grid, times)) + + FT = eltype(grid) + if isnothing(thermodynamics_parameters) + thermodynamics_parameters = AtmosphereThermodynamicsParameters(FT) + end + + atmosphere = PrescribedAtmosphere(grid, + clock, + velocities, + pressure, + tracers, + freshwater_flux, + auxiliary_freshwater_flux, + downwelling_radiation, + thermodynamics_parameters, + times, + convert(FT, surface_layer_height), + convert(FT, boundary_layer_height)) + update_state!(atmosphere) + + return atmosphere +end + +struct TwoBandDownwellingRadiation{SW, LW} + shortwave :: SW + longwave :: LW +end + +""" + TwoBandDownwellingRadiation(shortwave=nothing, longwave=nothing) + +Return a two-band model for downwelling radiation (split into a shortwave band +and a longwave band) that passes through the atmosphere and arrives at the surface of ocean +or sea ice. +""" +TwoBandDownwellingRadiation(; shortwave=nothing, longwave=nothing) = + TwoBandDownwellingRadiation(shortwave, longwave) + +Adapt.adapt_structure(to, tsdr::TwoBandDownwellingRadiation) = + TwoBandDownwellingRadiation(adapt(to, tsdr.shortwave), + adapt(to, tsdr.longwave)) diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl index 0a4dc2de0..42a4c2765 100644 --- a/src/ClimaOcean.jl +++ b/src/ClimaOcean.jl @@ -130,4 +130,7 @@ using PrecompileTools: @setup_workload, @compile_workload end end +function + + end # module diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl index 6ae8f2bf5..7e0af42f7 100644 --- a/src/OceanSeaIceModels/OceanSeaIceModels.jl +++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl @@ -83,6 +83,10 @@ const NoSeaIceModel = Union{OceanSeaIceModel{Nothing}, OceanSeaIceModel{<:Freezi ##### Some implementation ##### +import ClimaOcean: compute_net_sea_ice_fluxes!, + compute_net_ocean_fluxes!, + compute_net_atmosphere_fluxes! + # Atmosphere interface interpolate_atmosphere_state!(interfaces, atmosphere, coupled_model) = nothing @@ -111,7 +115,7 @@ include("freezing_limited_ocean_temperature.jl") # TODO: import this last include("PrescribedAtmospheres.jl") -using .PrescribedAtmospheres: +using ClimaOcean.AtmosphereSimulations: PrescribedAtmosphere, AtmosphereThermodynamicsParameters, TwoBandDownwellingRadiation