Rename the Breeze seam for consistency: base_pressure and ρE

ext/NumericalEarthBreezeExt/breeze_atmosphere_interface.jl

@@ -127,7 +127,7 @@ function NumericalEarth.EarthSystemModels.interpolate_state!(exchanger, exchange
     qᵛ = specific_humidity(atmosphere)
 
     # Breeze's diagnosed vapor mass fraction, not the scheme-dependent moisture prognostic;
-    # `dynamics_pressure` gives the per-column pressure (a single scalar `surface_pressure`
+    # `dynamics_pressure` gives the per-column pressure (a single scalar `base_pressure`
     # would bias fluxes over terrain).
     p = dynamics_pressure(atmosphere.dynamics)
 
@@ -160,12 +160,12 @@ function NumericalEarth.EarthSystemModels.InterfaceComputations.net_fluxes(atmos
         "flux field, so its surface stress cannot come from the coupler. Build the atmosphere " *
         "without `bottom_drag_coefficient` (Breeze `BulkDrag`) when coupling to land or ocean."))
 
-    ρe = thermodynamic_density(atmosphere.formulation).boundary_conditions.bottom.condition.condition
+    ρE = thermodynamic_density(atmosphere.formulation).boundary_conditions.bottom.condition.condition
 
     # Moisture flux field
     ρqᵛᵉ = atmosphere.moisture_density.boundary_conditions.bottom.condition
 
-    return (; ρu, ρv, ρe, ρqᵛᵉ)
+    return (; ρu, ρv, ρE, ρqᵛᵉ)
 end
 
 NumericalEarth.EarthSystemModels.InterfaceComputations.net_fluxes(atmos::BreezeAtmosphereSim) =
@@ -184,7 +184,7 @@ NumericalEarth.EarthSystemModels.InterfaceComputations.net_fluxes(atmos::BreezeA
         # interpolate stresses on variable's location
         net.ρu[i, j, 1]  = ℑxᶠᵃᵃ(i, j, 1, grid, ao_fluxes.x_momentum)
         net.ρv[i, j, 1]  = ℑyᵃᶠᵃ(i, j, 1, grid, ao_fluxes.y_momentum)
-        net.ρe[i, j, 1]  = Qc   # sensible heat only; latent heat handled by moisture flux
+        net.ρE[i, j, 1]  = Qc   # sensible heat only; latent heat handled by moisture flux
         net.ρqᵛᵉ[i, j, 1] = Fv
     end
 end

ext/NumericalEarthBreezeExt/breeze_atmosphere_simulation.jl

@@ -40,10 +40,10 @@ end
 
 """
     atmosphere_model(grid;
-                     surface_pressure = 101325,
+                     base_pressure = 101325,
                      potential_temperature = 285,
                      thermodynamic_constants = ThermodynamicConstants(eltype(grid)),
-                     dynamics = CompressibleDynamics(; base_pressure = surface_pressure,
+                     dynamics = CompressibleDynamics(; base_pressure,
                                                      reference_potential_temperature = potential_temperature),
                      microphysics = SaturationAdjustment(equilibrium = WarmPhaseEquilibrium()),
                      momentum_advection = WENO(order=9),
@@ -61,7 +61,7 @@ the role of [`ocean_simulation`](@ref)).
 
 When `initialize` (the default) and `dynamics isa CompressibleDynamics`, the returned model is
 set to a resting, hydrostatically balanced state at the reference `potential_temperature` and
-`surface_pressure`, so its density is valid for anything that divides by it (e.g. the MOST
+`base_pressure`, so its density is valid for anything that divides by it (e.g. the MOST
 surface-flux coupling). `AnelasticDynamics` needs no such step: its `reference_state` already
 prescribes a valid resting density, and zero prognostic perturbation is by construction the
 resting hydrostatically balanced state. Pass `initialize = false` when a caller derives the full
@@ -86,10 +86,10 @@ extension), which derives the lateral BCs and Davies relaxation from the parent
 in a `NestedModel`.
 """
 function NumericalEarth.Atmospheres.atmosphere_model(grid;
-                                                     surface_pressure = 101325,
+                                                     base_pressure = 101325,
                                                      potential_temperature = 285,
                                                      thermodynamic_constants = ThermodynamicConstants(eltype(grid)),
-                                                     dynamics = CompressibleDynamics(; base_pressure = surface_pressure,
+                                                     dynamics = CompressibleDynamics(; base_pressure,
                                                                                      reference_potential_temperature = potential_temperature),
                                                      microphysics = SaturationAdjustment(equilibrium = WarmPhaseEquilibrium()),
                                                      momentum_advection = Oceananigans.WENO(order=9),
@@ -111,12 +111,12 @@ function NumericalEarth.Atmospheres.atmosphere_model(grid;
     # Create 2D coupling-flux fields populated by the ESM coupler each step.
     ρτˣ = Field{Center, Center, Nothing}(grid)
     ρτʸ = Field{Center, Center, Nothing}(grid)
-    Jᵉ  = Field{Center, Center, Nothing}(grid)
+    Jᴱ  = Field{Center, Center, Nothing}(grid)
     Jᵛ  = Field{Center, Center, Nothing}(grid)
 
     moisture_key = moisture_prognostic_name(microphysics)
     moisture_bc = NamedTuple{tuple(moisture_key)}(tuple(FieldBoundaryConditions(bottom = FluxBoundaryCondition(Jᵛ))))
-    energy_bc = NamedTuple{(energy_bc_key(),)}((FieldBoundaryConditions(bottom = FluxBoundaryCondition(Jᵉ)),))
+    energy_bc = NamedTuple{(energy_bc_key(),)}((FieldBoundaryConditions(bottom = FluxBoundaryCondition(Jᴱ)),))
 
     momentum_bcs = (
         ρu = FieldBoundaryConditions(bottom = FluxBoundaryCondition(ρτˣ)),
@@ -150,7 +150,7 @@ function NumericalEarth.Atmospheres.atmosphere_model(grid;
     # applying this compressible-only set! corrupts it in place (Breeze's set_to_mean.jl copies
     # between mismatched-shape Fields), NaN-ing every field that divides by density.
     initialize && dynamics isa CompressibleDynamics &&
-        set!(model; θ = potential_temperature, ρ = HydrostaticallyBalancedDensity(; surface_pressure))
+        set!(model; θ = potential_temperature, ρ = HydrostaticallyBalancedDensity())
 
     return model
 end

ext/NumericalEarthBreezeExt/breeze_nested_atmosphere.jl

@@ -119,13 +119,13 @@ default_lid_depth(grid) = convert(eltype(grid), grid.Lz / 4)
 # Default child dynamics: compressible with split-explicit acoustic substepping, an `UpperSponge`
 # Rayleigh layer over the top `damping_depth` meters at `damping_rate`, and no divergence damping
 # (its (ρθ)′-proxy damper injects a spurious force on an unbalanced cold start). When given,
-# `surface_pressure`/`reference_potential_temperature` anchor the hydrostatic reference and the
+# `base_pressure`/`reference_potential_temperature` anchor the hydrostatic reference and the
 # perturbation-form pressure-gradient reference profile.
-function default_nested_dynamics(grid; surface_pressure, reference_potential_temperature, damping_rate, damping_depth)
+function default_nested_dynamics(grid; base_pressure, reference_potential_temperature, damping_rate, damping_depth)
     time_discretization = SplitExplicitTimeDiscretization(sponge = UpperSponge(; damping_rate, depth = damping_depth),
                                                           damping = NoDivergenceDamping())
     kw = (;)
-    isnothing(surface_pressure)                || (kw = merge(kw, (; base_pressure = surface_pressure)))
+    isnothing(base_pressure)                   || (kw = merge(kw, (; base_pressure)))
     isnothing(reference_potential_temperature) || (kw = merge(kw, (; reference_potential_temperature)))
     return CompressibleDynamics(time_discretization; kw...)
 end
@@ -192,7 +192,7 @@ Provides sensible, overridable physics defaults: `microphysics` (1-moment mixed-
 `CloudMicrophysics` is loaded), `momentum_advection = WENO(order=9)`, `coriolis = SphericalCoriolis()`,
 and a compressible split-explicit `dynamics` with an `UpperSponge` over the top `damping_depth` m at
 `damping_rate`; a matching ρw Rayleigh lid sponge (`Relaxation` toward zero) is added to `forcing`. Pass
-`surface_pressure`/`reference_potential_temperature` to anchor the default dynamics. Any
+`base_pressure`/`reference_potential_temperature` to anchor the default dynamics. Any
 `boundary_conditions`/`forcing` the caller passes are merged with the parent-derived ones (caller wins).
 
 When `bottom_drag_coefficient` is given — a constant drag coefficient or a `Breeze.PolynomialCoefficient`
@@ -221,7 +221,7 @@ function NumericalEarth.NestedModels.nested_atmosphere_model(parent_atmosphere::
     relaxation_mask = davies_relaxation_mask(child_grid, relaxation_width),
     sides = (:west, :east, :south, :north),
     thermodynamic_constants = ThermodynamicConstants(eltype(child_grid)),
-    surface_pressure = nothing,
+    base_pressure = nothing,
     reference_potential_temperature = nothing,
     terrain = nothing,
     terrain_blend_length = 60_000,   # meters; physical blend width → resolution-invariant slope
@@ -236,7 +236,7 @@ function NumericalEarth.NestedModels.nested_atmosphere_model(parent_atmosphere::
     coriolis = SphericalCoriolis(),
     damping_rate = 1/5,
     damping_depth = default_lid_depth(child_grid),
-    dynamics = default_nested_dynamics(child_grid; surface_pressure, reference_potential_temperature, damping_rate, damping_depth),
+    dynamics = default_nested_dynamics(child_grid; base_pressure, reference_potential_temperature, damping_rate, damping_depth),
     boundary_conditions = NamedTuple(),
     forcing = NamedTuple(),
     kw...)
@@ -368,7 +368,7 @@ Build the parent `PrescribedAtmosphere`, nest a Breeze child in it, and initiali
 `parent_dataset` at `first(dates)` — the returned model is ready to step. The parent spans
 `child_grid`'s bounding box padded by `parent_padding` (default `parent_dataset`'s
 `default_horizontal_padding`, margin for the lateral-BC interpolation stencils) at `dates`, on
-`parent_dataset`'s native grid. Unless given, the default dynamics' `surface_pressure` anchor is the domain-mean dataset
+`parent_dataset`'s native grid. Unless given, the default dynamics' `base_pressure` anchor is the domain-mean dataset
 mean-sea-level pressure over the child at `first(dates)`. When `bottom_drag_coefficient` is given,
 `drag_surface_temperature` defaults to the dataset's skin temperature at `first(dates)` regridded onto
 the child grid (a static snapshot, not the dataset's diurnal cycle). `balancer` controls the
@@ -380,7 +380,7 @@ function NumericalEarth.NestedModels.nested_atmosphere_model(child_grid, parent_
     dir = default_download_directory(parent_dataset),
     parent_padding = default_horizontal_padding(parent_dataset),
     parent_time_indices_in_memory = nothing,   # nothing ⇒ every date resident; ≥3 streams a moving window
-    surface_pressure = nothing,
+    base_pressure = nothing,
     bottom_drag_coefficient = nothing,
     drag_surface_temperature = nothing,
     balancer = true,
@@ -391,15 +391,15 @@ function NumericalEarth.NestedModels.nested_atmosphere_model(child_grid, parent_
                                              architecture = architecture(child_grid), dir,
                                              time_indices_in_memory = parent_time_indices_in_memory)
 
-    if isnothing(surface_pressure)
-        surface_pressure = mean_sea_level_pressure(parent_dataset, child_grid, first(dates), dir)
+    if isnothing(base_pressure)
+        base_pressure = mean_sea_level_pressure(parent_dataset, child_grid, first(dates), dir)
     end
 
     if !isnothing(bottom_drag_coefficient) && isnothing(drag_surface_temperature)
         drag_surface_temperature = dataset_skin_temperature(parent_dataset, child_grid, first(dates), dir)
     end
 
-    nested_model = NumericalEarth.NestedModels.nested_atmosphere_model(parent_atmosphere, child_grid; surface_pressure,
+    nested_model = NumericalEarth.NestedModels.nested_atmosphere_model(parent_atmosphere, child_grid; base_pressure,
                                                                        bottom_drag_coefficient, drag_surface_temperature, kw...)
     initialize_nested_child!(nested_model, parent_dataset, first(dates), dir; balancer)
     return nested_model

test/test_breeze_coupling.jl

@@ -426,3 +426,28 @@ end
         end
     end
 end
+
+# `base_pressure` is the datum at z = 0, which `reference_state` reduces to each column's surface.
+@testset "Cold start agrees with its own reference on a raised domain" begin
+    for arch in test_architectures
+        A = typeof(arch)
+
+        @testset "on $A" begin
+            p₀ = 101325
+            grid = RectilinearGrid(arch; size = (8, 20), halo = (5, 5),
+                                   x = (0, 10kilometers), z = (2kilometers, 6kilometers),
+                                   topology = (Periodic, Flat, Bounded))
+
+            model = atmosphere_model(grid; base_pressure = p₀)
+
+            p  = Array(interior(model.dynamics.pressure))
+            pᵣ = Array(interior(model.dynamics.reference_state.pressure))
+
+            @test all(p .≈ pᵣ)
+
+            # Both anchored at the datum would also agree, so pin that the reference really is
+            # reduced to the domain bottom: 2 km of hydrostatic descent is ~22 kPa.
+            @test maximum(p[:, :, 1]) < 0.9 * p₀
+        end
+    end
+end

test/test_nested_simulation.jl

@@ -812,7 +812,7 @@ end
     # Default microphysics: with CloudMicrophysics unloaded (as in the test env) this is the Breeze-native
     # `SaturationAdjustment(WarmPhaseEquilibrium())` — no extra dependency needed to exercise the seam.
     model = nested_atmosphere_model(parent, child_grid;
-                relaxation_rate = 1/300, relaxation_width = 3, surface_pressure = 1e5,
+                relaxation_rate = 1/300, relaxation_width = 3, base_pressure = 1e5,
                 coriolis = nothing, terrain = nothing, parent_condensates = nothing)
     ext.initialize_nested_child!(model, nothing, first(times), ""; balancer = false)
 
@@ -852,7 +852,7 @@ end
     set!(terrain, (λ, φ) -> 1000 * exp(-(λ^2 + (φ - 36.6)^2) / 0.25))
 
     model = nested_atmosphere_model(parent, child_grid; terrain, terrain_smoothing_passes = 0,
-                                    surface_pressure = 1e5, coriolis = nothing, parent_condensates = nothing)
+                                    base_pressure = 1e5, coriolis = nothing, parent_condensates = nothing)
     ext.initialize_nested_child!(model, nothing, 0.0, ""; balancer = false)
 
     expected_u = XFaceField(child_grid); set!(expected_u, (λ, φ, z) -> parent_u(z))