diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index 8a7884375..4cbe3f317 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/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/OceanSeaIceModels/PrescribedAtmospheres.jl b/src/AtmosphereSimulations/atmosphere_thermodynamics_parameters.jl similarity index 65% rename from src/OceanSeaIceModels/PrescribedAtmospheres.jl rename to src/AtmosphereSimulations/atmosphere_thermodynamics_parameters.jl index 49d72d32d..aa1ba4da0 100644 --- a/src/OceanSeaIceModels/PrescribedAtmospheres.jl +++ b/src/AtmosphereSimulations/atmosphere_thermodynamics_parameters.jl @@ -1,17 +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.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) @@ -40,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 ##### @@ -286,167 +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 - -@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 3759bca57..42a4c2765 100644 --- a/src/ClimaOcean.jl +++ b/src/ClimaOcean.jl @@ -52,6 +52,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} @@ -82,7 +92,8 @@ end function atmosphere_simulation 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") @@ -119,4 +130,7 @@ using PrecompileTools: @setup_workload, @compile_workload end end +function + + end # module 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_ocean_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_ocean_fluxes.jl index 2292fbf1f..7a5574aaf 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_ocean_fluxes.jl @@ -14,8 +14,10 @@ using ClimaOcean.OceanSeaIceModels: sea_ice_concentration, NoAtmosphereModel, No 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 @@ -34,7 +36,7 @@ computed_sea_ice_ocean_fluxes(::Nothing) = nothing # A generic ocean flux assembler for a coupled model with both an atmosphere and sea ice function compute_net_ocean_fluxes!(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 @@ -58,13 +60,14 @@ function compute_net_ocean_fluxes!(coupled_model, ocean) 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!, @@ -82,6 +85,8 @@ function compute_net_ocean_fluxes!(coupled_model, ocean) atmos_ocean_properties, ocean_properties) + fill_net_fluxes!(coupled_model.ocean, net_ocean_fluxes) + return nothing end 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 215a9bcaa..622277a0c 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, e end initialize!(exchanger::StateExchanger, ::Nothing) = nothing +initialize!(exchanger::StateExchanger, atmosphere) = initialize!(exchanger.atmosphere_exchanger, exchanger.exchange_grid, atmosphere) -function initialize!(exchanger::StateExchanger, atmosphere::PrescribedAtmosphere) +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 a5226076b..7e0af42f7 100644 --- a/src/OceanSeaIceModels/OceanSeaIceModels.jl +++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl @@ -27,7 +27,6 @@ using ClimaSeaIce: SeaIceModel using ClimaSeaIce.SeaIceThermodynamics: melting_temperature using ClimaOcean: stateindex - using KernelAbstractions: @kernel, @index using KernelAbstractions.Extras.LoopInfo: @unroll @@ -43,10 +42,10 @@ import Oceananigans.Simulations: reset!, initialize!, iteration import Oceananigans.TimeSteppers: time_step!, update_state!, time import Oceananigans.Utils: prettytime +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 = Oceananigans.defaults.gravitational_acceleration const default_freshwater_density = 1000 # kg m⁻³ @@ -55,6 +54,8 @@ 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 @@ -82,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 @@ -110,7 +115,7 @@ include("freezing_limited_ocean_temperature.jl") # TODO: import this last include("PrescribedAtmospheres.jl") -using .PrescribedAtmospheres: +using ClimaOcean.AtmosphereSimulations: PrescribedAtmosphere, AtmosphereThermodynamicsParameters, TwoBandDownwellingRadiation diff --git a/src/OceanSeaIceModels/ocean_sea_ice_model.jl b/src/OceanSeaIceModels/ocean_sea_ice_model.jl index 06f849793..5fc27c4b3 100644 --- a/src/OceanSeaIceModels/ocean_sea_ice_model.jl +++ b/src/OceanSeaIceModels/ocean_sea_ice_model.jl @@ -56,25 +56,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/ocean_simulation.jl b/src/OceanSimulations/ocean_simulation.jl index cd2094cae..af3a96eeb 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 new file mode 100644 index 000000000..70276cd83 --- /dev/null +++ b/src/OceanSimulations/prescribed_ocean.jl @@ -0,0 +1,83 @@ +using Oceananigans.Models: AbstractModel +using Oceananigans.Fields: ZeroField +using Oceananigans.OutputReaders: extract_field_time_series, update_field_time_series! + +import Oceananigans.TimeSteppers: time_step! + +import ClimaOcean: reference_density, heat_capacity +import Oceananigans.Architectures: on_architecture + +##### +##### A prescribed ocean... +##### + +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 + +""" + 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`: 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(grid; + velocities = default_ocean_velocities(grid), + tracers = default_ocean_tracers(grid), + clock = Clock{Float64}(time = 0), + reference_density = 1029, + heat_capacity = 3998) + + reference_density = convert(eltype(grid), reference_density) + heat_capacity = convert(eltype(grid), heat_capacity) + + return PrescribedOcean(grid, + clock, + 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 +##### + +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 + + return nothing +end + + +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 94% rename from src/SeaIceSimulations.jl rename to src/SeaIceSimulations/SeaIceSimulations.jl index 4b45f52fe..8342f82da 100644 --- a/src/SeaIceSimulations.jl +++ b/src/SeaIceSimulations/SeaIceSimulations.jl @@ -19,6 +19,11 @@ using ClimaSeaIce.Rheologies: IceStrength, ElastoViscoPlasticRheology using ClimaOcean.OceanSimulations: Default +import ClimaOcean: reference_density, 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 + g_Earth = Oceananigans.defaults.gravitational_acceleration Ω_Earth = Oceananigans.defaults.planet_rotation_rate diff --git a/test/test_ocean_sea_ice_model.jl b/test/test_ocean_sea_ice_model.jl index 6d8edcdbc..c0549f3ae 100644 --- a/test/test_ocean_sea_ice_model.jl +++ b/test/test_ocean_sea_ice_model.jl @@ -104,4 +104,23 @@ using ClimaSeaIce.Rheologies end 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