diff --git a/post_processing/post_processing_funcs.jl b/post_processing/post_processing_funcs.jl
index 28d1f5f410e..1b1eae488b5 100644
--- a/post_processing/post_processing_funcs.jl
+++ b/post_processing/post_processing_funcs.jl
@@ -228,7 +228,7 @@ function custom_postprocessing(sol, output_dir, p)
     anim = @animate for (Y, t) in zip(sol.u, sol.t)
         CA.set_precomputed_quantities!(Y, p, t) # sets ᶜts
         Plots.plot(
-            vec(TD.air_temperature.(thermo_params, p.ᶜts)),
+            vec(TD.air_temperature.(thermo_params, p.precomputed.ᶜts)),
             vec(Fields.coordinate_field(Y.c).z ./ 1000);
             xlabel = "T [K]",
             ylabel = "z [km]",
@@ -294,7 +294,7 @@ function postprocessing_plane(sol, output_dir, p)
     end
 
     gen_plot_plane(
-        Geometry.UVector.(p.ᶜu),
+        Geometry.UVector.(p.precomputed.ᶜu),
         "horz_velocity.png",
         "Horizontal Velocity",
         "u[m/s]",
@@ -302,7 +302,7 @@ function postprocessing_plane(sol, output_dir, p)
     )
 
     gen_plot_plane(
-        Geometry.WVector.(p.ᶜu),
+        Geometry.WVector.(p.precomputed.ᶜu),
         "vert_velocity.png",
         "Vertical Velocity",
         "w[m/s]",
@@ -310,7 +310,7 @@ function postprocessing_plane(sol, output_dir, p)
     )
 
     gen_plot_plane(
-        TD.virtual_pottemp.(thermo_params, p.ᶜts),
+        TD.virtual_pottemp.(thermo_params, p.precomputed.ᶜts),
         "virtual_pottemp.png",
         "Virtual Pottemp",
         "Theta[K]",
diff --git a/src/cache/cache.jl b/src/cache/cache.jl
index 4152af389df..e22dc629ae2 100644
--- a/src/cache/cache.jl
+++ b/src/cache/cache.jl
@@ -52,12 +52,7 @@ function build_cache(Y, atmos, params, surface_setup, simulation)
     net_energy_flux_toa = [Geometry.WVector(FT(0))]
     net_energy_flux_sfc = [Geometry.WVector(FT(0))]
 
-    default_cache = (;
-        is_init = Ref(true),
-        simulation,
-        atmos,
-        sfc_setup = surface_setup(params),
-        limiter,
+    core = (
         ᶜΦ,
         ᶠgradᵥ_ᶜΦ = ᶠgradᵥ.(ᶜΦ),
         ᶜρ_ref,
@@ -65,42 +60,64 @@ function build_cache(Y, atmos, params, surface_setup, simulation)
         ᶜT = similar(Y.c, FT),
         ᶜf,
         ∂ᶜK_∂ᶠu₃ = similar(Y.c, BidiagonalMatrixRow{Adjoint{FT, CT3{FT}}}),
+    )
+
+    sfc_setup = surface_setup(params)
+    scratch = temporary_quantities(Y, atmos)
+    env_thermo_quad = SGSQuadrature(FT)
+
+    precomputed = precomputed_quantities(Y, atmos)
+    precomputing_arguments =
+        (; atmos, core, params, sfc_setup, precomputed, scratch, simulation)
+
+    # Coupler compatibility
+    isnothing(precomputing_arguments.sfc_setup) &&
+        SurfaceConditions.set_dummy_surface_conditions!(precomputing_arguments)
+
+    set_precomputed_quantities!(Y, precomputing_arguments, FT(0))
+
+    radiation_args =
+        atmos.radiation_mode isa RRTMGPI.AbstractRRTMGPMode ?
+        (params, precomputed.ᶜp) : ()
+
+    hyperdiff = hyperdiffusion_cache(Y, atmos)
+    rayleigh_sponge = rayleigh_sponge_cache(Y, atmos)
+    viscous_sponge = viscous_sponge_cache(Y, atmos)
+    precipitation = precipitation_cache(Y, atmos)
+    subsidence = subsidence_cache(Y, atmos)
+    large_scale_advection = large_scale_advection_cache(Y, atmos)
+    edmf_coriolis = edmf_coriolis_cache(Y, atmos)
+    forcing = forcing_cache(Y, atmos)
+    non_orographic_gravity_wave = non_orographic_gravity_wave_cache(Y, atmos)
+    orographic_gravity_wave = orographic_gravity_wave_cache(Y, atmos)
+    radiation = radiation_model_cache(Y, atmos, radiation_args...)
+
+    p = (;
+        simulation,
+        atmos,
+        core,
+        sfc_setup,
         params,
         do_dss,
+        precomputed,
+        scratch,
+        hyperdiff,
+        rayleigh_sponge,
+        viscous_sponge,
+        precipitation,
+        subsidence,
+        large_scale_advection,
+        edmf_coriolis,
+        forcing,
+        non_orographic_gravity_wave,
+        orographic_gravity_wave,
+        radiation,
+        env_thermo_quad,
         ghost_buffer,
+        dt = simulation.dt,
         net_energy_flux_toa,
         net_energy_flux_sfc,
-        env_thermo_quad = SGSQuadrature(FT),
-        precomputed_quantities(Y, atmos)...,
-        scratch = temporary_quantities(Y, atmos),
         numerics,
     )
-    set_precomputed_quantities!(Y, default_cache, FT(0))
-    default_cache.is_init[] = false
-
-    radiation_args =
-        atmos.radiation_mode isa RRTMGPI.AbstractRRTMGPMode ?
-        (params, default_cache.ᶜp) : ()
-
-    return merge(
-        (;
-            hyperdiff = hyperdiffusion_cache(Y, atmos),
-            rayleigh_sponge = rayleigh_sponge_cache(Y, atmos),
-            viscous_sponge = viscous_sponge_cache(Y, atmos),
-            precipitation = precipitation_cache(Y, atmos),
-            large_scale_advection = large_scale_advection_cache(Y, atmos),
-            subsidence = subsidence_cache(Y, atmos),
-            edmf_coriolis = edmf_coriolis_cache(Y, atmos),
-            forcing = forcing_cache(Y, atmos),
-            non_orographic_gravity_wave = non_orographic_gravity_wave_cache(
-                Y,
-                atmos,
-            ),
-            orographic_gravity_wave = orographic_gravity_wave_cache(Y, atmos),
-            radiation = radiation_model_cache(Y, atmos, radiation_args...),
-        ),
-        (; Δt = simulation.dt),
-        (; atmos.turbconv_model),
-        default_cache,
-    )
+    return p
 end
diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl
index a78389c3801..828deb9a577 100644
--- a/src/cache/diagnostic_edmf_precomputed_quantities.jl
+++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl
@@ -105,13 +105,13 @@ function set_diagnostic_edmf_precomputed_quantities_bottom_bc!(Y, p, t)
     (; turbconv_model) = p.atmos
     FT = eltype(Y)
     n = n_mass_flux_subdomains(turbconv_model)
-
-    (; ᶜp, ᶜΦ, ᶠu³, ᶜh_tot, ᶜK) = p
-    (; q_tot) = p.ᶜspecific
+    (; ᶜΦ) = p.core
+    (; ᶜp, ᶠu³, ᶜh_tot, ᶜK) = p.precomputed
+    (; q_tot) = p.precomputed.ᶜspecific
     (; ustar, obukhov_length, buoyancy_flux, ρ_flux_h_tot, ρ_flux_q_tot) =
-        p.sfc_conditions
-    (; ᶜρaʲs, ᶠu³ʲs, ᶜuʲs, ᶜKʲs, ᶜh_totʲs, ᶜq_totʲs, ᶜtsʲs, ᶜρʲs) = p
-    (; ᶠu³⁰, ᶜK⁰, ᶜh_tot⁰) = p
+        p.precomputed.sfc_conditions
+    (; ᶜρaʲs, ᶠu³ʲs, ᶜKʲs, ᶜh_totʲs, ᶜq_totʲs, ᶜtsʲs, ᶜρʲs) = p.precomputed
+    (; ᶠu³⁰, ᶜK⁰, ᶜh_tot⁰) = p.precomputed
 
     thermo_params = CAP.thermodynamics_params(p.params)
 
@@ -231,9 +231,10 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
     ᶜz = Fields.coordinate_field(Y.c).z
     (; params) = p
     (; dt) = p.simulation
-    (; ᶜp, ᶜΦ, ᶜρ_ref, ᶠu³, ᶜts, ᶜh_tot, ᶜK) = p
-    (; q_tot) = p.ᶜspecific
-    (; buoyancy_flux) = p.sfc_conditions
+    (; ᶜΦ, ᶜρ_ref) = p.core
+    (; ᶜp, ᶠu³, ᶜts, ᶜh_tot, ᶜK) = p.precomputed
+    (; q_tot) = p.precomputed.ᶜspecific
+    (; buoyancy_flux) = p.precomputed.sfc_conditions
     (;
         ᶜρaʲs,
         ᶠu³ʲs,
@@ -247,8 +248,8 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
         ᶜnh_pressureʲs,
         ᶜS_q_totʲs,
         ᶜS_e_totʲs_helper,
-    ) = p
-    (; ᶠu³⁰, ᶜK⁰, ᶜh_tot⁰) = p
+    ) = p.precomputed
+    (; ᶠu³⁰, ᶜK⁰, ᶜh_tot⁰) = p.precomputed
     thermo_params = CAP.thermodynamics_params(params)
     microphys_params = CAP.microphysics_params(params)
 
@@ -617,8 +618,8 @@ Updates the top boundary condition of precomputed quantities stored in `p` for d
 """
 function set_diagnostic_edmf_precomputed_quantities_top_bc!(Y, p, t)
     n = n_mass_flux_subdomains(p.atmos.turbconv_model)
-    (; ᶜentrʲs, ᶜdetrʲs, ᶜS_q_totʲs, ᶜS_e_totʲs_helper) = p
-    (; ᶠu³⁰, ᶠu³ʲs, ᶜuʲs) = p
+    (; ᶜentrʲs, ᶜdetrʲs, ᶜS_q_totʲs, ᶜS_e_totʲs_helper) = p.precomputed
+    (; ᶠu³⁰, ᶠu³ʲs, ᶜuʲs) = p.precomputed
 
     # set values for the top level
     i_top = Spaces.nlevels(axes(Y.c))
@@ -664,13 +665,13 @@ function set_diagnostic_edmf_precomputed_quantities_env_closures!(Y, p, t)
     ᶜdz = Fields.Δz_field(axes(Y.c))
     (; params) = p
     (; dt) = p.simulation
-    (; ᶜp, ᶜu, ᶜts) = p
-    (; q_tot) = p.ᶜspecific
-    (; ustar, obukhov_length) = p.sfc_conditions
-    (; ᶜρaʲs, ᶠu³ʲs, ᶜdetrʲs) = p
-    (; ᶜtke⁰, ᶠu³⁰, ᶜS_q_tot⁰, ᶜu⁰) = p
-    (; ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜmixing_length) = p
-    (; ᶜK_h, ᶜK_u, ρatke_flux) = p
+    (; ᶜp, ᶜu, ᶜts) = p.precomputed
+    (; q_tot) = p.precomputed.ᶜspecific
+    (; ustar, obukhov_length) = p.precomputed.sfc_conditions
+    (; ᶜρaʲs, ᶠu³ʲs, ᶜdetrʲs) = p.precomputed
+    (; ᶜtke⁰, ᶠu³⁰, ᶜS_q_tot⁰, ᶜu⁰) = p.precomputed
+    (; ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜmixing_length) = p.precomputed
+    (; ᶜK_h, ᶜK_u, ρatke_flux) = p.precomputed
     thermo_params = CAP.thermodynamics_params(params)
     microphys_params = CAP.microphysics_params(params)
     ᶜlg = Fields.local_geometry_field(Y.c)
diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
index 0edba97307f..3624298c3b4 100644
--- a/src/cache/precomputed_quantities.jl
+++ b/src/cache/precomputed_quantities.jl
@@ -302,7 +302,8 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
     thermo_params = CAP.thermodynamics_params(p.params)
     n = n_mass_flux_subdomains(turbconv_model)
     thermo_args = (thermo_params, energy_form, moisture_model)
-    (; ᶜspecific, ᶜu, ᶠu³, ᶜK, ᶜts, ᶜp, ᶜΦ) = p
+    (; ᶜΦ) = p.core
+    (; ᶜspecific, ᶜu, ᶠu³, ᶜK, ᶜts, ᶜp) = p.precomputed
     ᶠuₕ³ = p.scratch.ᶠtemp_CT3
 
     @. ᶜspecific = specific_gs(Y.c)
@@ -333,12 +334,14 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
     @. ᶜp = TD.air_pressure(thermo_params, ᶜts)
 
     if energy_form isa TotalEnergy
-        (; ᶜh_tot) = p
+        (; ᶜh_tot) = p.precomputed
         @. ᶜh_tot =
             TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜspecific.e_tot)
     end
 
-    SurfaceConditions.update_surface_conditions!(Y, p, t)
+    if !isnothing(p.sfc_setup)
+        SurfaceConditions.update_surface_conditions!(Y, p, t)
+    end
 
     if turbconv_model isa PrognosticEDMFX
         set_prognostic_edmf_precomputed_quantities_environment!(Y, p, ᶠuₕ³, t)
@@ -365,10 +368,16 @@ values of the first updraft.
 function output_prognostic_sgs_quantities(Y, p, t)
     (; turbconv_model) = p.atmos
     thermo_params = CAP.thermodynamics_params(p.params)
-    (; ᶜρa⁰, ᶜρ⁰, ᶜtsʲs) = p
+    (; ᶜρa⁰, ᶜρ⁰, ᶜtsʲs) = p.precomputed
     ᶠuₕ³ = p.scratch.ᶠtemp_CT3
     set_ᶠuₕ³!(ᶠuₕ³, Y)
-    (ᶠu₃⁺, ᶜu⁺, ᶠu³⁺, ᶜK⁺) = similar.((p.ᶠu₃⁰, p.ᶜu⁰, p.ᶠu³⁰, p.ᶜK⁰))
+    (ᶠu₃⁺, ᶜu⁺, ᶠu³⁺, ᶜK⁺) =
+        similar.((
+            p.precomputed.ᶠu₃⁰,
+            p.precomputed.ᶜu⁰,
+            p.precomputed.ᶠu³⁰,
+            p.precomputed.ᶜK⁰,
+        ))
     set_sgs_ᶠu₃!(u₃⁺, ᶠu₃⁺, Y, turbconv_model)
     set_velocity_quantities!(ᶜu⁺, ᶠu³⁺, ᶜK⁺, ᶠu₃⁺, Y.c.uₕ, ᶠuₕ³)
     ᶜts⁺ = ᶜtsʲs.:1
@@ -385,8 +394,8 @@ values of the first updraft.
 """
 function output_diagnostic_sgs_quantities(Y, p, t)
     thermo_params = CAP.thermodynamics_params(p.params)
-    (; ᶜρaʲs, ᶜtsʲs) = p
-    ᶠu³⁺ = p.ᶠu³ʲs.:1
+    (; ᶜρaʲs, ᶜtsʲs) = p.precomputed
+    ᶠu³⁺ = p.precomputed.ᶠu³ʲs.:1
     ᶜu⁺ = @. (C123(Y.c.uₕ) + C123(ᶜinterp(ᶠu³⁺)))
     ᶜts⁺ = @. ᶜtsʲs.:1
     ᶜa⁺ = @. draft_area(ᶜρaʲs.:1, TD.air_density(thermo_params, ᶜts⁺))
diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl
index b9e336e3572..3ffa1b9f02a 100644
--- a/src/cache/prognostic_edmf_precomputed_quantities.jl
+++ b/src/cache/prognostic_edmf_precomputed_quantities.jl
@@ -14,8 +14,10 @@ function set_prognostic_edmf_precomputed_quantities_environment!(Y, p, ᶠuₕ³
 
     thermo_params = CAP.thermodynamics_params(p.params)
     (; turbconv_model) = p.atmos
-    (; ᶜp, ᶜΦ, ᶜh_tot) = p
-    (; ᶜtke⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜh_tot⁰, ᶜq_tot⁰) = p
+    (; ᶜΦ,) = p.core
+    (; ᶜp, ᶜh_tot) = p.precomputed
+    (; ᶜtke⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜh_tot⁰, ᶜq_tot⁰) =
+        p.precomputed
 
     @. ᶜρa⁰ = ρa⁰(Y.c)
     @. ᶜtke⁰ = divide_by_ρa(Y.c.sgs⁰.ρatke, ᶜρa⁰, 0, Y.c.ρ, turbconv_model)
@@ -59,9 +61,10 @@ function set_prognostic_edmf_precomputed_quantities_draft_and_bc!(Y, p, ᶠuₕ
     (; params) = p
     thermo_params = CAP.thermodynamics_params(params)
 
-    (; ᶜspecific, ᶜp, ᶜΦ, ᶜh_tot) = p
-    (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜtsʲs, ᶜρʲs) = p
-    (; ustar, obukhov_length, buoyancy_flux) = p.sfc_conditions
+    (; ᶜΦ,) = p.core
+    (; ᶜspecific, ᶜp, ᶜh_tot) = p.precomputed
+    (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜtsʲs, ᶜρʲs) = p.precomputed
+    (; ustar, obukhov_length, buoyancy_flux) = p.precomputed.sfc_conditions
 
     for j in 1:n
         ᶜuʲ = ᶜuʲs.:($j)
@@ -90,7 +93,8 @@ function set_prognostic_edmf_precomputed_quantities_draft_and_bc!(Y, p, ᶠuₕ
         )
         ᶜρ_int_val = Fields.field_values(Fields.level(Y.c.ρ, 1))
         ᶜp_int_val = Fields.field_values(Fields.level(ᶜp, 1))
-        (; ρ_flux_h_tot, ρ_flux_q_tot, ustar, obukhov_length) = p.sfc_conditions
+        (; ρ_flux_h_tot, ρ_flux_q_tot, ustar, obukhov_length) =
+            p.precomputed.sfc_conditions
         buoyancy_flux_val = Fields.field_values(buoyancy_flux)
         ρ_flux_h_tot_val = Fields.field_values(ρ_flux_h_tot)
         ρ_flux_q_tot_val = Fields.field_values(ρ_flux_q_tot)
@@ -168,8 +172,9 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
     FT = eltype(params)
     n = n_mass_flux_subdomains(turbconv_model)
 
-    (; ᶜu, ᶜp, ᶜρ_ref) = p
-    (; ᶜtke⁰, ᶜρa⁰, ᶜu⁰, ᶠu³⁰, ᶜts⁰, ᶜρ⁰, ᶜq_tot⁰) = p
+    (; ᶜρ_ref) = p.core
+    (; ᶜspecific, ᶜtke⁰, ᶜu, ᶜp, ᶜρa⁰, ᶜu⁰, ᶠu³⁰, ᶜts⁰, ᶜρ⁰, ᶜq_tot⁰) =
+        p.precomputed
     (;
         ᶜmixing_length,
         ᶜlinear_buoygrad,
@@ -177,9 +182,9 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
         ᶜK_u,
         ᶜK_h,
         ρatke_flux,
-    ) = p
-    (; ᶜuʲs, ᶠu³ʲs, ᶜtsʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs) = p
-    (; ustar, obukhov_length, buoyancy_flux) = p.sfc_conditions
+    ) = p.precomputed
+    (; ᶜuʲs, ᶜtsʲs, ᶠu³ʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs) = p.precomputed
+    (; ustar, obukhov_length, buoyancy_flux) = p.precomputed.sfc_conditions
 
     ᶜz = Fields.coordinate_field(Y.c).z
     z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl
index 368b163a37a..5806d564764 100644
--- a/src/callbacks/callbacks.jl
+++ b/src/callbacks/callbacks.jl
@@ -42,7 +42,8 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator)
     p = integrator.p
     t = integrator.t
 
-    (; ᶜts, sfc_conditions, params) = p
+    (; ᶜts, sfc_conditions) = p.precomputed
+    (; params) = p
     (; idealized_insolation, idealized_h2o, idealized_clouds) = p.radiation
     (; insolation_tuple, ᶠradiation_flux, radiation_model) = p.radiation
 
@@ -182,7 +183,7 @@ function common_diagnostics(p, ᶜu, ᶜts)
         temperature = TD.air_temperature.(thermo_params, ᶜts),
         potential_temperature = TD.dry_pottemp.(thermo_params, ᶜts),
         specific_enthalpy = TD.specific_enthalpy.(thermo_params, ᶜts),
-        buoyancy = CAP.grav(p.params) .* (p.ᶜρ_ref .- ᶜρ) ./ ᶜρ,
+        buoyancy = CAP.grav(p.params) .* (p.core.ᶜρ_ref .- ᶜρ) ./ ᶜρ,
         density = TD.air_density.(thermo_params, ᶜts),
     )
     if !(p.atmos.moisture_model isa DryModel)
@@ -196,7 +197,7 @@ function common_diagnostics(p, ᶜu, ᶜts)
             cloud_fraction_gm = get_cloud_fraction.(
                 thermo_params,
                 env_thermo_quad,
-                p.ᶜp,
+                p.precomputed.ᶜp,
                 ᶜts,
             ),
         )
@@ -219,7 +220,7 @@ NVTX.@annotate function compute_diagnostics(integrator)
     FT = eltype(params)
     thermo_params = CAP.thermodynamics_params(params)
 
-    (; ᶜu, ᶜK, ᶜts, ᶜp, sfc_conditions) = p
+    (; ᶜu, ᶜK, ᶜts, ᶜp, sfc_conditions) = p.precomputed
     dycore_diagnostic = (;
         common_diagnostics(p, ᶜu, ᶜts)...,
         pressure = ᶜp,
@@ -263,7 +264,7 @@ NVTX.@annotate function compute_diagnostics(integrator)
     end
 
     if p.atmos.turbconv_model isa PrognosticEDMFX
-        (; ᶜtke⁰, ᶜu⁰, ᶜts⁰, ᶜmixing_length) = p
+        (; ᶜtke⁰, ᶜu⁰, ᶜts⁰, ᶜmixing_length) = p.precomputed
         (; ᶜu⁺, ᶜts⁺, ᶜa⁺, ᶜa⁰) = output_prognostic_sgs_quantities(Y, p, t)
         env_diagnostics = (;
             common_diagnostics(p, ᶜu⁰, ᶜts⁰)...,
@@ -294,7 +295,7 @@ NVTX.@annotate function compute_diagnostics(integrator)
             ) .+ ᶜa⁺ .* cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺),
         )
     elseif p.atmos.turbconv_model isa DiagnosticEDMFX
-        (; ᶜtke⁰, ᶜmixing_length) = p
+        (; ᶜtke⁰, ᶜmixing_length) = p.precomputed
         (; ᶜu⁺, ᶜts⁺, ᶜa⁺) = output_diagnostic_sgs_quantities(Y, p, t)
         env_diagnostics = (;
             cloud_fraction = get_cloud_fraction.(
@@ -330,7 +331,7 @@ NVTX.@annotate function compute_diagnostics(integrator)
             Fields.level(Fields.local_geometry_field(Y.f), Fields.half)
         surface_ct3_unit =
             CT3.(unit_basis_vector_data.(CT3, sfc_local_geometry))
-        (; ρ_flux_uₕ, ρ_flux_h_tot) = p.sfc_conditions
+        (; ρ_flux_uₕ, ρ_flux_h_tot) = p.precomputed.sfc_conditions
         sfc_flux_momentum =
             Geometry.UVVector.(
                 adjoint.(ρ_flux_uₕ ./ Spaces.level(ᶠinterp.(Y.c.ρ), half)) .*
@@ -342,7 +343,7 @@ NVTX.@annotate function compute_diagnostics(integrator)
             sfc_flux_energy = dot.(ρ_flux_h_tot, surface_ct3_unit),
         )
         if :ρq_tot in propertynames(Y.c)
-            (; ρ_flux_q_tot) = p.sfc_conditions
+            (; ρ_flux_q_tot) = p.precomputed.sfc_conditions
             vert_diff_diagnostic = (;
                 vert_diff_diagnostic...,
                 sfc_evaporation = dot.(ρ_flux_q_tot, surface_ct3_unit),
diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl
index 7de8b7838bb..3595e8a5255 100644
--- a/src/diagnostics/core_diagnostics.jl
+++ b/src/diagnostics/core_diagnostics.jl
@@ -74,9 +74,9 @@ add_diagnostic_variable!(
     comments = "Eastward (zonal) wind component",
     compute! = (out, state, cache, time) -> begin
         if isnothing(out)
-            return copy(Geometry.UVector.(cache.ᶜu).components.data.:1)
+            return copy(Geometry.UVector.(cache.precomputed.ᶜu).components.data.:1)
         else
-            out .= Geometry.UVector.(cache.ᶜu).components.data.:1
+            out .= Geometry.UVector.(cache.precomputed.ᶜu).components.data.:1
         end
     end,
 )
@@ -92,9 +92,9 @@ add_diagnostic_variable!(
     comments = "Northward (meridional) wind component",
     compute! = (out, state, cache, time) -> begin
         if isnothing(out)
-            return copy(Geometry.VVector.(cache.ᶜu).components.data.:1)
+            return copy(Geometry.VVector.(cache.precomputed.ᶜu).components.data.:1)
         else
-            out .= Geometry.VVector.(cache.ᶜu).components.data.:1
+            out .= Geometry.VVector.(cache.precomputed.ᶜu).components.data.:1
         end
     end,
 )
@@ -113,9 +113,9 @@ add_diagnostic_variable!(
     comments = "Vertical wind component",
     compute! = (out, state, cache, time) -> begin
         if isnothing(out)
-            return copy(Geometry.WVector.(cache.ᶜu).components.data.:1)
+            return copy(Geometry.WVector.(cache.precomputed.ᶜu).components.data.:1)
         else
-            out .= Geometry.WVector.(cache.ᶜu).components.data.:1
+            out .= Geometry.WVector.(cache.precomputed.ᶜu).components.data.:1
         end
     end,
 )
@@ -131,9 +131,9 @@ add_diagnostic_variable!(
     compute! = (out, state, cache, time) -> begin
         thermo_params = CAP.thermodynamics_params(cache.params)
         if isnothing(out)
-            return TD.air_temperature.(thermo_params, cache.ᶜts)
+            return TD.air_temperature.(thermo_params, cache.precomputed.ᶜts)
         else
-            out .= TD.air_temperature.(thermo_params, cache.ᶜts)
+            out .= TD.air_temperature.(thermo_params, cache.precomputed.ᶜts)
         end
     end,
 )
@@ -149,9 +149,9 @@ add_diagnostic_variable!(
     compute! = (out, state, cache, time) -> begin
         thermo_params = CAP.thermodynamics_params(cache.params)
         if isnothing(out)
-            return TD.dry_pottemp.(thermo_params, cache.ᶜts)
+            return TD.dry_pottemp.(thermo_params, cache.precomputed.ᶜts)
         else
-            out .= TD.dry_pottemp.(thermo_params, cache.ᶜts)
+            out .= TD.dry_pottemp.(thermo_params, cache.precomputed.ᶜts)
         end
     end,
 )
@@ -166,9 +166,9 @@ add_diagnostic_variable!(
     compute! = (out, state, cache, time) -> begin
         thermo_params = CAP.thermodynamics_params(cache.params)
         if isnothing(out)
-            return TD.specific_enthalpy.(thermo_params, cache.ᶜts)
+            return TD.specific_enthalpy.(thermo_params, cache.precomputed.ᶜts)
         else
-            out .= TD.specific_enthalpy.(thermo_params, cache.ᶜts)
+            out .= TD.specific_enthalpy.(thermo_params, cache.precomputed.ᶜts)
         end
     end,
 )
@@ -182,9 +182,9 @@ add_diagnostic_variable!(
     units = "Pa",
     compute! = (out, state, cache, time) -> begin
         if isnothing(out)
-            return copy(cache.ᶜp)
+            return copy(cache.precomputed.ᶜp)
         else
-            out .= cache.ᶜp
+            out .= cache.precomputed.ᶜp
         end
     end,
 )
@@ -228,9 +228,9 @@ function compute_hur!(
 ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}}
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.relative_humidity.(thermo_params, cache.ᶜts)
+        return TD.relative_humidity.(thermo_params, cache.precomputed.ᶜts)
     else
-        out .= TD.relative_humidity.(thermo_params, cache.ᶜts)
+        out .= TD.relative_humidity.(thermo_params, cache.precomputed.ᶜts)
     end
 end
 
@@ -291,9 +291,13 @@ function compute_clw!(
 ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}}
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.liquid_specific_humidity.(thermo_params, cache.ᶜts)
+        return TD.liquid_specific_humidity.(
+            thermo_params,
+            cache.precomputed.ᶜts,
+        )
     else
-        out .= TD.liquid_specific_humidity.(thermo_params, cache.ᶜts)
+        out .=
+            TD.liquid_specific_humidity.(thermo_params, cache.precomputed.ᶜts)
     end
 end
 
@@ -303,9 +307,9 @@ add_diagnostic_variable!(
     standard_name = "mass_fraction_of_cloud_liquid_water_in_air",
     units = "kg kg^-1",
     comments = """
-    Includes both large-scale and convective cloud. 
-    This is calculated as the mass of cloud liquid water in the grid cell divided by 
-    the mass of air (including the water in all phases) in the grid cells. 
+    Includes both large-scale and convective cloud.
+    This is calculated as the mass of cloud liquid water in the grid cell divided by
+    the mass of air (including the water in all phases) in the grid cells.
     """,
     compute! = compute_clw!,
 )
@@ -327,9 +331,9 @@ function compute_cli!(
 ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}}
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.ice_specific_humidity.(thermo_params, cache.ᶜts)
+        return TD.ice_specific_humidity.(thermo_params, cache.precomputed.ᶜts)
     else
-        out .= TD.ice_specific_humidity.(thermo_params, cache.ᶜts)
+        out .= TD.ice_specific_humidity.(thermo_params, cache.precomputed.ᶜts)
     end
 end
 
@@ -339,8 +343,8 @@ add_diagnostic_variable!(
     standard_name = "mass_fraction_of_cloud_ice_in_air",
     units = "kg kg^-1",
     comments = """
-    Includes both large-scale and convective cloud. 
-    This is calculated as the mass of cloud ice in the grid cell divided by 
+    Includes both large-scale and convective cloud.
+    This is calculated as the mass of cloud ice in the grid cell divided by
     the mass of air (including the water in all phases) in the grid cell.
     """,
     compute! = compute_cli!,
@@ -365,11 +369,14 @@ function compute_hussfc!(
     if isnothing(out)
         return TD.total_specific_humidity.(
             thermo_params,
-            cache.sfc_conditions.ts,
+            cache.precomputed.sfc_conditions.ts,
         )
     else
         out .=
-            TD.total_specific_humidity.(thermo_params, cache.sfc_conditions.ts)
+            TD.total_specific_humidity.(
+                thermo_params,
+                cache.precomputed.sfc_conditions.ts,
+            )
     end
 end
 
@@ -394,9 +401,16 @@ add_diagnostic_variable!(
     compute! = (out, state, cache, time) -> begin
         thermo_params = CAP.thermodynamics_params(cache.params)
         if isnothing(out)
-            return TD.air_temperature.(thermo_params, cache.sfc_conditions.ts)
+            return TD.air_temperature.(
+                thermo_params,
+                cache.precomputed.sfc_conditions.ts,
+            )
         else
-            out .= TD.air_temperature.(thermo_params, cache.sfc_conditions.ts)
+            out .=
+                TD.air_temperature.(
+                    thermo_params,
+                    cache.precomputed.sfc_conditions.ts,
+                )
         end
     end,
 )
@@ -411,7 +425,7 @@ function compute_tau!(out, state, cache, component, energy_form::TotalEnergy)
     sfc_local_geometry =
         Fields.level(Fields.local_geometry_field(state.f), Fields.half)
     surface_ct3_unit = CT3.(unit_basis_vector_data.(CT3, sfc_local_geometry))
-    (; ρ_flux_uₕ) = cache.sfc_conditions
+    (; ρ_flux_uₕ) = cache.precomputed.sfc_conditions
 
     if isnothing(out)
         return getproperty(
@@ -464,7 +478,7 @@ compute_hfes!(_, _, _, _, energy_form::T) where {T} =
     error_diagnostic_variable("hfes", energy_form)
 
 function compute_hfes!(out, state, cache, time, energy_form::TotalEnergy)
-    (; ρ_flux_h_tot) = cache.sfc_conditions
+    (; ρ_flux_h_tot) = cache.precomputed.sfc_conditions
     sfc_local_geometry =
         Fields.level(Fields.local_geometry_field(state.f), Fields.half)
     surface_ct3_unit = CT3.(unit_basis_vector_data.(CT3, sfc_local_geometry))
@@ -513,7 +527,7 @@ function compute_evspsbl!(
     moisture_model::T,
     energy_form::TotalEnergy,
 ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}}
-    (; ρ_flux_q_tot) = cache.sfc_conditions
+    (; ρ_flux_q_tot) = cache.precomputed.sfc_conditions
     sfc_local_geometry =
         Fields.level(Fields.local_geometry_field(state.f), Fields.half)
     surface_ct3_unit = CT3.(unit_basis_vector_data.(CT3, sfc_local_geometry))
diff --git a/src/diagnostics/edmfx_diagnostics.jl b/src/diagnostics/edmfx_diagnostics.jl
index 54a2b076209..1c301876283 100644
--- a/src/diagnostics/edmfx_diagnostics.jl
+++ b/src/diagnostics/edmfx_diagnostics.jl
@@ -15,13 +15,13 @@ function compute_arup!(out, state, cache, time, turbconv_model::PrognosticEDMFX)
     if isnothing(out)
         return draft_area.(
             (state.c.sgsʲs.:1).ρa,
-            TD.air_density.(thermo_params, cache.ᶜtsʲs.:1),
+            TD.air_density.(thermo_params, cache.precomputed.ᶜtsʲs.:1),
         )
     else
         out .=
             draft_area.(
                 (state.c.sgsʲs.:1).ρa,
-                TD.air_density.(thermo_params, cache.ᶜtsʲs.:1),
+                TD.air_density.(thermo_params, cache.precomputed.ᶜtsʲs.:1),
             )
     end
 end
@@ -30,14 +30,14 @@ function compute_arup!(out, state, cache, time, turbconv_model::DiagnosticEDMFX)
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
         return draft_area.(
-            cache.ᶜρaʲs.:1,
-            TD.air_density.(thermo_params, cache.ᶜtsʲs.:1),
+            cache.precomputed.ᶜρaʲs.:1,
+            TD.air_density.(thermo_params, cache.precomputed.ᶜtsʲs.:1),
         )
     else
         out .=
             draft_area.(
-                cache.ᶜρaʲs.:1,
-                TD.air_density.(thermo_params, cache.ᶜtsʲs.:1),
+                cache.precomputed.ᶜρaʲs.:1,
+                TD.air_density.(thermo_params, cache.precomputed.ᶜtsʲs.:1),
             )
     end
 end
@@ -67,9 +67,9 @@ function compute_rhoaup!(
 )
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.air_density.(thermo_params, cache.ᶜtsʲs.:1)
+        return TD.air_density.(thermo_params, cache.precomputed.ᶜtsʲs.:1)
     else
-        out .= TD.air_density.(thermo_params, cache.ᶜtsʲs.:1)
+        out .= TD.air_density.(thermo_params, cache.precomputed.ᶜtsʲs.:1)
     end
 end
 
@@ -97,9 +97,11 @@ function compute_waup!(
     turbconv_model::Union{PrognosticEDMFX, DiagnosticEDMFX},
 )
     if isnothing(out)
-        return copy(Geometry.WVector.(cache.ᶜuʲs.:1).components.data.:1)
+        return copy(
+            Geometry.WVector.(cache.precomputed.ᶜuʲs.:1).components.data.:1,
+        )
     else
-        out .= Geometry.WVector.(cache.ᶜuʲs.:1).components.data.:1
+        out .= Geometry.WVector.(cache.precomputed.ᶜuʲs.:1).components.data.:1
     end
 end
 
@@ -128,9 +130,9 @@ function compute_taup!(
 )
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.air_temperature.(thermo_params, cache.ᶜtsʲs.:1)
+        return TD.air_temperature.(thermo_params, cache.precomputed.ᶜtsʲs.:1)
     else
-        out .= TD.air_temperature.(thermo_params, cache.ᶜtsʲs.:1)
+        out .= TD.air_temperature.(thermo_params, cache.precomputed.ᶜtsʲs.:1)
     end
 end
 
@@ -159,9 +161,9 @@ function compute_thetaaup!(
 )
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.dry_pottemp.(thermo_params, cache.ᶜtsʲs.:1)
+        return TD.dry_pottemp.(thermo_params, cache.precomputed.ᶜtsʲs.:1)
     else
-        out .= TD.dry_pottemp.(thermo_params, cache.ᶜtsʲs.:1)
+        out .= TD.dry_pottemp.(thermo_params, cache.precomputed.ᶜtsʲs.:1)
     end
 end
 
@@ -190,9 +192,9 @@ function compute_haup!(
 )
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.specific_enthalpy.(thermo_params, cache.ᶜtsʲs.:1)
+        return TD.specific_enthalpy.(thermo_params, cache.precomputed.ᶜtsʲs.:1)
     else
-        out .= TD.specific_enthalpy.(thermo_params, cache.ᶜtsʲs.:1)
+        out .= TD.specific_enthalpy.(thermo_params, cache.precomputed.ᶜtsʲs.:1)
     end
 end
 
@@ -236,9 +238,16 @@ function compute_husup!(
 )
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.total_specific_humidity.(thermo_params, cache.ᶜtsʲs.:1)
+        return TD.total_specific_humidity.(
+            thermo_params,
+            cache.precomputed.ᶜtsʲs.:1,
+        )
     else
-        out .= TD.total_specific_humidity.(thermo_params, cache.ᶜtsʲs.:1)
+        out .=
+            TD.total_specific_humidity.(
+                thermo_params,
+                cache.precomputed.ᶜtsʲs.:1,
+            )
     end
 end
 
@@ -282,9 +291,9 @@ function compute_hurup!(
 )
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.relative_humidity.(thermo_params, cache.ᶜtsʲs.:1)
+        return TD.relative_humidity.(thermo_params, cache.precomputed.ᶜtsʲs.:1)
     else
-        out .= TD.relative_humidity.(thermo_params, cache.ᶜtsʲs.:1)
+        out .= TD.relative_humidity.(thermo_params, cache.precomputed.ᶜtsʲs.:1)
     end
 end
 
@@ -328,9 +337,16 @@ function compute_clwup!(
 )
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.liquid_specific_humidity.(thermo_params, cache.ᶜtsʲs.:1)
+        return TD.liquid_specific_humidity.(
+            thermo_params,
+            cache.precomputed.ᶜtsʲs.:1,
+        )
     else
-        out .= TD.liquid_specific_humidity.(thermo_params, cache.ᶜtsʲs.:1)
+        out .=
+            TD.liquid_specific_humidity.(
+                thermo_params,
+                cache.precomputed.ᶜtsʲs.:1,
+            )
     end
 end
 
@@ -339,7 +355,7 @@ add_diagnostic_variable!(
     long_name = "Updraft Mass Fraction of Cloud Liquid Water",
     units = "kg kg^-1",
     comments = """
-    This is calculated as the mass of cloud liquid water in the first updraft divided by 
+    This is calculated as the mass of cloud liquid water in the first updraft divided by
     the mass of air (including the water in all phases) in the first updraft.
     """,
     compute! = compute_clwup!,
@@ -377,9 +393,13 @@ function compute_cliup!(
 )
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.ice_specific_humidity.(thermo_params, cache.ᶜtsʲs.:1)
+        return TD.ice_specific_humidity.(
+            thermo_params,
+            cache.precomputed.ᶜtsʲs.:1,
+        )
     else
-        out .= TD.ice_specific_humidity.(thermo_params, cache.ᶜtsʲs.:1)
+        out .=
+            TD.ice_specific_humidity.(thermo_params, cache.precomputed.ᶜtsʲs.:1)
     end
 end
 
@@ -388,7 +408,7 @@ add_diagnostic_variable!(
     long_name = "Updraft Mass Fraction of Cloud Ice",
     units = "kg kg^-1",
     comments = """
-    This is calculated as the mass of cloud ice in the first updraft divided by 
+    This is calculated as the mass of cloud ice in the first updraft divided by
     the mass of air (including the water in all phases) in the first updraft.
     """,
     compute! = compute_cliup!,
@@ -406,12 +426,15 @@ function compute_aren!(out, state, cache, time, turbconv_model::PrognosticEDMFX)
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
         return draft_area.(
-            cache.ᶜρa⁰,
-            TD.air_density.(thermo_params, cache.ᶜts⁰),
+            cache.precomputed.ᶜρa⁰,
+            TD.air_density.(thermo_params, cache.precomputed.ᶜts⁰),
         )
     else
         out .=
-            draft_area.(cache.ᶜρa⁰, TD.air_density.(thermo_params, cache.ᶜts⁰))
+            draft_area.(
+                cache.precomputed.ᶜρa⁰,
+                TD.air_density.(thermo_params, cache.precomputed.ᶜts⁰),
+            )
     end
 end
 
@@ -419,14 +442,14 @@ function compute_aren!(out, state, cache, time, turbconv_model::DiagnosticEDMFX)
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
         return draft_area.(
-            cache.ᶜρaʲs.:1,
-            TD.air_density.(thermo_params, cache.ᶜtsʲs.:1),
+            cache.precomputed.ᶜρaʲs.:1,
+            TD.air_density.(thermo_params, cache.precomputed.ᶜtsʲs.:1),
         )
     else
         out .=
             draft_area.(
-                cache.ᶜρaʲs.:1,
-                TD.air_density.(thermo_params, cache.ᶜtsʲs.:1),
+                cache.precomputed.ᶜρaʲs.:1,
+                TD.air_density.(thermo_params, cache.precomputed.ᶜtsʲs.:1),
             )
     end
 end
@@ -455,9 +478,9 @@ function compute_rhoaen!(
 )
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.air_density.(thermo_params, cache.ᶜts⁰)
+        return TD.air_density.(thermo_params, cache.precomputed.ᶜts⁰)
     else
-        out .= TD.air_density.(thermo_params, cache.ᶜts⁰)
+        out .= TD.air_density.(thermo_params, cache.precomputed.ᶜts⁰)
     end
 end
 
@@ -484,9 +507,9 @@ function compute_waen!(
     turbconv_model::Union{PrognosticEDMFX, DiagnosticEDMFX},
 )
     if isnothing(out)
-        return copy(Geometry.WVector.(cache.ᶜu⁰).components.data.:1)
+        return copy(Geometry.WVector.(cache.precomputed.ᶜu⁰).components.data.:1)
     else
-        out .= Geometry.WVector.(cache.ᶜu⁰).components.data.:1
+        out .= Geometry.WVector.(cache.precomputed.ᶜu⁰).components.data.:1
     end
 end
 
@@ -509,9 +532,9 @@ compute_taen!(_, _, _, _, turbconv_model::T) where {T} =
 function compute_taen!(out, state, cache, time, turbconv_model::PrognosticEDMFX)
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.air_temperature.(thermo_params, cache.ᶜts⁰)
+        return TD.air_temperature.(thermo_params, cache.precomputed.ᶜts⁰)
     else
-        out .= TD.air_temperature.(thermo_params, cache.ᶜts⁰)
+        out .= TD.air_temperature.(thermo_params, cache.precomputed.ᶜts⁰)
     end
 end
 
@@ -539,9 +562,9 @@ function compute_thetaaen!(
 )
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.dry_pottemp.(thermo_params, cache.ᶜts⁰)
+        return TD.dry_pottemp.(thermo_params, cache.precomputed.ᶜts⁰)
     else
-        out .= TD.dry_pottemp.(thermo_params, cache.ᶜts⁰)
+        out .= TD.dry_pottemp.(thermo_params, cache.precomputed.ᶜts⁰)
     end
 end
 
@@ -563,9 +586,9 @@ compute_haen!(_, _, _, _, turbconv_model::T) where {T} =
 function compute_haen!(out, state, cache, time, turbconv_model::PrognosticEDMFX)
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.dry_pottemp.(thermo_params, cache.ᶜts⁰)
+        return TD.dry_pottemp.(thermo_params, cache.precomputed.ᶜts⁰)
     else
-        out .= TD.dry_pottemp.(thermo_params, cache.ᶜts⁰)
+        out .= TD.dry_pottemp.(thermo_params, cache.precomputed.ᶜts⁰)
     end
 end
 
@@ -608,9 +631,13 @@ function compute_husen!(
 )
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.total_specific_humidity.(thermo_params, cache.ᶜts⁰)
+        return TD.total_specific_humidity.(
+            thermo_params,
+            cache.precomputed.ᶜts⁰,
+        )
     else
-        out .= TD.total_specific_humidity.(thermo_params, cache.ᶜts⁰)
+        out .=
+            TD.total_specific_humidity.(thermo_params, cache.precomputed.ᶜts⁰)
     end
 end
 
@@ -653,9 +680,9 @@ function compute_huren!(
 )
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.relative_humidity.(thermo_params, cache.ᶜts⁰)
+        return TD.relative_humidity.(thermo_params, cache.precomputed.ᶜts⁰)
     else
-        out .= TD.relative_humidity.(thermo_params, cache.ᶜts⁰)
+        out .= TD.relative_humidity.(thermo_params, cache.precomputed.ᶜts⁰)
     end
 end
 
@@ -698,9 +725,13 @@ function compute_clwen!(
 )
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.liquid_specific_humidity.(thermo_params, cache.ᶜts⁰)
+        return TD.liquid_specific_humidity.(
+            thermo_params,
+            cache.precomputed.ᶜts⁰,
+        )
     else
-        out .= TD.liquid_specific_humidity.(thermo_params, cache.ᶜts⁰)
+        out .=
+            TD.liquid_specific_humidity.(thermo_params, cache.precomputed.ᶜts⁰)
     end
 end
 
@@ -709,7 +740,7 @@ add_diagnostic_variable!(
     long_name = "Envrionment Mass Fraction of Cloud Liquid Water",
     units = "kg kg^-1",
     comments = """
-    This is calculated as the mass of cloud liquid water in the environment divided by 
+    This is calculated as the mass of cloud liquid water in the environment divided by
     the mass of air (including the water in all phases) in the environment.
     """,
     compute! = compute_clwen!,
@@ -747,9 +778,9 @@ function compute_clien!(
 )
     thermo_params = CAP.thermodynamics_params(cache.params)
     if isnothing(out)
-        return TD.ice_specific_humidity.(thermo_params, cache.ᶜts⁰)
+        return TD.ice_specific_humidity.(thermo_params, cache.precomputed.ᶜts⁰)
     else
-        out .= TD.ice_specific_humidity.(thermo_params, cache.ᶜts⁰)
+        out .= TD.ice_specific_humidity.(thermo_params, cache.precomputed.ᶜts⁰)
     end
 end
 
@@ -758,7 +789,7 @@ add_diagnostic_variable!(
     long_name = "Environment Mass Fraction of Cloud Ice",
     units = "kg kg^-1",
     comments = """
-    This is calculated as the mass of cloud ice in the environment divided by 
+    This is calculated as the mass of cloud ice in the environment divided by
     the mass of air (including the water in all phases) in the environment.
     """,
     compute! = compute_clien!,
@@ -780,9 +811,9 @@ function compute_lmix!(
     turbconv_model::Union{PrognosticEDMFX, DiagnosticEDMFX},
 )
     if isnothing(out)
-        return copy(cache.ᶜmixing_length)
+        return copy(cache.precomputed.ᶜmixing_length)
     else
-        out .= cache.ᶜmixing_length
+        out .= cache.precomputed.ᶜmixing_length
     end
 end
 
@@ -809,9 +840,9 @@ function compute_tke!(
     turbconv_model::Union{PrognosticEDMFX, DiagnosticEDMFX},
 )
     if isnothing(out)
-        return copy(cache.ᶜtke⁰)
+        return copy(cache.precomputed.ᶜtke⁰)
     else
-        out .= cache.ᶜtke⁰
+        out .= cache.precomputed.ᶜtke⁰
     end
 end
 
diff --git a/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl b/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl
index 408b16ec00b..c71e424ae30 100644
--- a/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl
+++ b/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl
@@ -113,7 +113,9 @@ function non_orographic_gravity_wave_tendency!(
     ::NonOrographyGravityWave,
 )
     #unpack
-    (; ᶜts, ᶜT, params) = p
+    (; ᶜT,) = p.core
+    (; ᶜts) = p.precomputed
+    (; params) = p
     (; ᶜdTdz, ᶜbuoyancy_frequency) = p.non_orographic_gravity_wave
     (; model_config) = p.atmos
     (;
@@ -166,7 +168,7 @@ function non_orographic_gravity_wave_tendency!(
             parent(damp_level[colidx]) .= length(parent(ᶜz[colidx]))
         end
     elseif model_config isa SphericalModel
-        (; ᶜp) = p
+        (; ᶜp) = p.precomputed
         # source level: the index of the highest level whose pressure is higher than source pressure
         source_level = similar(Fields.level(Y.c.ρ, 1))
         Fields.bycolumn(axes(ᶜρ)) do colidx
diff --git a/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave.jl b/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave.jl
index 202ef7a2c4b..23051648f79 100644
--- a/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave.jl
+++ b/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave.jl
@@ -75,7 +75,9 @@ function orographic_gravity_wave_cache(Y, ogw::OrographicGravityWave)
 end
 
 function orographic_gravity_wave_tendency!(Yₜ, Y, p, t, ::OrographicGravityWave)
-    (; params, ᶜts, ᶜT, ᶜp) = p
+    (; ᶜT) = p.core
+    (; params) = p
+    (; ᶜts, ᶜp) = p.precomputed
     (; ᶜdTdz) = p.orographic_gravity_wave
     (;
         topo_k_pbl,
diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl
index 42b09859fe5..d970e0b710a 100644
--- a/src/parameterized_tendencies/microphysics/precipitation.jl
+++ b/src/parameterized_tendencies/microphysics/precipitation.jl
@@ -35,7 +35,8 @@ function precipitation_cache(Y, precip_model::Microphysics0Moment)
 end
 
 function compute_precipitation_cache!(Y, p, colidx, ::Microphysics0Moment, _)
-    (; ᶜts, params) = p
+    (; params) = p
+    (; ᶜts) = p.precomputed
     (; ᶜS_ρq_tot) = p.precipitation
     cm_params = CAP.microphysics_params(params)
     thermo_params = CAP.thermodynamics_params(params)
@@ -57,7 +58,7 @@ function compute_precipitation_cache!(
     # For environment we multiply by grid mean ρ and not byᶜρa⁰
     # I.e. assuming a⁰=1
 
-    (; ᶜS_q_tot⁰, ᶜS_q_totʲs, ᶜρaʲs) = p
+    (; ᶜS_q_tot⁰, ᶜS_q_totʲs, ᶜρaʲs) = p.precomputed
     (; ᶜS_ρq_tot) = p.precipitation
     n = n_mass_flux_subdomains(p.atmos.turbconv_model)
     ρ = Y.c.ρ
@@ -76,7 +77,10 @@ function precipitation_tendency!(
     colidx,
     precip_model::Microphysics0Moment,
 )
-    (; ᶜts, ᶜΦ, ᶜT, params, turbconv_model) = p # assume ᶜts has been updated
+    (; ᶜT, ᶜΦ) = p.core
+    (; ᶜts,) = p.precomputed  # assume ᶜts has been updated
+    (; params) = p
+    (; turbconv_model) = p.atmos
     (;
         ᶜ3d_rain,
         ᶜ3d_snow,
@@ -118,11 +122,11 @@ function precipitation_tendency!(
         if turbconv_model isa DiagnosticEDMFX && !isnothing(Yₜ)
             @. Yₜ.c.ρe_tot[colidx] +=
                 sum(
-                    p.ᶜS_q_totʲs[colidx] *
-                    p.ᶜρaʲs[colidx] *
-                    p.ᶜS_e_totʲs_helper[colidx],
+                    p.precomputed.ᶜS_q_totʲs[colidx] *
+                    p.precomputed.ᶜρaʲs[colidx] *
+                    p.precomputed.ᶜS_e_totʲs_helper[colidx],
                 ) +
-                p.ᶜS_q_tot⁰[colidx] *
+                p.precomputed.ᶜS_q_tot⁰[colidx] *
                 Y.c.ρ[colidx] *
                 e_tot_0M_precipitation_sources_helper(
                     thermo_params,
diff --git a/src/parameterized_tendencies/radiation/held_suarez.jl b/src/parameterized_tendencies/radiation/held_suarez.jl
index c56c50b6fb5..03f8309c23a 100644
--- a/src/parameterized_tendencies/radiation/held_suarez.jl
+++ b/src/parameterized_tendencies/radiation/held_suarez.jl
@@ -46,7 +46,8 @@ function held_suarez_ΔT_y_T_equator(
 end
 
 function forcing_tendency!(Yₜ, Y, p, t, colidx, ::HeldSuarezForcing)
-    (; sfc_conditions, ᶜp, params) = p
+    (; params) = p
+    (; ᶜp, sfc_conditions) = p.precomputed
     (; ᶜσ, ᶜheight_factor, ᶜΔρT, ᶜφ) = p.forcing
 
     # TODO: Don't need to enforce FT here, it should be done at param creation.
diff --git a/src/parameterized_tendencies/radiation/radiation.jl b/src/parameterized_tendencies/radiation/radiation.jl
index 2672f326981..1c23fb91195 100644
--- a/src/parameterized_tendencies/radiation/radiation.jl
+++ b/src/parameterized_tendencies/radiation/radiation.jl
@@ -265,7 +265,8 @@ function radiation_tendency!(
     @assert !(p.atmos.moisture_model isa DryModel)
     @assert p.atmos.energy_form isa TotalEnergy
 
-    (; params, ᶜspecific, ᶜts) = p
+    (; params) = p
+    (; ᶜspecific, ᶜts) = p.precomputed
     (; ᶜκρq, ∫_0_∞_κρq, ᶠ∫_0_z_κρq, isoline_z_ρ_q, ᶠradiation_flux) =
         p.radiation
     thermo_params = CAP.thermodynamics_params(params)
@@ -355,7 +356,7 @@ function radiation_tendency!(
     thermo_params = CAP.thermodynamics_params(params)
     ᶜdTdt_rad = p.radiation.ᶜdTdt_rad[colidx]
     ᶜρ = Y.c.ρ[colidx]
-    ᶜts_gm = p.ᶜts[colidx]
+    ᶜts_gm = p.precomputed.ᶜts[colidx]
     zc = Fields.coordinate_field(axes(ᶜρ)).z
     @. ᶜdTdt_rad = rad(FT(t), zc)
     @. Yₜ.c.ρe_tot[colidx] += ᶜρ * TD.cv_m(thermo_params, ᶜts_gm) * ᶜdTdt_rad
diff --git a/src/parameterized_tendencies/sponge/viscous_sponge.jl b/src/parameterized_tendencies/sponge/viscous_sponge.jl
index 16c9e2517d6..97279acfce3 100644
--- a/src/parameterized_tendencies/sponge/viscous_sponge.jl
+++ b/src/parameterized_tendencies/sponge/viscous_sponge.jl
@@ -30,7 +30,7 @@ add_viscous_sponge_energy_tendency!(Yₜ, Y, p, t) =
 
 function add_viscous_sponge_energy_tendency!(Yₜ, Y, p, t, ::TotalEnergy)
     (; ᶜβ_viscous) = p.viscous_sponge
-    (; ᶜh_tot) = p
+    (; ᶜh_tot) = p.precomputed
     ᶜρ = Y.c.ρ
     @. Yₜ.c.ρe_tot += ᶜβ_viscous * wdivₕ(ᶜρ * gradₕ(ᶜh_tot))
 end
@@ -50,9 +50,7 @@ end
 function viscous_sponge_tendency!(Yₜ, Y, p, t, ::ViscousSponge)
     (; ᶜβ_viscous, ᶠβ_viscous) = p.viscous_sponge
     ᶜuₕ = Y.c.uₕ
-
     add_viscous_sponge_energy_tendency!(Yₜ, Y, p, t)
-
     @. Yₜ.c.uₕ +=
         ᶜβ_viscous * (
             wgradₕ(divₕ(ᶜuₕ)) - Geometry.project(
diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl
index c3a500edcc4..2d290f30bfd 100644
--- a/src/prognostic_equations/advection.jl
+++ b/src/prognostic_equations/advection.jl
@@ -8,12 +8,13 @@ import ClimaCore.Geometry as Geometry
 
 NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t)
     n = n_mass_flux_subdomains(p.atmos.turbconv_model)
-    (; ᶜu, ᶜK, ᶜp, ᶜΦ, ᶜp_ref) = p
+    (; ᶜΦ, ᶜp_ref) = p.core
+    (; ᶜu, ᶜK, ᶜp) = p.precomputed
     if p.atmos.turbconv_model isa AbstractEDMF
-        (; ᶜu⁰) = p
+        (; ᶜu⁰) = p.precomputed
     end
     if p.atmos.turbconv_model isa PrognosticEDMFX
-        (; ᶜuʲs) = p
+        (; ᶜuʲs) = p.precomputed
     end
 
     @. Yₜ.c.ρ -= wdivₕ(Y.c.ρ * ᶜu)
@@ -26,7 +27,7 @@ NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t)
     if :ρθ in propertynames(Y.c)
         @. Yₜ.c.ρθ -= wdivₕ(Y.c.ρθ * ᶜu)
     elseif :ρe_tot in propertynames(Y.c)
-        (; ᶜh_tot) = p
+        (; ᶜh_tot) = p.precomputed
         @. Yₜ.c.ρe_tot -= wdivₕ(Y.c.ρ * ᶜh_tot * ᶜu)
     end
 
@@ -49,9 +50,9 @@ end
 
 NVTX.@annotate function horizontal_tracer_advection_tendency!(Yₜ, Y, p, t)
     n = n_mass_flux_subdomains(p.atmos.turbconv_model)
-    (; ᶜu) = p
+    (; ᶜu) = p.precomputed
     if p.atmos.turbconv_model isa PrognosticEDMFX
-        (; ᶜuʲs) = p
+        (; ᶜuʲs) = p.precomputed
     end
 
     for ρχ_name in filter(is_tracer_var, propertynames(Y.c))
@@ -76,14 +77,17 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     point_type = eltype(Fields.coordinate_field(Y.c))
     (; dt) = p.simulation
     ᶜJ = Fields.local_geometry_field(Y.c).J
-    (; ᶜu, ᶠu³, ᶜK, ᶜf) = p
+    (; ᶜf) = p.core
+    (; ᶜu, ᶠu³, ᶜK) = p.precomputed
     (; edmfx_upwinding) = n > 0 || advect_tke ? p.atmos.numerics : all_nothing
-    (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜρʲs) = n > 0 ? p : all_nothing
-    (; ᶜp, ᶜp_ref, ᶜρ_ref, ᶠgradᵥ_ᶜΦ) = n > 0 ? p : all_nothing
-    (; ᶠu³⁰) = advect_tke ? p : all_nothing
-    ᶜρa⁰ = advect_tke ? (n > 0 ? p.ᶜρa⁰ : Y.c.ρ) : nothing
-    ᶜρ⁰ = advect_tke ? (n > 0 ? p.ᶜρ⁰ : Y.c.ρ) : nothing
-    ᶜtke⁰ = advect_tke ? p.ᶜtke⁰ : nothing
+    (; ᶜp, ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜρʲs) = n > 0 ? p.precomputed : all_nothing
+    (; ᶜp_ref, ᶜρ_ref, ᶠgradᵥ_ᶜΦ) = n > 0 ? p.core : all_nothing
+    (; ᶠu³⁰) = advect_tke ? p.precomputed : all_nothing
+    ᶜρa⁰ = advect_tke ? (n > 0 ? p.precomputed.ᶜρa⁰ : Y.c.ρ) : nothing
+    ᶜρ⁰ = advect_tke ? (n > 0 ? p.precomputed.ᶜρ⁰ : Y.c.ρ) : nothing
+    ᶜtke⁰ =
+        advect_tke ?
+        (n > 0 ? p.precomputed.ᶜspecific⁰.tke : p.precomputed.ᶜtke⁰) : nothing
     ᶜa_scalar = p.scratch.ᶜtemp_scalar
     ᶜω³ = p.scratch.ᶜtemp_CT3
     ᶠω¹² = p.scratch.ᶠtemp_CT12
diff --git a/src/prognostic_equations/edmfx_closures.jl b/src/prognostic_equations/edmfx_closures.jl
index 5e1869d093d..57ae0fa753a 100644
--- a/src/prognostic_equations/edmfx_closures.jl
+++ b/src/prognostic_equations/edmfx_closures.jl
@@ -106,7 +106,9 @@ function edmfx_nh_pressure_tendency!(
 )
 
     n = n_mass_flux_subdomains(turbconv_model)
-    (; params, ᶜρʲs, ᶠgradᵥ_ᶜΦ, ᶠu₃⁰) = p
+    (; params) = p
+    (; ᶠgradᵥ_ᶜΦ) = p.core
+    (; ᶜρʲs, ᶠu₃⁰) = p.precomputed
     FT = eltype(Y)
     ᶜz = Fields.coordinate_field(Y.c).z
     z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl
index e968841bb64..e923160dffc 100644
--- a/src/prognostic_equations/edmfx_entr_detr.jl
+++ b/src/prognostic_equations/edmfx_entr_detr.jl
@@ -327,8 +327,8 @@ function edmfx_entr_detr_tendency!(
 )
 
     n = n_mass_flux_subdomains(turbconv_model)
-    (; ᶜentrʲs, ᶜdetrʲs) = p
-    (; ᶜq_tot⁰, ᶜh_tot⁰, ᶠu₃⁰) = p
+    (; ᶜentrʲs, ᶜdetrʲs) = p.precomputed
+    (; ᶜq_tot⁰, ᶜh_tot⁰, ᶠu₃⁰) = p.precomputed
 
     for j in 1:n
 
diff --git a/src/prognostic_equations/edmfx_sgs_flux.jl b/src/prognostic_equations/edmfx_sgs_flux.jl
index bb159e5ddd1..1337079b490 100644
--- a/src/prognostic_equations/edmfx_sgs_flux.jl
+++ b/src/prognostic_equations/edmfx_sgs_flux.jl
@@ -15,9 +15,9 @@ function edmfx_sgs_mass_flux_tendency!(
 
     n = n_mass_flux_subdomains(turbconv_model)
     (; edmfx_upwinding) = p.atmos.numerics
-    (; ᶠu³, ᶜh_tot, ᶜspecific) = p
-    (; ᶠu³ʲs) = p
-    (; ᶜρa⁰, ᶠu³⁰, ᶜh_tot⁰, ᶜq_tot⁰) = p
+    (; ᶠu³, ᶜh_tot, ᶜspecific) = p.precomputed
+    (; ᶠu³ʲs) = p.precomputed
+    (; ᶜρa⁰, ᶠu³⁰, ᶜh_tot⁰, ᶜq_tot⁰) = p.precomputed
     (; dt) = p.simulation
     ᶜJ = Fields.local_geometry_field(Y.c).J
 
@@ -95,11 +95,10 @@ function edmfx_sgs_mass_flux_tendency!(
     turbconv_model::DiagnosticEDMFX,
 )
 
-    FT = Spaces.undertype(axes(Y.c))
     n = n_mass_flux_subdomains(turbconv_model)
     (; edmfx_upwinding) = p.atmos.numerics
-    (; ᶠu³, ᶜh_tot, ᶜspecific) = p
-    (; ᶜρaʲs, ᶠu³ʲs, ᶜh_totʲs, ᶜq_totʲs) = p
+    (; ᶠu³, ᶜh_tot, ᶜspecific) = p.precomputed
+    (; ᶜρaʲs, ᶠu³ʲs, ᶜh_totʲs, ᶜq_totʲs) = p.precomputed
     (; dt) = p.simulation
     ᶜJ = Fields.local_geometry_field(Y.c).J
 
@@ -158,9 +157,9 @@ function edmfx_sgs_diffusive_flux_tendency!(
 )
 
     FT = Spaces.undertype(axes(Y.c))
-    (; sfc_conditions) = p
-    (; ᶜρa⁰, ᶜu⁰, ᶜh_tot⁰, ᶜq_tot⁰) = p
-    (; ᶜK_u, ᶜK_h) = p
+    (; sfc_conditions) = p.precomputed
+    (; ᶜρa⁰, ᶜu⁰, ᶜh_tot⁰, ᶜq_tot⁰) = p.precomputed
+    (; ᶜK_u, ᶜK_h) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F()
 
     if p.atmos.edmfx_sgs_diffusive_flux
@@ -219,9 +218,9 @@ function edmfx_sgs_diffusive_flux_tendency!(
 )
 
     FT = Spaces.undertype(axes(Y.c))
-    (; sfc_conditions) = p
-    (; ᶜu, ᶜh_tot, ᶜspecific) = p
-    (; ᶜK_u, ᶜK_h) = p
+    (; sfc_conditions) = p.precomputed
+    (; ᶜu, ᶜh_tot, ᶜspecific) = p.precomputed
+    (; ᶜK_u, ᶜK_h) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F()
 
     if p.atmos.edmfx_sgs_diffusive_flux
diff --git a/src/prognostic_equations/edmfx_tke.jl b/src/prognostic_equations/edmfx_tke.jl
index 9ab75ac970a..a39039b0b5c 100644
--- a/src/prognostic_equations/edmfx_tke.jl
+++ b/src/prognostic_equations/edmfx_tke.jl
@@ -18,10 +18,11 @@ function edmfx_tke_tendency!(
     (; params) = p
     turbconv_params = CAP.turbconv_params(params)
     c_d = TCP.tke_diss_coeff(turbconv_params)
-    (; ᶜentrʲs, ᶜdetrʲs, ᶠu³ʲs) = p
-    (; ᶠu³⁰, ᶜstrain_rate_norm, ᶜlinear_buoygrad, ᶜtke⁰, ᶜmixing_length) = p
-    (; ᶜK_u, ᶜK_h, ρatke_flux) = p
-    ᶜρa⁰ = turbconv_model isa PrognosticEDMFX ? p.ᶜρa⁰ : Y.c.ρ
+    (; ᶜentrʲs, ᶜdetrʲs, ᶠu³ʲs) = p.precomputed
+    (; ᶠu³⁰, ᶜstrain_rate_norm, ᶜlinear_buoygrad, ᶜtke⁰, ᶜmixing_length) =
+        p.precomputed
+    (; ᶜK_u, ᶜK_h, ρatke_flux) = p.precomputed
+    ᶜρa⁰ = turbconv_model isa PrognosticEDMFX ? p.precomputed.ᶜρa⁰ : Y.c.ρ
     ᶠgradᵥ = Operators.GradientC2F()
 
     ᶠρaK_u = p.scratch.ᶠtemp_scalar
diff --git a/src/prognostic_equations/forcing/large_scale_advection.jl b/src/prognostic_equations/forcing/large_scale_advection.jl
index 35a0809680d..beebd235a3c 100644
--- a/src/prognostic_equations/forcing/large_scale_advection.jl
+++ b/src/prognostic_equations/forcing/large_scale_advection.jl
@@ -27,7 +27,7 @@ function large_scale_advection_tendency!(
     (; prof_dTdt, prof_dqtdt) = ls_adv
 
     thermo_params = CAP.thermodynamics_params(p.params)
-    ᶜts = p.ᶜts[colidx]
+    ᶜts = p.precomputed.ᶜts[colidx]
     ᶜdqtdt_hadv = p.large_scale_advection.ᶜdqtdt_hadv[colidx]
     ᶜdTdt_hadv = p.large_scale_advection.ᶜdTdt_hadv[colidx]
     z = Fields.coordinate_field(axes(ᶜdqtdt_hadv)).z
diff --git a/src/prognostic_equations/forcing/subsidence.jl b/src/prognostic_equations/forcing/subsidence.jl
index 5fa3f9f8d7a..44f256662a8 100644
--- a/src/prognostic_equations/forcing/subsidence.jl
+++ b/src/prognostic_equations/forcing/subsidence.jl
@@ -40,8 +40,8 @@ function subsidence_tendency!(Yₜ, Y, p, t, colidx, ::Subsidence)
     ᶜ∇MSE_gm = p.subsidence.ᶜ∇MSE_gm[colidx]
     ᶜsubsidence = p.subsidence.ᶜsubsidence[colidx]
     ᶜ∇q_tot_gm = p.subsidence.ᶜ∇q_tot_gm[colidx]
-    ᶜK = p.ᶜK[colidx]
-    ᶜh_tot = p.ᶜh_tot[colidx]
+    ᶜK = p.precomputed.ᶜK[colidx]
+    ᶜh_tot = p.precomputed.ᶜh_tot[colidx]
     ᶜMSE_gm_toa = p.subsidence.ᶜMSE_gm_toa[colidx]
     ᶜq_tot_gm_toa = p.subsidence.ᶜq_tot_gm_toa[colidx]
 
diff --git a/src/prognostic_equations/hyperdiffusion.jl b/src/prognostic_equations/hyperdiffusion.jl
index 040a90c4e1d..1e06f37eaea 100644
--- a/src/prognostic_equations/hyperdiffusion.jl
+++ b/src/prognostic_equations/hyperdiffusion.jl
@@ -64,14 +64,15 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t)
     diffuse_tke = use_prognostic_tke(turbconv_model)
     ᶜJ = Fields.local_geometry_field(Y.c).J
     point_type = eltype(Fields.coordinate_field(Y.c))
-    (; do_dss, ᶜp, ᶜspecific) = p
+    (; do_dss) = p
+    (; ᶜp, ᶜspecific) = p.precomputed
     (; ᶜ∇²u, ᶜ∇²specific_energy) = p.hyperdiff
     if turbconv_model isa PrognosticEDMFX
-        (; ᶜρa⁰, ᶜtke⁰) = p
+        (; ᶜρa⁰, ᶜtke⁰) = p.precomputed
         (; ᶜ∇²tke⁰, ᶜ∇²uₕʲs, ᶜ∇²uᵥʲs, ᶜ∇²uʲs, ᶜ∇²h_totʲs) = p.hyperdiff
     end
     if turbconv_model isa DiagnosticEDMFX
-        (; ᶜtke⁰) = p
+        (; ᶜtke⁰) = p.precomputed
         (; ᶜ∇²tke⁰) = p.hyperdiff
     end
     if do_dss
@@ -79,7 +80,9 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t)
     end
 
     # Grid scale hyperdiffusion
-    @. ᶜ∇²u = C123(wgradₕ(divₕ(p.ᶜu))) - C123(wcurlₕ(C123(curlₕ(p.ᶜu))))
+    @. ᶜ∇²u =
+        C123(wgradₕ(divₕ(p.precomputed.ᶜu))) -
+        C123(wcurlₕ(C123(curlₕ(p.precomputed.ᶜu))))
 
     if :θ in propertynames(ᶜspecific)
         @. ᶜ∇²specific_energy = wdivₕ(gradₕ(ᶜspecific.θ))
@@ -94,8 +97,8 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t)
     if turbconv_model isa PrognosticEDMFX
         for j in 1:n
             @. ᶜ∇²uʲs.:($$j) =
-                C123(wgradₕ(divₕ(p.ᶜuʲs.:($$j)))) -
-                C123(wcurlₕ(C123(curlₕ(p.ᶜuʲs.:($$j)))))
+                C123(wgradₕ(divₕ(p.precomputed.ᶜuʲs.:($$j)))) -
+                C123(wcurlₕ(C123(curlₕ(p.precomputed.ᶜuʲs.:($$j)))))
             @. ᶜ∇²h_totʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).h_tot))
         end
     end
@@ -173,7 +176,8 @@ NVTX.@annotate function tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t)
     (; κ₄) = hyperdiff
     n = n_mass_flux_subdomains(turbconv_model)
 
-    (; do_dss, ᶜspecific) = p
+    (; ᶜspecific) = p.precomputed
+    (; do_dss) = p
     (; ᶜ∇²specific_tracers) = p.hyperdiff
     if turbconv_model isa PrognosticEDMFX
         (; ᶜ∇²q_totʲs) = p.hyperdiff
diff --git a/src/prognostic_equations/implicit/implicit_solver.jl b/src/prognostic_equations/implicit/implicit_solver.jl
index 8448ec25ec9..e7c89ae21ca 100644
--- a/src/prognostic_equations/implicit/implicit_solver.jl
+++ b/src/prognostic_equations/implicit/implicit_solver.jl
@@ -208,19 +208,19 @@ function Wfact!(A, Y, p, dtγ, t)
         # Remove unnecessary values from p to avoid allocations in bycolumn.
         (; energy_form, rayleigh_sponge) = p.atmos
         p′ = (;
-            p.ᶜspecific,
-            p.ᶠu³,
-            p.ᶜK,
-            p.ᶜp,
-            p.∂ᶜK_∂ᶠu₃,
-            p.ᶜΦ,
-            p.ᶠgradᵥ_ᶜΦ,
-            p.ᶜρ_ref,
-            p.ᶜp_ref,
+            p.precomputed.ᶜspecific,
+            p.precomputed.ᶠu³,
+            p.precomputed.ᶜK,
+            p.precomputed.ᶜp,
+            p.core.∂ᶜK_∂ᶠu₃,
+            p.core.ᶜΦ,
+            p.core.ᶠgradᵥ_ᶜΦ,
+            p.core.ᶜρ_ref,
+            p.core.ᶜp_ref,
             p.scratch.ᶜtemp_scalar,
             p.params,
             p.atmos,
-            (energy_form isa TotalEnergy ? (; p.ᶜh_tot) : (;))...,
+            (energy_form isa TotalEnergy ? (; p.precomputed.ᶜh_tot) : (;))...,
             (
                 rayleigh_sponge isa RayleighSponge ?
                 (; p.rayleigh_sponge.ᶠβ_rayleigh_w) : (;)
diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl
index 4b61b505eea..323332a983b 100644
--- a/src/prognostic_equations/implicit/implicit_tendency.jl
+++ b/src/prognostic_equations/implicit/implicit_tendency.jl
@@ -62,7 +62,8 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx)
     (; dt) = p.simulation
     n = n_mass_flux_subdomains(turbconv_model)
     ᶜJ = Fields.local_geometry_field(Y.c).J
-    (; ᶜspecific, ᶠu³, ᶜp, ᶠgradᵥ_ᶜΦ, ᶜρ_ref, ᶜp_ref) = p
+    (; ᶠgradᵥ_ᶜΦ, ᶜρ_ref, ᶜp_ref) = p.core
+    (; ᶜspecific, ᶠu³, ᶜp) = p.precomputed
 
     ᶜ1 = p.scratch.ᶜtemp_scalar
     @. ᶜ1[colidx] = one(Y.c.ρ[colidx])
@@ -77,7 +78,7 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx)
     )
 
     if :ρe_tot in propertynames(Yₜ.c)
-        (; ᶜh_tot) = p
+        (; ᶜh_tot) = p.precomputed
         vertical_transport!(
             Yₜ.c.ρe_tot[colidx],
             ᶜJ[colidx],
diff --git a/src/prognostic_equations/remaining_tendency.jl b/src/prognostic_equations/remaining_tendency.jl
index 35ab1e19147..366071981de 100644
--- a/src/prognostic_equations/remaining_tendency.jl
+++ b/src/prognostic_equations/remaining_tendency.jl
@@ -51,8 +51,8 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
             colidx,
             p.atmos.turbconv_model,
         )
-        edmfx_nh_pressure_tendency!(Yₜ, Y, p, t, colidx, p.turbconv_model)
-        edmfx_tke_tendency!(Yₜ, Y, p, t, colidx, p.turbconv_model)
+        edmfx_nh_pressure_tendency!(Yₜ, Y, p, t, colidx, p.atmos.turbconv_model)
+        edmfx_tke_tendency!(Yₜ, Y, p, t, colidx, p.atmos.turbconv_model)
         precipitation_tendency!(Yₜ, Y, p, t, colidx, p.atmos.precip_model)
 
         # NOTE: All ρa tendencies should be applied before calling this function
diff --git a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl b/src/prognostic_equations/vertical_diffusion_boundary_layer.jl
index dc68cc32800..92cbee09ed2 100644
--- a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl
+++ b/src/prognostic_equations/vertical_diffusion_boundary_layer.jl
@@ -77,7 +77,7 @@ function vertical_diffusion_boundary_layer_tendency!(
 )
     ᶜρ = Y.c.ρ
     FT = Spaces.undertype(axes(ᶜρ))
-    (; ᶜp, sfc_conditions, ᶜspecific) = p # assume ᶜts and ᶜp have been updated
+    (; ᶜp, ᶜspecific, sfc_conditions) = p.precomputed # assume ᶜts and ᶜp have been updated
     (; C_E) = p.atmos.vert_diff
 
     ᶠgradᵥ = Operators.GradientC2F() # apply BCs to ᶜdivᵥ, which wraps ᶠgradᵥ
@@ -105,7 +105,7 @@ function vertical_diffusion_boundary_layer_tendency!(
     end
 
     if :ρe_tot in propertynames(Y.c)
-        (; ᶜh_tot) = p
+        (; ᶜh_tot) = p.precomputed
 
         ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
             top = Operators.SetValue(C3(FT(0))),
diff --git a/src/surface_conditions/surface_conditions.jl b/src/surface_conditions/surface_conditions.jl
index 9f2d04f66cf..7ff2723ce85 100644
--- a/src/surface_conditions/surface_conditions.jl
+++ b/src/surface_conditions/surface_conditions.jl
@@ -1,7 +1,7 @@
 """
     update_surface_conditions!(Y, p, t)
 
-Updates the value of `p.sfc_conditions` based on the current state `Y` and time
+Updates the value of `p.precomputed.sfc_conditions` based on the current state `Y` and time
 `t`. This function will only update the surface conditions if the surface_setup
 is not a PrescribedSurface.
 """
@@ -10,16 +10,13 @@ function update_surface_conditions!(Y, p, t)
     # Need to extract the field values so that we can do
     # a DataLayout broadcast rather than a Field broadcast
     # because we are mixing surface and interior fields
-    if isnothing(p.sfc_setup)
-        p.is_init[] && set_dummy_surface_conditions!(p)
-        return
-    end
     sfc_local_geometry_values = Fields.field_values(
         Fields.level(Fields.local_geometry_field(Y.f), Fields.half),
     )
     int_local_geometry_values =
         Fields.field_values(Fields.level(Fields.local_geometry_field(Y.c), 1))
-    (; ᶜts, ᶜu, sfc_conditions, params, sfc_setup, atmos) = p
+    (; ᶜts, ᶜu, sfc_conditions) = p.precomputed
+    (; params, sfc_setup, atmos) = p
     int_ts_values = Fields.field_values(Fields.level(ᶜts, 1))
     int_u_values = Fields.field_values(Fields.level(ᶜu, 1))
     int_z_values =
@@ -74,9 +71,10 @@ surface_state(
 # conditions, but without throwing an error during the computation of
 # precomputed quantities for diagnostic EDMF due to uninitialized surface
 # conditions.
-# TODO: Refactor the surface conditions API to avoid needing to do this. 
+# TODO: Refactor the surface conditions API to avoid needing to do this.
 function set_dummy_surface_conditions!(p)
-    (; sfc_conditions, params, atmos) = p
+    (; params, atmos) = p
+    (; sfc_conditions) = p.precomputed
     FT = eltype(params)
     thermo_params = CAP.thermodynamics_params(params)
     @. sfc_conditions.ustar = FT(0.2)
@@ -101,13 +99,14 @@ end
 """
     set_surface_conditions!(p, surface_conditions, surface_ts)
 
-Sets `p.sfc_conditions` according to `surface_conditions` and `surface_ts`,
+Sets `p.precomputed.sfc_conditions` according to `surface_conditions` and `surface_ts`,
 which are `Field`s of `SurfaceFluxes.SurfaceFluxConditions` and `Thermodynamics.ThermodynamicState`s
 This functions needs to be called by the coupler whenever either field changes
 to ensure that the simulation is properly updated.
 """
 function set_surface_conditions!(p, surface_conditions, surface_ts)
-    (; sfc_conditions, params, atmos) = p
+    (; params, atmos) = p
+    (; sfc_conditions,) = p.precomputed
     (; ᶠtemp_scalar) = p.scratch
 
     FT = eltype(params)
diff --git a/src/utils/debug_utils.jl b/src/utils/debug_utils.jl
index 44b0e74a804..0eb5812e2cb 100644
--- a/src/utils/debug_utils.jl
+++ b/src/utils/debug_utils.jl
@@ -109,7 +109,6 @@ feature with threading disabled.
 function precomputed_quantities!(Y, p, t, colidx)
     ᶜuₕ = Y.c.uₕ
     ᶠu₃ = Y.f.u₃
-    (; ᶜu_bar, ᶜK, ᶜts, ᶜp, params, thermo_dispatcher) = p
 
     @. ᶜu_bar[colidx] = C123(ᶜuₕ[colidx]) + C123(ᶜinterp(ᶠu₃[colidx]))
     @. ᶜK[colidx] = norm_sqr(ᶜu_bar[colidx]) / 2
diff --git a/src/utils/discrete_hydrostatic_balance.jl b/src/utils/discrete_hydrostatic_balance.jl
index 57a45f47d1d..87064718043 100644
--- a/src/utils/discrete_hydrostatic_balance.jl
+++ b/src/utils/discrete_hydrostatic_balance.jl
@@ -6,7 +6,7 @@ import .InitialConditions as ICs
 
 """
     set_discrete_hydrostatic_balanced_state!(Y, p)
-Modify the energy variable in state `Y` given Y and the cache `p` so that 
+Modify the energy variable in state `Y` given Y and the cache `p` so that
 `Y` is in discrete hydrostatic balance.
 """
 function set_discrete_hydrostatic_balanced_state!(Y, p)
@@ -14,20 +14,25 @@ function set_discrete_hydrostatic_balanced_state!(Y, p)
     ᶠgradᵥ_ᶜp = similar(Y.f.u₃)
     Fields.bycolumn(axes(Y.c.ρ)) do colidx
         set_discrete_hydrostatic_balanced_pressure!(
-            p.ᶜp,
+            p.precomputed.ᶜp,
             ᶠgradᵥ_ᶜp,
             Y.c.ρ,
-            p.ᶠgradᵥ_ᶜΦ,
+            p.core.ᶠgradᵥ_ᶜΦ,
             FT(CAP.MSLP(p.params)),
             colidx,
         )
     end
     thermo_params = CAP.thermodynamics_params(p.params)
     if p.atmos.moisture_model isa DryModel
-        @. p.ᶜts = TD.PhaseDry_ρp(thermo_params, Y.c.ρ, p.ᶜp)
+        @. p.precomputed.ᶜts =
+            TD.PhaseDry_ρp(thermo_params, Y.c.ρ, p.precomputed.ᶜp)
     elseif p.atmos.moisture_model isa EquilMoistModel
-        @. p.ᶜts =
-            TD.PhaseEquil_ρpq(thermo_params, Y.c.ρ, p.ᶜp, Y.c.ρq_tot / Y.c.ρ)
+        @. p.precomputed.ᶜts = TD.PhaseEquil_ρpq(
+            thermo_params,
+            Y.c.ρ,
+            p.precomputed.ᶜp,
+            Y.c.ρq_tot / Y.c.ρ,
+        )
     else
         error("Unsupported moisture model")
     end
@@ -39,7 +44,7 @@ function set_discrete_hydrostatic_balanced_state!(Y, p)
         ICs.energy_variables(
             ls(
                 p.params,
-                p.ᶜts,
+                p.precomputed.ᶜts,
                 ᶜlocal_geometry,
                 Geometry.UVWVector(Y.c.uₕ) +
                 Geometry.UVWVector(ᶜinterp(Y.f.u₃)),
@@ -51,7 +56,7 @@ end
 
 """
     set_discrete_hydrostatic_balanced_pressure!(ᶜp, ᶠgradᵥ_ᶜp, ᶜρ, ᶠgradᵥ_ᶜΦ, p1, colidx)
-Construct discrete hydrostatic balanced pressure `ᶜp` from density `ᶜρ`, 
+Construct discrete hydrostatic balanced pressure `ᶜp` from density `ᶜρ`,
 potential energy gradient `ᶠgradᵥ_ᶜΦ`, and surface pressure `p1`.
 
 Yₜ.f.u₃ = 0 ==>
diff --git a/test/coupler_compatibility.jl b/test/coupler_compatibility.jl
index 6737105fad3..56836196f02 100644
--- a/test/coupler_compatibility.jl
+++ b/test/coupler_compatibility.jl
@@ -61,11 +61,13 @@ const T2 = 290
     # temperature to T1 and then to T2.
     @. sfc_setup.T = FT(T1)
     CA.set_precomputed_quantities!(Y, p_overwritten, t)
-    sfc_T = @. TD.air_temperature(thermo_params, p.sfc_conditions.ts)
+    sfc_T =
+        @. TD.air_temperature(thermo_params, p.precomputed.sfc_conditions.ts)
     @test all(isequal(T1), parent(sfc_T))
     @. sfc_setup.T = FT(T2)
     CA.set_precomputed_quantities!(Y, p_overwritten, t)
-    sfc_T = @. TD.air_temperature(thermo_params, p.sfc_conditions.ts)
+    sfc_T =
+        @. TD.air_temperature(thermo_params, p.precomputed.sfc_conditions.ts)
     @test all(isequal(T2), parent(sfc_T))
 end
 
@@ -92,7 +94,8 @@ end
     # Define a function for updating these fields, given the temperature at the
     # surface and the current cache p.
     function update_surface_fields!(sfc_ts, sfc_conditions, surface_T, p)
-        (; ᶜts, ᶜu, params) = p
+        (; params) = p
+        (; ᶜts, ᶜu) = p.precomputed
         thermo_params = CAP.thermodynamics_params(params)
         surface_params = CAP.surface_fluxes_params(params)
 
@@ -158,11 +161,13 @@ end
     # temperature to T1 and then to T2.
     update_surface_fields!(sfc_ts, sfc_conditions, FT(T1), p)
     CA.SurfaceConditions.set_surface_conditions!(p, sfc_conditions, sfc_ts)
-    sfc_T = @. TD.air_temperature(thermo_params, p.sfc_conditions.ts)
+    sfc_T =
+        @. TD.air_temperature(thermo_params, p.precomputed.sfc_conditions.ts)
     @test all(isequal(T1), parent(sfc_T))
     update_surface_fields!(sfc_ts, sfc_conditions, FT(T2), p)
     CA.SurfaceConditions.set_surface_conditions!(p, sfc_conditions, sfc_ts)
-    sfc_T = @. TD.air_temperature(thermo_params, p.sfc_conditions.ts)
+    sfc_T =
+        @. TD.air_temperature(thermo_params, p.precomputed.sfc_conditions.ts)
     @test all(isequal(T2), parent(sfc_T))
 end