diff --git a/docs/make.jl b/docs/make.jl
index bee842b8..c50a4a0e 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -39,6 +39,7 @@ examples = [
     Example("Stationary parcel model", "stationary_parcel_model", true),
     Example("Acoustic wave in shear layer", "acoustic_wave", true),
     Example("Cloud formation in prescribed updraft", "kinematic_driver", true),
+    Example("Splitting supercell", "splitting_supercell", false),
 ]
 
 # Filter out long-running example if necessary
diff --git a/docs/src/breeze.bib b/docs/src/breeze.bib
index 1208c84a..4171b20a 100644
--- a/docs/src/breeze.bib
+++ b/docs/src/breeze.bib
@@ -322,3 +322,24 @@ @article{Shu1988Efficient
   doi = {10.1016/0021-9991(88)90177-5},
   author = {Chi-Wang Shu and Stanley Osher},
 }
+
+@article{KlempEtAl2015,
+  author = {Klemp, Joseph B. and Skamarock, William C. and Park, Sang-Hun},
+  title = {Idealized global nonhydrostatic atmospheric test cases on a reduced-radius sphere},
+  journal = {Journal of Advances in Modeling Earth Systems},
+  volume = {7},
+  number = {3},
+  pages = {1155--1177},
+  year = {2015},
+  doi = {10.1002/2015MS000435}
+}
+
+@article{Zarzycki2019,
+  author = {Zarzycki, Colin M. and Jablonowski, Christiane and Kent, James and Lauritzen, Peter H. and Nair, Ramachandran and Reed, Kevin A. and Ullrich, Paul A. and Hall, David M. and Dazlich, Don and Heber, Ross and Achatz, Ulrich and Butter, Tobias and Galewsky, Joseph and Goodman, Justin and Klein, Rupert and Lemarié, Florian and Malardel, Sylvie and Rauscher, Sara A. and Schar, Christoph and Sprenger, Michael and Taylor, Mark A. and Vogl, Christian and Wan, Hui and Williamson, David L.},
+  title = {{DCMIP2016}: the splitting supercell test case},
+  journal = {Geoscientific Model Development},
+  volume = {12},
+  pages = {879--892},
+  year = {2019},
+  doi = {10.5194/gmd-12-879-2019}
+}
diff --git a/examples/splitting_supercell.jl b/examples/splitting_supercell.jl
new file mode 100644
index 00000000..06842a40
--- /dev/null
+++ b/examples/splitting_supercell.jl
@@ -0,0 +1,371 @@
+# # Splitting supercell
+#
+# This example simulates the development of a splitting supercell thunderstorm, following the
+# idealized test case described by [KlempEtAl2015](@citet) and the DCMIP2016
+# supercell intercomparison by [Zarzycki2019](@citet). This benchmark evaluates the model's
+# ability to capture deep moist convection with warm-rain microphysics and strong updrafts.
+#
+# For microphysics we use the Kessler scheme, which includes prognostic cloud water
+# and rain water with autoconversion, accretion, rain evaporation, and sedimentation processes.
+# This is the same scheme used in the DCMIP2016 supercell intercomparison [Zarzycki2019](@cite).
+#
+# ## Physical setup
+#
+# The simulation initializes a conditionally unstable atmosphere with a warm bubble perturbation
+# that triggers deep convection. The environment includes:
+# - A realistic tropospheric potential temperature profile with a tropopause at 12 km
+# - Moisture that decreases with height, with relative humidity dropping above the tropopause
+# - Wind shear in the lower 5 km to promote storm rotation and supercell development
+#
+# ### Potential temperature profile
+#
+# The background potential temperature follows a piecewise profile
+# (Equation 14 in [KlempEtAl2015](@citet)):
+#
+# ```math
+# θ(z) = \begin{cases}
+#     θ_0 + (θ_{\rm tr} - θ_0) \left(\frac{z}{z_{\rm tr}}\right)^{5/4} & z \leq z_{\rm tr} \\
+#     θ_{\rm tr} \exp\left(\frac{g}{c_p^d T_{\rm tr}} (z - z_{\rm tr})\right) & z > z_{\rm tr}
+# \end{cases}
+# ```
+#
+# where ``θ_0 = 300 \, {\rm K}`` is the surface potential temperature,
+# ``θ_{\rm tr} = 343 \, {\rm K}`` is the tropopause potential temperature,
+# ``z_{\rm tr} = 12 \, {\rm km}`` is the tropopause height, and
+# ``T_{\rm tr} = 213 \, {\rm K}`` is the tropopause temperature.
+#
+# ### Warm bubble perturbation
+#
+# A localized warm bubble triggers convection (Equations 17–18 in [KlempEtAl2015](@citet)):
+#
+# ```math
+# θ'(x, y, z) = \begin{cases}
+#     Δθ \cos^2\left(\frac{\pi}{2} R\right) & R < 1 \\
+#     0 & R \geq 1
+# \end{cases}
+# ```
+#
+# where ``R = \sqrt{(r/r_h)^2 + ((z-z_c)/r_z)^2}`` is the normalized radius,
+# ``r = \sqrt{(x-x_c)^2 + (y-y_c)^2}`` is the horizontal distance from the bubble center,
+# ``Δθ = 3 \, {\rm K}`` is the perturbation amplitude, ``r_h = 10 \, {\rm km}`` is the
+# horizontal radius, and ``r_z = 1.5 \, {\rm km}`` is the vertical radius.
+#
+# ### Wind shear profile
+#
+# The zonal wind increases linearly with height up to the shear layer ``z_s = 5 \, {\rm km}``,
+# with a smooth transition zone, providing the environmental shear necessary for supercell
+# development and mesocyclone formation (Equations 15-16 in [KlempEtAl2015](@citet)).
+
+using Breeze
+using Breeze: DCMIP2016KesslerMicrophysics, TetensFormula
+using Oceananigans: Oceananigans
+using Oceananigans.Units
+using Oceananigans.Grids: znodes
+
+using CairoMakie
+using CUDA
+using Printf
+
+# ## Domain and grid
+#
+# The domain is 168 km × 168 km × 20 km with 168 × 168 × 40 grid points, giving
+# 1 km horizontal resolution and 500 m vertical resolution. The grid uses periodic
+# lateral boundary conditions and bounded top/bottom boundaries.
+
+Oceananigans.defaults.FloatType = Float32
+
+Nx, Ny, Nz = 168, 168, 40
+Lx, Ly, Lz = 168kilometers, 168kilometers, 20kilometers
+
+grid = RectilinearGrid(GPU(),
+                       size = (Nx, Ny, Nz),
+                       x = (0, Lx),
+                       y = (0, Ly),
+                       z = (0, Lz),
+                       halo = (5, 5, 5),
+                       topology = (Periodic, Periodic, Bounded))
+
+# ## Reference state and dynamics
+#
+# We define the anelastic reference state with surface pressure ``p_0 = 1000 \, {\rm hPa}``
+# and reference potential temperature ``θ_0 = 300 \, {\rm K}``.
+
+constants = ThermodynamicConstants(saturation_vapor_pressure = TetensFormula())
+
+reference_state = ReferenceState(grid, constants,
+                                 surface_pressure = 100000,
+                                 potential_temperature = 300)
+
+dynamics = AnelasticDynamics(reference_state)
+
+# ## Background atmosphere profiles
+#
+# The atmospheric stratification parameters define the troposphere-stratosphere transition.
+
+θ₀ = 300       # K - surface potential temperature
+θᵖ = 343       # K - tropopause potential temperature
+zᵖ = 12000     # m - tropopause height
+Tᵖ = 213       # K - tropopause temperature
+nothing #hide
+
+# Wind shear parameters control the low-level environmental wind profile:
+
+zˢ = 5kilometers  # m - shear layer height
+uˢ = 30           # m/s - maximum shear wind speed
+uᶜ = 15           # m/s - storm motion (Galilean translation speed)
+nothing #hide
+
+# Extract thermodynamic constants for profile calculations:
+
+g = constants.gravitational_acceleration
+cᵖᵈ = constants.dry_air.heat_capacity
+nothing #hide
+
+# Background potential temperature profile (Equation 14 in [KlempEtAl2015](@citet)):
+
+function θ_background(z)
+    θᵗ = θ₀ + (θᵖ - θ₀) * (z / zᵖ)^(5/4)
+    θˢ = θᵖ * exp(g / (cᵖᵈ * Tᵖ) * (z - zᵖ))
+    return (z <= zᵖ) * θᵗ + (z > zᵖ) * θˢ
+end
+
+# Relative humidity profile (decreases with height, 25% above tropopause):
+
+ℋ_background(z) = (1 - 3/4 * (z / zᵖ)^(5/4)) * (z <= zᵖ) + 1/4 * (z > zᵖ)
+
+# Zonal wind profile with linear shear below ``zˢ`` and smooth transition (Equations 15-16):
+
+function u_background(z)
+    uˡ = uˢ * (z / zˢ) - uᶜ
+    uᵗ = (-4/5 + 3 * (z / zˢ) - 5/4 * (z / zˢ)^2) * uˢ - uᶜ
+    uᵘ = uˢ - uᶜ
+    return (z < (zˢ - 1000)) * uˡ +
+           (abs(z - zˢ) <= 1000) * uᵗ +
+           (z > (zˢ + 1000)) * uᵘ
+end
+
+# ## Warm bubble perturbation
+#
+# The warm bubble parameters following Equations 17–18 in [KlempEtAl2015](@citet):
+
+Δθ = 3              # K - perturbation amplitude
+rᵇʰ = 10kilometers  # m - bubble horizontal radius
+rᵇᵛ = 1500          # m - bubble vertical radius
+zᵇ = 1500           # m - bubble center height
+xᵇ = Lx / 2         # m - bubble center x-coordinate
+yᵇ = Ly / 2         # m - bubble center y-coordinate
+nothing #hide
+
+# The total initial potential temperature combines the background profile with the
+# cosine-squared warm bubble perturbation:
+
+function θᵢ(x, y, z)
+    θ̄ = θ_background(z)
+    r = sqrt((x - xᵇ)^2 + (y - yᵇ)^2)
+    R = sqrt((r / rᵇʰ)^2 + ((z - zᵇ) / rᵇᵛ)^2)
+    θ′ = ifelse(R < 1, Δθ * cos((π / 2) * R)^2, 0.0)
+    return θ̄ + θ′
+end
+
+uᵢ(x, y, z) = u_background(z)
+
+# ## Visualization of initial conditions and warm bubble perturbation
+#
+# We visualize the background potential temperature, relative humidity, and wind shear profiles
+# that define the environmental stratification:
+
+θ_profile = set!(Field{Nothing, Nothing, Center}(grid), z -> θ_background(z))
+ℋ_profile = set!(Field{Nothing, Nothing, Center}(grid), z -> ℋ_background(z) * 100)
+u_profile = set!(Field{Nothing, Nothing, Center}(grid), z -> u_background(z))
+
+fig = Figure(size=(1000, 400), fontsize=14)
+
+axθ = Axis(fig[1, 1], xlabel="θ (K)", ylabel="z (km)", title="Potential temperature")
+lines!(axθ, θ_profile, linewidth=2, color=:magenta)
+hlines!(axθ, [zᵖ / 1000], color=:gray, linestyle=:dash)
+
+axℋ = Axis(fig[1, 2], xlabel="ℋ (%)", ylabel="z (km)", title="Relative humidity")
+lines!(axℋ, ℋ_profile, linewidth=2, color=:dodgerblue)
+hlines!(axℋ, [zᵖ / 1000], color=:gray, linestyle=:dash)
+
+axu = Axis(fig[1, 3], xlabel="u (m/s)", ylabel="z (km)", title="Wind profile")
+lines!(axu, u_profile, linewidth=2, color=:orangered)
+hlines!(axu, [zˢ / 1000], color=:gray, linestyle=:dash)
+vlines!(axu, [0], color=:black, linestyle=:dot)
+
+save("supercell_initial_conditions.png", fig) #src
+fig
+
+# Visualize the warm bubble perturbation on a vertical slice through the domain center:
+
+θ′_slice = set!(Field{Center, Nothing, Center}(grid), (x, z) -> θᵢ(x, yᵇ, z) - θ_background(z))
+
+fig = Figure(size=(700, 400), fontsize=14)
+ax = Axis(fig[1, 1], xlabel="x (km)", ylabel="z (km)",
+          title="Warm bubble perturbation θ′")
+
+hm = heatmap!(ax, θ′_slice, colormap=:thermal, colorrange=(0, Δθ))
+Colorbar(fig[1, 2], hm, label="θ′ (K)")
+
+save("supercell_warm_bubble.png", fig) #src
+fig
+
+# ## Model setup
+#
+# We use the DCMIP2016 Kessler microphysics scheme with high-order WENO advection.
+# The Kessler scheme includes prognostic cloud water and rain water with autoconversion,
+# accretion, rain evaporation, and sedimentation processes.
+
+microphysics = DCMIP2016KesslerMicrophysics()
+advection = WENO(order=9, minimum_buffer_upwind_order=3)
+
+model = AtmosphereModel(grid; dynamics, microphysics, advection, thermodynamic_constants=constants)
+
+# ## Model initialization
+#
+# We initialize the model with the previously described initial conditions, including a warm-bubble perturbation.
+# We precompute the RH field to ensure GPU compatibility.
+
+ℋᵢ = set!(CenterField(grid), (x, y, z) -> ℋ_background(z))
+
+set!(model, θ=θᵢ, ℋ=ℋᵢ, u=uᵢ)
+
+# ## Simulation
+#
+# Run for 2 hours with adaptive time stepping (CFL = 0.7):
+
+simulation = Simulation(model; Δt=2, stop_time=2hours)
+conjure_time_step_wizard!(simulation, cfl=0.7)
+
+# ## Output and progress
+#
+# We set up callbacks to monitor simulation health and collect diagnostics.
+# The maximum vertical velocity is tracked during the simulation to avoid
+# saving large 3D datasets.
+
+θˡⁱ = liquid_ice_potential_temperature(model)
+qᶜˡ = model.microphysical_fields.qᶜˡ
+qʳ = model.microphysical_fields.qʳ
+qᵛ = model.microphysical_fields.qᵛ
+u, v, w = model.velocities
+
+wall_clock = Ref(time_ns())
+
+function progress(sim)
+    elapsed = 1e-9 * (time_ns() - wall_clock[])
+
+    msg = @sprintf("Iter: %d, t: %s, Δt: %s, wall time: %s, max|u|: %.2f m/s, max w: %.2f m/s, min w: %.2f m/s",
+                   iteration(sim), prettytime(sim), prettytime(sim.Δt), prettytime(elapsed),
+                   maximum(abs, u), maximum(w), minimum(w))
+
+    msg *= @sprintf(", max(qᵛ): %.2e, max(qᶜˡ): %.2e, max(qʳ): %.2e",
+                    maximum(qᵛ), maximum(qᶜˡ), maximum(qʳ))
+    @info msg
+
+    return nothing
+end
+
+add_callback!(simulation, progress, IterationInterval(100))
+
+# Collect maximum vertical velocity time series during simulation:
+
+max_w_ts = []
+max_w_times = []
+
+function collect_max_w(sim)
+    push!(max_w_times, time(sim))
+    push!(max_w_ts, maximum(w))
+    return nothing
+end
+
+add_callback!(simulation, collect_max_w, TimeInterval(1minutes))
+
+# Save horizontal slices at z ≈ 5 km for animation:
+
+z = znodes(grid, Center())
+k_5km = searchsortedfirst(z, 5000)
+@info "Saving xy slices at z = $(z[k_5km]) m (k = $k_5km)"
+
+slice_outputs = (
+    wxy = view(w, :, :, k_5km),
+    qʳxy = view(qʳ, :, :, k_5km),
+    qᶜˡxy = view(qᶜˡ, :, :, k_5km),
+)
+
+slices_filename = "splitting_supercell_slices.jld2"
+simulation.output_writers[:slices] = JLD2Writer(model, slice_outputs; filename=slices_filename,
+                                                including = [:grid],
+                                                schedule = TimeInterval(2minutes),
+                                                overwrite_existing = true)
+
+run!(simulation)
+
+# ## Animation: horizontal slices at z ≈ 5 km
+#
+# We create a 3-panel animation showing the storm structure at mid-levels:
+# - Vertical velocity ``w``: reveals the updraft/downdraft structure
+# - Cloud liquid ``qᶜˡ``: shows the cloud boundaries
+# - Rain ``qʳ``: indicates precipitation regions
+#
+# The simulated supercell exhibits splitting behavior, with the initial storm
+# dividing into right-moving and left-moving cells, consistent with the
+# DCMIP2016 intercomparison results [Zarzycki2019](@cite).
+
+wxy_ts = FieldTimeSeries(slices_filename, "wxy")
+qʳxy_ts = FieldTimeSeries(slices_filename, "qʳxy")
+qᶜˡxy_ts = FieldTimeSeries(slices_filename, "qᶜˡxy")
+
+times = wxy_ts.times
+Nt = length(times)
+
+wlim = maximum(abs, wxy_ts) / 2
+qʳlim = maximum(qʳxy_ts) / 4
+qᶜˡlim = maximum(qᶜˡxy_ts) / 4
+
+fig = Figure(size=(900, 400), fontsize=12)
+
+axw = Axis(fig[1, 1], aspect=1, xlabel="x (m)", ylabel="y (m)", title="w (m/s)")
+axqᶜˡ = Axis(fig[1, 2], aspect=1, xlabel="x (m)", ylabel="y (m)", title="qᶜˡ (kg/kg)")
+axqʳ = Axis(fig[1, 3], aspect=1, xlabel="x (m)", ylabel="y (m)", title="qʳ (kg/kg)")
+
+n = Observable(1)
+wxy_n = @lift wxy_ts[$n]
+qᶜˡxy_n = @lift qᶜˡxy_ts[$n]
+qʳxy_n = @lift qʳxy_ts[$n]
+title = @lift "Splitting supercell at z ≈ 5 km, t = " * prettytime(times[$n])
+
+hmw = heatmap!(axw, wxy_n, colormap=:balance, colorrange=(-wlim, wlim))
+hmqᶜˡ = heatmap!(axqᶜˡ, qᶜˡxy_n, colormap=:dense, colorrange=(0, qᶜˡlim))
+hmqʳ = heatmap!(axqʳ, qʳxy_n, colormap=:amp, colorrange=(0, qʳlim))
+
+Colorbar(fig[2, 1], hmw, vertical=false)
+Colorbar(fig[2, 2], hmqᶜˡ, vertical=false)
+Colorbar(fig[2, 3], hmqʳ, vertical=false)
+
+fig[0, :] = Label(fig, title, fontsize=14, tellwidth=false)
+
+CairoMakie.record(fig, "splitting_supercell_slices.mp4", 1:Nt, framerate=10) do nn
+    n[] = nn
+end
+nothing #hide
+
+# ![](splitting_supercell_slices.mp4)
+
+# ## Results: maximum vertical velocity time series
+#
+# The maximum updraft velocity is a key diagnostic for supercell intensity.
+# Strong supercells typically develop updrafts exceeding 30–50 m/s.
+#
+# Our simulated storm intensity is notably stronger than the DCMIP2016 intercomparison
+# results reported by [Zarzycki2019](@citet). One explanation is that
+# no explicit numerical diffusion is applied in this simulation. As noted by
+# [KlempEtAl2015](@citet), the simulated storm intensity and structure
+# are highly sensitive to numerical diffusion.
+
+fig = Figure(size=(700, 400), fontsize=14)
+ax = Axis(fig[1, 1], xlabel="Time (s)", ylabel="Maximum w (m/s)", title="Maximum Vertical Velocity",
+          xticks=0:1800:7200)
+lines!(ax, max_w_times, max_w_ts, linewidth=2)
+
+save("supercell_max_w.png", fig) #src
+fig
diff --git a/examples/supercell.jl b/examples/supercell.jl
deleted file mode 100644
index 5a77bd6e..00000000
--- a/examples/supercell.jl
+++ /dev/null
@@ -1,423 +0,0 @@
-# # Supercell thunderstorm
-#
-# This example simulates the development of a splitting supercell thunderstorm, following the idealized
-# test case described by [KlempEtAl2015](@citet) ("Idealized global nonhydrostatic atmospheric
-# test cases on a reduced-radius sphere"). This benchmark evaluates the model's ability to
-# capture deep moist convection with cloud microphysics, and strong updrafts.
-#
-# ## Physical setup
-#
-# The simulation initializes a conditionally unstable atmosphere with a warm bubble perturbation
-# that triggers deep convection. The environment includes:
-# - A realistic tropospheric potential temperature profile with a tropopause at 12 km
-# - Moisture that decreases with height, with relative humidity dropping above the tropopause
-# - Wind shear in the lower 5 km to promote storm rotation and supercell development
-#
-# ### Potential temperature profile
-#
-# The background potential temperature follows a piecewise profile:
-#
-# ```math
-# θ^{\rm bg}(z) = \begin{cases}
-#     θ_0 + (θ_{\rm tr} - θ_0) \left(\frac{z}{z_{\rm tr}}\right)^{5/4} & z \leq z_{\rm tr} \\
-#     θ_{\rm tr} \exp\left(\frac{g}{c_p^d T_{\rm tr}} (z - z_{\rm tr})\right) & z > z_{\rm tr}
-# \end{cases}
-# ```
-#
-# where ``θ_0 = 300 \, {\rm K}`` is the surface potential temperature,
-# ``θ_{\rm tr} = 343 \, {\rm K}`` is the tropopause potential temperature,
-# ``z_{\rm tr} = 12 \, {\rm km}`` is the tropopause height, and
-# ``T_{\rm tr} = 213 \, {\rm K}`` is the tropopause temperature.
-#
-# ### Warm bubble perturbation
-#
-# A localized warm bubble triggers convection (Equations 17–18 in [KlempEtAl2015](@cite)):
-#
-# ```math
-# θ'(x, y, z) = \begin{cases}
-#     Δθ \cos^2\left(\frac{\pi}{2} R\right) & R < 1 \\
-#     0 & R \geq 1
-# \end{cases}
-# ```
-#
-# where ``R = \sqrt{(r/r_h)^2 + ((z-z_c)/r_z)^2}`` is the normalized radius,
-# ``r = \sqrt{(x-x_c)^2 + (y-y_c)^2}`` is the horizontal distance from the bubble center,
-# ``Δθ = 3 \, {\rm K}`` is the perturbation amplitude, ``r_h = 10 \, {\rm km}`` is the
-# horizontal radius, and ``r_z = 1.5 \, {\rm km}`` is the vertical radius.
-#
-# ### Wind shear profile
-#
-# The zonal wind increases linearly with height up to the shear layer ``z_s = 5 \, {\rm km}``,
-# with a smooth transition zone, providing the environmental shear necessary for supercell
-# development and mesocyclone formation.
-
-using Breeze
-using Oceananigans.Units
-using Statistics
-using Printf
-using CairoMakie
-
-using Oceananigans.Grids: znode
-using Oceananigans: Center, Face
-using Oceananigans.Operators: Δzᶜᶜᶜ, Δzᶜᶜᶠ, ℑzᵃᵃᶠ
-using Breeze.Thermodynamics: dry_air_gas_constant
-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
-# 1 km horizontal resolution and 500 m vertical resolution. The grid uses periodic
-# lateral boundary conditions and bounded top/bottom boundaries.
-
-Nx, Ny, Nz = 168, 168, 40
-Lx, Ly, Lz = 168kilometers, 168kilometers, 20kilometers
-
-grid = RectilinearGrid(GPU(),
-                       size = (Nx, Ny, Nz),
-                       x = (0, Lx),
-                       y = (0, Ly),
-                       z = (0, Lz),
-                       halo = (5, 5, 5),
-                       topology = (Periodic, Periodic, Bounded))
-
-# ## Thermodynamic parameters
-#
-# We define the reference state with surface pressure ``p_0 = 1000 \, {\rm hPa}`` and
-# reference potential temperature ``θ_0 = 300 \, {\rm K}``:
-
-p₀ = 100000  # Pa - surface pressure
-θ₀ = 300     # K - reference potential temperature
-constants = ThermodynamicConstants()
-reference_state = ReferenceState(grid, constants, surface_pressure=p₀, potential_temperature=θ₀)
-dynamics = AnelasticDynamics(reference_state)
-
-# ## Background atmosphere profiles
-#
-# The atmospheric stratification parameters define the troposphere-stratosphere transition:
-
-θₜᵣ = 343           # K - tropopause potential temperature
-zₜᵣ = 12000         # m - tropopause height
-Tₜᵣ = 213           # K - tropopause temperature
-
-# Wind shear parameters control the low-level environmental wind profile:
-
-zₛ = 5kilometers    # m - shear layer height
-uₛ = 30             # m/s - maximum shear wind speed
-u_c = 15            # m/s - storm motion (Galilean translation speed)
-
-# Extract thermodynamic constants for profile calculations:
-
-g = constants.gravitational_acceleration
-cᵖᵈ = constants.dry_air.heat_capacity
-Rᵈ = dry_air_gas_constant(constants)
-
-# Background potential temperature profile (Equation 14 in [KlempEtAl2015](@cite)):
-
-function θᵇᵍ(x, y, z)
-    θ_troposphere = θ₀ + (θₜᵣ - θ₀) * (z / zₜᵣ)^(5/4)
-    θ_stratosphere = θₜᵣ * exp(g / (cᵖᵈ * Tₜᵣ) * (z - zₜᵣ))
-    return (z <= zₜᵣ) * θ_troposphere + (z > zₜᵣ) * θ_stratosphere
-end
-
-# Relative humidity profile (decreases with height, 25% above tropopause):
-
-RHᵇᵍ(z) = (1 - 3/4 * (z / zₜᵣ)^(5/4)) * (z <= zₜᵣ) + 1/4 * (z > zₜᵣ)
-
-# ## Plot: Initial thermodynamic profiles
-#
-# Visualize the background potential temperature and relative humidity profiles:
-
-z_plot = range(0, Lz, length=200)
-θ_profile = [θᵇᵍ(0, 0, z) for z in z_plot]
-RH_profile = [RHᵇᵍ(z) * 100 for z in z_plot]  # Convert to percentage
-
-fig_thermo = Figure(size=(900, 600))
-
-ax_theta = Axis(fig_thermo[1, 1],
-                xlabel = "Potential temperature θ (K)",
-                ylabel = "Height (m)",
-                title = "Background Potential Temperature")
-lines!(ax_theta, θ_profile, collect(z_plot), linewidth = 2, color = :red)
-
-ax_rh = Axis(fig_thermo[1, 2],
-             xlabel = "Relative humidity (%)",
-             ylabel = "Height (m)",
-             title = "Relative Humidity Profile")
-lines!(ax_rh, RH_profile, collect(z_plot), linewidth = 2, color = :blue)
-
-save("supercell_thermo_profiles.png", fig_thermo) #src
-fig_thermo
-
-# Zonal wind profile with linear shear below ``z_s`` and smooth transition (Equation 15-16):
-
-function uᵢ(x, y, z)
-    u_lower = uₛ * (z / zₛ) - u_c
-    u_transition = (-4/5 + 3 * (z / zₛ) - 5/4 * (z / zₛ)^2) * uₛ - u_c
-    u_upper = uₛ - u_c
-    return (z < (zₛ - 1000)) * u_lower +
-           (abs(z - zₛ) <= 1000) * u_transition +
-           (z > (zₛ + 1000)) * u_upper
-end
-
-# Meridional wind is zero (no cross-flow):
-
-vᵢ(x, y, z) = 0.0
-
-# ## Plot: Background wind profile
-#
-# Visualize the environmental wind shear that promotes supercell rotation:
-
-u_profile = [uᵢ(0, 0, z) for z in z_plot]
-v_profile = [vᵢ(0, 0, z) for z in z_plot]
-
-fig_wind = Figure(size=(500, 600))
-ax_wind = Axis(fig_wind[1, 1],
-               xlabel = "Wind speed (m/s)",
-               ylabel = "Height (m)",
-               title = "Background Wind Profile")
-
-lines!(ax_wind, u_profile, collect(z_plot), label = "u (zonal)", linewidth = 2)
-lines!(ax_wind, v_profile, collect(z_plot), label = "v (meridional)", linewidth = 2, linestyle = :dash)
-axislegend(ax_wind, position = :rb)
-
-save("supercell_wind_profile.png", fig_wind) #src
-fig_wind
-
-# ## Warm bubble initial perturbation
-#
-# The warm bubble parameters following Equations 17–18 in [KlempEtAl2015](@cite):
-
-Δθ = 3              # K - perturbation amplitude
-rₕ = 10kilometers   # m - bubble horizontal radius
-rᵥ = 1500           # m - bubble vertical radius
-zᵦ = 1500           # m - bubble center height
-xᵦ = Lx / 2         # m - bubble center x-coordinate
-yᵦ = Ly / 2         # m - bubble center y-coordinate
-
-# The total initial potential temperature combines the background profile with the
-# cosine-squared warm bubble perturbation:
-
-function θᵢ(x, y, z)
-    θ_base = θᵇᵍ(x, y, z)
-    r = sqrt((x - xᵦ)^2 + (y - yᵦ)^2)
-    R = sqrt((r / rₕ)^2 + ((z - zᵦ) / rᵥ)^2)
-    θ_pert = ifelse(R < 1, Δθ * cos((π / 2) * R)^2, 0.0)
-    return θ_base + θ_pert
-end
-
-# ## Plot: Warm bubble perturbation
-#
-# Visualize the warm bubble perturbation on a vertical slice through the domain center:
-
-x_slice = range(0, Lx, length=200)
-z_slice = range(0, Lz, length=200)  # Focus on lower atmosphere where bubble is located
-
-θ_pert_slice = [θᵢ(x, yᵦ, z) - θᵇᵍ(x, yᵦ, z) for x in x_slice, z in z_slice]
-
-fig_bubble = Figure(size=(800, 400))
-ax_bubble = Axis(fig_bubble[1, 1],
-                 xlabel = "x (km)",
-                 ylabel = "Height (m)",
-                 title = "Warm Bubble Perturbation θ' (K) at y = Ly/2")
-
-hm = heatmap!(ax_bubble, collect(x_slice) ./ 1000, collect(z_slice), θ_pert_slice,
-              colormap = :thermal, colorrange = (0, Δθ))
-Colorbar(fig_bubble[1, 2], hm, label = "θ' (K)")
-
-save("supercell_warm_bubble.png", fig_bubble) #src
-fig_bubble
-
-# ## Model initialization
-#
-# Create the atmosphere model with one-moment cloud microphysics from CloudMicrophysics.jl,
-# high-order WENO advection, and anisotropic minimum dissipation turbulence closure:
-
-microphysics = BreezeOneMomentCloudMicrophysics()
-advection = WENO(order=9, minimum_buffer_upwind_order=3)
-closure = AnisotropicMinimumDissipation()
-model = AtmosphereModel(grid; dynamics, closure, microphysics, advection)
-
-# Initialize with background potential temperature to compute hydrostatic pressure:
-
-set!(model, θ = θᵇᵍ)
-
-# ## Water vapor initialization
-#
-# Compute the initial water vapor mixing ratio from the saturation mixing ratio
-# and relative humidity profile. The saturation mixing ratio uses the Tetens formula:
-#
-# ```math
-# q_v^* = \frac{380}{p} \exp\left(17.27 \frac{T - 273}{T - 36}\right)
-# ```
-
-ph = Breeze.AtmosphereModels.compute_hydrostatic_pressure!(CenterField(grid), model)
-T = model.temperature
-
-qᵛᵢ = Field{Center, Center, Center}(grid)
-
-# Transfer to CPU for scalar indexing (required when using GPU arrays):
-
-ph_host = Array(parent(ph))
-T_host = Array(parent(T))
-qᵛᵢ_host = similar(ph_host)
-
-# Compute initial water vapor from relative humidity and saturation mixing ratio:
-
-for k in axes(qᵛᵢ_host, 3), j in axes(qᵛᵢ_host, 2), i in axes(qᵛᵢ_host, 1)
-    z = znode(i, j, k, grid, Center(), Center(), Center())
-    T_eq = @inbounds T_host[i, j, k]
-    p_eq = @inbounds ph_host[i, j, k]
-    qᵛ_sat = 380 / p_eq * exp(17.27 * ((T_eq - 273) / (T_eq - 36)))
-    @inbounds qᵛᵢ_host[i, j, k] = RHᵇᵍ(z) * qᵛ_sat
-end
-
-copyto!(parent(qᵛᵢ), qᵛᵢ_host)
-
-# Set the full initial conditions (water vapor, potential temperature with bubble, and wind):
-
-set!(model, qᵗ = qᵛᵢ, θ = θᵢ, u = uᵢ)
-
-# Compute potential temperature perturbation for diagnostics:
-
-θ = Breeze.AtmosphereModels.liquid_ice_potential_temperature(model)
-θᵇᵍf = CenterField(grid)
-set!(θᵇᵍf, (x, y, z) -> θᵇᵍ(x, y, z))
-θ′ = θ - θᵇᵍf
-
-# Extract microphysical fields for output:
-
-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)
-
-# 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ᶜⁱ
-
-    ρe = Breeze.AtmosphereModels.static_energy_density(sim.model)
-    ρemean = mean(ρe)
-
-    msg = @sprintf("Iter: %d, t: %s, Δt: %s, mean(ρe): %.6e J/kg, max|u|: %.5f m/s, max w: %.5f m/s, min w: %.5f m/s",
-                   iteration(sim), prettytime(sim), prettytime(sim.Δt), ρemean,
-                   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ᶜⁱ))
-    @info msg
-    return nothing
-end
-
-add_callback!(simulation, progress, IterationInterval(100))
-
-# ## Output
-#
-# Save full 3D fields for post-processing analysis:
-
-outputs = merge(model.velocities, model.tracers, (; θ, θ′, qᶜˡ, qᶜⁱ, qᵛ))
-
-filename = "supercell.jld2"
-ow = JLD2Writer(model, outputs; filename,
-                including = [:grid],
-                schedule = TimeInterval(1minutes),
-                overwrite_existing = true)
-simulation.output_writers[:jld2] = ow
-
-run!(simulation)
-
-# ## Results: maximum vertical velocity time series
-#
-# The maximum updraft velocity is a key diagnostic for supercell intensity.
-# Strong supercells typically develop updrafts exceeding 30–50 m/s:
-
-wt = FieldTimeSeries(filename, "w")
-times = wt.times
-max_w = [maximum(wt[n]) for n in 1:length(times)]
-
-fig = Figure()
-ax = Axis(fig[1, 1],
-          xlabel = "Time [s]",
-          ylabel = "Max w [m/s]",
-          title = "Maximum Vertical Velocity",
-          xticks = 0:900:maximum(times))
-lines!(ax, times, max_w)
-
-save("max_w_timeseries.png", fig) #src
-fig
-
-# ## Animation: horizontal slices at 5 km
-#
-# 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
-
-wxy_ts = FieldTimeSeries("supercell_slices.jld2", "wxy")
-qʳxy_ts = FieldTimeSeries("supercell_slices.jld2", "qʳxy")
-qᶜˡxy_ts = FieldTimeSeries("supercell_slices.jld2", "qᶜˡxy")
-
-times = wxy_ts.times
-Nt = length(times)
-
-# Set color limits for visualization:
-
-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
-
-# 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[3, 1], xlabel="x (m)", ylabel="y (m)", title="Rain water qʳ")
-
-# Set up observables for animation:
-
-n = Observable(1)
-wxy_n = @lift wxy_ts[$n]
-qᶜˡxy_n = @lift qᶜˡxy_ts[$n]
-qʳxy_n = @lift qʳxy_ts[$n]
-title_text = @lift "Supercell: Horizontal slices at z ≈ 5 km, t = " * prettytime(times[$n])
-
-# Create heatmaps and colorbars:
-
-hmw = heatmap!(axw, wxy_n, colormap=:balance, colorrange=(-wlim, wlim))
-hmqᶜˡ = heatmap!(axqᶜˡ, qᶜˡxy_n, colormap=:viridis, colorrange=(0, qᶜˡlim))
-hmqʳ = heatmap!(axqʳ, qʳxy_n, colormap=:viridis, colorrange=(0, qʳlim))
-
-Colorbar(slices_fig[2, 1], hmw, label="w (m/s)", vertical=false, flipaxis=false)
-Colorbar(slices_fig[2, 2], hmqᶜˡ, label="qᶜˡ (kg/kg)", vertical=false, flipaxis=false)
-Colorbar(slices_fig[4, 1], hmqʳ, label="qʳ (kg/kg)", vertical=false, flipaxis=false)
-
-slices_fig[0, :] = Label(slices_fig, title_text, fontsize=18, tellwidth=false)
-
-# Record the animation:
-
-CairoMakie.record(slices_fig, "supercell_horizontal_5km.mp4", 1:Nt, framerate=10) do nn
-    n[] = nn
-end
-
-@info "Animation saved to supercell_horizontal_5km.mp4"
-
-# ![](supercell_horizontal_5km.mp4)