Introduce an interface to allow attaching a prognostic atmosphere

src/OceanSeaIceModels/Atmospheres/Atmospheres.jl

@@ -0,0 +1,20 @@
+module Atmospheres
+
+export PrescribedAtmosphere
+
+using KernelAbstractions: @kernel, @index
+using Oceananigans.Utils: KernelParameters
+using Oceananigans.Grids: grid_name, architecture
+using Oceananigans.Utils: prettysummary
+using Oceananigans.Fields: Center
+using Oceananigans.OutputReaders: FieldTimeSeries, update_field_time_series!, extract_field_time_series
+
+include("atmospheric_parameters.jl")
+include("prescribed_atmospheres.jl")
+
+# Need to be extended by atmospheric models
+surface_layer_height(atmos::PrescribedAtmosphere)      = atmos.reference_height
+boundary_layer_height(atmos::PrescribedAtmosphere)     = atmos.boundary_layer_height
+thermodynamics_parameters(atmos::PrescribedAtmosphere) = atmos.thermodynamics_parameters
+
+end

src/OceanSeaIceModels/PrescribedAtmospheres.jl → src/OceanSeaIceModels/Atmospheres/atmospheric_parameters.jl

@@ -1,10 +1,3 @@
-module PrescribedAtmospheres
-
-using Oceananigans.Grids: grid_name
-using Oceananigans.Utils: prettysummary
-using Oceananigans.Fields: Center
-using Oceananigans.OutputReaders: FieldTimeSeries, update_field_time_series!, extract_field_time_series
-
 using Adapt
 using Thermodynamics.Parameters: AbstractThermodynamicsParameters
 
@@ -38,10 +31,6 @@ import Thermodynamics.Parameters:
     pow_icenuc      # "Power parameter" that controls liquid/ice condensate partitioning
                     # during partial ice nucleation
 
-import ..OceanSeaIceModels:
-    downwelling_radiation,
-    freshwater_flux
-
 #####
 ##### Atmospheric thermodynamics parameters
 #####
@@ -282,147 +271,3 @@ const PATP = PrescribedAtmosphereThermodynamicsParameters
 @inline cp_d(p::PATP)           = R_d(p) / kappa_d(p)
 @inline cv_d(p::PATP)           = cp_d(p) - R_d(p)
 @inline cv_v(p::PATP)           = cp_v(p) - R_v(p)
-
-#####
-##### Prescribed atmosphere (as opposed to dynamically evolving / prognostic)
-#####
-
-struct PrescribedAtmosphere{FT, M, G, U, P, C, F, I, R, TP, TI}
-    grid :: G
-    metadata :: M
-    velocities :: U
-    pressure :: P
-    tracers :: C
-    freshwater_flux :: F
-    auxiliary_freshwater_flux :: I
-    downwelling_radiation :: R
-    thermodynamics_parameters :: TP
-    times :: TI
-    reference_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, "├── reference_height: ", prettysummary(pa.reference_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
-
-function default_atmosphere_pressure(grid, times)
-    pa = FieldTimeSeries{Center, Center, Nothing}(grid, times)
-    parent(pa) .= 101325
-    return pa
-end
-
-"""
-    PrescribedAtmosphere(grid, times;
-                         metadata = nothing,
-                         reference_height = 10, # meters
-                         boundary_layer_height = 600 # meters,
-                         thermodynamics_parameters = PrescribedAtmosphereThermodynamicsParameters(FT),
-                              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;
-                              metadata = nothing,  
-                              reference_height = convert(eltype(grid), 10),
-                              boundary_layer_height = convert(eltype(grid), 600),
-                              thermodynamics_parameters = nothing,
-                              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 = PrescribedAtmosphereThermodynamicsParameters(FT)
-    end
-
-    return PrescribedAtmosphere(grid,
-                                metadata,
-                                velocities,
-                                pressure,
-                                tracers,
-                                freshwater_flux,
-                                auxiliary_freshwater_flux,
-                                downwelling_radiation,
-                                thermodynamics_parameters,
-                                times,
-                                convert(FT, reference_height),
-                                convert(FT, boundary_layer_height))
-end
-
-update_model_field_time_series!(::Nothing, time) = nothing
-
-function update_model_field_time_series!(atmos::PrescribedAtmosphere, time)
-    ftses = extract_field_time_series(atmos)
-    for fts in ftses
-        update_field_time_series!(fts, time)
-    end
-
-    return nothing
-end
-
-struct TwoBandDownwellingRadiation{SW, LW}
-    shortwave :: SW
-    longwave :: LW
-end
-
-"""
-    TwoBandDownwellingRadiation(shortwave=nothing, longwave=nothing)
-
-Return a two-band model for downwelling radiation (split in 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
-