Add a prescribed ocean to the src

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...
 #####

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

src/OceanSeaIceModels/PrescribedAtmospheres.jl → 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