diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl
index 416d29f62b..81e40db889 100644
--- a/src/ClimaAtmos.jl
+++ b/src/ClimaAtmos.jl
@@ -21,6 +21,7 @@ include(joinpath("utils", "read_gcm_driven_scm_data.jl"))
 include(joinpath("utils", "AtmosArtifacts.jl"))
 import .AtmosArtifacts as AA
 
+include(joinpath("core", "core_quantities.jl"))
 include(
     joinpath("parameterized_tendencies", "radiation", "radiation_utilities.jl"),
 )
diff --git a/src/cache/cache.jl b/src/cache/cache.jl
index abf60dd48d..94c5973111 100644
--- a/src/cache/cache.jl
+++ b/src/cache/cache.jl
@@ -33,7 +33,7 @@ struct AtmosCache{
     """ClimaAtmosParameters that have to be used"""
     params::CAP
 
-    """Variables that are used generally, such as ᶜΦ"""
+    """Variables that are used generally"""
     core::COR
 
     """Used by update_surface_conditions! in set_precomputed_quantities! and coupler"""
@@ -89,10 +89,8 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names)
     ᶜcoord = Fields.local_geometry_field(Y.c).coordinates
     ᶠcoord = Fields.local_geometry_field(Y.f).coordinates
     grav = FT(CAP.grav(params))
-    ᶜΦ = grav .* ᶜcoord.z
-    ᶠΦ = grav .* ᶠcoord.z
-
-    (; ᶜf³, ᶠf¹²) = compute_coriolis(ᶜcoord, ᶠcoord, params)
+    ᶜz = ᶜcoord.z
+    ᶠz = ᶠcoord.z
 
     ghost_buffer =
         !do_dss(axes(Y.c)) ? (;) :
@@ -121,11 +119,8 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names)
         Fields.level(Fields.local_geometry_field(Y.f), Fields.half)
 
     core = (
-        ᶜΦ,
-        ᶠgradᵥ_ᶜΦ = ᶠgradᵥ.(ᶜΦ),
-        ᶜgradᵥ_ᶠΦ = ᶜgradᵥ.(ᶠΦ),
-        ᶜf³,
-        ᶠf¹²,
+        ᶠgradᵥ_ᶜΦ = ᶠgradᵥ.(Φ.(grav, ᶜz)),
+        ᶜgradᵥ_ᶠΦ = ᶜgradᵥ.(Φ.(grav, ᶠz)),
         # Used by diagnostics such as hfres, evspblw
         surface_ct3_unit = CT3.(
             unit_basis_vector_data.(CT3, sfc_local_geometry)
@@ -185,32 +180,3 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names)
 
     return AtmosCache{map(typeof, args)...}(args...)
 end
-
-
-function compute_coriolis(ᶜcoord, ᶠcoord, params)
-    if eltype(ᶜcoord) <: Geometry.LatLongZPoint
-        Ω = CAP.Omega(params)
-        global_geom = Spaces.global_geometry(axes(ᶜcoord))
-        if global_geom isa Geometry.DeepSphericalGlobalGeometry
-            @info "using deep atmosphere"
-            coriolis_deep(coord::Geometry.LatLongZPoint) = Geometry.LocalVector(
-                Geometry.Cartesian123Vector(zero(Ω), zero(Ω), 2 * Ω),
-                global_geom,
-                coord,
-            )
-            ᶜf³ = @. CT3(CT123(coriolis_deep(ᶜcoord)))
-            ᶠf¹² = @. CT12(CT123(coriolis_deep(ᶠcoord)))
-        else
-            coriolis_shallow(coord::Geometry.LatLongZPoint) =
-                Geometry.WVector(2 * Ω * sind(coord.lat))
-            ᶜf³ = @. CT3(coriolis_shallow(ᶜcoord))
-            ᶠf¹² = nothing
-        end
-    else
-        f = CAP.f_plane_coriolis_frequency(params)
-        coriolis_f_plane(coord) = Geometry.WVector(f)
-        ᶜf³ = @. CT3(coriolis_f_plane(ᶜcoord))
-        ᶠf¹² = nothing
-    end
-    return (; ᶜf³, ᶠf¹²)
-end
diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl
index 59d88ba0c3..41e554cc06 100644
--- a/src/cache/diagnostic_edmf_precomputed_quantities.jl
+++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl
@@ -34,13 +34,14 @@ NVTX.@annotate function set_diagnostic_edmfx_draft_quantities_level!(
     mse_level,
     q_tot_level,
     p_level,
-    Φ_level,
+    z_level,
 )
     FT = eltype(thermo_params)
+    grav = TDP.grav(thermo_params)
     @. ts_level = TD.PhaseEquil_phq(
         thermo_params,
         p_level,
-        mse_level - Φ_level,
+        mse_level - Φ(grav, z_level),
         q_tot_level,
         8,
         FT(0.0003),
@@ -92,7 +93,6 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!(
     (; turbconv_model) = p.atmos
     FT = eltype(Y)
     n = n_mass_flux_subdomains(turbconv_model)
-    (; ᶜΦ) = 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) =
@@ -112,7 +112,6 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!(
     q_tot_int_level = Fields.field_values(Fields.level(q_tot, 1))
 
     p_int_level = Fields.field_values(Fields.level(ᶜp, 1))
-    Φ_int_level = Fields.field_values(Fields.level(ᶜΦ, 1))
 
     local_geometry_int_level =
         Fields.field_values(Fields.level(Fields.local_geometry_field(Y.c), 1))
@@ -187,7 +186,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!(
             mseʲ_int_level,
             q_totʲ_int_level,
             p_int_level,
-            Φ_int_level,
+            z_int_level,
         )
         @. ρaʲ_int_level = ρʲ_int_level * FT(turbconv_params.surface_area)
     end
@@ -296,7 +295,6 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
     ᶜdz = Fields.Δz_field(axes(Y.c))
     (; params) = p
     (; dt) = p
-    (; ᶜΦ) = p.core
     (; ᶜp, ᶠu³, ᶜts, ᶜh_tot, ᶜK) = p.precomputed
     (; q_tot) = p.precomputed.ᶜspecific
     (;
@@ -324,13 +322,17 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
     microphys_1m_params = CAP.microphysics_1m_params(params)
     turbconv_params = CAP.turbconv_params(params)
 
-    ᶠΦ = p.scratch.ᶠtemp_scalar
-    @. ᶠΦ = CAP.grav(params) * ᶠz
+    g = CAP.grav(params)
+    ᶜcoords = Fields.coordinate_field(Y.c)
+    ᶜz = Fields.coordinate_field(Y.c).z
+    ᶠz = Fields.coordinate_field(Y.f).z
+    global_geom = Spaces.global_geometry(axes(ᶜcoords))
+
     ᶜ∇Φ³ = p.scratch.ᶜtemp_CT3
-    @. ᶜ∇Φ³ = CT3(ᶜgradᵥ(ᶠΦ))
-    @. ᶜ∇Φ³ += CT3(gradₕ(ᶜΦ))
+    @. ᶜ∇Φ³ = CT3(ᶜgradᵥ(Φ(g, ᶠz)))
+    @. ᶜ∇Φ³ += CT3(gradₕ(Φ(g, ᶜz)))
     ᶜ∇Φ₃ = p.scratch.ᶜtemp_C3
-    @. ᶜ∇Φ₃ = ᶜgradᵥ(ᶠΦ)
+    @. ᶜ∇Φ₃ = ᶜgradᵥ(Φ(g, ᶠz))
 
     z_sfc_halflevel =
         Fields.field_values(Fields.level(Fields.coordinate_field(Y.f).z, half))
@@ -344,7 +346,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
         h_tot_level = Fields.field_values(Fields.level(ᶜh_tot, i))
         q_tot_level = Fields.field_values(Fields.level(q_tot, i))
         p_level = Fields.field_values(Fields.level(ᶜp, i))
-        Φ_level = Fields.field_values(Fields.level(ᶜΦ, i))
+        z_level = Fields.field_values(Fields.level(ᶜz, i))
         local_geometry_level = Fields.field_values(
             Fields.level(Fields.local_geometry_field(Y.c), i),
         )
@@ -355,7 +357,6 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
         end_index = fieldcount(eltype(∂x∂ξ_level)) # This will be 4 in 2D and 9 in 3D.
         ∂x³∂ξ³_level = ∂x∂ξ_level.:($end_index)
 
-        Φ_prev_level = Fields.field_values(Fields.level(ᶜΦ, i - 1))
         ∇Φ³_prev_level = Fields.field_values(Fields.level(ᶜ∇Φ³, i - 1))
         ∇Φ³_data_prev_level = ∇Φ³_prev_level.components.data.:1
         ∇Φ₃_prev_level = Fields.field_values(Fields.level(ᶜ∇Φ₃, i - 1))
@@ -555,7 +556,6 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                     q_rai_prev_level,
                     q_sno_prev_level,
                     tsʲ_prev_level,
-                    Φ_prev_level,
                     dt,
                     microphys_1m_params,
                     thermo_params,
@@ -702,7 +702,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                             e_tot_0M_precipitation_sources_helper(
                                 thermo_params,
                                 tsʲ_prev_level,
-                                Φ_prev_level,
+                                Φ(g, z_prev_level),
                             )
                         )
                     )
@@ -807,7 +807,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                 mseʲ_level,
                 q_totʲ_level,
                 p_level,
-                Φ_level,
+                z_level,
             )
         end
         ρaʲs_level = Fields.field_values(Fields.level(ᶜρaʲs, i))
@@ -1072,7 +1072,6 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita
         ᶜqᵣ,
         ᶜqₛ,
         ᶜts,
-        p.core.ᶜΦ,
         p.dt,
         microphys_1m_params,
         thermo_params,
diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
index f7fa5668e0..9e4a0fc4ec 100644
--- a/src/cache/precomputed_quantities.jl
+++ b/src/cache/precomputed_quantities.jl
@@ -465,9 +465,10 @@ 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, moisture_model)
-    (; ᶜΦ) = p.core
     (; ᶜspecific, ᶜu, ᶠu³, ᶜK, ᶜts, ᶜp) = p.precomputed
     ᶠuₕ³ = p.scratch.ᶠtemp_CT3
+    ᶜz = Fields.coordinate_field(axes(Y.c)).z
+    grav = TDP.grav(thermo_params)
 
     @. ᶜspecific = specific_gs(Y.c)
     set_ᶠuₕ³!(ᶠuₕ³, Y)
@@ -493,7 +494,7 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
         # @. ᶜK += Y.c.sgs⁰.ρatke / Y.c.ρ
         # TODO: We should think more about these increments before we use them.
     end
-    @. ᶜts = ts_gs(thermo_args..., ᶜspecific, ᶜK, ᶜΦ, Y.c.ρ)
+    @. ᶜts = ts_gs(thermo_args..., ᶜspecific, ᶜK, Φ(grav, ᶜz), Y.c.ρ)
     @. ᶜp = TD.air_pressure(thermo_params, ᶜts)
 
     if turbconv_model isa AbstractEDMF
diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl
index cf6b889f41..a93eec06a3 100644
--- a/src/cache/prognostic_edmf_precomputed_quantities.jl
+++ b/src/cache/prognostic_edmf_precomputed_quantities.jl
@@ -20,7 +20,8 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_environment!(
 
     thermo_params = CAP.thermodynamics_params(p.params)
     (; turbconv_model) = p.atmos
-    (; ᶜΦ,) = p.core
+    ᶜz = Fields.coordinate_field(axes(Y.c)).z
+    grav = TDP.grav(thermo_params)
     (; ᶜp, ᶜh_tot, ᶜK) = p.precomputed
     (; ᶜtke⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜmse⁰, ᶜq_tot⁰) =
         p.precomputed
@@ -44,7 +45,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_environment!(
     set_sgs_ᶠu₃!(u₃⁰, ᶠu₃⁰, Y, turbconv_model)
     set_velocity_quantities!(ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶠu₃⁰, Y.c.uₕ, ᶠuₕ³)
     # @. ᶜK⁰ += ᶜtke⁰
-    @. ᶜts⁰ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmse⁰ - ᶜΦ, ᶜq_tot⁰)
+    @. ᶜts⁰ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmse⁰ - Φ(grav, ᶜz), ᶜq_tot⁰)
     @. ᶜρ⁰ = TD.air_density(thermo_params, ᶜts⁰)
     return nothing
 end
@@ -72,7 +73,8 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft_and_bc!
     thermo_params = CAP.thermodynamics_params(params)
     turbconv_params = CAP.turbconv_params(params)
 
-    (; ᶜΦ,) = p.core
+    grav = TDP.grav(thermo_params)
+    ᶜz = Fields.coordinate_field(Y.c).z
     (; ᶜspecific, ᶜp, ᶜh_tot, ᶜK) = p.precomputed
     (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶠKᵥʲs, ᶜtsʲs, ᶜρʲs) = p.precomputed
     (; ustar, obukhov_length, buoyancy_flux) = p.precomputed.sfc_conditions
@@ -90,7 +92,8 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft_and_bc!
 
         set_velocity_quantities!(ᶜuʲ, ᶠu³ʲ, ᶜKʲ, ᶠu₃ʲ, Y.c.uₕ, ᶠuₕ³)
         @. ᶠKᵥʲ = (adjoint(CT3(ᶠu₃ʲ)) * ᶠu₃ʲ) / 2
-        @. ᶜtsʲ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmseʲ - ᶜΦ, ᶜq_totʲ)
+        @. ᶜtsʲ =
+            TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmseʲ - Φ(grav, ᶜz), ᶜq_totʲ)
         @. ᶜρʲ = TD.air_density(thermo_params, ᶜtsʲ)
 
         # EDMFX boundary condition:
@@ -153,12 +156,11 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft_and_bc!
         )
 
         # Then overwrite the prognostic variables at first inetrior point.
-        ᶜΦ_int_val = Fields.field_values(Fields.level(ᶜΦ, 1))
         ᶜtsʲ_int_val = Fields.field_values(Fields.level(ᶜtsʲ, 1))
         @. ᶜtsʲ_int_val = TD.PhaseEquil_phq(
             thermo_params,
             ᶜp_int_val,
-            ᶜmseʲ_int_val - ᶜΦ_int_val,
+            ᶜmseʲ_int_val - Φ(grav, ᶜz_int_val),
             ᶜq_totʲ_int_val,
         )
         sgsʲs_ρ_int_val = Fields.field_values(Fields.level(ᶜρʲs.:($j), 1))
@@ -424,7 +426,6 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation
     @assert !(p.atmos.moisture_model isa DryModel)
 
     (; params, dt) = p
-    (; ᶜΦ,) = p.core
     thp = CAP.thermodynamics_params(params)
     cmp = CAP.microphysics_1m_params(params)
 
@@ -451,7 +452,6 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation
             ᶜqᵣ,
             ᶜqₛ,
             ᶜtsʲs.:($j),
-            ᶜΦ,
             dt,
             cmp,
             thp,
@@ -470,7 +470,6 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation
         ᶜqᵣ,
         ᶜqₛ,
         ᶜts⁰,
-        ᶜΦ,
         dt,
         cmp,
         thp,
diff --git a/src/core/core_quantities.jl b/src/core/core_quantities.jl
new file mode 100644
index 0000000000..e09b55a21b
--- /dev/null
+++ b/src/core/core_quantities.jl
@@ -0,0 +1,64 @@
+import ClimaCore.Geometry
+
+f³(lg, global_geom, params) = f³(lg, lg.coordinates, global_geom, params)
+
+f¹²(lg, global_geom, params) = f¹²(lg, lg.coordinates, global_geom, params)
+
+function f³(
+    lg,
+    coord::Geometry.LatLongZPoint,
+    global_geom::Geometry.DeepSphericalGlobalGeometry,
+    params,
+)
+    Ω = CAP.Omega(params)
+    CT3(
+        CT123(
+            Geometry.LocalVector(
+                Geometry.Cartesian123Vector(zero(Ω), zero(Ω), 2 * Ω),
+                global_geom,
+                coord,
+            ),
+            lg,
+        ),
+        lg,
+    )
+end
+
+function f¹²(
+    lg,
+    coord::Geometry.LatLongZPoint,
+    global_geom::Geometry.DeepSphericalGlobalGeometry,
+    params,
+)
+    Ω = CAP.Omega(params)
+    CT12(
+        CT123(
+            Geometry.LocalVector(
+                Geometry.Cartesian123Vector(zero(Ω), zero(Ω), 2 * Ω),
+                global_geom,
+                coord,
+            ),
+            lg,
+        ),
+        lg,
+    )
+end
+
+# Shallow sphere
+function f³(lg, coord::Geometry.LatLongZPoint, global_geom, params)
+    Ω = CAP.Omega(params)
+    CT3(Geometry.WVector(2 * Ω * sind(coord.lat), lg), lg)
+end
+
+# Shallow cartesian
+function f³(lg, coord, global_geom, params)
+    f = CAP.f_plane_coriolis_frequency(params)
+    CT3(Geometry.WVector(f, lg), lg)
+end
+
+f¹²(lg, coord::Geometry.LatLongZPoint, global_geom, params) =
+    error("Not supported for $coord, $global_geom.")
+f¹²(lg, coord, global_geom, p) =
+    error("Not supported for $coord, $global_geom.")
+
+Φ(g, z) = g * z
diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl
index cec21cb153..c4c035f394 100644
--- a/src/diagnostics/Diagnostics.jl
+++ b/src/diagnostics/Diagnostics.jl
@@ -12,6 +12,7 @@ import ClimaCore.Utilities: half
 import Thermodynamics as TD
 
 import ..AtmosModel
+import ..Φ
 import ..AtmosCallback
 import ..EveryNSteps
 
diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl
index 4ab948a469..82d9a73589 100644
--- a/src/diagnostics/core_diagnostics.jl
+++ b/src/diagnostics/core_diagnostics.jl
@@ -220,9 +220,16 @@ add_diagnostic_variable!(
     units = "m",
     compute! = (out, state, cache, time) -> begin
         if isnothing(out)
-            return cache.core.ᶜΦ ./ CAP.grav(cache.params)
+            return Φ.(
+                CAP.grav(cache.params),
+                Fields.coordinate_field(axes(state.c.ρ)).z,
+            ) ./ CAP.grav(cache.params)
         else
-            out .= cache.core.ᶜΦ ./ CAP.grav(cache.params)
+            out .=
+                Φ.(
+                    CAP.grav(cache.params),
+                    Fields.coordinate_field(axes(state.c.ρ)).z,
+                ) ./ CAP.grav(cache.params)
         end
     end,
 )
@@ -1025,6 +1032,8 @@ add_diagnostic_variable!(
 ###
 function compute_dsevi!(out, state, cache, time)
     thermo_params = CAP.thermodynamics_params(cache.params)
+    z = Fields.coordinate_field(axes(state.c.ρ)).z
+    grav = CAP.grav(cache.params)
     if isnothing(out)
         out = zeros(axes(Fields.level(state.f, half)))
         cp = CAP.cp_d(cache.params)
@@ -1032,7 +1041,7 @@ function compute_dsevi!(out, state, cache, time)
         @. dse =
             state.c.ρ * (
                 cp * TD.air_temperature(thermo_params, cache.precomputed.ᶜts) +
-                cache.core.ᶜΦ
+                Φ(grav, z)
             )
         Operators.column_integral_definite!(out, dse)
         return out
@@ -1042,7 +1051,7 @@ function compute_dsevi!(out, state, cache, time)
         @. dse =
             state.c.ρ * (
                 cp * TD.air_temperature(thermo_params, cache.precomputed.ᶜts) +
-                cache.core.ᶜΦ
+                Φ(grav, z)
             )
         Operators.column_integral_definite!(out, dse)
     end
diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl
index 341a3a26bb..af21d0faca 100644
--- a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl
+++ b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl
@@ -92,7 +92,7 @@ function q_tot_precipitation_sources(
 end
 
 """
-    e_tot_0M_precipitation_sources_helper(thp, ts, Φ)
+    e_tot_0M_precipitation_sources_helper(thp, ts, g, z)
 
  - thp - set with thermodynamics parameters
  - ts - thermodynamic state (see td package for details)
@@ -109,14 +109,13 @@ function e_tot_0M_precipitation_sources_helper(thp, ts, Φ)
 end
 
 """
-    compute_precipitation_sources!(Sᵖ, Sᵖ_snow, Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ, ρ, qᵣ, qₛ, ts, Φ, dt, mp, thp)
+    compute_precipitation_sources!(Sᵖ, Sᵖ_snow, Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ, ρ, qᵣ, qₛ, ts, dt, mp, thp)
 
  - Sᵖ, Sᵖ_snow - temporary containters to help compute precipitation source terms
  - Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ - cached storage for precipitation source terms
  - ρ - air density
  - qᵣ, qₛ - precipitation (rain and snow) specific humidity
  - ts - thermodynamic state (see td package for details)
- - Φ - geopotential
  - dt - model time step
  - thp, cmp - structs with thermodynamic and microphysics parameters
 
@@ -136,11 +135,12 @@ function compute_precipitation_sources!(
     qᵣ,
     qₛ,
     ts,
-    Φ,
     dt,
     mp,
     thp,
 )
+    g = TDP.grav(thp)
+    z = Fields.coordinate_field(axes(Sᵖ)).z
     FT = eltype(thp)
     # @. Sqₜᵖ = FT(0) should work after fixing
     # https://github.com/CliMA/ClimaCore.jl/issues/1786
@@ -159,7 +159,7 @@ function compute_precipitation_sources!(
     @. Sᵖ = min(limit(qₗ(thp, ts), dt, 5), Sᵖ)
     @. Sqₜᵖ -= Sᵖ
     @. Sqᵣᵖ += Sᵖ
-    @. Seₜᵖ -= Sᵖ * (Iₗ(thp, ts) + Φ)
+    @. Seₜᵖ -= Sᵖ * (Iₗ(thp, ts) + Φ(g, z))
 
     # snow autoconversion assuming no supersaturation: q_ice -> q_snow
     @. Sᵖ = min(
@@ -168,7 +168,7 @@ function compute_precipitation_sources!(
     )
     @. Sqₜᵖ -= Sᵖ
     @. Sqₛᵖ += Sᵖ
-    @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ)
+    @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ(g, z))
 
     # accretion: q_liq + q_rain -> q_rain
     @. Sᵖ = min(
@@ -177,7 +177,7 @@ function compute_precipitation_sources!(
     )
     @. Sqₜᵖ -= Sᵖ
     @. Sqᵣᵖ += Sᵖ
-    @. Seₜᵖ -= Sᵖ * (Iₗ(thp, ts) + Φ)
+    @. Seₜᵖ -= Sᵖ * (Iₗ(thp, ts) + Φ(g, z))
 
     # accretion: q_ice + q_snow -> q_snow
     @. Sᵖ = min(
@@ -186,7 +186,7 @@ function compute_precipitation_sources!(
     )
     @. Sqₜᵖ -= Sᵖ
     @. Sqₛᵖ += Sᵖ
-    @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ)
+    @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ(g, z))
 
     # accretion: q_liq + q_sno -> q_sno or q_rai
     # sink of cloud water via accretion cloud water + snow
@@ -207,8 +207,8 @@ function compute_precipitation_sources!(
     @. Sqᵣᵖ += ifelse(Tₐ(thp, ts) < mp.ps.T_freeze, FT(0), Sᵖ - Sᵖ_snow)
     @. Seₜᵖ -= ifelse(
         Tₐ(thp, ts) < mp.ps.T_freeze,
-        Sᵖ * (Iᵢ(thp, ts) + Φ),
-        Sᵖ * (Iₗ(thp, ts) + Φ) - Sᵖ_snow * (Iₗ(thp, ts) - Iᵢ(thp, ts)),
+        Sᵖ * (Iᵢ(thp, ts) + Φ(g, z)),
+        Sᵖ * (Iₗ(thp, ts) + Φ(g, z)) - Sᵖ_snow * (Iₗ(thp, ts) - Iᵢ(thp, ts)),
     )
 
     # accretion: q_ice + q_rai -> q_sno
@@ -218,7 +218,7 @@ function compute_precipitation_sources!(
     )
     @. Sqₜᵖ -= Sᵖ
     @. Sqₛᵖ += Sᵖ
-    @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ)
+    @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ(g, z))
     # sink of rain via accretion cloud ice - rain
     @. Sᵖ = min(
         limit(qᵣ, dt, 5),
@@ -285,14 +285,13 @@ function compute_precipitation_heating!(
     @. ᶜSeₜᵖ -= dot(ᶜ∇T, (ᶜu - C123(Geometry.WVector(ᶜwₛ)))) * cᵥᵢ(thp) * ᶜqₛ
 end
 """
-    compute_precipitation_sinks!(Sᵖ, Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ, ρ, qᵣ, qₛ, ts, Φ, dt, mp, thp)
+    compute_precipitation_sinks!(Sᵖ, Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ, ρ, qᵣ, qₛ, ts, dt, mp, thp)
 
  - Sᵖ - a temporary containter to help compute precipitation source terms
  - Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ - cached storage for precipitation source terms
  - ρ - air density
  - qᵣ, qₛ - precipitation (rain and snow) specific humidities
  - ts - thermodynamic state (see td package for details)
- - Φ - geopotential
  - dt - model time step
  - thp, cmp - structs with thermodynamic and microphysics parameters
 
@@ -311,7 +310,6 @@ function compute_precipitation_sinks!(
     qᵣ,
     qₛ,
     ts,
-    Φ,
     dt,
     mp,
     thp,
@@ -319,6 +317,8 @@ function compute_precipitation_sinks!(
     FT = eltype(Sqₜᵖ)
     sps = (mp.ps, mp.tv.snow, mp.aps, thp)
     rps = (mp.pr, mp.tv.rain, mp.aps, thp)
+    g = TDP.grav(thp)
+    z = Fields.coordinate_field(axes(ρ)).z
 
     #! format: off
     # evaporation: q_rai -> q_vap
@@ -328,7 +328,7 @@ function compute_precipitation_sinks!(
     )
     @. Sqₜᵖ -= Sᵖ
     @. Sqᵣᵖ += Sᵖ
-    @. Seₜᵖ -= Sᵖ * (Iₗ(thp, ts) + Φ)
+    @. Seₜᵖ -= Sᵖ * (Iₗ(thp, ts) + Φ(g, z))
 
     # melting: q_sno -> q_rai
     @. Sᵖ = min(
@@ -348,6 +348,6 @@ function compute_precipitation_sinks!(
     )
     @. Sqₜᵖ -= Sᵖ
     @. Sqₛᵖ += Sᵖ
-    @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ)
+    @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ(g, z))
     #! format: on
 end
diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl
index 86606a5253..681b8245c8 100644
--- a/src/parameterized_tendencies/microphysics/precipitation.jl
+++ b/src/parameterized_tendencies/microphysics/precipitation.jl
@@ -46,7 +46,6 @@ function compute_precipitation_cache!(Y, p, ::Microphysics0Moment, _)
     (; params, dt) = p
     (; ᶜts) = p.precomputed
     (; ᶜS_ρq_tot, ᶜS_ρe_tot) = p.precipitation
-    (; ᶜΦ) = p.core
     cm_params = CAP.microphysics_0m_params(params)
     thermo_params = CAP.thermodynamics_params(params)
     @. ᶜS_ρq_tot =
@@ -58,9 +57,11 @@ function compute_precipitation_cache!(Y, p, ::Microphysics0Moment, _)
             Y.c.ρq_tot / Y.c.ρ,
             ᶜts,
         )
+    grav = TDP.grav(thermo_params)
+    ᶜz = Fields.coordinate_field(axes(Y.c)).z
     @. ᶜS_ρe_tot =
         ᶜS_ρq_tot *
-        e_tot_0M_precipitation_sources_helper(thermo_params, ᶜts, ᶜΦ)
+        e_tot_0M_precipitation_sources_helper(thermo_params, ᶜts, Φ(grav, ᶜz))
 end
 function compute_precipitation_cache!(
     Y,
@@ -70,11 +71,12 @@ function compute_precipitation_cache!(
 )
     # For environment we multiply by grid mean ρ and not byᶜρa⁰
     # assuming a⁰=1
-    (; ᶜΦ) = p.core
     (; ᶜSqₜᵖ⁰, ᶜSqₜᵖʲs, ᶜρaʲs) = p.precomputed
     (; ᶜS_ρq_tot, ᶜS_ρe_tot) = p.precipitation
     (; ᶜts, ᶜtsʲs) = p.precomputed
     thermo_params = CAP.thermodynamics_params(p.params)
+    grav = TDP.grav(thermo_params)
+    ᶜz = Fields.coordinate_field(axes(Y.c)).z
 
     n = n_mass_flux_subdomains(p.atmos.turbconv_model)
     ρ = Y.c.ρ
@@ -83,7 +85,7 @@ function compute_precipitation_cache!(
     @. ᶜS_ρe_tot =
         ᶜSqₜᵖ⁰ *
         ρ *
-        e_tot_0M_precipitation_sources_helper(thermo_params, ᶜts, ᶜΦ)
+        e_tot_0M_precipitation_sources_helper(thermo_params, ᶜts, Φ(grav, ᶜz))
     for j in 1:n
         @. ᶜS_ρq_tot += ᶜSqₜᵖʲs.:($$j) * ᶜρaʲs.:($$j)
         @. ᶜS_ρe_tot +=
@@ -92,7 +94,7 @@ function compute_precipitation_cache!(
             e_tot_0M_precipitation_sources_helper(
                 thermo_params,
                 ᶜtsʲs.:($$j),
-                ᶜΦ,
+                Φ(grav, ᶜz),
             )
     end
 end
@@ -102,19 +104,20 @@ function compute_precipitation_cache!(
     ::Microphysics0Moment,
     ::PrognosticEDMFX,
 )
-    (; ᶜΦ) = p.core
     (; ᶜSqₜᵖ⁰, ᶜSqₜᵖʲs, ᶜρa⁰) = p.precomputed
     (; ᶜS_ρq_tot, ᶜS_ρe_tot) = p.precipitation
     (; ᶜts⁰, ᶜtsʲs) = p.precomputed
     thermo_params = CAP.thermodynamics_params(p.params)
 
     n = n_mass_flux_subdomains(p.atmos.turbconv_model)
+    grav = TDP.grav(thermo_params)
+    ᶜz = Fields.coordinate_field(axes(Y.c)).z
 
     @. ᶜS_ρq_tot = ᶜSqₜᵖ⁰ * ᶜρa⁰
     @. ᶜS_ρe_tot =
         ᶜSqₜᵖ⁰ *
         ᶜρa⁰ *
-        e_tot_0M_precipitation_sources_helper(thermo_params, ᶜts⁰, ᶜΦ)
+        e_tot_0M_precipitation_sources_helper(thermo_params, ᶜts⁰, Φ(grav, ᶜz))
     for j in 1:n
         @. ᶜS_ρq_tot += ᶜSqₜᵖʲs.:($$j) * Y.c.sgsʲs.:($$j).ρa
         @. ᶜS_ρe_tot +=
@@ -123,7 +126,7 @@ function compute_precipitation_cache!(
             e_tot_0M_precipitation_sources_helper(
                 thermo_params,
                 ᶜtsʲs.:($$j),
-                ᶜΦ,
+                Φ(grav, ᶜz),
             )
     end
 end
@@ -197,8 +200,10 @@ function compute_precipitation_cache!(Y, p, ::Microphysics1Moment, _)
     FT = Spaces.undertype(axes(Y.c))
     (; dt) = p
     (; ᶜts, ᶜqᵣ, ᶜqₛ, ᶜwᵣ, ᶜwₛ, ᶜu) = p.precomputed
-    (; ᶜΦ) = p.core
     (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation
+    thermo_params = CAP.thermodynamics_params(p.params)
+    grav = TDP.grav(thermo_params)
+    ᶜz = Fields.coordinate_field(axes(Y.c)).z
 
     ᶜSᵖ = p.scratch.ᶜtemp_scalar
     ᶜSᵖ_snow = p.scratch.ᶜtemp_scalar_2
@@ -221,7 +226,6 @@ function compute_precipitation_cache!(Y, p, ::Microphysics1Moment, _)
         ᶜqᵣ,
         ᶜqₛ,
         ᶜts,
-        ᶜΦ,
         dt,
         cmp,
         thp,
@@ -239,7 +243,6 @@ function compute_precipitation_cache!(Y, p, ::Microphysics1Moment, _)
         ᶜqᵣ,
         ᶜqₛ,
         ᶜts,
-        ᶜΦ,
         dt,
         cmp,
         thp,
@@ -255,8 +258,10 @@ function compute_precipitation_cache!(
 )
     FT = Spaces.undertype(axes(Y.c))
     (; dt) = p
+    thermo_params = CAP.thermodynamics_params(p.params)
+    grav = TDP.grav(thermo_params)
+    ᶜz = Fields.coordinate_field(axes(Y.c)).z
     (; ᶜts, ᶜqᵣ, ᶜqₛ, ᶜwᵣ, ᶜwₛ, ᶜu) = p.precomputed
-    (; ᶜΦ) = p.core
     # Grid mean precipitation sinks
     (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation
     # additional scratch storage
@@ -285,7 +290,6 @@ function compute_precipitation_cache!(
         ᶜqᵣ,
         ᶜqₛ,
         ᶜts,
-        ᶜΦ,
         dt,
         cmp,
         thp,
diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl
index 44dcbf6ba5..adde170128 100644
--- a/src/prognostic_equations/advection.jl
+++ b/src/prognostic_equations/advection.jl
@@ -8,7 +8,6 @@ import ClimaCore.Geometry as Geometry
 
 NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t)
     n = n_mass_flux_subdomains(p.atmos.turbconv_model)
-    (; ᶜΦ) = p.core
     (; ᶜu, ᶜK, ᶜp) = p.precomputed
     if p.atmos.turbconv_model isa AbstractEDMF
         (; ᶜu⁰) = p.precomputed
@@ -39,7 +38,11 @@ NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t)
         @. Yₜ.c.sgs⁰.ρatke -= wdivₕ(Y.c.sgs⁰.ρatke * ᶜu⁰)
     end
 
-    @. Yₜ.c.uₕ -= C12(gradₕ(ᶜp) / Y.c.ρ + gradₕ(ᶜK + ᶜΦ))
+    (; params) = p
+    thermo_params = CAP.thermodynamics_params(params)
+    grav = TDP.grav(thermo_params)
+    ᶜz = Fields.coordinate_field(axes(Y.c)).z
+    @. Yₜ.c.uₕ -= C12(gradₕ(ᶜp) / Y.c.ρ + gradₕ(ᶜK + Φ(grav, ᶜz)))
     # Without the C12(), the right-hand side would be a C1 or C2 in 2D space.
     return nothing
 end
@@ -66,20 +69,28 @@ NVTX.@annotate function horizontal_tracer_advection_tendency!(Yₜ, Y, p, t)
     return nothing
 end
 
+# TODO: move to ClimaCore:
+import ClimaCore.Geometry: CartesianGlobalGeometry
+Base.broadcastable(x::CartesianGlobalGeometry) = tuple(x)
+
 NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     (; turbconv_model) = p.atmos
+    (; params) = p
     n = n_prognostic_mass_flux_subdomains(turbconv_model)
     advect_tke = use_prognostic_tke(turbconv_model)
     point_type = eltype(Fields.coordinate_field(Y.c))
     (; dt) = p
-    ᶜJ = Fields.local_geometry_field(Y.c).J
-    (; ᶜf³, ᶠf¹², ᶜΦ) = p.core
+    ᶜlocal_geometry = Fields.local_geometry_field(Y.c)
+    ᶠlocal_geometry = Fields.local_geometry_field(Y.f)
+    ᶜJ = ᶜlocal_geometry.J
     (; ᶜu, ᶠu³, ᶜK) = p.precomputed
     (; edmfx_upwinding) = n > 0 || advect_tke ? p.atmos.numerics : all_nothing
     (; ᶜuʲs, ᶜKʲs, ᶠKᵥʲs) = n > 0 ? p.precomputed : all_nothing
     (; ᶠu³⁰) = advect_tke ? p.precomputed : all_nothing
     (; energy_upwinding, tracer_upwinding) = p.atmos.numerics
     (; ᶜspecific) = p.precomputed
+    ᶜcoords = Fields.coordinate_field(Y.c)
+    global_geom = Spaces.global_geometry(axes(ᶜcoords))
 
     ᶜρa⁰ = advect_tke ? (n > 0 ? p.precomputed.ᶜρa⁰ : Y.c.ρ) : nothing
     ᶜρ⁰ = advect_tke ? (n > 0 ? p.precomputed.ᶜρ⁰ : Y.c.ρ) : nothing
@@ -129,11 +140,11 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
         end
     end
 
-    if isnothing(ᶠf¹²)
+    if !(global_geom isa Geometry.DeepSphericalGlobalGeometry)
         # shallow atmosphere
         @. Yₜ.c.uₕ -=
             ᶜinterp(ᶠω¹² × (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³)) / (Y.c.ρ * ᶜJ) +
-            (ᶜf³ + ᶜω³) × CT12(ᶜu)
+            (f³(ᶜlocal_geometry, global_geom, params) + ᶜω³) × CT12(ᶜu)
         @. Yₜ.f.u₃ -= ᶠω¹² × ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK)
         for j in 1:n
             @. Yₜ.f.sgsʲs.:($$j).u₃ -=
@@ -143,12 +154,18 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     else
         # deep atmosphere
         @. Yₜ.c.uₕ -=
-            ᶜinterp((ᶠf¹² + ᶠω¹²) × (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³)) /
-            (Y.c.ρ * ᶜJ) + (ᶜf³ + ᶜω³) × CT12(ᶜu)
-        @. Yₜ.f.u₃ -= (ᶠf¹² + ᶠω¹²) × ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK)
+            ᶜinterp(
+                (f¹²(ᶠlocal_geometry, global_geom, params) + ᶠω¹²) ×
+                (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³),
+            ) / (Y.c.ρ * ᶜJ) +
+            (f³(ᶜlocal_geometry, global_geom, params) + ᶜω³) × CT12(ᶜu)
+        @. Yₜ.f.u₃ -=
+            (f¹²(ᶠlocal_geometry, global_geom, params) + ᶠω¹²) ×
+            ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK)
         for j in 1:n
             @. Yₜ.f.sgsʲs.:($$j).u₃ -=
-                (ᶠf¹² + ᶠω¹²ʲs.:($$j)) × ᶠinterp(CT12(ᶜuʲs.:($$j))) +
+                (f¹²(ᶠlocal_geometry, global_geom, params) + ᶠω¹²ʲs.:($$j)) ×
+                ᶠinterp(CT12(ᶜuʲs.:($$j))) +
                 ᶠgradᵥ(ᶜKʲs.:($$j) - ᶜinterp(ᶠKᵥʲs.:($$j)))
         end
     end
diff --git a/src/prognostic_equations/edmfx_precipitation.jl b/src/prognostic_equations/edmfx_precipitation.jl
index 52a697d161..5cad00374a 100644
--- a/src/prognostic_equations/edmfx_precipitation.jl
+++ b/src/prognostic_equations/edmfx_precipitation.jl
@@ -16,7 +16,8 @@ function edmfx_precipitation_tendency!(
     n = n_mass_flux_subdomains(turbconv_model)
     (; ᶜSqₜᵖʲs, ᶜtsʲs) = p.precomputed
     thermo_params = CAP.thermodynamics_params(p.params)
-    (; ᶜΦ) = p.core
+    grav = TDP.grav(thermo_params)
+    ᶜz = Fields.coordinate_field(axes(Y.c)).z
 
     for j in 1:n
 
@@ -27,7 +28,7 @@ function edmfx_precipitation_tendency!(
                 e_tot_0M_precipitation_sources_helper(
                     thermo_params,
                     ᶜtsʲs.:($$j),
-                    ᶜΦ,
+                    Φ(grav, ᶜz),
                 ) - TD.internal_energy(thermo_params, ᶜtsʲs.:($$j))
             )
 
@@ -48,7 +49,6 @@ function edmfx_precipitation_tendency!(
     n = n_mass_flux_subdomains(turbconv_model)
     (; ᶜSeₜᵖʲs, ᶜSqₜᵖʲs, ᶜtsʲs) = p.precomputed
     thp = CAP.thermodynamics_params(p.params)
-    (; ᶜΦ) = p.core
 
     for j in 1:n
 
diff --git a/src/prognostic_equations/implicit/implicit_solver.jl b/src/prognostic_equations/implicit/implicit_solver.jl
index d331f47a56..821439e01a 100644
--- a/src/prognostic_equations/implicit/implicit_solver.jl
+++ b/src/prognostic_equations/implicit/implicit_solver.jl
@@ -434,7 +434,6 @@ NVTX.@annotate function Wfact!(A, Y, p, dtγ, t)
                 p.precomputed.bdmr,
             ) : (;)
         )...,
-        p.core.ᶜΦ,
         p.core.ᶠgradᵥ_ᶜΦ,
         p.scratch.ᶜtemp_scalar,
         p.scratch.ᶜtemp_C3,
@@ -463,7 +462,7 @@ end
 
 function update_implicit_equation_jacobian!(A, Y, p, dtγ)
     (; matrix, diffusion_flag, sgs_advection_flag, topography_flag) = A
-    (; ᶜspecific, ᶜK, ᶜts, ᶜp, ᶜΦ, ᶠgradᵥ_ᶜΦ) = p
+    (; ᶜspecific, ᶜK, ᶜts, ᶜp, ᶠgradᵥ_ᶜΦ) = p
     (;
         ᶜtemp_C3,
         ∂ᶜK_∂ᶜuₕ,
@@ -490,6 +489,8 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
     # minus ∂e_int_∂q_tot
     ∂e_int_∂q_tot = T_0 * (Δcv_v - R_d) - FT(CAP.e_int_v0(params))
     thermo_params = CAP.thermodynamics_params(params)
+    grav = TDP.grav(thermo_params)
+    ᶜz = Fields.coordinate_field(axes(Y.c)).z
 
     ᶜρ = Y.c.ρ
     ᶜuₕ = Y.c.uₕ
@@ -550,7 +551,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
     @. ∂ᶠu₃_err_∂ᶜρ =
         dtγ * (
             ᶠp_grad_matrix ⋅
-            DiagonalMatrixRow(ᶜkappa_m * (T_0 * cp_d - ᶜK - ᶜΦ)) +
+            DiagonalMatrixRow(ᶜkappa_m * (T_0 * cp_d - ᶜK - Φ(grav, ᶜz))) +
             DiagonalMatrixRow(ᶠgradᵥ(ᶜp) / abs2(ᶠinterp(ᶜρ))) ⋅
             ᶠinterp_matrix()
         )