diff --git a/examples/supercell.jl b/examples/supercell.jl index 89a901db..c4f445ba 100644 --- a/examples/supercell.jl +++ b/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,9 +286,10 @@ 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 @@ -303,7 +297,7 @@ qᵛ = model.microphysical_fields.qᵛ # 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: @@ -311,7 +305,7 @@ 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: diff --git a/src/AtmosphereModels/atmosphere_model.jl b/src/AtmosphereModels/atmosphere_model.jl index 83e899cd..8d125e5d 100644 --- a/src/AtmosphereModels/atmosphere_model.jl +++ b/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...) diff --git a/src/Breeze.jl b/src/Breeze.jl index 92fc486e..067cce50 100644 --- a/src/Breeze.jl +++ b/src/Breeze.jl @@ -49,6 +49,8 @@ export RelativeHumidity, RelativeHumidityField, BulkMicrophysics, + KesslerMicrophysics, + compute_hydrostatic_pressure!, NonEquilibriumCloudFormation, # BoundaryConditions diff --git a/src/Microphysics/Microphysics.jl b/src/Microphysics/Microphysics.jl index f778733b..fe454d75 100644 --- a/src/Microphysics/Microphysics.jl +++ b/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 diff --git a/src/Microphysics/kessler_microphysics.jl b/src/Microphysics/kessler_microphysics.jl new file mode 100644 index 00000000..8516717d --- /dev/null +++ b/src/Microphysics/kessler_microphysics.jl @@ -0,0 +1,708 @@ +""" +Kessler warm-rain bulk microphysics scheme. + +A "warm-rain" (Kessler-type) bulk microphysics scheme with water vapor, cloud liquid, and rain. + +Breeze uses mass fractions (q = mass_species / mass_total), while Kessler formulas use +mixing ratios (r = mass_species / mass_dry_air). Conversion: +- r = q / (1 - qᵗ) where qᵗ is total moisture mass fraction +- q = r * (1 - qᵗ) + +Prognostic variables (in Breeze mass fraction form): +- qᶜˡ: cloud liquid water mass fraction +- qʳ: rain water mass fraction + +Diagnostic variable: +- qᵛ: water vapor mass fraction = qᵗ - qᶜˡ - qʳ (from Breeze's total moisture qᵗ) + +Reference: Kessler (1969), "On the Distribution and Continuity of Water Substance in Atmospheric Circulations" +""" + +using Oceananigans: Oceananigans, CenterField, Field, Center, Face +using Oceananigans.BoundaryConditions: FieldBoundaryConditions +using Oceananigans.Fields: ZFaceField, ZeroField +using Oceananigans.Operators: Δzᶜᶜᶜ +using DocStringExtensions: TYPEDSIGNATURES + +import ..AtmosphereModels: + prognostic_field_names, + materialize_microphysical_fields, + microphysical_velocities, + compute_moisture_fractions, + microphysical_tendency, + update_microphysical_fields!, + precipitation_rate, + surface_precipitation_flux, + maybe_adjust_thermodynamic_state + +using ..Thermodynamics: + MoistureMassFractions, + PlanarLiquidSurface, + saturation_specific_humidity, + temperature, + density, + liquid_latent_heat, + mixture_heat_capacity, + total_specific_moisture, + exner_function + +##### +##### Kessler microphysics struct +##### + +""" +$(TYPEDSIGNATURES) + +Kessler warm-rain microphysics scheme with cloud liquid and rain. + +# Fields +- `autoconversion_rate`: Rate constant for autoconversion (cloud → rain), k₁ [s⁻¹]. Default: 0.001 s⁻¹ +- `autoconversion_threshold`: Cloud water threshold for autoconversion, a [kg kg⁻¹]. Default: 0.001 kg kg⁻¹ +- `accretion_rate`: Rate constant for accretion (collection of cloud by rain), k₂ [s⁻¹]. Default: 2.2 s⁻¹ + +Note: The reference density ρ₀ for terminal velocity is obtained from Breeze's reference state +(ρᵣ[i,j,1]) rather than being stored as a parameter. +""" +struct KesslerMicrophysics{FT} + autoconversion_rate :: FT # k₁ [s⁻¹] + autoconversion_threshold :: FT # a [kg kg⁻¹] + accretion_rate :: FT # k₂ [s⁻¹] +end + +Base.summary(::KesslerMicrophysics) = "KesslerMicrophysics" + +function Base.show(io::IO, km::KesslerMicrophysics{FT}) where FT + print(io, "KesslerMicrophysics{$FT}:\n", + "├── autoconversion_rate: ", km.autoconversion_rate, " s⁻¹\n", + "├── autoconversion_threshold: ", km.autoconversion_threshold, " kg kg⁻¹\n", + "└── accretion_rate: ", km.accretion_rate, " s⁻¹") +end + +""" +$(TYPEDSIGNATURES) + +Construct a `KesslerMicrophysics` scheme with default parameters from Kessler (1969). + +# Arguments +- `FT`: Float type to use (defaults to `Oceananigans.defaults.FloatType`) + +# Keyword Arguments +- `autoconversion_rate`: Rate constant k₁ [s⁻¹]. Default: 0.001 s⁻¹ +- `autoconversion_threshold`: Cloud water threshold a [kg kg⁻¹]. Default: 0.001 kg kg⁻¹ +- `accretion_rate`: Rate constant k₂ [s⁻¹]. Default: 2.2 s⁻¹ +""" +function KesslerMicrophysics(FT::DataType = Oceananigans.defaults.FloatType; + autoconversion_rate = 0.001, + autoconversion_threshold = 0.001, + accretion_rate = 2.2) + + return KesslerMicrophysics{FT}(convert(FT, autoconversion_rate), + convert(FT, autoconversion_threshold), + convert(FT, accretion_rate)) +end + +const KM = KesslerMicrophysics + +##### +##### Mass fraction ↔ mixing ratio conversion +##### + +""" +Convert mass fraction q to mixing ratio r. + +r = q / (1 - qᵗ) + +where qᵗ is total moisture mass fraction and (1 - qᵗ) is dry air mass fraction. +""" +@inline function mass_fraction_to_mixing_ratio(q, qᵗ) + qᵈ = 1 - qᵗ # dry air mass fraction + return q / qᵈ +end + +""" +Convert mixing ratio r to mass fraction q. + +q = r * (1 - qᵗ) + +where qᵗ is total moisture mass fraction and (1 - qᵗ) is dry air mass fraction. +Also used to convert mixing ratio tendencies to mass fraction tendencies. +""" +@inline function mixing_ratio_to_mass_fraction(r, qᵗ) + qᵈ = 1 - qᵗ # dry air mass fraction + return r * qᵈ +end + +##### +##### Microphysics interface implementation +##### + +# Only cloud liquid and rain are prognostic; vapor is diagnosed from qᵗ +prognostic_field_names(::KM) = (:ρqᶜˡ, :ρqʳ) + +function materialize_microphysical_fields(km::KM, grid, boundary_conditions) + # Prognostic fields (density-weighted mass fractions) + ρqᶜˡ = CenterField(grid; boundary_conditions=boundary_conditions.ρqᶜˡ) + ρqʳ = CenterField(grid; boundary_conditions=boundary_conditions.ρqʳ) + + # Diagnostic fields (mass fractions) + qᵛ = CenterField(grid) + qᶜˡ = CenterField(grid) + qʳ = CenterField(grid) + + # Cached microphysics rates (computed once per timestep in update_microphysical_fields!) + # These are tendencies in mixing ratio space [kg kg⁻¹ s⁻¹] + Cₖ = CenterField(grid) # Condensation rate + Eₖ = CenterField(grid) # Cloud evaporation rate + Aₖ = CenterField(grid) # Autoconversion rate + Kₖ = CenterField(grid) # Accretion rate + Eʳ = CenterField(grid) # Rain evaporation rate + + # Rain terminal velocity (negative = downward) + # bottom = nothing ensures the kernel-set value is preserved during fill_halo_regions! + wʳ_bcs = FieldBoundaryConditions(grid, (Center(), Center(), Face()); bottom=nothing) + wʳ = ZFaceField(grid; boundary_conditions=wʳ_bcs) + + # Surface precipitation rate (2D field, m/s) + # This is the volume flux of rain at the surface: wʳ * qʳ (positive = precipitation out of domain) + precipitation_rate = Field{Center, Center, Nothing}(grid) + + return (; ρqᶜˡ, ρqʳ, qᵛ, qᶜˡ, qʳ, Cₖ, Eₖ, Aₖ, Kₖ, Eʳ, wʳ, precipitation_rate) +end + +@inline function update_microphysical_fields!(μ, km::KM, i, j, k, grid, ρ, 𝒰, constants) + FT = eltype(grid) + @inbounds begin + # Get total moisture from thermodynamic state + # In the moisture_mass_fractions, vapor contains qᵛ and liquid contains total condensate (qᶜˡ + qʳ) + # But we need to separate qᶜˡ and qʳ from prognostic fields + qᶜˡ = μ.ρqᶜˡ[i, j, k] / ρ + qʳ = μ.ρqʳ[i, j, k] / ρ + + # Vapor is diagnosed: qᵛ = qᵗ - qᶜˡ - qʳ + # where qᵗ = total moisture from Breeze's prognostic ρqᵗ + qᵗ = total_specific_moisture(𝒰) + qᵛ = max(zero(qᵗ), qᵗ - qᶜˡ - qʳ) + + # Update diagnostic fields + μ.qᵛ[i, j, k] = qᵛ + μ.qᶜˡ[i, j, k] = qᶜˡ + μ.qʳ[i, j, k] = qʳ + + # Compute and cache microphysics rates (once per timestep) + T = temperature(𝒰, constants) + + # Convert mass fractions to mixing ratios for Kessler formulas + rᵛ = mass_fraction_to_mixing_ratio(qᵛ, qᵗ) + rᶜˡ = mass_fraction_to_mixing_ratio(qᶜˡ, qᵗ) + rʳ = mass_fraction_to_mixing_ratio(qʳ, qᵗ) + + # Saturation: compute in mixing ratio space + qᵛ⁺ = saturation_specific_humidity(T, ρ, constants, PlanarLiquidSurface()) + rᵛ⁺ = mass_fraction_to_mixing_ratio(qᵛ⁺, qᵗ) + + # Latent heat and heat capacity + L = liquid_latent_heat(T, constants) + q = MoistureMassFractions(qᵛ, qᶜˡ + qʳ) + cₚ = mixture_heat_capacity(q, constants) + + # Compute all rates in mixing ratio space + D = condensation_denominator(T, rᵛ⁺, L, cₚ) + Cₖ_val = condensation_rate(rᵛ, rᵛ⁺, D) + Eₖ_val = cloud_evaporation_rate(rᵛ, rᶜˡ, rᵛ⁺, D) + Aₖ_val = autoconversion_rate(rᶜˡ, km) + Kₖ_val = accretion_rate(rᶜˡ, rʳ, km) + + # Compute rain evaporation rate, coupled to cloud evaporation + # Following Fortran DCMIP2016: rain evaporation is limited by the + # remaining saturation deficit after cloud evaporation + Eʳ_uncoupled = rain_evaporation_rate(ρ, rᵛ, rʳ, rᵛ⁺) + + # Net cloud condensate change: Cₖ - Eₖ (positive = condensation) + # The saturation deficit "used up" by cloud evaporation is Eₖ + # Remaining deficit available for rain evaporation: max(subsaturation - Eₖ/D, 0) + # But simpler: limit rain evaporation to max(-Cₖ + Eₖ - rᶜˡ, 0) based on available deficit + # Following Fortran: ern = min(ern, max(-prod - qc, 0), qr) + # where prod = (qv - qvs)/D is net condensation (negative when subsaturated) + # Here: prod ≈ Cₖ - Eₖ (net condensation rate) + # Remaining deficit after cloud processes: max(-(Cₖ - Eₖ) - rᶜˡ, 0) doesn't quite work + # Simpler interpretation: rain evaporation limited by remaining subsaturation after cloud evaporation + # If cloud fully evaporates, remaining deficit = subsaturation - rᶜˡ + remaining_deficit = max(zero(FT), (rᵛ⁺ - rᵛ) - Eₖ_val) + Eʳ_val = min(Eʳ_uncoupled, remaining_deficit, rʳ) + + # Store rates for use in microphysical_tendency + μ.Cₖ[i, j, k] = Cₖ_val + μ.Eₖ[i, j, k] = Eₖ_val + μ.Aₖ[i, j, k] = Aₖ_val + μ.Kₖ[i, j, k] = Kₖ_val + μ.Eʳ[i, j, k] = Eʳ_val + + # Compute terminal velocity at face k using surface density as reference + # Following Fortran DCMIP2016: ρ₀ = ρ(1) (surface density from column) + # Since we don't have direct access to dynamics here, use a fixed standard value + # TODO: Pass reference_density field to update_microphysical_fields! for proper ρ₀ + ρ₀ = convert(FT, 1.225) # Standard sea-level air density [kg/m³] + wₜ = rain_terminal_velocity(ρ, rʳ, ρ₀) + wʳ = -wₜ # Negative = downward + μ.wʳ[i, j, k] = wʳ + + # Compute surface precipitation rate at k=1 only (2D field) + # precipitation_rate = -wʳ * qʳ [m/s] (positive = precipitation falling out) + if k == 1 + μ.precipitation_rate[i, j, 1] = wₜ * qʳ + end + end + return nothing +end + +@inline function compute_moisture_fractions(i, j, k, grid, ::KM, ρ, qᵗ, μ) + @inbounds begin + qᶜˡ = μ.ρqᶜˡ[i, j, k] / ρ + qʳ = μ.ρqʳ[i, j, k] / ρ + end + # Vapor is diagnosed from total moisture + qᵛ = max(zero(qᵗ), qᵗ - qᶜˡ - qʳ) + + # Rain is counted as liquid in the liquid-ice potential temperature definition + # Total liquid for θˡⁱ = cloud liquid + rain + return MoistureMassFractions(qᵛ, qᶜˡ + qʳ) +end + +# No saturation adjustment for explicit Kessler scheme +@inline maybe_adjust_thermodynamic_state(i, j, k, 𝒰, ::KM, ρ, μ, qᵗ, constants) = 𝒰 + +##### +##### Terminal velocity for rain sedimentation +##### + +""" +$(TYPEDSIGNATURES) + +Compute the terminal fall speed of rain droplets [m s⁻¹]. + +The terminal velocity is given by (following the DCMIP2016 Fortran Kessler reference): + +```math +wₜ = 36.34 (0.001 ρ rʳ)^{0.1364} (ρ₀ / ρ)^{1/2} +``` + +where ρ is air density [kg m⁻³], rʳ is rain mixing ratio [kg kg⁻¹], and ρ₀ is reference +surface density [kg m⁻³]. + +Note: The original formula gives velocity in cm s⁻¹ with coefficient 3634. +Here we use 36.34 m s⁻¹ for SI units. +""" +@inline function rain_terminal_velocity(ρ, rʳ, ρ₀) + FT = typeof(ρ) + # Match Fortran: r(k) = 0.001 * rho(k) is used inside (qr * r)^0.1364. + ρrʳ = convert(FT, 0.001) * ρ * max(zero(FT), rʳ) + + # Avoid issues when there's no rain + ρrʳ <= zero(FT) && return zero(FT) + + # Coefficient 36.34 m/s (converted from 3634 cm/s) + # rhalf = sqrt(ρ₀/ρ) as in Fortran reference + wₜ = convert(FT, 36.34) * ρrʳ^convert(FT, 0.1364) * sqrt(ρ₀ / ρ) + + return wₜ +end + +""" +$(TYPEDSIGNATURES) + +Compute the sedimentation flux for rain at level k. + +Uses upstream differencing following the Fortran Kessler reference: +```math +\\text{sed}_k = \\frac{(ρ r^r w_t)_{k+1} - (ρ r^r w_t)_k}{ρ_k Δz_k} +``` + +At the top boundary (k = Nz), uses: +```math +\\text{sed}_{Nz} = -\\frac{r^r_{Nz} \\cdot w_{t,Nz}}{0.5 \\cdot Δz_{Nz}} +``` + +At the bottom boundary (k = 1), rain falling out is removed (precip). +""" +@inline function sedimentation_tendency(i, j, k, grid, ρᵣ, μ) + FT = eltype(grid) + Nz = size(grid, 3) + + # Get Δz at this level + Δz = Δzᶜᶜᶜ(i, j, k, grid) + + @inbounds begin + # Column densities (use reference-state profile to access k+1 in a local kernel) + ρ_k = ρᵣ[i, j, k] + ρ₀ = ρᵣ[i, j, 1] + + # Current level moisture: convert mass fractions -> mixing ratios (no q≈r shortcut) + qʳ_k = μ.qʳ[i, j, k] + qᵛ_k = μ.qᵛ[i, j, k] + qᶜˡ_k = μ.qᶜˡ[i, j, k] + qᵗ_k = min(qᵛ_k + qᶜˡ_k + qʳ_k, one(FT) - eps(one(FT))) + rʳ_k = mass_fraction_to_mixing_ratio(qʳ_k, qᵗ_k) + + wₜ_k = rain_terminal_velocity(ρ_k, rʳ_k, ρ₀) + + if k == Nz + # Top boundary: no flux from above, only outflow + # sed = -qr * vt / (0.5 * Δz) following Fortran + Δz_half = Δz / 2 + sed = -rʳ_k * wₜ_k / Δz_half + else + # Interior: Fortran-style flux divergence normalized by local density (ρ_k) + ρ_kp1 = ρᵣ[i, j, k+1] + + qʳ_kp1 = μ.qʳ[i, j, k+1] + qᵛ_kp1 = μ.qᵛ[i, j, k+1] + qᶜˡ_kp1 = μ.qᶜˡ[i, j, k+1] + qᵗ_kp1 = min(qᵛ_kp1 + qᶜˡ_kp1 + qʳ_kp1, one(FT) - eps(one(FT))) + rʳ_kp1 = mass_fraction_to_mixing_ratio(qʳ_kp1, qᵗ_kp1) + + wₜ_kp1 = rain_terminal_velocity(ρ_kp1, rʳ_kp1, ρ₀) + + F_kp1 = ρ_kp1 * rʳ_kp1 * wₜ_kp1 + F_k = ρ_k * rʳ_k * wₜ_k + sed = (F_kp1 - F_k) / (ρ_k * Δz) + end + + # At bottom (k=1), rain that would fall below is removed (precipitation) + # This is handled by the flux divergence naturally - flux out at bottom + # is not balanced by flux from below + end + + return sed +end + +""" +$(TYPEDSIGNATURES) + +Return the microphysical velocities for the Kessler scheme. + +For rain (`ρqʳ`), returns the terminal velocity field `wʳ` so that Breeze's +advection machinery handles sedimentation. Cloud liquid has no sedimentation velocity. +""" +@inline function microphysical_velocities(::KM, μ, ::Val{:ρqʳ}) + wʳ = μ.wʳ + return (; u = ZeroField(), v = ZeroField(), w = wʳ) +end +@inline microphysical_velocities(::KM, μ, ::Val{:ρqᶜˡ}) = nothing +@inline microphysical_velocities(::KM, μ, name) = nothing + +##### +##### Source term calculations (in mixing ratio space) +##### + +""" +$(TYPEDSIGNATURES) + +Compute the denominator D for condensation/evaporation rate. + +This follows Klemp & Wilhelmson (1978) eq. 3.10 and the DCMIP Kessler implementation. +The formula derives from the Tetens saturation vapor pressure approximation. + +```math +D = 1 + \\frac{rᵛ⁺ \\cdot 4093 \\cdot L}{cₚ (T - 36)^2} +``` + +where T is temperature in **Kelvin**. The constant 36 K comes from the Tetens formula: +in Celsius, the denominator is (Tc + 237.3), and converting to Kelvin gives +(T - 273.15 + 237.3) = (T - 35.85) ≈ (T - 36). +""" +@inline function condensation_denominator(T, rᵛ⁺, L, cₚ) + FT = typeof(T) + return one(FT) + rᵛ⁺ * convert(FT, 4093) * L / (cₚ * (T - convert(FT, 36))^2) +end + +""" +$(TYPEDSIGNATURES) + +Compute condensation rate Cₖ [kg kg⁻¹ s⁻¹] in mixing ratio space. + +If supersaturated (rᵛ > rᵛ⁺): Cₖ = (rᵛ - rᵛ⁺) / D +Otherwise: Cₖ = 0 +""" +@inline function condensation_rate(rᵛ, rᵛ⁺, D) + FT = typeof(rᵛ) + return rᵛ > rᵛ⁺ ? (rᵛ - rᵛ⁺) / D : zero(FT) +end + +""" +$(TYPEDSIGNATURES) + +Compute cloud evaporation rate Eₖ [kg kg⁻¹ s⁻¹] in mixing ratio space. + +If subsaturated (rᵛ < rᵛ⁺): Eₖ = min(rᶜˡ, (rᵛ⁺ - rᵛ) / D) +Otherwise: Eₖ = 0 + +The evaporation is limited by available cloud water. +""" +@inline function cloud_evaporation_rate(rᵛ, rᶜˡ, rᵛ⁺, D) + FT = typeof(rᵛ) + if rᵛ < rᵛ⁺ + return min(rᶜˡ, (rᵛ⁺ - rᵛ) / D) + else + return zero(FT) + end +end + +""" +$(TYPEDSIGNATURES) + +Compute autoconversion rate Aₖ [kg kg⁻¹ s⁻¹] in mixing ratio space. + +```math +Aₖ = \\max(0, k₁ (rᶜˡ - a)) +``` +""" +@inline function autoconversion_rate(rᶜˡ, km::KM) + FT = typeof(rᶜˡ) + k₁ = km.autoconversion_rate + a = km.autoconversion_threshold + return max(zero(FT), k₁ * (rᶜˡ - a)) +end + +""" +$(TYPEDSIGNATURES) + +Compute accretion rate Kₖ [kg kg⁻¹ s⁻¹] in mixing ratio space. + +```math +Kₖ = k₂ rᶜˡ rʳ^{0.875} +``` +""" +@inline function accretion_rate(rᶜˡ, rʳ, km::KM) + FT = typeof(rᶜˡ) + k₂ = km.accretion_rate + rʳ_safe = max(zero(FT), rʳ) + return k₂ * rᶜˡ * rʳ_safe^convert(FT, 0.875) +end + +""" +$(TYPEDSIGNATURES) + +Compute rain evaporation rate Eʳ [kg kg⁻¹ s⁻¹] in mixing ratio space. + +```math +Eʳ = \\frac{(1 - rᵛ/rᵛ⁺) C (ρ rʳ)^{0.525}}{ρ (5.4 \\times 10^5 + 2.55 \\times 10^6 / (ρ rᵛ⁺))} +``` + +where the ventilation factor is: +```math +C = 1.6 + 124.9 (ρ rʳ)^{0.2046} +``` +""" +@inline function rain_evaporation_rate(ρ, rᵛ, rʳ, rᵛ⁺) + FT = typeof(ρ) + + # No evaporation if saturated or supersaturated + rᵛ >= rᵛ⁺ && return zero(FT) + + # No evaporation if no rain + rʳ <= zero(FT) && return zero(FT) + + ρrʳ = ρ * rʳ + ρrᵛ⁺ = ρ * rᵛ⁺ + + # Ventilation factor + C = convert(FT, 1.6) + convert(FT, 124.9) * ρrʳ^convert(FT, 0.2046) + + # Subsaturation factor + subsaturation = one(FT) - rᵛ / rᵛ⁺ + + # Denominator + denom = convert(FT, 5.4e5) + convert(FT, 2.55e6) / ρrᵛ⁺ + + # Rain evaporation rate + Eʳ = subsaturation * C * ρrʳ^convert(FT, 0.525) / (ρ * denom) + + # Limit by available rain + return min(Eʳ, rʳ) +end + +##### +##### Microphysical tendencies +##### + +""" +$(TYPEDSIGNATURES) + +Compute the tendency for cloud liquid density (ρqᶜˡ). + +The rates Cₖ, Eₖ, Aₖ, Kₖ are computed once per timestep in `update_microphysical_fields!` +and cached in the microphysical fields. + +```math +\\frac{∂(ρqᶜˡ)}{∂t} = ρ \\cdot (1 - qᵗ) \\cdot (Cₖ - Eₖ - Aₖ - Kₖ) +``` + +where the rates Cₖ, Eₖ, Aₖ, Kₖ are in mixing ratio space. +""" +@inline function microphysical_tendency(i, j, k, grid, km::KM, ::Val{:ρqᶜˡ}, ρᵣ, μ, 𝒰, constants) + # Get thermodynamic quantities + ρ = density(𝒰, constants) + qᵗ = total_specific_moisture(𝒰) + + # Get cached rates (computed in update_microphysical_fields!) + @inbounds begin + Cₖ = μ.Cₖ[i, j, k] + Eₖ = μ.Eₖ[i, j, k] + Aₖ = μ.Aₖ[i, j, k] + Kₖ = μ.Kₖ[i, j, k] + end + + # Tendency in mixing ratio space: drᶜˡ/dt = Cₖ - Eₖ - Aₖ - Kₖ + drᶜˡdt = Cₖ - Eₖ - Aₖ - Kₖ + + # Convert to mass fraction tendency + dqᶜˡdt = mixing_ratio_to_mass_fraction(drᶜˡdt, qᵗ) + + return ρ * dqᶜˡdt +end + +""" +$(TYPEDSIGNATURES) + +Compute the tendency for rain density (ρqʳ). + +The rates Aₖ, Kₖ, Eʳ are computed once per timestep in `update_microphysical_fields!` +and cached in the microphysical fields. + +**Sedimentation** is handled by Breeze's advection machinery via `microphysical_velocities`, +which adds the terminal velocity `wʳ` to the rain tracer advection. + +```math +\\frac{∂(ρqʳ)}{∂t} = ρ \\cdot (1 - qᵗ) \\cdot (Aₖ + Kₖ - Eʳ) +``` +""" +@inline function microphysical_tendency(i, j, k, grid, km::KM, ::Val{:ρqʳ}, ρᵣ, μ, 𝒰, constants) + # Get thermodynamic quantities + ρ = density(𝒰, constants) + qᵗ = total_specific_moisture(𝒰) + + # Get cached rates (computed in update_microphysical_fields!) + @inbounds begin + Aₖ = μ.Aₖ[i, j, k] + Kₖ = μ.Kₖ[i, j, k] + Eʳ = μ.Eʳ[i, j, k] + end + + # Tendency in mixing ratio space: drʳ/dt = Aₖ + Kₖ - Eʳ + # Note: sedimentation is handled via microphysical_velocities, not here + drʳdt = Aₖ + Kₖ - Eʳ + + # Convert to mass fraction tendency + dqʳdt = mixing_ratio_to_mass_fraction(drʳdt, qᵗ) + + return ρ * dqʳdt +end + +# Default: no tendency for other variables +# Phase changes (condensation/evaporation of cloud and rain) conserve liquid-ice potential temperature by design. +@inline microphysical_tendency(i, j, k, grid, ::KM, name, ρ, μ, 𝒰, constants) = zero(grid) + +""" +$(TYPEDSIGNATURES) + +Compute the tendency for potential temperature density (ρθˡⁱ) due to rain sedimentation. + +Sedimentation removes rain water (qʳ) from the control volume. Since liquid-ice potential +temperature is defined as θˡⁱ ≈ θ - L/(cₚ Π) qˡ, removing liquid water (qˡ) while +maintaining air temperature (θ) results in an increase in θˡⁱ. + +```math +\\frac{∂(ρθ)}{∂t} = -ρ \\cdot \\frac{L}{cₚ Π} \\cdot \\left(\\frac{∂q_r}{∂t}\\right)_{sed} +``` +""" +@inline function microphysical_tendency(i, j, k, grid, km::KM, ::Val{:ρθ}, ρᵣ, μ, 𝒰, constants) + # Get thermodynamic quantities + ρ = density(𝒰, constants) + T = temperature(𝒰, constants) + qᵗ = total_specific_moisture(𝒰) + + # Get moisture fractions for heat capacity calculation + @inbounds qᵛ = μ.qᵛ[i, j, k] + @inbounds qᶜˡ = μ.qᶜˡ[i, j, k] + @inbounds qʳ = μ.qʳ[i, j, k] + q = MoistureMassFractions(qᵛ, qᶜˡ + qʳ) + + # Calculate sedimentation tendency for rain (in mixing ratio space) + drʳdt_sed = sedimentation_tendency(i, j, k, grid, ρᵣ, μ) + + # Convert to mass fraction tendency + dqʳdt_sed = mixing_ratio_to_mass_fraction(drʳdt_sed, qᵗ) + + # Latent heat and heat capacity + L = liquid_latent_heat(T, constants) + cₚ = mixture_heat_capacity(q, constants) + + # Exner function + Π = exner_function(𝒰, constants) + + # Tendency for θ_li + # Note: dqʳdt_sed is negative for removal. The term -L/(Cp*Π) * dq makes the result positive. + dθdt = -L / (cₚ * Π) * dqʳdt_sed + + return ρ * dθdt +end + +##### +##### Precipitation rate diagnostics +##### + +""" +$(TYPEDSIGNATURES) + +Return the precipitation rate field for the Kessler scheme. + +For `phase = :liquid`, returns the pre-computed `precipitation_rate` 2D field +from `model.microphysical_fields`, which represents the surface precipitation rate [m/s]. + +For `phase = :ice`, returns `nothing` (Kessler is a warm-rain scheme). +""" +precipitation_rate(model, ::KM, ::Val{:liquid}) = model.microphysical_fields.precipitation_rate +precipitation_rate(model, ::KM, ::Val{:ice}) = nothing + +""" +$(TYPEDSIGNATURES) + +Return the surface precipitation flux for the Kessler scheme. + +The surface precipitation flux is `|wʳ| * ρqʳ` at k=1 (bottom face), representing +the rate at which rain mass leaves the domain through the bottom boundary. + +Units: kg/m²/s (positive = downward, out of domain) +""" +function surface_precipitation_flux(model, ::KM) + grid = model.grid + μ = model.microphysical_fields + ρᵣ = model.formulation.reference_state.density + kernel = KesslerSurfacePrecipitationFluxKernel(μ.wʳ, μ.ρqʳ, ρᵣ) + op = KernelFunctionOperation{Center, Center, Nothing}(kernel, grid) + return Field(op) +end + +using Oceananigans.AbstractOperations: KernelFunctionOperation +using Adapt: Adapt, adapt + +struct KesslerSurfacePrecipitationFluxKernel{W, R, D} + terminal_velocity :: W + rain_density :: R + reference_density :: D +end + +Adapt.adapt_structure(to, k::KesslerSurfacePrecipitationFluxKernel) = + KesslerSurfacePrecipitationFluxKernel(adapt(to, k.terminal_velocity), + adapt(to, k.rain_density), + adapt(to, k.reference_density)) + +@inline function (kernel::KesslerSurfacePrecipitationFluxKernel)(i, j, k_idx, grid) + # Flux at bottom face (k=1), ignore k_idx since this is a 2D field + # wʳ < 0 (downward), so -wʳ * ρqʳ > 0 represents flux out of domain + @inbounds wʳ = kernel.terminal_velocity[i, j, 1] + @inbounds ρqʳ = kernel.rain_density[i, j, 1] + + # Return positive flux for rain leaving domain (downward) + return -wʳ * ρqʳ +end diff --git a/test/kessler_microphysics.jl b/test/kessler_microphysics.jl new file mode 100644 index 00000000..5b989db6 --- /dev/null +++ b/test/kessler_microphysics.jl @@ -0,0 +1,194 @@ +using Breeze +using Breeze.Microphysics: KesslerMicrophysics, prognostic_field_names +using Test + +# Import helper functions directly from the module +const BM = Breeze.Microphysics + +@testset "KesslerMicrophysics construction" begin + @testset "Default construction [$(FT)]" for FT in (Float32, Float64) + km = KesslerMicrophysics(FT) + + @test km isa KesslerMicrophysics{FT} + @test km.autoconversion_rate == FT(0.001) + @test km.autoconversion_threshold == FT(0.001) + @test km.accretion_rate == FT(2.2) + end + + @testset "Custom parameters [$(FT)]" for FT in (Float32, Float64) + km = KesslerMicrophysics(FT; + autoconversion_rate = 0.002, + autoconversion_threshold = 0.0005, + accretion_rate = 3.0 + ) + + @test km.autoconversion_rate == FT(0.002) + @test km.autoconversion_threshold == FT(0.0005) + @test km.accretion_rate == FT(3.0) + end +end + +@testset "KesslerMicrophysics interface" begin + @testset "Prognostic field names" begin + km = KesslerMicrophysics() + names = prognostic_field_names(km) + # Only cloud liquid and rain are prognostic; vapor is diagnosed from qᵗ + @test names == (:ρqᶜˡ, :ρqʳ) + end +end + +@testset "Mass fraction ↔ mixing ratio conversion" begin + @testset "mass_fraction_to_mixing_ratio" begin + # Simple case: no moisture → division by 1 + @test BM.mass_fraction_to_mixing_ratio(0.01, 0.0) ≈ 0.01 + + # With 2% total moisture → division by 0.98 + qᵗ = 0.02 + q = 0.01 + r = BM.mass_fraction_to_mixing_ratio(q, qᵗ) + @test r ≈ q / (1 - qᵗ) + end + + @testset "mixing_ratio_to_mass_fraction" begin + qᵗ = 0.02 + r = 0.001 # mixing ratio (or tendency) + q = BM.mixing_ratio_to_mass_fraction(r, qᵗ) + @test q ≈ r * (1 - qᵗ) + end +end + +@testset "Source term calculations" begin + @testset "Autoconversion rate" begin + km = KesslerMicrophysics(Float64) + + # Below threshold: no autoconversion + @test BM.autoconversion_rate(0.0005, km) == 0.0 + + # At threshold: no autoconversion + @test BM.autoconversion_rate(0.001, km) == 0.0 + + # Above threshold: positive autoconversion + Aₖ = BM.autoconversion_rate(0.002, km) + @test Aₖ ≈ 0.001 * (0.002 - 0.001) atol=1e-10 + end + + @testset "Accretion rate" begin + km = KesslerMicrophysics(Float64) + + # No cloud water: no accretion + @test BM.accretion_rate(0.0, 0.001, km) == 0.0 + + # No rain: no accretion + @test BM.accretion_rate(0.001, 0.0, km) == 0.0 + + # Both present: positive accretion + qˡ = 0.001 + qʳ = 0.001 + Kₖ = BM.accretion_rate(qˡ, qʳ, km) + @test Kₖ ≈ 2.2 * qˡ * qʳ^0.875 atol=1e-10 + end + + @testset "Condensation rate" begin + # Supersaturated: condensation occurs + qᵛ = 0.012 + qᵛ⁺ = 0.010 + D = 1.5 + Cₖ = BM.condensation_rate(qᵛ, qᵛ⁺, D) + @test Cₖ ≈ (qᵛ - qᵛ⁺) / D atol=1e-10 + + # Subsaturated: no condensation + Cₖ = BM.condensation_rate(0.008, 0.010, D) + @test Cₖ == 0.0 + end + + @testset "Cloud evaporation rate" begin + D = 1.5 + + # Supersaturated: no evaporation + Eₖ = BM.cloud_evaporation_rate(0.012, 0.001, 0.010, D) + @test Eₖ == 0.0 + + # Subsaturated with cloud water: evaporation occurs + qᵛ = 0.008 + qˡ = 0.002 + qᵛ⁺ = 0.010 + Eₖ = BM.cloud_evaporation_rate(qᵛ, qˡ, qᵛ⁺, D) + expected = min(qˡ, (qᵛ⁺ - qᵛ) / D) + @test Eₖ ≈ expected atol=1e-10 + + # Subsaturated with limited cloud water + qˡ_small = 0.0001 + Eₖ = BM.cloud_evaporation_rate(qᵛ, qˡ_small, qᵛ⁺, D) + @test Eₖ == qˡ_small # Limited by available cloud water + end + + @testset "Rain evaporation rate" begin + # Saturated: no evaporation + Eʳ = BM.rain_evaporation_rate(1.0, 0.010, 0.001, 0.010) + @test Eʳ == 0.0 + + # Supersaturated: no evaporation + Eʳ = BM.rain_evaporation_rate(1.0, 0.012, 0.001, 0.010) + @test Eʳ == 0.0 + + # No rain: no evaporation + Eʳ = BM.rain_evaporation_rate(1.0, 0.008, 0.0, 0.010) + @test Eʳ == 0.0 + + # Subsaturated with rain: evaporation occurs + Eʳ = BM.rain_evaporation_rate(1.0, 0.008, 0.001, 0.010) + @test Eʳ > 0.0 + end + + @testset "Terminal velocity" begin + # Reference density at surface (would come from ρᵣ[1,1,1] in practice) + ρ₀ = 1.0 + + # No rain: zero terminal velocity + wₜ = BM.rain_terminal_velocity(1.0, 0.0, ρ₀) + @test wₜ == 0.0 + + # With rain: positive terminal velocity + wₜ = BM.rain_terminal_velocity(1.0, 0.001, ρ₀) + @test wₜ > 0.0 + + # Higher density decreases terminal velocity + wₜ_low_ρ = BM.rain_terminal_velocity(0.5, 0.001, ρ₀) + wₜ_high_ρ = BM.rain_terminal_velocity(1.5, 0.001, ρ₀) + @test wₜ_low_ρ > wₜ_high_ρ + end +end + +@testset "Mass conservation" begin + # The sum of vapor, cloud, and rain tendencies should be zero + # (neglecting sedimentation which is handled separately) + km = KesslerMicrophysics(Float64) + + # Test parameters + ρ = 1.0 + qᵛ = 0.012 + qˡ = 0.001 + qʳ = 0.0005 + qᵛ⁺ = 0.010 + T = 288.0 + L = 2.5e6 + cₚ = 1005.0 + + D = BM.condensation_denominator(T, qᵛ⁺, L, cₚ) + + # Compute all rates + Cₖ = BM.condensation_rate(qᵛ, qᵛ⁺, D) + Eₖ = BM.cloud_evaporation_rate(qᵛ, qˡ, qᵛ⁺, D) + Aₖ = BM.autoconversion_rate(qˡ, km) + Kₖ = BM.accretion_rate(qˡ, qʳ, km) + Eʳ = BM.rain_evaporation_rate(ρ, qᵛ, qʳ, qᵛ⁺) + + # Tendencies (without density factor) + dqᵛ_dt = -Cₖ + Eₖ + Eʳ + dqˡ_dt = Cₖ - Eₖ - Aₖ - Kₖ + dqʳ_dt = Aₖ + Kₖ - Eʳ + + # Total water tendency should be zero (mass conservation) + total_tendency = dqᵛ_dt + dqˡ_dt + dqʳ_dt + @test abs(total_tendency) < 1e-15 +end