Clean implementation of Kessler microphysics

examples/supercell.jl

@@ -61,16 +61,9 @@ using Oceananigans.Grids: znode
 using Oceananigans: Center, Face
 using Oceananigans.Operators: Δzᶜᶜᶜ, Δzᶜᶜᶠ, ℑzᵃᵃᶠ
 using Breeze.Thermodynamics: dry_air_gas_constant
+using Breeze.Microphysics: KesslerMicrophysics
 using CUDA
 
-using CloudMicrophysics
-import Breeze: Breeze
-
-# Access extension module and define aliases to avoid namespace conflicts:
-
-const BreezeCloudMicrophysicsExt = Base.get_extension(Breeze, :BreezeCloudMicrophysicsExt)
-const BreezeOneMomentCloudMicrophysics = BreezeCloudMicrophysicsExt.OneMomentCloudMicrophysics
-
 # ## Grid configuration
 #
 # The domain is 168 km × 168 km × 20 km with 168 × 168 × 40 grid points, giving
@@ -237,10 +230,10 @@ fig_bubble
 
 # ## Model initialization
 #
-# Create the atmosphere model with one-moment cloud microphysics from CloudMicrophysics.jl,
+# Create the atmosphere model with Kessler warm-rain microphysics,
 # high-order WENO advection, and anisotropic minimum dissipation turbulence closure:
 
-microphysics = BreezeOneMomentCloudMicrophysics()
+microphysics = KesslerMicrophysics()
 advection = WENO(order=9, minimum_buffer_upwind_order=3)
 closure = AnisotropicMinimumDissipation()
 model = AtmosphereModel(grid; dynamics, closure, microphysics, advection)
@@ -293,25 +286,26 @@ set!(θᵇᵍf, (x, y, z) -> θᵇᵍ(x, y, z))
 θ′ = θ - θᵇᵍf
 
 # Extract microphysical fields for output:
+# Kessler scheme has cloud liquid (qᶜˡ), rain (qʳ), and diagnosed vapor (qᵛ)
 
 qᶜˡ = model.microphysical_fields.qᶜˡ
-qᶜⁱ = model.microphysical_fields.qᶜⁱ
+qʳ = model.microphysical_fields.qʳ
 qᵛ = model.microphysical_fields.qᵛ
 
 # ## Simulation
 #
 # Run for 2 hours with adaptive time stepping (CFL = 0.7) starting from Δt = 2 s:
 
 simulation = Simulation(model; Δt=2, stop_time=2hours)
-conjure_time_step_wizard!(simulation, cfl=0.7)
+#conjure_time_step_wizard!(simulation, cfl=0.7)
 
 # Progress callback to monitor simulation health:
 
 function progress(sim)
     u, v, w = sim.model.velocities
     qᵛ = model.microphysical_fields.qᵛ
     qᶜˡ = model.microphysical_fields.qᶜˡ
-    qᶜⁱ = model.microphysical_fields.qᶜⁱ
+    qʳ = model.microphysical_fields.qʳ
 
     ρe = Breeze.AtmosphereModels.static_energy_density(sim.model)
     ρemean = mean(ρe)
@@ -321,8 +315,8 @@ function progress(sim)
                    maximum(abs, u), maximum(w), minimum(w))
     @info msg
 
-    msg *= @sprintf(", max(qᵛ): %.5e, max(qᶜˡ): %.5e, max(qᶜⁱ): %.5e",
-                    maximum(qᵛ), maximum(qᶜˡ), maximum(qᶜⁱ))
+    msg *= @sprintf(", max(qᵛ): %.5e, max(qᶜˡ): %.5e, max(qʳ): %.5e",
+                    maximum(qᵛ), maximum(qᶜˡ), maximum(qʳ))
     @info msg
     return nothing
 end
@@ -333,7 +327,7 @@ add_callback!(simulation, progress, IterationInterval(100))
 #
 # Save full 3D fields for post-processing analysis:
 
-outputs = merge(model.velocities, model.tracers, (; θ, θ′, qᶜˡ, qᶜⁱ, qᵛ))
+outputs = merge(model.velocities, model.tracers, (; θ, θ′, qᶜˡ, qʳ, qᵛ))
 
 filename = "supercell.jld2"
 ow = JLD2Writer(model, outputs; filename,
@@ -366,8 +360,8 @@ save("max_w_timeseries.png", fig)
 #
 # We create a 3-panel animation showing the storm structure at mid-levels (z ≈ 5 km):
 # - Vertical velocity ``w``: reveals the updraft/downdraft structure
-# - Cloud water ``q^{cl}``: shows the cloud boundaries
-# - Rain water ``q^r``: indicates precipitation regions
+# - Cloud liquid water ``qᶜˡ``: shows the cloud boundaries
+# - Rain water ``qʳ``: indicates precipitation regions
 
 wxy_ts = FieldTimeSeries("supercell_slices.jld2", "wxy")
 qʳxy_ts = FieldTimeSeries("supercell_slices.jld2", "qʳxy")
@@ -380,14 +374,14 @@ Nt = length(times)
 
 wlim = 25       # m/s - vertical velocity range
 qʳlim = 0.01    # kg/kg - rain water range
-qᶜˡlim = 0.001  # kg/kg - cloud water range
+qᶜˡlim = 0.001   # kg/kg - cloud liquid water range
 
 # Create the figure with 3 panels:
 
 slices_fig = Figure(size=(900, 1000), fontsize=14)
 
 axw = Axis(slices_fig[1, 1], xlabel="x (m)", ylabel="y (m)", title="Vertical velocity w")
-axqᶜˡ = Axis(slices_fig[1, 2], xlabel="x (m)", ylabel="y (m)", title="Cloud water qᶜˡ")
+axqᶜˡ = Axis(slices_fig[1, 2], xlabel="x (m)", ylabel="y (m)", title="Cloud liquid water qᶜˡ")
 axqʳ = Axis(slices_fig[3, 1], xlabel="x (m)", ylabel="y (m)", title="Rain water qʳ")
 
 # Set up observables for animation:

src/AtmosphereModels/atmosphere_model.jl

@@ -189,7 +189,7 @@ function AtmosphereModel(grid;
                                        grid, coriolis, density,
                                        velocities, dynamics, formulation, specific_moisture)
 
-    # Include thermodynamic density (ρe or ρθ), ρqᵗ plus user tracers for closure field construction
+    # Include thermodynamic density (ρe or ρθ), ρqᵗ, microphysical prognostic fields, plus user tracers
     closure_thermo_name = thermodynamic_density_name(formulation)
     microphysical_names = prognostic_field_names(microphysics)
     scalar_names = tuple(closure_thermo_name, :ρqᵗ, microphysical_names..., tracer_names...)

src/Breeze.jl

@@ -49,6 +49,8 @@ export
     RelativeHumidity,
     RelativeHumidityField,
     BulkMicrophysics,
+    KesslerMicrophysics,
+    compute_hydrostatic_pressure!,
     NonEquilibriumCloudFormation,
 
     # BoundaryConditions

src/Microphysics/Microphysics.jl

@@ -14,13 +14,18 @@ export
     SaturationSpecificHumidity,
     SaturationSpecificHumidityField,
     RelativeHumidity,
-    RelativeHumidityField
+    RelativeHumidityField,
+    KesslerMicrophysics
 
-using ..AtmosphereModels: AtmosphereModels, compute_moisture_fractions,
-    materialize_microphysical_fields, update_microphysical_fields!
+import ..AtmosphereModels: compute_moisture_fractions,
+    materialize_microphysical_fields, update_microphysical_fields!,
+    maybe_adjust_thermodynamic_state
+
+using ..AtmosphereModels: AtmosphereModels
 
 include("saturation_adjustment.jl")
 include("bulk_microphysics.jl")
+include("kessler_microphysics.jl")
 include("microphysics_diagnostics.jl")
 
 end # module Microphysics