This commit adds a new diagnostic which computes the convective

available potential energy based on the buoyancy difference between
a parcel lifted adiabatically from the surface and the envrionment.
CAPE is useful for diagnosing the potential for extreme weather
like thunderstorms, tornados, etc. It is accessed by including `cape`
in the diagnostics of any config file. Typically cape is computed
by integrating between the LFC and the LNB, however here we simply
integrate the positive buoyancy between the surface and 20km as this
was significantly simpler to implement in ClimaCore. Future work
should revisit this.

Author: Julians42

NEWS.md

@@ -4,6 +4,9 @@ ClimaAtmos.jl Release Notes
 main
 -------
 
+### Add diagnostic for CAPE
+PR [#3820](https://github.com/CliMA/ClimaAtmos.jl/pull/3820) adds support for computing convective available potential energy (CAPE), or the vertical integral of the buoyancy differential between a parcel lifted from the surface and the environment. Exemplified in the TRMM deep convection case.
+
 v0.30.2
 -------
 

config/model_configs/prognostic_edmfx_trmm_column_0M.yml

@@ -30,11 +30,10 @@ z_stretch: false
 dz_bottom: 30
 dt: 150secs
 t_end: 6hours
-dt_save_state_to_disk: 10mins
 toml: [toml/prognostic_edmfx_implicit_scm_calibrated_5_cases_shallow_deep_v1.toml]
 netcdf_interpolation_num_points: [8, 8, 82]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr, cape]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins

src/diagnostics/Diagnostics.jl

@@ -11,6 +11,9 @@ import ClimaCore:
 import ClimaCore.Utilities: half
 import Thermodynamics as TD
 
+# compute lazily to reduce allocations
+import ..lazy
+
 import ..AtmosModel
 import ..AtmosCallback
 import ..EveryNSteps

src/diagnostics/core_diagnostics.jl

@@ -1349,3 +1349,62 @@ add_diagnostic_variable!(
         end
     end,
 )
+
+###
+# Convective Available Potential Energy (2d)
+###
+function compute_cape!(out, state, cache, time)
+    thermo_params = lazy.(CAP.thermodynamics_params(cache.params))
+    g = lazy.(TD.Parameters.grav(thermo_params))
+
+    # Get surface parcel properties
+    surface_thermal_state = lazy.(cache.precomputed.sfc_conditions.ts)
+    surface_q =
+        lazy.(TD.total_specific_humidity.(thermo_params, surface_thermal_state))
+    surface_θ_liq_ice =
+        lazy.(TD.liquid_ice_pottemp.(thermo_params, surface_thermal_state))
+
+    # Create parcel thermodynamic states at each level based on energy & moisture at surface
+    parcel_ts_moist =
+        lazy.(
+            TD.PhaseEquil_pθq.(
+                thermo_params,
+                cache.precomputed.ᶜp,
+                surface_θ_liq_ice,
+                surface_q,
+            )
+        )
+
+    # Calculate virtual temperatures for parcel & environment
+    parcel_Tv = lazy.(TD.virtual_temperature.(thermo_params, parcel_ts_moist))
+    env_Tv =
+        lazy.(TD.virtual_temperature.(thermo_params, cache.precomputed.ᶜts))
+
+    # Calculate buoyancy from the difference in virtual temperatures
+    ᶜbuoyancy = cache.scratch.ᶜtemp_scalar
+    ᶜbuoyancy .= g .* (parcel_Tv .- env_Tv) ./ env_Tv
+
+    # restrict to tropospheric buoyancy (generously below 20km) TODO: integrate from LFC to LNB 
+    FT = Spaces.undertype(axes(ᶜbuoyancy))
+    ᶜbuoyancy .=
+        ᶜbuoyancy .*
+        ifelse.(
+            Fields.coordinate_field(state.c.ρ).z .< FT(20000.0),
+            FT(1.0),
+            FT(0.0),
+        )
+
+    # Integrate positive buoyancy to get CAPE
+    out′ = isnothing(out) ? zeros(axes(Fields.level(state.f, half))) : out
+    Operators.column_integral_definite!(out′, lazy.(max.(ᶜbuoyancy, 0)))
+    return out′
+end
+
+add_diagnostic_variable!(
+    short_name = "cape",
+    long_name = "Convective Available Potential Energy",
+    standard_name = "convective_available_potential_energy",
+    units = "J kg^-1",
+    comments = "Energy available to a parcel lifted moist adiabatically from the surface. We assume fully reversible phase changes and no precipitation.",
+    compute! = compute_cape!,
+)

src/diagnostics/default_diagnostics.jl

@@ -209,6 +209,7 @@ function default_diagnostics(
         "lwp",
         "clwvi",
         "clivi",
+        "cape",
     ]
     average_func = frequency_averages(duration)
     return [average_func(moist_diagnostics...; output_writer, start_date)...]