Merge pull request #3820 from CliMA/j/cape

Add CAPE diagnostic

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)...]