Revisit flux climatology

src/DataWrangling/DataWrangling.jl

@@ -20,7 +20,7 @@ using KernelAbstractions: @kernel, @index
 using Oceananigans.DistributedComputations
 using Adapt
 
-const dataset_defaults = Oceananigans.Defaults(FloatType=Float32)
+import Oceananigans.Fields: set!
 
 #####
 ##### Downloading utilities

src/DataWrangling/metadata.jl

@@ -286,59 +286,3 @@ function compute_native_date_range(native_dates, start_date, end_date)
 
     return native_dates[start_idx:end_idx]
 end
-
-"""
-    vertical_interfaces(metadata)
-
-Return an array with the vertical interfaces (``z``-faces) of the dataset
-that `metadata` corresponds to.
-"""
-vertical_interfaces(metadata::Metadata{V}) where V =
-    error("vertical_interfaces not implemented for $V")
-
-variable_is_three_dimensional(metadata::Metadata{V}) where V =
-    error("variable_is_three_dimensional not implemented for $V")
-
-function dataset_latitude_extent end
-
-"""
-    empty_field(metadata::Metadata;
-                architecture = CPU(),
-                horizontal_halo = (7, 7))
-
-Return an empty field of `metadata` on `architecture` and with `horizontal_halo`s.
-"""
-function empty_field(metadata::Metadata, FT=dataset_defaults.FloatType;
-                     architecture = CPU(),
-                     horizontal_halo = (7, 7))
-
-    Nx, Ny, Nz, _ = size(metadata)
-    loc = location(metadata)
-    longitude = (0, 360)
-    latitude = dataset_latitude_extent(metadata)
-    TX, TY = (Periodic, Bounded)
-
-    if variable_is_three_dimensional(metadata)
-        TZ = Bounded
-        LZ = Center
-        z = vertical_interfaces(metadata)
-        halo = (horizontal_halo..., 3)
-        sz = (Nx, Ny, Nz)
-    else # the variable is two-dimensional
-        TZ = Flat
-        LZ = Nothing
-        z = nothing
-        halo = horizontal_halo
-        sz = (Nx, Ny)
-    end
-
-    grid = LatitudeLongitudeGrid(architecture, FT; halo, longitude, latitude, z,
-                                 size = sz,
-                                 topology = (TX, TY, TZ))
-    return Field{loc...}(grid)
-end
-
-struct Celsius end
-struct Kelvin end
-
-function dataset_temperature_units end

src/OceanSeaIceModels/InterfaceComputations/radiation.jl

@@ -110,6 +110,5 @@ end
 
 # Inside the solver we lose both spatial and temporal information, but the
 # radiative properties have already been computed correctly
-@inline net_downwelling_radiation(Qs, Qℓ, α, ϵ) = - (1 - α) * Qs - ϵ * Qℓ
-
-@inline upwelling_radiation(T, σ, ϵ) = σ * ϵ * T^4
+@inline net_absorbed_interface_radiation(Qs, Qℓ, α, ϵ) = - (1 - α) * Qs - ϵ * Qℓ
+@inline emitted_longwave_radiation(T, σ, ϵ) = σ * ϵ * T^4

src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl

@@ -87,9 +87,11 @@ Keyword Arguments
 function SimilarityTheoryFluxes(FT::DataType = Oceananigans.defaults.FloatType;
                                 von_karman_constant = 0.4,
                                 turbulent_prandtl_number = 1,
-                                gustiness_parameter = 6.5,
-                                stability_functions = edson_stability_functions(FT),
-                                roughness_lengths = default_roughness_lengths(FT),
+                                gustiness_parameter = 1,
+                                stability_functions = atmosphere_ocean_stability_functions(FT),
+                                momentum_roughness_length = MomentumRoughnessLength(FT),
+                                temperature_roughness_length = ScalarRoughnessLength(FT),
+                                water_vapor_roughness_length = ScalarRoughnessLength(FT),
                                 similarity_form = LogarithmicSimilarityProfile(),
                                 solver_stop_criteria = nothing,
                                 solver_tolerance = 1e-8,
@@ -328,16 +330,16 @@ stability function for _stable_ or _unstable_ atmospheric conditions, respective
 """
 @kwdef struct EdsonMomentumStabilityFunction{FT} <: AbstractStabilityFunction
     ζmax :: FT = 50.0
-    Aˢ   :: FT = 0.35
-    Bˢ   :: FT = 0.7
-    Cˢ   :: FT = 0.75
-    Dˢ   :: FT = 5 / 0.35
-    Aᵘ   :: FT = 15.0
-    Bᵘ   :: FT = 2.0
-    Cᵘ   :: FT = π/2
-    Dᵘ   :: FT = 10.15
-    Eᵘ   :: FT = 3.0
-    Fᵘ   :: FT = π / sqrt(3)
+    A⁺   :: FT = 0.35
+    B⁺   :: FT = 0.7
+    C⁺   :: FT = 0.75
+    D⁺   :: FT = 5/0.35
+    A⁻   :: FT = 15.0
+    B⁻   :: FT = 2.0
+    C⁻   :: FT = π/2
+    D⁻   :: FT = 10.15
+    E⁻   :: FT = 3.0
+    F⁻   :: FT = π / sqrt(3)
 end
 
 @inline function stability_profile(ψ::EdsonMomentumStabilityFunction, ζ)
@@ -355,7 +357,7 @@ end
 
     ζ⁻ = min(zero(ζ), ζ)
     ζ⁺ = max(zero(ζ), ζ)
-    dζ = min(ζmax, Aˢ * ζ⁺) # ok
+    dζ = min(ζmax, A⁺ * ζ⁺)
 
     # Stability parameter for _stable_ atmospheric conditions
     ψ⁺ = - B⁺ * ζ⁺ - C⁺ * (ζ⁺ - D⁺) * exp(- dζ) - C⁺ * D⁺