diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 4087c552338..d1223e50184 100755
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -311,13 +311,13 @@ steps:
     steps:
 
       - label: ":balloon: Single column EDMF Bomex"
-        command: "julia --color=yes --project=examples examples/hybrid/driver.jl --config column --FLOAT_TYPE Float64 --hyperdiff false --moist equil --turbconv edmf --turbconv_case Bomex --split_ode false --ode_algo ODE.Euler --dt_save_to_sol 5mins --z_elem 60 --z_stretch false --z_max 3e3 --job_id edmf_bomex --dt 6secs --t_end 6hours --regression_test true --anelastic_dycore true --apply_limiter false"
+        command: "julia --color=yes --project=examples examples/hybrid/driver.jl --config column --FLOAT_TYPE Float64 --hyperdiff false --moist equil --turbconv edmf --turbconv_case Bomex --split_ode false --ode_algo ODE.Euler --dt_save_to_sol 5mins --z_elem 60 --z_stretch false --z_max 3e3 --job_id edmf_bomex --dt 6secs --t_end 6hours --regression_test true --anelastic_dycore true --apply_limiter false --debugging_tc true --dt_save_to_disk 6hours"
         artifact_paths: "edmf_bomex/*"
         agents:
           slurm_mem: 20GB
 
       - label: ":balloon: Compressible single column EDMF Bomex"
-        command: "julia --color=yes --project=examples examples/hybrid/driver.jl --config column --FLOAT_TYPE Float64 --hyperdiff false --moist equil --turbconv edmf --turbconv_case Bomex --dt_save_to_sol 5mins --z_elem 60 --z_stretch false --z_max 3e3 --job_id compressible_edmf_bomex --dt 1secs --t_end 6hours --regression_test true --apply_limiter false"
+        command: "julia --color=yes --project=examples examples/hybrid/driver.jl --config column --FLOAT_TYPE Float64 --hyperdiff false --moist equil --turbconv edmf --turbconv_case Bomex --dt_save_to_sol 5mins --z_elem 60 --z_stretch false --z_max 3e3 --job_id compressible_edmf_bomex --dt 1secs --t_end 6hours --regression_test true --apply_limiter false --debugging_tc true --dt_save_to_disk 6hours"
         artifact_paths: "compressible_edmf_bomex/*"
         agents:
           slurm_mem: 20GB
diff --git a/examples/hybrid/callbacks.jl b/examples/hybrid/callbacks.jl
index 291940c87f5..31154a6b23c 100644
--- a/examples/hybrid/callbacks.jl
+++ b/examples/hybrid/callbacks.jl
@@ -151,6 +151,7 @@ function turb_conv_affect_filter!(integrator)
 end
 
 function save_to_disk_func(integrator)
+
     (; t, u, p) = integrator
     (; output_dir) = p.simulation
     Y = u
@@ -160,6 +161,7 @@ function save_to_disk_func(integrator)
     else
         (; ᶜts, ᶜp, params, ᶜK) = p
     end
+
     thermo_params = CAP.thermodynamics_params(params)
     cm_params = CAP.microphysics_params(params)
 
@@ -189,6 +191,7 @@ function save_to_disk_func(integrator)
 
     # cloudwater (liquid and ice), watervapor and RH for moist simulation
     if :ρq_tot in propertynames(Y.c)
+
         ᶜq = @. TD.PhasePartition(thermo_params, ᶜts)
         ᶜcloud_liquid = @. ᶜq.liq
         ᶜcloud_ice = @. ᶜq.ice
@@ -230,6 +233,53 @@ function save_to_disk_func(integrator)
         moist_diagnostic = NamedTuple()
     end
 
+    if :edmf_cache in propertynames(p) && p.simulation.is_debugging_tc
+
+        tc_cent(p) = p.edmf_cache.aux.cent.turbconv
+        tc_face(p) = p.edmf_cache.aux.face.turbconv
+        turbulence_convection_diagnostic = (;
+            bulk_up_area = tc_cent(p).bulk.area,
+            bulk_up_h_tot = tc_cent(p).bulk.h_tot,
+            bulk_up_buoyancy = tc_cent(p).bulk.buoy,
+            bulk_up_q_tot = tc_cent(p).bulk.q_tot,
+            bulk_up_q_liq = tc_cent(p).bulk.q_liq,
+            bulk_up_q_ice = tc_cent(p).bulk.q_ice,
+            bulk_up_temperature = tc_cent(p).bulk.T,
+            bulk_up_cloud_fraction = tc_cent(p).bulk.cloud_fraction,
+            bulk_up_e_tot_tendency_precip_formation = tc_cent(p).bulk.e_tot_tendency_precip_formation,
+            bulk_up_qt_tendency_precip_formation = tc_cent(p).bulk.qt_tendency_precip_formation,
+            env_w = tc_cent(p).en.w,
+            env_area = tc_cent(p).en.area,
+            env_q_tot = tc_cent(p).en.q_tot,
+            env_q_liq = tc_cent(p).en.q_liq,
+            env_q_ice = tc_cent(p).en.q_ice,
+            env_theta_liq_ice = tc_cent(p).en.θ_liq_ice,
+            env_theta_virt = tc_cent(p).en.θ_virt,
+            env_theta_dry = tc_cent(p).en.θ_dry,
+            env_e_tot = tc_cent(p).en.e_tot,
+            env_e_kin = tc_cent(p).en.e_kin,
+            env_h_tot = tc_cent(p).en.h_tot,
+            env_RH = tc_cent(p).en.RH,
+            env_s = tc_cent(p).en.s,
+            env_temperature = tc_cent(p).en.T,
+            env_buoyancy = tc_cent(p).en.buoy,
+            env_cloud_fraction = tc_cent(p).en.cloud_fraction,
+            env_TKE = tc_cent(p).en.tke,
+            env_Hvar = tc_cent(p).en.Hvar,
+            env_QTvar = tc_cent(p).en.QTvar,
+            env_HQTcov = tc_cent(p).en.HQTcov,
+            env_e_tot_tendency_precip_formation = tc_cent(p).en.e_tot_tendency_precip_formation,
+            env_qt_tendency_precip_formation = tc_cent(p).en.qt_tendency_precip_formation,
+            env_Hvar_rain_dt = tc_cent(p).en.Hvar_rain_dt,
+            env_QTvar_rain_dt = tc_cent(p).en.QTvar_rain_dt,
+            env_HQTcov_rain_dt = tc_cent(p).en.HQTcov_rain_dt,
+            face_bulk_w = tc_face(p).bulk.w,
+            face_env_w = tc_face(p).en.w,
+        )
+    else
+        turbulence_convection_diagnostic = NamedTuple()
+    end
+
     if vert_diff
         (; dif_flux_uₕ, dif_flux_energy, dif_flux_ρq_tot) = p
         vert_diff_diagnostic = (;
@@ -290,6 +340,7 @@ function save_to_disk_func(integrator)
         vert_diff_diagnostic,
         rad_diagnostic,
         rad_clear_diagnostic,
+        turbulence_convection_diagnostic,
     )
 
     day = floor(Int, t / (60 * 60 * 24))
diff --git a/examples/hybrid/cli_options.jl b/examples/hybrid/cli_options.jl
index 46e4943080c..1d9d8e7e8fc 100644
--- a/examples/hybrid/cli_options.jl
+++ b/examples/hybrid/cli_options.jl
@@ -218,6 +218,10 @@ function parse_commandline()
         help = "Apply a horizontal limiter to every tracer [`true` (default), `false`]"
         arg_type = Bool
         default = true
+        "--debugging_tc"
+        help = "Save most of the tc aux state to HDF5 file [`false` (default), `true`]"
+        arg_type = Bool
+        default = false
     end
     parsed_args = ArgParse.parse_args(ARGS, s)
     return (s, parsed_args)
diff --git a/examples/hybrid/define_tc_quicklook_profiles.jl b/examples/hybrid/define_tc_quicklook_profiles.jl
new file mode 100644
index 00000000000..df6b9e934eb
--- /dev/null
+++ b/examples/hybrid/define_tc_quicklook_profiles.jl
@@ -0,0 +1,119 @@
+import ClimaCore: Fields, InputOutput
+using Plots
+
+function plot_tc_profiles(folder, hdf5_filename)
+
+    input_filename = joinpath(folder, hdf5_filename)
+    output_filename = joinpath(folder, "final_profiles.png")
+
+    reader = InputOutput.HDF5Reader(input_filename)
+    Y = InputOutput.read_field(reader, "Y")
+    D = InputOutput.read_field(reader, "diagnostics")
+
+    zc = Fields.coordinate_field(Y.c).z
+    zf = Fields.coordinate_field(Y.f).z
+
+    args =
+        (; tickfontsize = 13, guidefontsize = 16, legendfontsize = 10, lw = 3)
+
+    p1 = plot(
+        parent(D.bulk_up_area)[:],
+        parent(zc)[:],
+        label = "up area";
+        args...,
+    )
+    p1 = plot!(parent(D.env_area)[:], parent(zc)[:], label = "en area"; args...)
+    p2 = plot(
+        parent(D.bulk_up_q_tot)[:],
+        parent(zc)[:],
+        label = "up qt";
+        args...,
+    )
+    p3 = plot(
+        parent(D.bulk_up_q_liq)[:],
+        parent(zc)[:],
+        label = "up ql";
+        args...,
+    )
+    p4 = plot(
+        parent(D.bulk_up_q_ice)[:],
+        parent(zc)[:],
+        label = "up qi";
+        args...,
+    )
+    p5 = plot(parent(D.face_bulk_w)[:], parent(zf)[:], label = "up w"; args...)
+    p5 = plot!(parent(D.face_env_w)[:], parent(zf)[:], label = "en w"; args...)
+    p6 = plot(parent(D.env_q_tot)[:], parent(zc)[:], label = "en qt"; args...)
+    p7 = plot(parent(D.env_q_liq)[:], parent(zc)[:], label = "en ql"; args...)
+    p8 = plot(parent(D.env_q_ice)[:], parent(zc)[:], label = "en qi"; args...)
+    p9 = plot(
+        parent(D.bulk_up_buoyancy)[:],
+        parent(zc)[:],
+        label = "up buoy";
+        args...,
+    )
+    p9 = plot!(
+        parent(D.env_buoyancy)[:],
+        parent(zc)[:],
+        label = "en buoy";
+        args...,
+    )
+    p10 = plot(
+        parent(D.bulk_up_temperature)[:],
+        parent(zc)[:],
+        label = "up T";
+        args...,
+    )
+    p10 = plot!(
+        parent(D.env_temperature)[:],
+        parent(zc)[:],
+        label = "en T";
+        args...,
+    )
+    p11 = plot(
+        parent(D.bulk_up_cloud_fraction)[:],
+        parent(zc)[:],
+        label = "up CF";
+        args...,
+    )
+    p11 = plot!(
+        parent(D.env_cloud_fraction)[:],
+        parent(zc)[:],
+        label = "env CF";
+        args...,
+    )
+    p12 = plot(parent(D.env_RH)[:], parent(zc)[:], label = "en RH"; args...)
+    p13 = plot(parent(D.env_TKE)[:], parent(zc)[:], label = "en TKE"; args...)
+    p14 = plot(parent(D.env_Hvar)[:], parent(zc)[:], label = "en Hvar"; args...)
+    p15 =
+        plot(parent(D.env_QTvar)[:], parent(zc)[:], label = "en QTvar"; args...)
+    p16 = plot(
+        parent(D.env_HQTcov)[:],
+        parent(zc)[:],
+        label = "en HQTcov";
+        args...,
+    )
+
+    p = plot(
+        p1,
+        p2,
+        p3,
+        p4,
+        p5,
+        p6,
+        p7,
+        p8,
+        p9,
+        p10,
+        p11,
+        p12,
+        p13,
+        p14,
+        p15,
+        p16,
+        size = (2400.0, 1500.0),
+        bottom_margin = 20.0 * Plots.PlotMeasures.px,
+        left_margin = 20.0 * Plots.PlotMeasures.px,
+    )
+    png(p, output_filename)
+end
diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl
index c15e45d4e64..03dc53c5383 100644
--- a/examples/hybrid/driver.jl
+++ b/examples/hybrid/driver.jl
@@ -319,6 +319,14 @@ if !simulation.is_distributed && parsed_args["post_process"]
     end
 end
 
+if parsed_args["debugging_tc"]
+    include(joinpath(@__DIR__, "define_tc_quicklook_profiles.jl"))
+    plot_tc_profiles(
+        simulation.output_dir,
+        "day0."*string(Int(simulation.t_end))*".hdf5"
+    )
+end
+
 if parsed_args["regression_test"]
     # Test results against main branch
     include(
diff --git a/examples/hybrid/types.jl b/examples/hybrid/types.jl
index 48af42eb587..586412f88ec 100644
--- a/examples/hybrid/types.jl
+++ b/examples/hybrid/types.jl
@@ -151,6 +151,7 @@ function get_simulation(::Type{FT}, parsed_args) where {FT}
 
     sim = (;
         is_distributed = haskey(ENV, "CLIMACORE_DISTRIBUTED"),
+        is_debugging_tc = parsed_args["debugging_tc"],
         output_dir,
         restart = haskey(ENV, "RESTART_FILE"),
         job_id,