diff --git a/docs/src/PlotReferenceStates.jl b/docs/src/PlotReferenceStates.jl
index 8ab1f427ae..9e063bf2e9 100644
--- a/docs/src/PlotReferenceStates.jl
+++ b/docs/src/PlotReferenceStates.jl
@@ -24,7 +24,7 @@ function export_ref_profile(case_name::String)
     case = Cases.get_case(namelist)
     surf_ref_state = Cases.surface_ref_state(case, param_set, namelist)
 
-    aux_vars(FT) = (; ref_state = (ρ0 = FT(0), α0 = FT(0), p0 = FT(0)))
+    aux_vars(FT) = (; ref_state = (ρ = FT(0), p = FT(0)))
     aux_cent_fields = TC.FieldFromNamedTuple(TC.center_space(grid), aux_vars(FT))
     aux_face_fields = TC.FieldFromNamedTuple(TC.face_space(grid), aux_vars(FT))
 
@@ -35,39 +35,26 @@ function export_ref_profile(case_name::String)
 
     zc = vec(grid.zc)
     zf = vec(grid.zf)
-    ρ0_c = vec(aux.cent.ref_state.ρ0)
-    p0_c = vec(aux.cent.ref_state.p0)
-    α0_c = vec(aux.cent.ref_state.α0)
-    ρ0_f = vec(aux.face.ref_state.ρ0)
-    p0_f = vec(aux.face.ref_state.p0)
-    α0_f = vec(aux.face.ref_state.α0)
+    ρ_c = vec(aux.cent.ref_state.ρ)
+    p_c = vec(aux.cent.ref_state.p)
+    ρ_f = vec(aux.face.ref_state.ρ)
+    p_f = vec(aux.face.ref_state.p)
 
-    p1 = Plots.plot(ρ0_c, zc ./ 1000; label = "centers")
-    Plots.plot!(ρ0_f, zf ./ 1000; label = "faces")
+    p1 = Plots.plot(ρ_c, zc ./ 1000; label = "centers")
+    Plots.plot!(ρ_f, zf ./ 1000; label = "faces")
     Plots.plot!(size = (1000, 400))
     Plots.plot!(margin = 5 * Plots.mm)
     Plots.xlabel!("ρ_0")
     Plots.ylabel!("z (km)")
     Plots.title!("ρ_0")
 
-    p2 = Plots.plot(p0_c ./ 1000, zc ./ 1000; label = "centers")
-    Plots.plot!(p0_f ./ 1000, zf ./ 1000; label = "faces")
+    p2 = Plots.plot(p_c ./ 1000, zc ./ 1000; label = "centers")
+    Plots.plot!(p_f ./ 1000, zf ./ 1000; label = "faces")
     Plots.plot!(size = (1000, 400))
     Plots.plot!(margin = 5 * Plots.mm)
     Plots.xlabel!("p_0 (kPa)")
     Plots.ylabel!("z (km)")
     Plots.title!("p_0 (kPa)")
-
-    p3 = Plots.plot(α0_c, zc ./ 1000; label = "centers")
-    Plots.plot!(α0_f, zf ./ 1000; label = "faces")
-    Plots.plot!(size = (1000, 400))
-    Plots.plot!(margin = 5 * Plots.mm)
-    Plots.xlabel!("α_0")
-    Plots.ylabel!("z (km)")
-    Plots.title!("α_0")
-    Plots.plot(p1, p2, p3; layout = (1, 3))
-    Plots.savefig("$case_name.svg")
-
 end
 
 Logging.with_logger(Logging.NullLogger()) do # silence output
diff --git a/driver/Cases.jl b/driver/Cases.jl
index 7738d464dc..69452beccb 100644
--- a/driver/Cases.jl
+++ b/driver/Cases.jl
@@ -182,7 +182,7 @@ end
 function initialize_profiles(self::CasesBase{Soares}, grid::Grid, param_set, state)
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
-    ρ0_c = TC.center_ref_state(state).ρ0
+    ρ_c = prog_gm.ρ
     FT = eltype(grid)
     prof_q_tot = APL.Soares_q_tot(FT)
     prof_θ_liq_ice = APL.Soares_θ_liq_ice(FT)
@@ -192,26 +192,26 @@ function initialize_profiles(self::CasesBase{Soares}, grid::Grid, param_set, sta
     @inbounds for k in real_center_indices(grid)
         z = grid.zc[k]
         aux_gm.q_tot[k] = prof_q_tot(z)
-        prog_gm.ρq_tot[k] = ρ0_c[k] * aux_gm.q_tot[k]
+        prog_gm.ρq_tot[k] = ρ_c[k] * aux_gm.q_tot[k]
         aux_gm.θ_liq_ice[k] = prof_θ_liq_ice(z)
-        prog_gm.ρθ_liq_ice[k] = ρ0_c[k] * aux_gm.θ_liq_ice[k]
+        prog_gm.ρθ_liq_ice[k] = ρ_c[k] * aux_gm.θ_liq_ice[k]
         prog_gm.u[k] = prof_u(z)
         aux_gm.tke[k] = prof_tke(z)
     end
 end
 
 function surface_params(case::Soares, surf_ref_state, param_set; Ri_bulk_crit)
-    p0_f_surf = TD.air_pressure(param_set, surf_ref_state)
-    ρ0_f_surf = TD.air_density(param_set, surf_ref_state)
-    FT = eltype(p0_f_surf)
+    p_f_surf = TD.air_pressure(param_set, surf_ref_state)
+    ρ_f_surf = TD.air_density(param_set, surf_ref_state)
+    FT = eltype(p_f_surf)
     zrough = 0.16 #1.0e-4 0.16 is the value specified in the Nieuwstadt paper.
     Tsurface = 300.0
     qsurface = 5.0e-3
     θ_flux = 6.0e-2
     qt_flux = 2.5e-5
-    ts = TD.PhaseEquil_pTq(param_set, p0_f_surf, Tsurface, qsurface)
-    lhf = qt_flux * ρ0_f_surf * TD.latent_heat_vapor(param_set, ts)
-    shf = θ_flux * TD.cp_m(param_set, ts) * ρ0_f_surf
+    ts = TD.PhaseEquil_pTq(param_set, p_f_surf, Tsurface, qsurface)
+    lhf = qt_flux * ρ_f_surf * TD.latent_heat_vapor(param_set, ts)
+    shf = θ_flux * TD.cp_m(param_set, ts) * ρ_f_surf
     ustar = 0.28 # just to initilize grid mean covariances
     kwargs = (; zrough, Tsurface, qsurface, shf, lhf, ustar, Ri_bulk_crit)
     return TC.FixedSurfaceFlux(FT, TC.VariableFrictionVelocity; kwargs...)
@@ -232,7 +232,7 @@ end
 function initialize_profiles(self::CasesBase{Nieuwstadt}, grid::Grid, param_set, state)
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
-    ρ0_c = TC.center_ref_state(state).ρ0
+    ρ_c = prog_gm.ρ
 
     FT = eltype(grid)
     prof_θ_liq_ice = APL.Nieuwstadt_θ_liq_ice(FT)
@@ -242,23 +242,23 @@ function initialize_profiles(self::CasesBase{Nieuwstadt}, grid::Grid, param_set,
     @inbounds for k in real_center_indices(grid)
         z = grid.zc[k]
         aux_gm.θ_liq_ice[k] = prof_θ_liq_ice(z)
-        prog_gm.ρθ_liq_ice[k] = ρ0_c[k] * aux_gm.θ_liq_ice[k]
+        prog_gm.ρθ_liq_ice[k] = ρ_c[k] * aux_gm.θ_liq_ice[k]
         prog_gm.u[k] = prof_u(z)
         aux_gm.tke[k] = prof_tke(z)
     end
 end
 
 function surface_params(case::Nieuwstadt, surf_ref_state, param_set; Ri_bulk_crit)
-    p0_f_surf = TD.air_pressure(param_set, surf_ref_state)
-    ρ0_f_surf = TD.air_density(param_set, surf_ref_state)
-    FT = eltype(p0_f_surf)
+    p_f_surf = TD.air_pressure(param_set, surf_ref_state)
+    ρ_f_surf = TD.air_density(param_set, surf_ref_state)
+    FT = eltype(p_f_surf)
     zrough = 0.16 #1.0e-4 0.16 is the value specified in the Nieuwstadt paper.
     Tsurface = 300.0
     qsurface = 0.0
     θ_flux = 6.0e-2
     lhf = 0.0 # It would be 0.0 if we follow Nieuwstadt.
-    ts = TD.PhaseEquil_pTq(param_set, p0_f_surf, Tsurface, qsurface)
-    shf = θ_flux * TD.cp_m(param_set, ts) * ρ0_f_surf
+    ts = TD.PhaseEquil_pTq(param_set, p_f_surf, Tsurface, qsurface)
+    shf = θ_flux * TD.cp_m(param_set, ts) * ρ_f_surf
     ustar = 0.28 # just to initilize grid mean covariances
     kwargs = (; zrough, Tsurface, qsurface, shf, lhf, ustar, Ri_bulk_crit)
     return TC.FixedSurfaceFlux(FT, TC.VariableFrictionVelocity; kwargs...)
@@ -281,7 +281,7 @@ end
 function initialize_profiles(self::CasesBase{Bomex}, grid::Grid, param_set, state)
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
-    ρ0_c = TC.center_ref_state(state).ρ0
+    ρ_c = prog_gm.ρ
 
     FT = eltype(grid)
     prof_q_tot = APL.Bomex_q_tot(FT)
@@ -292,27 +292,27 @@ function initialize_profiles(self::CasesBase{Bomex}, grid::Grid, param_set, stat
     @inbounds for k in real_center_indices(grid)
         z = grid.zc[k]
         aux_gm.θ_liq_ice[k] = prof_θ_liq_ice(z)
-        prog_gm.ρθ_liq_ice[k] = ρ0_c[k] * aux_gm.θ_liq_ice[k]
+        prog_gm.ρθ_liq_ice[k] = ρ_c[k] * aux_gm.θ_liq_ice[k]
         aux_gm.q_tot[k] = prof_q_tot(z)
-        prog_gm.ρq_tot[k] = ρ0_c[k] * aux_gm.q_tot[k]
+        prog_gm.ρq_tot[k] = ρ_c[k] * aux_gm.q_tot[k]
         prog_gm.u[k] = prof_u(z)
         aux_gm.tke[k] = prof_tke(z)
     end
 end
 
 function surface_params(case::Bomex, surf_ref_state, param_set; Ri_bulk_crit)
-    p0_f_surf = TD.air_pressure(param_set, surf_ref_state)
-    ρ0_f_surf = TD.air_density(param_set, surf_ref_state)
-    FT = eltype(p0_f_surf)
+    p_f_surf = TD.air_pressure(param_set, surf_ref_state)
+    ρ_f_surf = TD.air_density(param_set, surf_ref_state)
+    FT = eltype(p_f_surf)
     zrough = 1.0e-4
     qsurface = 22.45e-3 # kg/kg
     θ_surface = 299.1
     θ_flux = 8.0e-3
     qt_flux = 5.2e-5
-    ts = TD.PhaseEquil_pθq(param_set, p0_f_surf, θ_surface, qsurface)
+    ts = TD.PhaseEquil_pθq(param_set, p_f_surf, θ_surface, qsurface)
     Tsurface = TD.air_temperature(param_set, ts)
-    lhf = qt_flux * ρ0_f_surf * TD.latent_heat_vapor(param_set, ts)
-    shf = θ_flux * TD.cp_m(param_set, ts) * ρ0_f_surf
+    lhf = qt_flux * ρ_f_surf * TD.latent_heat_vapor(param_set, ts)
+    shf = θ_flux * TD.cp_m(param_set, ts) * ρ_f_surf
     ustar = 0.28 # m/s
     kwargs = (; zrough, Tsurface, qsurface, shf, lhf, ustar, Ri_bulk_crit)
     return TC.FixedSurfaceFlux(FT, TC.FixedFrictionVelocity; kwargs...)
@@ -320,10 +320,10 @@ end
 
 function initialize_forcing(self::CasesBase{Bomex}, grid::Grid, state, param_set)
     initialize(self.Fo, grid, state)
-    p0_c = TC.center_ref_state(state).p0
     prog_gm = TC.center_prog_grid_mean(state)
     aux_gm = TC.center_aux_grid_mean(state)
     ts_gm = aux_gm.ts
+    p_c = aux_gm.p
 
     FT = eltype(grid)
     prof_ug = APL.Bomex_geostrophic_u(FT)
@@ -362,7 +362,7 @@ end
 function initialize_profiles(self::CasesBase{life_cycle_Tan2018}, grid::Grid, param_set, state)
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
-    ρ0_c = TC.center_ref_state(state).ρ0
+    ρ_c = prog_gm.ρ
 
     FT = eltype(grid)
     prof_q_tot = APL.LifeCycleTan2018_q_tot(FT)
@@ -373,9 +373,9 @@ function initialize_profiles(self::CasesBase{life_cycle_Tan2018}, grid::Grid, pa
     @inbounds for k in real_center_indices(grid)
         z = grid.zc[k]
         aux_gm.θ_liq_ice[k] = prof_θ_liq_ice(z)
-        prog_gm.ρθ_liq_ice[k] = ρ0_c[k] * aux_gm.θ_liq_ice[k]
+        prog_gm.ρθ_liq_ice[k] = ρ_c[k] * aux_gm.θ_liq_ice[k]
         aux_gm.q_tot[k] = prof_q_tot(z)
-        prog_gm.ρq_tot[k] = ρ0_c[k] * aux_gm.q_tot[k]
+        prog_gm.ρq_tot[k] = ρ_c[k] * aux_gm.q_tot[k]
         prog_gm.u[k] = prof_u(z)
         aux_gm.tke[k] = prof_tke(z)
     end
@@ -391,18 +391,18 @@ function life_cycle_buoyancy_flux(param_set, weight = 1)
 end
 
 function surface_params(case::life_cycle_Tan2018, surf_ref_state, param_set; Ri_bulk_crit)
-    p0_f_surf = TD.air_pressure(param_set, surf_ref_state)
-    ρ0_f_surf = TD.air_density(param_set, surf_ref_state)
-    FT = eltype(p0_f_surf)
+    p_f_surf = TD.air_pressure(param_set, surf_ref_state)
+    ρ_f_surf = TD.air_density(param_set, surf_ref_state)
+    FT = eltype(p_f_surf)
     zrough = 1.0e-4 # not actually used, but initialized to reasonable value
     qsurface = 22.45e-3 # kg/kg
     θ_surface = 299.1
     θ_flux = 8.0e-3
     qt_flux = 5.2e-5
-    ts = TD.PhaseEquil_pθq(param_set, p0_f_surf, θ_surface, qsurface)
+    ts = TD.PhaseEquil_pθq(param_set, p_f_surf, θ_surface, qsurface)
     Tsurface = TD.air_temperature(param_set, ts)
-    lhf0 = qt_flux * ρ0_f_surf * TD.latent_heat_vapor(param_set, ts)
-    shf0 = θ_flux * TD.cp_m(param_set, ts) * ρ0_f_surf
+    lhf0 = qt_flux * ρ_f_surf * TD.latent_heat_vapor(param_set, ts)
+    shf0 = θ_flux * TD.cp_m(param_set, ts) * ρ_f_surf
 
     weight_factor(t) = 0.01 + 0.99 * (cos(2.0 * π * t / 3600.0) + 1.0) / 2.0
     weight = 1.0
@@ -416,9 +416,9 @@ end
 
 function initialize_forcing(self::CasesBase{life_cycle_Tan2018}, grid::Grid, state, param_set)
     initialize(self.Fo, grid, state)
-    p0_c = TC.center_ref_state(state).p0
     prog_gm = TC.center_prog_grid_mean(state)
     aux_gm = TC.center_aux_grid_mean(state)
+    p_c = aux_gm.p
     ts_gm = aux_gm.ts
 
     FT = eltype(grid)
@@ -468,8 +468,8 @@ function initialize_profiles(self::CasesBase{Rico}, grid::Grid, param_set, state
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
     aux_tc = TC.center_aux_turbconv(state)
-    p0 = TC.center_ref_state(state).p0
-    ρ0_c = TC.center_ref_state(state).ρ0
+    p = aux_gm.p
+    ρ_c = prog_gm.ρ
 
     FT = eltype(grid)
     prof_u = APL.Rico_u(FT)
@@ -482,16 +482,16 @@ function initialize_profiles(self::CasesBase{Rico}, grid::Grid, param_set, state
         prog_gm.u[k] = prof_u(z)
         prog_gm.v[k] = prof_v(z)
         aux_gm.θ_liq_ice[k] = prof_θ_liq_ice(z)
-        prog_gm.ρθ_liq_ice[k] = ρ0_c[k] * aux_gm.θ_liq_ice[k]
+        prog_gm.ρθ_liq_ice[k] = ρ_c[k] * aux_gm.θ_liq_ice[k]
         aux_gm.q_tot[k] = prof_q_tot(z)
-        prog_gm.ρq_tot[k] = ρ0_c[k] * aux_gm.q_tot[k]
+        prog_gm.ρq_tot[k] = ρ_c[k] * aux_gm.q_tot[k]
     end
 
     # Need to get θ_virt
     @inbounds for k in real_center_indices(grid)
         # Thermo state field cache is not yet
         # defined, so we can't use it yet.
-        ts = TD.PhaseEquil_pθq(param_set, p0[k], aux_gm.θ_liq_ice[k], aux_gm.q_tot[k])
+        ts = TD.PhaseEquil_pθq(param_set, p[k], aux_gm.θ_liq_ice[k], aux_gm.q_tot[k])
         aux_gm.θ_virt[k] = TD.virtual_pottemp(param_set, ts)
     end
     zi = 0.6 * get_inversion(grid, state, param_set, 0.2)
@@ -507,8 +507,8 @@ function initialize_profiles(self::CasesBase{Rico}, grid::Grid, param_set, state
 end
 
 function surface_params(case::Rico, surf_ref_state, param_set; kwargs...)
-    p0_f_surf = TD.air_pressure(param_set, surf_ref_state)
-    FT = eltype(p0_f_surf)
+    p_f_surf = TD.air_pressure(param_set, surf_ref_state)
+    FT = eltype(p_f_surf)
 
     zrough = 0.00015
     cm0 = 0.001229
@@ -522,7 +522,7 @@ function surface_params(case::Rico, surf_ref_state, param_set; kwargs...)
     Tsurface = 299.8
 
     # For Rico we provide values of transfer coefficients
-    ts = TD.PhaseEquil_pTq(param_set, p0_f_surf, Tsurface, FT(0)) # TODO: is this correct?
+    ts = TD.PhaseEquil_pTq(param_set, p_f_surf, Tsurface, FT(0)) # TODO: is this correct?
     qsurface = TD.q_vap_saturation(param_set, ts)
     kwargs = (; zrough, Tsurface, qsurface, cm, ch)
     return TC.FixedSurfaceCoeffs(FT; kwargs...)
@@ -530,10 +530,10 @@ end
 
 function initialize_forcing(self::CasesBase{Rico}, grid::Grid, state, param_set)
     initialize(self.Fo, grid, state)
-    p0_c = TC.center_ref_state(state).p0
     prog_gm = TC.center_prog_grid_mean(state)
     aux_gm = TC.center_aux_grid_mean(state)
     ts_gm = aux_gm.ts
+    p_c = aux_gm.p
 
     FT = eltype(grid)
     prof_ug = APL.Rico_geostrophic_ug(FT)
@@ -565,16 +565,16 @@ ForcingBase(case::TRMM_LBA, param_set::APS; kwargs...) =
 function surface_ref_state(::TRMM_LBA, param_set::APS, namelist)
     molmass_ratio = CPP.molmass_ratio(param_set)
     Pg = 991.3 * 100  #Pressure at ground
-    Tg = 296.85   # surface values for reference state (RS) which outputs p0 rho0 alpha0
+    Tg = 296.85   # surface values for reference state (RS) which outputs p rho0 alpha0
     pvg = TD.saturation_vapor_pressure(param_set, Tg, TD.Liquid())
     qtg = (1 / molmass_ratio) * pvg / (Pg - pvg) #Total water mixing ratio at surface
     return TD.PhaseEquil_pTq(param_set, Pg, Tg, qtg)
 end
 function initialize_profiles(self::CasesBase{TRMM_LBA}, grid::Grid, param_set, state)
-    p0 = TC.center_ref_state(state).p0
-    ρ0_c = TC.center_ref_state(state).ρ0
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
+    ρ_c = prog_gm.ρ
+    p = aux_gm.p
 
     FT = eltype(grid)
     # Get profiles from AtmosphericProfilesLibrary.jl
@@ -601,21 +601,21 @@ function initialize_profiles(self::CasesBase{TRMM_LBA}, grid::Grid, param_set, s
         denom = (prof_p(z) - pv_star + (1 / molmass_ratio) * pv_star * RH / 100.0)
         qv_star = pv_star * (1 / molmass_ratio) / denom
         aux_gm.q_tot[k] = qv_star * RH / 100
-        prog_gm.ρq_tot[k] = ρ0_c[k] * aux_gm.q_tot[k]
+        prog_gm.ρq_tot[k] = ρ_c[k] * aux_gm.q_tot[k]
         phase_part = TD.PhasePartition(aux_gm.q_tot[k], 0.0, 0.0) # initial state is not saturated
-        aux_gm.θ_liq_ice[k] = TD.liquid_ice_pottemp_given_pressure(param_set, aux_gm.T[k], p0[k], phase_part)
-        prog_gm.ρθ_liq_ice[k] = ρ0_c[k] * aux_gm.θ_liq_ice[k]
+        aux_gm.θ_liq_ice[k] = TD.liquid_ice_pottemp_given_pressure(param_set, aux_gm.T[k], p[k], phase_part)
+        prog_gm.ρθ_liq_ice[k] = ρ_c[k] * aux_gm.θ_liq_ice[k]
         aux_gm.tke[k] = prof_tke(z)
     end
 end
 
 function surface_params(case::TRMM_LBA, surf_ref_state, param_set; Ri_bulk_crit)
-    p0_f_surf = TD.air_pressure(param_set, surf_ref_state)
-    FT = eltype(p0_f_surf)
+    p_f_surf = TD.air_pressure(param_set, surf_ref_state)
+    FT = eltype(p_f_surf)
     # zrough = 1.0e-4 # not actually used, but initialized to reasonable value
     qsurface = 22.45e-3 # kg/kg
     θ_surface = (273.15 + 23)
-    ts = TD.PhaseEquil_pθq(param_set, p0_f_surf, θ_surface, qsurface)
+    ts = TD.PhaseEquil_pθq(param_set, p_f_surf, θ_surface, qsurface)
     Tsurface = TD.air_temperature(param_set, ts)
     ustar = 0.28 # this is taken from Bomex -- better option is to approximate from LES tke above the surface
     lhf = t -> 554.0 * max(0, cos(π / 2 * ((5.25 * 3600.0 - t) / 5.25 / 3600.0)))^1.3
@@ -647,17 +647,17 @@ ForcingBase(case::ARM_SGP, param_set::APS; kwargs...) =
 
 function surface_ref_state(::ARM_SGP, param_set::APS, namelist)
     Pg = 970.0 * 100 #Pressure at ground
-    Tg = 299.0   # surface values for reference state (RS) which outputs p0 rho0 alpha0
+    Tg = 299.0   # surface values for reference state (RS) which outputs p rho0 alpha0
     qtg = 15.2 / 1000 #Total water mixing ratio at surface
     return TD.PhaseEquil_pTq(param_set, Pg, Tg, qtg)
 end
 
 function initialize_profiles(self::CasesBase{ARM_SGP}, grid::Grid, param_set, state)
     # ARM_SGP inputs
-    p0 = TC.center_ref_state(state).p0
     prog_gm = TC.center_prog_grid_mean(state)
     aux_gm = TC.center_aux_grid_mean(state)
-    ρ0_c = TC.center_ref_state(state).ρ0
+    ρ_c = prog_gm.ρ
+    p = aux_gm.p
 
     FT = eltype(grid)
     prof_u = APL.ARM_SGP_u(FT)
@@ -669,23 +669,23 @@ function initialize_profiles(self::CasesBase{ARM_SGP}, grid::Grid, param_set, st
         z = grid.zc[k]
         # TODO figure out how to use ts here
         phase_part = TD.PhasePartition(aux_gm.q_tot[k], aux_gm.q_liq[k], 0.0)
-        Π = TD.exner_given_pressure(param_set, p0[k], phase_part)
+        Π = TD.exner_given_pressure(param_set, p[k], phase_part)
         prog_gm.u[k] = prof_u(z)
         aux_gm.q_tot[k] = prof_q_tot(z)
-        prog_gm.ρq_tot[k] = ρ0_c[k] * aux_gm.q_tot[k]
+        prog_gm.ρq_tot[k] = ρ_c[k] * aux_gm.q_tot[k]
         aux_gm.T[k] = prof_θ_liq_ice(z) * Π
-        aux_gm.θ_liq_ice[k] = TD.liquid_ice_pottemp_given_pressure(param_set, aux_gm.T[k], p0[k], phase_part)
-        prog_gm.ρθ_liq_ice[k] = ρ0_c[k] * aux_gm.θ_liq_ice[k]
+        aux_gm.θ_liq_ice[k] = TD.liquid_ice_pottemp_given_pressure(param_set, aux_gm.T[k], p[k], phase_part)
+        prog_gm.ρθ_liq_ice[k] = ρ_c[k] * aux_gm.θ_liq_ice[k]
         aux_gm.tke[k] = prof_tke(z)
     end
 end
 
 function surface_params(case::ARM_SGP, surf_ref_state, param_set; Ri_bulk_crit)
-    p0_f_surf = TD.air_pressure(param_set, surf_ref_state)
-    FT = eltype(p0_f_surf)
+    p_f_surf = TD.air_pressure(param_set, surf_ref_state)
+    FT = eltype(p_f_surf)
     qsurface = 15.2e-3 # kg/kg
     θ_surface = 299.0
-    ts = TD.PhaseEquil_pθq(param_set, p0_f_surf, θ_surface, qsurface)
+    ts = TD.PhaseEquil_pθq(param_set, p_f_surf, θ_surface, qsurface)
     Tsurface = TD.air_temperature(param_set, ts)
     ustar = 0.28 # this is taken from Bomex -- better option is to approximate from LES tke above the surface
 
@@ -710,9 +710,9 @@ end
 
 function update_forcing(self::CasesBase{ARM_SGP}, grid, state, t::Real, param_set)
     aux_gm = TC.center_aux_grid_mean(state)
-    p0_c = TC.center_ref_state(state).p0
     ts_gm = TC.center_aux_grid_mean(state).ts
     prog_gm = TC.center_prog_grid_mean(state)
+    p_c = prog_gm.ρ
     FT = eltype(grid)
     @inbounds for k in real_center_indices(grid)
         Π = TD.exner(param_set, ts_gm[k])
@@ -731,33 +731,33 @@ ForcingBase(case::GATE_III, param_set::APS; kwargs...) = ForcingBase(get_forcing
 
 function surface_ref_state(::GATE_III, param_set::APS, namelist)
     Pg = 1013.0 * 100  #Pressure at ground
-    Tg = 299.184   # surface values for reference state (RS) which outputs p0 rho0 alpha0
+    Tg = 299.184   # surface values for reference state (RS) which outputs p rho0 alpha0
     qtg = 16.5 / 1000 #Total water mixing ratio at surface
     return TD.PhaseEquil_pTq(param_set, Pg, Tg, qtg)
 end
 
 function initialize_profiles(self::CasesBase{GATE_III}, grid::Grid, param_set, state)
     FT = eltype(grid)
-    p0 = TC.center_ref_state(state).p0
-    ρ0_c = TC.center_ref_state(state).ρ0
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
+    p = aux_gm.p
+    ρ_c = prog_gm.ρ
     @inbounds for k in real_center_indices(grid)
         z = grid.zc[k]
         aux_gm.q_tot[k] = APL.GATE_III_q_tot(FT)(z)
-        prog_gm.ρq_tot[k] = ρ0_c[k] * aux_gm.q_tot[k]
+        prog_gm.ρq_tot[k] = ρ_c[k] * aux_gm.q_tot[k]
         aux_gm.T[k] = APL.GATE_III_T(FT)(z)
         prog_gm.u[k] = APL.GATE_III_u(FT)(z)
-        ts = TD.PhaseEquil_pTq(param_set, p0[k], aux_gm.T[k], aux_gm.q_tot[k])
+        ts = TD.PhaseEquil_pTq(param_set, p[k], aux_gm.T[k], aux_gm.q_tot[k])
         aux_gm.θ_liq_ice[k] = TD.liquid_ice_pottemp(param_set, ts)
-        prog_gm.ρθ_liq_ice[k] = ρ0_c[k] * aux_gm.θ_liq_ice[k]
+        prog_gm.ρθ_liq_ice[k] = ρ_c[k] * aux_gm.θ_liq_ice[k]
         aux_gm.tke[k] = APL.GATE_III_tke(FT)(z)
     end
 end
 
 function surface_params(case::GATE_III, surf_ref_state, param_set; kwargs...)
-    p0_f_surf = TD.air_pressure(param_set, surf_ref_state)
-    FT = eltype(p0_f_surf)
+    p_f_surf = TD.air_pressure(param_set, surf_ref_state)
+    FT = eltype(p_f_surf)
 
     qsurface = 16.5 / 1000.0 # kg/kg
     cm = zc_surf -> 0.0012
@@ -766,7 +766,7 @@ function surface_params(case::GATE_III, surf_ref_state, param_set; kwargs...)
     Tsurface = 299.184
 
     # For GATE_III we provide values of transfer coefficients
-    ts = TD.PhaseEquil_pθq(param_set, p0_f_surf, Tsurface, qsurface)
+    ts = TD.PhaseEquil_pθq(param_set, p_f_surf, Tsurface, qsurface)
     qsurface = TD.q_vap_saturation(param_set, ts)
     kwargs = (; Tsurface, qsurface, cm, ch)
     return TC.FixedSurfaceCoeffs(FT; kwargs...)
@@ -797,17 +797,17 @@ end
 
 function initialize_profiles(self::CasesBase{DYCOMS_RF01}, grid::Grid, param_set, state)
     FT = eltype(grid)
-    p0 = TC.center_ref_state(state).p0
-    ρ0_c = TC.center_ref_state(state).ρ0
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
+    ρ_c = prog_gm.ρ
+    p = aux_gm.p
     @inbounds for k in real_center_indices(grid)
         # thetal profile as defined in DYCOMS
         z = grid.zc[k]
         aux_gm.q_tot[k] = APL.Dycoms_RF01_q_tot(FT)(z)
-        prog_gm.ρq_tot[k] = ρ0_c[k] * aux_gm.q_tot[k]
+        prog_gm.ρq_tot[k] = ρ_c[k] * aux_gm.q_tot[k]
         aux_gm.θ_liq_ice[k] = APL.Dycoms_RF01_θ_liq_ice(FT)(z)
-        prog_gm.ρθ_liq_ice[k] = ρ0_c[k] * aux_gm.θ_liq_ice[k]
+        prog_gm.ρθ_liq_ice[k] = ρ_c[k] * aux_gm.θ_liq_ice[k]
         # velocity profile (geostrophic)
         prog_gm.u[k] = APL.Dycoms_RF01_u0(FT)(z)
         prog_gm.v[k] = APL.Dycoms_RF01_v0(FT)(z)
@@ -882,17 +882,17 @@ end
 
 function initialize_profiles(self::CasesBase{DYCOMS_RF02}, grid::Grid, param_set, state)
     FT = eltype(grid)
-    p0 = TC.center_ref_state(state).p0
-    ρ0_c = TC.center_ref_state(state).ρ0
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
+    p = aux_gm.p
+    ρ_c = prog_gm.ρ
     @inbounds for k in real_center_indices(grid)
         # θ_liq_ice profile as defined in DYCOM RF02
         z = grid.zc[k]
         aux_gm.q_tot[k] = APL.Dycoms_RF02_q_tot(FT)(z)
-        prog_gm.ρq_tot[k] = ρ0_c[k] * aux_gm.q_tot[k]
+        prog_gm.ρq_tot[k] = ρ_c[k] * aux_gm.q_tot[k]
         aux_gm.θ_liq_ice[k] = APL.Dycoms_RF02_θ_liq_ice(FT)(z)
-        prog_gm.ρθ_liq_ice[k] = ρ0_c[k] * aux_gm.θ_liq_ice[k]
+        prog_gm.ρθ_liq_ice[k] = ρ_c[k] * aux_gm.θ_liq_ice[k]
         # velocity profile
         prog_gm.v[k] = APL.Dycoms_RF02_v(FT)(z)
         prog_gm.u[k] = APL.Dycoms_RF02_u(FT)(z)
@@ -973,7 +973,7 @@ end
 function initialize_profiles(self::CasesBase{GABLS}, grid::Grid, param_set, state)
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
-    ρ0_c = TC.center_ref_state(state).ρ0
+    ρ_c = prog_gm.ρ
     FT = eltype(grid)
 
     @inbounds for k in real_center_indices(grid)
@@ -982,9 +982,9 @@ function initialize_profiles(self::CasesBase{GABLS}, grid::Grid, param_set, stat
         prog_gm.u[k] = APL.GABLS_u(FT)(z)
         prog_gm.v[k] = APL.GABLS_v(FT)(z)
         aux_gm.θ_liq_ice[k] = APL.GABLS_θ_liq_ice(FT)(z)
-        prog_gm.ρθ_liq_ice[k] = ρ0_c[k] * aux_gm.θ_liq_ice[k]
+        prog_gm.ρθ_liq_ice[k] = ρ_c[k] * aux_gm.θ_liq_ice[k]
         aux_gm.q_tot[k] = APL.GABLS_q_tot(FT)(z)
-        prog_gm.ρq_tot[k] = ρ0_c[k] * aux_gm.q_tot[k]
+        prog_gm.ρq_tot[k] = ρ_c[k] * aux_gm.q_tot[k]
         aux_gm.tke[k] = APL.GABLS_tke(FT)(z)
         aux_gm.Hvar[k] = aux_gm.tke[k]
     end
@@ -1031,23 +1031,23 @@ end
 function surface_ref_state(::SP, param_set::APS, namelist)
     Pg = 1.0e5  #Pressure at ground
     Tg = 300.0  #Temperature at ground
-    qtg = 1.0e-4   #Total water mixing ratio at TC. if set to 0, alpha0, rho0, p0 are NaN.
+    qtg = 1.0e-4   #Total water mixing ratio at TC. if set to 0, alpha0, rho0, p are NaN.
     return TD.PhaseEquil_pTq(param_set, Pg, Tg, qtg)
 end
 
 function initialize_profiles(self::CasesBase{SP}, grid::Grid, param_set, state)
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
-    ρ0_c = TC.center_ref_state(state).ρ0
+    ρ_c = prog_gm.ρ
     FT = eltype(grid)
     @inbounds for k in real_center_indices(grid)
         z = grid.zc[k]
         prog_gm.u[k] = APL.SP_u(FT)(z)
         prog_gm.v[k] = APL.SP_v(FT)(z)
         aux_gm.θ_liq_ice[k] = APL.SP_θ_liq_ice(FT)(z)
-        prog_gm.ρθ_liq_ice[k] = ρ0_c[k] * aux_gm.θ_liq_ice[k]
+        prog_gm.ρθ_liq_ice[k] = ρ_c[k] * aux_gm.θ_liq_ice[k]
         aux_gm.q_tot[k] = APL.SP_q_tot(FT)(z)
-        prog_gm.ρq_tot[k] = ρ0_c[k] * aux_gm.q_tot[k]
+        prog_gm.ρq_tot[k] = ρ_c[k] * aux_gm.q_tot[k]
         aux_gm.tke[k] = APL.SP_tke(FT)(z)
     end
 end
@@ -1078,13 +1078,13 @@ function initialize_profiles(self::CasesBase{DryBubble}, grid::Grid, param_set,
 
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
-    ρ0_c = TC.center_ref_state(state).ρ0
+    ρ_c = prog_gm.ρ
     # initialize Grid Mean Profiles of thetali and qt
     zc_in = grid.zc
     FT = eltype(grid)
     prof_θ_liq_ice = APL.DryBubble_θ_liq_ice(FT)
     aux_gm.θ_liq_ice .= prof_θ_liq_ice.(zc_in)
-    @. prog_gm.ρθ_liq_ice = ρ0_c * aux_gm.θ_liq_ice
+    @. prog_gm.ρθ_liq_ice = ρ_c * aux_gm.θ_liq_ice
     parent(prog_gm.u) .= 0.01
     parent(prog_gm.ρq_tot) .= 0
     parent(aux_gm.q_tot) .= 0
@@ -1153,7 +1153,7 @@ function initialize_profiles(self::CasesBase{LES_driven_SCM}, grid::Grid, param_
 
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
-    ρ0_c = TC.center_ref_state(state).ρ0
+    ρ_c = prog_gm.ρ
     NC.Dataset(self.LESDat.les_filename, "r") do data
         t = data.group["profiles"]["t"][:]
         # define time interval
@@ -1172,8 +1172,8 @@ function initialize_profiles(self::CasesBase{LES_driven_SCM}, grid::Grid, param_
         parent(aux_gm.q_tot) .= pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "qt_mean", imin, imax))
         parent(prog_gm.u) .= pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "u_mean", imin, imax))
         parent(prog_gm.v) .= pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "v_mean", imin, imax))
-        @. prog_gm.ρθ_liq_ice = ρ0_c * aux_gm.θ_liq_ice
-        @. prog_gm.ρq_tot .= ρ0_c * aux_gm.q_tot
+        @. prog_gm.ρθ_liq_ice = ρ_c * aux_gm.θ_liq_ice
+        @. prog_gm.ρq_tot .= ρ_c * aux_gm.q_tot
     end
     @inbounds for k in real_center_indices(grid)
         z = grid.zc[k]
diff --git a/driver/Radiation.jl b/driver/Radiation.jl
index 6ea28a2cf8..ea4221b1fe 100644
--- a/driver/Radiation.jl
+++ b/driver/Radiation.jl
@@ -6,12 +6,12 @@ see eq. 3 in Stevens et. al. 2005 DYCOMS paper
 """
 function update_radiation(self::TC.RadiationBase{TC.RadiationDYCOMS_RF01}, grid, state, t::Real, param_set)
     cp_d = CPP.cp_d(param_set)
-    ρ0_f = TC.face_ref_state(state).ρ0
-    ρ0_c = TC.center_ref_state(state).ρ0
     aux_gm = TC.center_aux_grid_mean(state)
     aux_gm_f = TC.face_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
     q_tot_f = TC.face_aux_turbconv(state).ϕ_temporary
+    ρ_f = aux_gm_f.ρ
+    ρ_c = prog_gm.ρ
     # find zi (level of 8.0 g/kg isoline of qt)
     # TODO: report bug: zi and ρ_i are not initialized
     zi = 0
@@ -25,12 +25,12 @@ function update_radiation(self::TC.RadiationBase{TC.RadiationDYCOMS_RF01}, grid,
             idx_zi = k
             # will be used at cell faces
             zi = grid.zf[k]
-            ρ_i = ρ0_f[k]
+            ρ_i = ρ_f[k]
             break
         end
     end
 
-    ρ_z = Dierckx.Spline1D(vec(grid.zc), vec(ρ0_c); k = 1)
+    ρ_z = Dierckx.Spline1D(vec(grid.zc), vec(ρ_c); k = 1)
     q_liq_z = Dierckx.Spline1D(vec(grid.zc), vec(aux_gm.q_liq); k = 1)
 
     integrand(ρq_l, params, z) = params.κ * ρ_z(z) * q_liq_z(z)
@@ -61,7 +61,7 @@ function update_radiation(self::TC.RadiationBase{TC.RadiationDYCOMS_RF01}, grid,
 
     ∇c = CCO.DivergenceF2C()
     wvec = CC.Geometry.WVector
-    @. aux_gm.dTdt_rad = -∇c(wvec(aux_gm_f.f_rad)) / ρ0_c / cp_d
+    @. aux_gm.dTdt_rad = -∇c(wvec(aux_gm_f.f_rad)) / ρ_c / cp_d
 
     return
 end
diff --git a/driver/Surface.jl b/driver/Surface.jl
index 9e9aa44946..19c264a320 100644
--- a/driver/Surface.jl
+++ b/driver/Surface.jl
@@ -19,7 +19,8 @@ function get_surface(
     z_in = grid.zc[kc_surf].z
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
-    p0_f_surf = TC.face_ref_state(state).p0[kf_surf]
+    aux_gm_f = TC.face_aux_grid_mean(state)
+    p_f_surf = aux_gm_f.p[kf_surf]
     u_gm_surf = prog_gm.u[kc_surf]
     v_gm_surf = prog_gm.v[kc_surf]
     Tsurface = TC.surface_temperature(surf_params, t)
@@ -29,7 +30,7 @@ function get_surface(
     Ri_bulk_crit = surf_params.Ri_bulk_crit
     zrough = surf_params.zrough
 
-    ts_sfc = TD.PhaseEquil_pTq(param_set, p0_f_surf, Tsurface, qsurface)
+    ts_sfc = TD.PhaseEquil_pTq(param_set, p_f_surf, Tsurface, qsurface)
     ts_in = aux_gm.ts[kc_surf]
     universal_func = UF.Businger()
     scheme = SF.FVScheme()
@@ -85,13 +86,14 @@ function get_surface(
     FT = eltype(grid)
     kc_surf = TC.kc_surface(grid)
     kf_surf = TC.kf_surface(grid)
-    p0_f_surf = TC.face_ref_state(state).p0[kf_surf]
+    aux_gm_f = TC.face_aux_grid_mean(state)
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
-    u_gm_surf = prog_gm.u[kc_surf]
-    v_gm_surf = prog_gm.v[kc_surf]
     Tsurface = TC.surface_temperature(surf_params, t)
     qsurface = TC.surface_q_tot(surf_params, t)
+    p_f_surf = aux_gm_f.p[kf_surf]
+    u_gm_surf = prog_gm.u[kc_surf]
+    v_gm_surf = prog_gm.v[kc_surf]
     zrough = surf_params.zrough
     zc_surf = grid.zc[kc_surf]
     cm = surf_params.cm(zc_surf)
@@ -102,7 +104,7 @@ function get_surface(
     scheme = SF.FVScheme()
     z_sfc = FT(0)
     z_in = grid.zc[kc_surf].z
-    ts_sfc = TD.PhaseEquil_pθq(param_set, p0_f_surf, Tsurface, qsurface)
+    ts_sfc = TD.PhaseEquil_pθq(param_set, p_f_surf, Tsurface, qsurface)
     ts_in = aux_gm.ts[kc_surf]
     u_sfc = SA.SVector{2, FT}(0, 0)
     u_in = SA.SVector{2, FT}(u_gm_surf, v_gm_surf)
@@ -145,10 +147,11 @@ function get_surface(
     FT = eltype(grid)
     z_sfc = FT(0)
     z_in = grid.zc[kc_surf].z
-    p0_f_surf = TC.face_ref_state(state).p0[kf_surf]
-    p0_c_surf = TC.center_ref_state(state).p0[kc_surf]
     prog_gm = TC.center_prog_grid_mean(state)
     aux_gm = TC.center_aux_grid_mean(state)
+    aux_gm_f = TC.face_aux_grid_mean(state)
+    p_c_surf = aux_gm.p[kf_surf]
+    p_f_surf = aux_gm_f.p[kf_surf]
     u_gm_surf = prog_gm.u[kc_surf]
     v_gm_surf = prog_gm.v[kc_surf]
     q_tot_gm_surf = aux_gm.q_tot[kc_surf]
@@ -162,8 +165,8 @@ function get_surface(
 
     universal_func = UF.Businger()
     scheme = SF.FVScheme()
-    ts_sfc = TD.PhaseEquil_pTq(param_set, p0_f_surf, Tsurface, qsurface)
-    ts_in = TD.PhaseEquil_pθq(param_set, p0_c_surf, θ_liq_ice_gm_surf, qsurface)
+    ts_sfc = TD.PhaseEquil_pTq(param_set, p_f_surf, Tsurface, qsurface)
+    ts_in = TD.PhaseEquil_pθq(param_set, p_c_surf, θ_liq_ice_gm_surf, qsurface)
 
     u_sfc = SA.SVector{2, FT}(0, 0)
     u_in = SA.SVector{2, FT}(u_gm_surf, v_gm_surf)
@@ -205,10 +208,10 @@ function get_surface(
     FT = eltype(grid)
     z_sfc = FT(0)
     z_in = grid.zc[kc_surf].z
-    p0_f_surf = TC.face_ref_state(state).p0[kf_surf]
-    p0_c_surf = TC.center_ref_state(state).p0[kc_surf]
     prog_gm = TC.center_prog_grid_mean(state)
     aux_gm = TC.center_aux_grid_mean(state)
+    p_c_surf = aux_gm.p
+    p_f_surf = aux_gm_f.p
     u_gm_surf = prog_gm.u[kc_surf]
     v_gm_surf = prog_gm.v[kc_surf]
     q_tot_gm_surf = aux_gm.q_tot[kc_surf]
@@ -219,13 +222,13 @@ function get_surface(
 
     phase_part = TD.PhasePartition(q_tot_gm_surf, 0.0, 0.0)
     pvg = TD.saturation_vapor_pressure(param_set, TD.PhaseEquil, Tsurface)
-    qsurface = TD.q_vap_saturation_from_density(param_set, Tsurface, ρ0_f_surf, pvg)
-    θ_star = TD.liquid_ice_pottemp_given_pressure(param_set, Tsurface, p0_f_surf, phase_part)
+    qsurface = TD.q_vap_saturation_from_density(param_set, Tsurface, ρ_f_surf, pvg)
+    θ_star = TD.liquid_ice_pottemp_given_pressure(param_set, Tsurface, p_f_surf, phase_part)
 
     universal_func = UF.Businger()
     scheme = SF.FVScheme()
-    ts_sfc = TD.PhaseEquil_pθq(param_set, p0_f_surf, θ_star, qsurface)
-    ts_in = TD.PhaseEquil_pθq(param_set, p0_c_surf, θ_liq_ice_gm_surf, qsurface)
+    ts_sfc = TD.PhaseEquil_pθq(param_set, p_f_surf, θ_star, qsurface)
+    ts_in = TD.PhaseEquil_pθq(param_set, p_c_surf, θ_liq_ice_gm_surf, qsurface)
 
     u_sfc = SA.SVector{2, FT}(0, 0)
     u_in = SA.SVector{2, FT}(u_gm_surf, v_gm_surf)
diff --git a/driver/compute_diagnostics.jl b/driver/compute_diagnostics.jl
index 4a748b5603..989196d3ed 100644
--- a/driver/compute_diagnostics.jl
+++ b/driver/compute_diagnostics.jl
@@ -109,8 +109,6 @@ function compute_diagnostics!(
 ) where {D <: CC.Fields.FieldVector}
     FT = eltype(grid)
     N_up = TC.n_updrafts(edmf)
-    ρ0_c = TC.center_ref_state(state).ρ0
-    p0_c = TC.center_ref_state(state).p0
     aux_gm = TC.center_aux_grid_mean(state)
     aux_en = TC.center_aux_environment(state)
     aux_up = TC.center_aux_updrafts(state)
@@ -126,6 +124,8 @@ function compute_diagnostics!(
     prog_gm = TC.center_prog_grid_mean(state)
     diag_tc = center_diagnostics_turbconv(diagnostics)
     diag_tc_f = face_diagnostics_turbconv(diagnostics)
+    ρ_c = prog_gm.ρ
+    p_c = aux_gm.p
 
     diag_tc_svpc = svpc_diagnostics_turbconv(diagnostics)
     diag_svpc = svpc_diagnostics_grid_mean(diagnostics)
@@ -134,8 +134,8 @@ function compute_diagnostics!(
 
     @inbounds for k in TC.real_center_indices(grid)
         aux_gm.s[k] = TD.specific_entropy(param_set, ts_gm[k])
-        aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ0_c[k]
-        aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ0_c[k]
+        aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ_c[k]
+        aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k]
         aux_en.s[k] = TD.specific_entropy(param_set, ts_en[k])
     end
     @inbounds for k in TC.real_center_indices(grid)
@@ -148,7 +148,7 @@ function compute_diagnostics!(
                 end
                 ts_up = TC.thermo_state_pθq(
                     param_set,
-                    p0_c[k],
+                    p_c[k],
                     aux_up[i].θ_liq_ice[k],
                     aux_up[i].q_tot[k],
                     thermo_args...,
@@ -311,13 +311,13 @@ function compute_diagnostics!(
 
     TC.update_cloud_frac(edmf, grid, state)
 
-    diag_svpc.lwp_mean[cent] = sum(ρ0_c .* aux_gm.q_liq)
-    diag_svpc.iwp_mean[cent] = sum(ρ0_c .* aux_gm.q_ice)
-    diag_svpc.rwp_mean[cent] = sum(ρ0_c .* prog_pr.q_rai)
-    diag_svpc.swp_mean[cent] = sum(ρ0_c .* prog_pr.q_sno)
+    diag_svpc.lwp_mean[cent] = sum(ρ_c .* aux_gm.q_liq)
+    diag_svpc.iwp_mean[cent] = sum(ρ_c .* aux_gm.q_ice)
+    diag_svpc.rwp_mean[cent] = sum(ρ_c .* prog_pr.q_rai)
+    diag_svpc.swp_mean[cent] = sum(ρ_c .* prog_pr.q_sno)
 
-    diag_tc_svpc.env_lwp[cent] = sum(ρ0_c .* aux_en.q_liq .* aux_en.area)
-    diag_tc_svpc.env_iwp[cent] = sum(ρ0_c .* aux_en.q_ice .* aux_en.area)
+    diag_tc_svpc.env_lwp[cent] = sum(ρ_c .* aux_en.q_liq .* aux_en.area)
+    diag_tc_svpc.env_iwp[cent] = sum(ρ_c .* aux_en.q_ice .* aux_en.area)
 
     #TODO - change to rain rate that depends on rain model choice
 
@@ -325,13 +325,13 @@ function compute_diagnostics!(
     rho_cloud_liq = 1e3
     if (precip_model isa TC.CutoffPrecipitation)
         f =
-            (aux_en.qt_tendency_precip_formation .+ aux_bulk.qt_tendency_precip_formation) .* ρ0_c ./
-            TC.rho_cloud_liq .* 3.6 .* 1e6
+            (aux_en.qt_tendency_precip_formation .+ aux_bulk.qt_tendency_precip_formation) .* ρ_c ./ TC.rho_cloud_liq .*
+            3.6 .* 1e6
         diag_svpc.cutoff_precipitation_rate[cent] = sum(f)
     end
 
-    lwp = sum(i -> sum(ρ0_c .* aux_up[i].q_liq .* aux_up[i].area .* (aux_up[i].area .> 1e-3)), 1:N_up)
-    iwp = sum(i -> sum(ρ0_c .* aux_up[i].q_ice .* aux_up[i].area .* (aux_up[i].area .> 1e-3)), 1:N_up)
+    lwp = sum(i -> sum(ρ_c .* aux_up[i].q_liq .* aux_up[i].area .* (aux_up[i].area .> 1e-3)), 1:N_up)
+    iwp = sum(i -> sum(ρ_c .* aux_up[i].q_ice .* aux_up[i].area .* (aux_up[i].area .> 1e-3)), 1:N_up)
     plume_scale_height = map(1:N_up) do i
         TC.compute_plume_scale_height(grid, state, param_set, i)
     end
diff --git a/driver/dycore.jl b/driver/dycore.jl
index 3ad45dafe4..51f1fafe49 100644
--- a/driver/dycore.jl
+++ b/driver/dycore.jl
@@ -23,9 +23,6 @@ const APS = CLIMAParameters.AbstractEarthParameterSet
 
 ##### Auxiliary fields
 
-# Face & Center
-aux_vars_ref_state(FT) = (; ref_state = (ρ0 = FT(0), α0 = FT(0), p0 = FT(0)))
-
 # Center only
 cent_aux_vars_gm_moisture(FT, ::TC.NonEquilibriumMoisture) = (;
     ∇q_liq_gm = FT(0),
@@ -80,9 +77,9 @@ cent_aux_vars_gm(FT, edmf) = (;
     Ri = FT(0),
     θ_liq_ice = FT(0),
     q_tot = FT(0),
+    p = FT(0),
 )
-cent_aux_vars(FT, edmf) =
-    (; aux_vars_ref_state(FT)..., cent_aux_vars_gm(FT, edmf)..., TC.cent_aux_vars_edmf(FT, edmf)...)
+cent_aux_vars(FT, edmf) = (; cent_aux_vars_gm(FT, edmf)..., TC.cent_aux_vars_edmf(FT, edmf)...)
 
 # Face only
 face_aux_vars_gm_moisture(FT, ::TC.NonEquilibriumMoisture) = (; sgs_flux_q_liq = FT(0), sgs_flux_q_ice = FT(0))
@@ -97,9 +94,10 @@ face_aux_vars_gm(FT, edmf) = (;
     face_aux_vars_gm_moisture(FT, edmf.moisture_model)...,
     sgs_flux_u = FT(0),
     sgs_flux_v = FT(0),
+    p = FT(0),
+    ρ = FT(0),
 )
-face_aux_vars(FT, edmf) =
-    (; aux_vars_ref_state(FT)..., face_aux_vars_gm(FT, edmf)..., TC.face_aux_vars_edmf(FT, edmf)...)
+face_aux_vars(FT, edmf) = (; face_aux_vars_gm(FT, edmf)..., TC.face_aux_vars_edmf(FT, edmf)...)
 
 ##### Diagnostic fields
 
@@ -138,6 +136,7 @@ cent_prognostic_vars(::Type{FT}, local_geometry, edmf) where {FT} =
 cent_prognostic_vars_gm_moisture(::Type{FT}, ::TC.NonEquilibriumMoisture) where {FT} = (; q_liq = FT(0), q_ice = FT(0))
 cent_prognostic_vars_gm_moisture(::Type{FT}, ::TC.EquilibriumMoisture) where {FT} = NamedTuple()
 cent_prognostic_vars_gm(::Type{FT}, local_geometry, edmf) where {FT} = (;
+    ρ = FT(0),
     u = FT(0),
     v = FT(0),
     ρθ_liq_ice = FT(0),
@@ -185,12 +184,13 @@ The reference profiles, given
 function compute_ref_state!(state, grid::TC.Grid, param_set::PS; ts_g) where {PS}
 
     FT = eltype(grid)
-    p0_c = TC.center_ref_state(state).p0
-    ρ0_c = TC.center_ref_state(state).ρ0
-    α0_c = TC.center_ref_state(state).α0
-    p0_f = TC.face_ref_state(state).p0
-    ρ0_f = TC.face_ref_state(state).ρ0
-    α0_f = TC.face_ref_state(state).α0
+    aux_gm = TC.center_aux_grid_mean(state)
+    aux_gm_f = TC.face_aux_grid_mean(state)
+    prog_gm = TC.center_prog_grid_mean(state)
+    p_c = aux_gm.p
+    ρ_c = prog_gm.ρ
+    p_f = aux_gm_f.p
+    ρ_f = aux_gm_f.ρ
 
     qtg = TD.total_specific_humidity(param_set, ts_g)
     θ_liq_ice_g = TD.liquid_ice_pottemp(param_set, ts_g)
@@ -216,25 +216,22 @@ function compute_ref_state!(state, grid::TC.Grid, param_set::PS; ts_g) where {PS
     prob = ODE.ODEProblem(rhs, logp, z_span)
     sol = ODE.solve(prob, ODE.Tsit5(), reltol = 1e-12, abstol = 1e-12)
 
-    parent(p0_f) .= sol.(vec(grid.zf))
-    parent(p0_c) .= sol.(vec(grid.zc))
+    parent(p_f) .= sol.(vec(grid.zf))
+    parent(p_c) .= sol.(vec(grid.zc))
 
-    p0_f .= exp.(p0_f)
-    p0_c .= exp.(p0_c)
+    p_f .= exp.(p_f)
+    p_c .= exp.(p_c)
 
     # Compute reference state thermodynamic profiles
     @inbounds for k in TC.real_center_indices(grid)
-        ts = TD.PhaseEquil_pθq(param_set, p0_c[k], θ_liq_ice_g, qtg)
-        α0_c[k] = TD.specific_volume(param_set, ts)
+        ts = TD.PhaseEquil_pθq(param_set, p_c[k], θ_liq_ice_g, qtg)
+        ρ_c[k] = TD.air_density(param_set, ts)
     end
 
     @inbounds for k in TC.real_face_indices(grid)
-        ts = TD.PhaseEquil_pθq(param_set, p0_f[k], θ_liq_ice_g, qtg)
-        α0_f[k] = TD.specific_volume(param_set, ts)
+        ts = TD.PhaseEquil_pθq(param_set, p_f[k], θ_liq_ice_g, qtg)
+        ρ_f[k] = TD.air_density(param_set, ts)
     end
-
-    ρ0_f .= 1 ./ α0_f
-    ρ0_c .= 1 ./ α0_c
     return nothing
 end
 
@@ -242,11 +239,11 @@ function set_thermo_state!(state, grid, moisture_model, param_set)
     ts_gm = TC.center_aux_grid_mean(state).ts
     prog_gm = TC.center_prog_grid_mean(state)
     aux_gm = TC.center_aux_grid_mean(state)
-    p0_c = TC.center_ref_state(state).p0
-    ρ0_c = TC.center_ref_state(state).ρ0
+    p_c = aux_gm.p
+    ρ_c = prog_gm.ρ
     @inbounds for k in TC.real_center_indices(grid)
-        aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ0_c[k]
-        aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ0_c[k]
+        aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ_c[k]
+        aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k]
         thermo_args = if moisture_model isa TC.EquilibriumMoisture
             ()
         elseif moisture_model isa TC.NonEquilibriumMoisture
@@ -254,25 +251,25 @@ function set_thermo_state!(state, grid, moisture_model, param_set)
         else
             error("Something went wrong. The moisture_model options are equilibrium or nonequilibrium")
         end
-        ts_gm[k] = TC.thermo_state_pθq(param_set, p0_c[k], aux_gm.θ_liq_ice[k], aux_gm.q_tot[k], thermo_args...)
+        ts_gm[k] = TC.thermo_state_pθq(param_set, p_c[k], aux_gm.θ_liq_ice[k], aux_gm.q_tot[k], thermo_args...)
     end
     return nothing
 end
 
 function assign_thermo_aux!(state, grid, moisture_model, param_set)
-    ρ0_c = TC.center_ref_state(state).ρ0
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
     ts_gm = TC.center_aux_grid_mean(state).ts
+    ρ_c = prog_gm.ρ
     @inbounds for k in TC.real_center_indices(grid)
-        aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ0_c[k]
-        aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ0_c[k]
+        aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ_c[k]
+        aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k]
         ts = ts_gm[k]
         aux_gm.q_liq[k] = TD.liquid_specific_humidity(param_set, ts)
         aux_gm.q_ice[k] = TD.ice_specific_humidity(param_set, ts)
         aux_gm.T[k] = TD.air_temperature(param_set, ts)
         ρ = TD.air_density(param_set, ts)
-        aux_gm.buoy[k] = TC.buoyancy_c(param_set, ρ0_c[k], ρ)
+        aux_gm.buoy[k] = TC.buoyancy_c(param_set, ρ_c[k], ρ)
         aux_gm.RH[k] = TD.relative_humidity(param_set, ts)
     end
     return
@@ -367,10 +364,9 @@ function compute_gm_tendencies!(
     aux_en_f = TC.face_aux_environment(state)
     aux_up = TC.center_aux_updrafts(state)
     aux_bulk = TC.center_aux_bulk(state)
-    ρ0_f = TC.face_ref_state(state).ρ0
-    p0_c = TC.center_ref_state(state).p0
-    α0_c = TC.center_ref_state(state).α0
-    ρ0_c = TC.center_ref_state(state).ρ0
+    ρ_f = aux_gm_f.ρ
+    p_c = aux_gm.p
+    ρ_c = prog_gm.ρ
     aux_tc = TC.center_aux_turbconv(state)
     ts_gm = TC.center_aux_grid_mean(state).ts
 
@@ -404,20 +400,20 @@ function compute_gm_tendencies!(
             tendencies_gm.v[k] += force.coriolis_param * (aux_gm.ug[k] - prog_gm.u[k])
         end
         # LS Subsidence
-        tendencies_gm.ρθ_liq_ice[k] -= ρ0_c[k] * ∇θ_liq_ice_gm[k] * aux_gm.subsidence[k]
-        tendencies_gm.ρq_tot[k] -= ρ0_c[k] * ∇q_tot_gm[k] * aux_gm.subsidence[k]
+        tendencies_gm.ρθ_liq_ice[k] -= ρ_c[k] * ∇θ_liq_ice_gm[k] * aux_gm.subsidence[k]
+        tendencies_gm.ρq_tot[k] -= ρ_c[k] * ∇q_tot_gm[k] * aux_gm.subsidence[k]
         if edmf.moisture_model isa TC.NonEquilibriumMoisture
             tendencies_gm.q_liq[k] -= ∇q_liq_gm[k] * aux_gm.subsidence[k]
             tendencies_gm.q_ice[k] -= ∇q_ice_gm[k] * aux_gm.subsidence[k]
         end
         # Radiation
         if TC.rad_type(radiation) <: Union{TC.RadiationDYCOMS_RF01, TC.RadiationLES, TC.RadiationTRMM_LBA}
-            tendencies_gm.ρθ_liq_ice[k] += ρ0_c[k] * aux_gm.dTdt_rad[k] / Π
+            tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * aux_gm.dTdt_rad[k] / Π
         end
         # LS advection
-        tendencies_gm.ρq_tot[k] += ρ0_c[k] * aux_gm.dqtdt_hadv[k]
+        tendencies_gm.ρq_tot[k] += ρ_c[k] * aux_gm.dqtdt_hadv[k]
         if !(TC.force_type(force) <: TC.ForcingDYCOMS_RF01)
-            tendencies_gm.ρθ_liq_ice[k] += ρ0_c[k] * aux_gm.dTdt_hadv[k] / Π
+            tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * aux_gm.dTdt_hadv[k] / Π
         end
         if edmf.moisture_model isa TC.NonEquilibriumMoisture
             tendencies_gm.q_liq[k] += aux_gm.dqldt[k]
@@ -439,8 +435,8 @@ function compute_gm_tendencies!(
                 gm_q_ice_nudge_k = Γᵣ * (aux_gm.qi_nudge[k] - prog_gm.q_ice[k])
             end
 
-            tendencies_gm.ρθ_liq_ice[k] += ρ0_c[k] * (gm_H_nudge_k + H_fluc)
-            tendencies_gm.ρq_tot[k] += ρ0_c[k] * (gm_q_tot_nudge_k + aux_gm.dqtdt_fluc[k])
+            tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * (gm_H_nudge_k + H_fluc)
+            tendencies_gm.ρq_tot[k] += ρ_c[k] * (gm_q_tot_nudge_k + aux_gm.dqtdt_fluc[k])
             tendencies_gm.u[k] += gm_U_nudge_k
             tendencies_gm.v[k] += gm_V_nudge_k
             if edmf.moisture_model isa TC.NonEquilibriumMoisture
@@ -451,13 +447,13 @@ function compute_gm_tendencies!(
 
         # Apply precipitation tendencies
         tendencies_gm.ρq_tot[k] +=
-            ρ0_c[k] * (
+            ρ_c[k] * (
                 aux_bulk.qt_tendency_precip_formation[k] +
                 aux_en.qt_tendency_precip_formation[k] +
                 aux_tc.qt_tendency_precip_sinks[k]
             )
         tendencies_gm.ρθ_liq_ice[k] +=
-            ρ0_c[k] * (
+            ρ_c[k] * (
                 aux_bulk.θ_liq_ice_tendency_precip_formation[k] +
                 aux_en.θ_liq_ice_tendency_precip_formation[k] +
                 aux_tc.θ_liq_ice_tendency_precip_sinks[k]
@@ -489,8 +485,8 @@ function compute_gm_tendencies!(
 
     @. tends_ρθ_liq_ice += -∇sgs(wvec(sgs_flux_θ_liq_ice))
     @. tends_ρq_tot += -∇sgs(wvec(sgs_flux_q_tot))
-    @. tends_u += -α0_c * ∇sgs(wvec(sgs_flux_u))
-    @. tends_v += -α0_c * ∇sgs(wvec(sgs_flux_v))
+    @. tends_u += -∇sgs(wvec(sgs_flux_u)) / ρ_c
+    @. tends_v += -∇sgs(wvec(sgs_flux_v)) / ρ_c
 
     if edmf.moisture_model isa TC.NonEquilibriumMoisture
         sgs_flux_q_liq = aux_gm_f.sgs_flux_q_liq
@@ -499,8 +495,8 @@ function compute_gm_tendencies!(
         sgs_flux_q_ice[kf_surf] = surf.ρq_ice_flux
         tends_q_liq = tendencies_gm.q_liq
         tends_q_ice = tendencies_gm.q_ice
-        @. tends_q_liq += -α0_c * ∇sgs(wvec(sgs_flux_q_liq))
-        @. tends_q_ice += -α0_c * ∇sgs(wvec(sgs_flux_q_ice))
+        @. tends_q_liq += -∇sgs(wvec(sgs_flux_q_liq)) / ρ_c
+        @. tends_q_ice += -∇sgs(wvec(sgs_flux_q_ice)) / ρ_c
     end
 
     return nothing
diff --git a/driver/initial_conditions.jl b/driver/initial_conditions.jl
index 00fa02df8b..6eaf5aa7c2 100644
--- a/driver/initial_conditions.jl
+++ b/driver/initial_conditions.jl
@@ -14,7 +14,7 @@ function initialize_edmf(edmf::TC.EDMFModel, grid::TC.Grid, state::TC.State, cas
     aux_gm = TC.center_aux_grid_mean(state)
     ts_gm = aux_gm.ts
     prog_gm = TC.center_prog_grid_mean(state)
-    p0_c = TC.center_ref_state(state).p0
+    p_c = prog_gm.ρ
     parent(aux_tc.prandtl_nvec) .= edmf.prandtl_number
     @inbounds for k in TC.real_center_indices(grid)
         aux_gm.θ_virt[k] = TD.virtual_pottemp(param_set, ts_gm[k])
@@ -35,7 +35,8 @@ function initialize_covariance(edmf::TC.EDMFModel, grid::TC.Grid, state::TC.Stat
     aux_gm = TC.center_aux_grid_mean(state)
     prog_en = TC.center_prog_environment(state)
     aux_en = TC.center_aux_environment(state)
-    ρ0_c = TC.center_ref_state(state).ρ0
+    prog_gm = TC.center_prog_grid_mean(state)
+    ρ_c = prog_gm.ρ
     aux_bulk = TC.center_aux_bulk(state)
     ae = 1 .- aux_bulk.area # area of environment
 
@@ -44,11 +45,11 @@ function initialize_covariance(edmf::TC.EDMFModel, grid::TC.Grid, state::TC.Stat
     aux_en.QTvar .= aux_gm.QTvar
     aux_en.HQTcov .= aux_gm.HQTcov
 
-    prog_en.ρatke .= aux_en.tke .* ρ0_c .* ae
+    prog_en.ρatke .= aux_en.tke .* ρ_c .* ae
     if edmf.thermo_covariance_model isa TC.PrognosticThermoCovariances
-        prog_en.ρaHvar .= aux_gm.Hvar .* ρ0_c .* ae
-        prog_en.ρaQTvar .= aux_gm.QTvar .* ρ0_c .* ae
-        prog_en.ρaHQTcov .= aux_gm.HQTcov .* ρ0_c .* ae
+        prog_en.ρaHvar .= aux_gm.Hvar .* ρ_c .* ae
+        prog_en.ρaQTvar .= aux_gm.QTvar .* ρ_c .* ae
+        prog_en.ρaHQTcov .= aux_gm.HQTcov .* ρ_c .* ae
     end
     return
 end
@@ -63,7 +64,7 @@ function initialize_updrafts(edmf, grid, state, surf)
     aux_gm = TC.center_aux_grid_mean(state)
     prog_up = TC.center_prog_updrafts(state)
     prog_up_f = TC.face_prog_updrafts(state)
-    ρ0_c = TC.center_ref_state(state).ρ0
+    ρ_c = prog_gm.ρ
     @inbounds for i in 1:N_up
         @inbounds for k in TC.real_face_indices(grid)
             aux_up_f[i].w[k] = 0
@@ -91,7 +92,7 @@ function initialize_updrafts(edmf, grid, state, surf)
 
         a_surf = TC.area_surface_bc(surf, edmf, i)
         aux_up[i].area[kc_surf] = a_surf
-        prog_up[i].ρarea[kc_surf] = ρ0_c[kc_surf] * a_surf
+        prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_surf
     end
     return
 end
@@ -104,10 +105,11 @@ function initialize_updrafts_DryBubble(edmf, grid, state)
     aux_up = TC.center_aux_updrafts(state)
     aux_up_f = TC.face_aux_updrafts(state)
     aux_gm = TC.center_aux_grid_mean(state)
+    aux_gm_f = TC.face_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
     prog_up = TC.center_prog_updrafts(state)
-    ρ_0_c = TC.center_ref_state(state).ρ0
-    ρ_0_f = TC.face_ref_state(state).ρ0
+    ρ_0_c = prog_gm.ρ
+    ρ_0_f = aux_gm_f.ρ
     N_up = TC.n_updrafts(edmf)
     FT = eltype(grid)
     z_in = APL.DryBubble_updrafts_z(FT)
diff --git a/driver/main.jl b/driver/main.jl
index ae641a88b4..05b510746d 100644
--- a/driver/main.jl
+++ b/driver/main.jl
@@ -138,12 +138,10 @@ function Simulation1d(namelist)
     if !skip_io
         NC.Dataset(Stats.nc_filename, "a") do ds
             group = "reference"
-            add_write_field(ds, "ρ0_f", vec(TC.face_ref_state(state).ρ0), group, ("zf",))
-            add_write_field(ds, "ρ0_c", vec(TC.center_ref_state(state).ρ0), group, ("zc",))
-            add_write_field(ds, "p0_f", vec(TC.face_ref_state(state).p0), group, ("zf",))
-            add_write_field(ds, "p0_c", vec(TC.center_ref_state(state).p0), group, ("zc",))
-            add_write_field(ds, "α0_f", vec(TC.face_ref_state(state).α0), group, ("zf",))
-            add_write_field(ds, "α0_c", vec(TC.center_ref_state(state).α0), group, ("zc",))
+            add_write_field(ds, "ρ_f", vec(TC.face_aux_grid_mean(state).ρ), group, ("zf",))
+            add_write_field(ds, "ρ_c", vec(TC.center_prog_grid_mean(state).ρ), group, ("zc",))
+            add_write_field(ds, "p_f", vec(TC.face_aux_grid_mean(state).p), group, ("zf",))
+            add_write_field(ds, "p_c", vec(TC.center_aux_grid_mean(state).p), group, ("zc",))
         end
     end
 
diff --git a/post_processing/mse_tables.jl b/post_processing/mse_tables.jl
index aafa6299d4..d8b8723e4a 100644
--- a/post_processing/mse_tables.jl
+++ b/post_processing/mse_tables.jl
@@ -5,182 +5,173 @@
 all_best_mse = OrderedCollections.OrderedDict()
 #
 all_best_mse["ARM_SGP"] = OrderedCollections.OrderedDict()
-all_best_mse["ARM_SGP"]["qt_mean"] = 0.252123418092592
-all_best_mse["ARM_SGP"]["updraft_area"] = 327.0886467601207
-all_best_mse["ARM_SGP"]["updraft_w"] = 153.60644689673893
-all_best_mse["ARM_SGP"]["updraft_qt"] = 30.636706048970172
-all_best_mse["ARM_SGP"]["updraft_thetal"] = 172.01251579139668
+all_best_mse["ARM_SGP"]["qt_mean"] = 0.2521234075111047
+all_best_mse["ARM_SGP"]["updraft_area"] = 327.08864677891694
+all_best_mse["ARM_SGP"]["updraft_w"] = 153.60643978772632
+all_best_mse["ARM_SGP"]["updraft_qt"] = 30.63670615757075
+all_best_mse["ARM_SGP"]["updraft_thetal"] = 172.0125157915087
 all_best_mse["ARM_SGP"]["u_mean"] = 1.3235797273549681e-5
-all_best_mse["ARM_SGP"]["tke_mean"] = 1101.5753412241338
-all_best_mse["ARM_SGP"]["temperature_mean"] = 0.00012552649995235397
-all_best_mse["ARM_SGP"]["ql_mean"] = 211.93337065688647
+all_best_mse["ARM_SGP"]["tke_mean"] = 1101.5753414987716
+all_best_mse["ARM_SGP"]["temperature_mean"] = 0.00012552650044046315
+all_best_mse["ARM_SGP"]["ql_mean"] = 211.9333773104895
 all_best_mse["ARM_SGP"]["qi_mean"] = "NA"
-all_best_mse["ARM_SGP"]["thetal_mean"] = 0.0001181997362832056
-all_best_mse["ARM_SGP"]["Hvar_mean"] = 2195.6266107140377
-all_best_mse["ARM_SGP"]["QTvar_mean"] = 1309.9340334273952
+all_best_mse["ARM_SGP"]["thetal_mean"] = 0.00011819973689302238
+all_best_mse["ARM_SGP"]["Hvar_mean"] = 2195.6037821937093
+all_best_mse["ARM_SGP"]["QTvar_mean"] = 1309.9233327883517
 #
 all_best_mse["Bomex"] = OrderedCollections.OrderedDict()
-all_best_mse["Bomex"]["qt_mean"] = 0.10118029045976014
-all_best_mse["Bomex"]["updraft_area"] = 127.60478829277722
-all_best_mse["Bomex"]["updraft_w"] = 18.17115270923453
-all_best_mse["Bomex"]["updraft_qt"] = 6.876484718878337
-all_best_mse["Bomex"]["updraft_thetal"] = 69.78940065177694
-all_best_mse["Bomex"]["v_mean"] = 62.75101975458986
-all_best_mse["Bomex"]["u_mean"] = 0.25586605778026955
-all_best_mse["Bomex"]["tke_mean"] = 72.80350572399242
-all_best_mse["Bomex"]["temperature_mean"] = 4.138849941184466e-5
-all_best_mse["Bomex"]["ql_mean"] = 9.497944951749115
+all_best_mse["Bomex"]["qt_mean"] = 0.10178076257542494
+all_best_mse["Bomex"]["updraft_area"] = 127.61793250808256
+all_best_mse["Bomex"]["updraft_w"] = 17.947263820879698
+all_best_mse["Bomex"]["updraft_qt"] = 7.132188703920948
+all_best_mse["Bomex"]["updraft_thetal"] = 69.79490167272438
+all_best_mse["Bomex"]["v_mean"] = 62.767660503421865
+all_best_mse["Bomex"]["u_mean"] = 0.25600294047426886
+all_best_mse["Bomex"]["tke_mean"] = 72.73304427841286
+all_best_mse["Bomex"]["temperature_mean"] = 4.158889742893274e-5
+all_best_mse["Bomex"]["ql_mean"] = 9.385139174268497
 all_best_mse["Bomex"]["qi_mean"] = "NA"
-all_best_mse["Bomex"]["thetal_mean"] = 4.195534239162112e-5
-all_best_mse["Bomex"]["Hvar_mean"] = 842.4281340649104
-all_best_mse["Bomex"]["QTvar_mean"] = 313.0900983916904
+all_best_mse["Bomex"]["thetal_mean"] = 4.215083961811883e-5
+all_best_mse["Bomex"]["Hvar_mean"] = 2574.2571663224408
+all_best_mse["Bomex"]["QTvar_mean"] = 918.4100774052899
 #
 all_best_mse["DryBubble"] = OrderedCollections.OrderedDict()
-all_best_mse["DryBubble"]["updraft_area"] = 0.0
-all_best_mse["DryBubble"]["updraft_w"] = 0.0
-all_best_mse["DryBubble"]["updraft_thetal"] = 0.0
-all_best_mse["DryBubble"]["u_mean"] = 0.0
-all_best_mse["DryBubble"]["tke_mean"] = 0.0
-all_best_mse["DryBubble"]["temperature_mean"] = 0.0
-all_best_mse["DryBubble"]["thetal_mean"] = 0.0
-all_best_mse["DryBubble"]["Hvar_mean"] = 78.03740725881899
+all_best_mse["DryBubble"]["updraft_area"] = 2.15111668772059e-22
+all_best_mse["DryBubble"]["updraft_w"] = 6.700778035066527e-22
+all_best_mse["DryBubble"]["updraft_thetal"] = 5.597000941990776e-29
+all_best_mse["DryBubble"]["u_mean"] = 6.339213832488136e-22
+all_best_mse["DryBubble"]["tke_mean"] = 5.951184326624882e-22
+all_best_mse["DryBubble"]["temperature_mean"] = 2.1690026983354042e-29
+all_best_mse["DryBubble"]["thetal_mean"] = 2.1681633083031133e-29
+all_best_mse["DryBubble"]["Hvar_mean"] = 1.7445289662960883e-21
 #
 all_best_mse["DYCOMS_RF01"] = OrderedCollections.OrderedDict()
-all_best_mse["DYCOMS_RF01"]["qt_mean"] = 0.018639912355757106
-all_best_mse["DYCOMS_RF01"]["ql_mean"] = 5.334510672241666
-all_best_mse["DYCOMS_RF01"]["updraft_area"] = 29.951787390089795
-all_best_mse["DYCOMS_RF01"]["updraft_w"] = 4.2898983129371375
-all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 2.141607127622009
-all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 46.18824292060282
-all_best_mse["DYCOMS_RF01"]["v_mean"] = 0.0027609152547554807
-all_best_mse["DYCOMS_RF01"]["u_mean"] = 0.08875417460646276
-all_best_mse["DYCOMS_RF01"]["tke_mean"] = 21.124989822145707
-all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 7.13572405899015e-5
-all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 7.174079092534151e-5
-all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1287.4864641621884
-all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 513.0844832478505
+all_best_mse["DYCOMS_RF01"]["qt_mean"] = 0.018639850553854188
+all_best_mse["DYCOMS_RF01"]["ql_mean"] = 5.3343483702028145
+all_best_mse["DYCOMS_RF01"]["updraft_area"] = 29.951759526804484
+all_best_mse["DYCOMS_RF01"]["updraft_w"] = 4.289896238855623
+all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 2.141605193888237
+all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 46.18824291535536
+all_best_mse["DYCOMS_RF01"]["v_mean"] = 0.0027609455255945695
+all_best_mse["DYCOMS_RF01"]["u_mean"] = 0.08875426568437153
+all_best_mse["DYCOMS_RF01"]["tke_mean"] = 21.125128725885098
+all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 7.135752046133808e-5
+all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 7.174107117755792e-5
+all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1287.4864356667742
+all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 513.0844941712962
 #
 all_best_mse["DYCOMS_RF02"] = OrderedCollections.OrderedDict()
-all_best_mse["DYCOMS_RF02"]["qt_mean"] = 0.0646435532921293
-all_best_mse["DYCOMS_RF02"]["ql_mean"] = 3.649974400406772
-all_best_mse["DYCOMS_RF02"]["qr_mean"] = 7.627111872119355
-all_best_mse["DYCOMS_RF02"]["updraft_area"] = 27.75311549797335
-all_best_mse["DYCOMS_RF02"]["updraft_w"] = 6.7101305450505135
-all_best_mse["DYCOMS_RF02"]["updraft_qt"] = 5.062597050045206
-all_best_mse["DYCOMS_RF02"]["updraft_thetal"] = 40.54966454415253
-all_best_mse["DYCOMS_RF02"]["v_mean"] = 42.87028853678261
-all_best_mse["DYCOMS_RF02"]["u_mean"] = 20.17393420636257
-all_best_mse["DYCOMS_RF02"]["tke_mean"] = 21.409525043454675
-all_best_mse["DYCOMS_RF02"]["temperature_mean"] = 2.7353260020219525e-5
-all_best_mse["DYCOMS_RF02"]["thetal_mean"] = 3.0094035941549328e-5
-all_best_mse["DYCOMS_RF02"]["Hvar_mean"] = 1465.0320663329937
-all_best_mse["DYCOMS_RF02"]["QTvar_mean"] = 514.6532270001834
-#
-all_best_mse["GABLS"] = OrderedCollections.OrderedDict()
-all_best_mse["GABLS"]["updraft_thetal"] = 0.0
-all_best_mse["GABLS"]["v_mean"] = 0.0
-all_best_mse["GABLS"]["u_mean"] = 0.0
-all_best_mse["GABLS"]["tke_mean"] = 0.0
-all_best_mse["GABLS"]["temperature_mean"] = 0.0
-all_best_mse["GABLS"]["thetal_mean"] = 0.0
-all_best_mse["GABLS"]["Hvar_mean"] = 0.18056270269686944
+all_best_mse["DYCOMS_RF02"]["qt_mean"] = 0.06464355295106711
+all_best_mse["DYCOMS_RF02"]["ql_mean"] = 3.6499743731380554
+all_best_mse["DYCOMS_RF02"]["qr_mean"] = 7.627111527506284
+all_best_mse["DYCOMS_RF02"]["updraft_area"] = 27.753115447091403
+all_best_mse["DYCOMS_RF02"]["updraft_w"] = 6.710130415020492
+all_best_mse["DYCOMS_RF02"]["updraft_qt"] = 5.062597048091387
+all_best_mse["DYCOMS_RF02"]["updraft_thetal"] = 40.54966454414522
+all_best_mse["DYCOMS_RF02"]["v_mean"] = 42.870288545493935
+all_best_mse["DYCOMS_RF02"]["u_mean"] = 20.17393420753407
+all_best_mse["DYCOMS_RF02"]["tke_mean"] = 21.409525076941335
+all_best_mse["DYCOMS_RF02"]["temperature_mean"] = 2.735325975217763e-5
+all_best_mse["DYCOMS_RF02"]["thetal_mean"] = 3.0094035826352356e-5
+all_best_mse["DYCOMS_RF02"]["Hvar_mean"] = 1465.0320663302696
+all_best_mse["DYCOMS_RF02"]["QTvar_mean"] = 514.6532269648812
 #
 all_best_mse["GATE_III"] = OrderedCollections.OrderedDict()
-all_best_mse["GATE_III"]["qt_mean"] = 1.3588998119672754e-11
-all_best_mse["GATE_III"]["updraft_area"] = 1.643225468965924e-8
-all_best_mse["GATE_III"]["updraft_w"] = 6.968244880700734e-12
-all_best_mse["GATE_III"]["updraft_qt"] = 5.754070370979365e-13
-all_best_mse["GATE_III"]["updraft_thetal"] = 1.20844940388337e-16
-all_best_mse["GATE_III"]["u_mean"] = 2.697130991236307e-14
-all_best_mse["GATE_III"]["tke_mean"] = 3.785033065600374e-12
-all_best_mse["GATE_III"]["temperature_mean"] = 9.97144496437381e-16
-all_best_mse["GATE_III"]["ql_mean"] = 2.2266770145966954e-10
-all_best_mse["GATE_III"]["qi_mean"] = 1.1842087295930235e-13
-all_best_mse["GATE_III"]["qr_mean"] = 7.065296147098808e-12
-all_best_mse["GATE_III"]["qs_mean"] = 8.685885117244753e-15
-all_best_mse["GATE_III"]["thetal_mean"] = 5.178875959127746e-16
-all_best_mse["GATE_III"]["Hvar_mean"] = 0.054763947080733774
-all_best_mse["GATE_III"]["QTvar_mean"] = 0.11541911365297315
+all_best_mse["GATE_III"]["qt_mean"] = 2.751763072097516e-17
+all_best_mse["GATE_III"]["updraft_area"] = 1.0067003077503561e-12
+all_best_mse["GATE_III"]["updraft_w"] = 1.9866903035375464e-15
+all_best_mse["GATE_III"]["updraft_qt"] = 9.943454893675091e-17
+all_best_mse["GATE_III"]["updraft_thetal"] = 3.976136212397589e-20
+all_best_mse["GATE_III"]["u_mean"] = 1.7786046247604032e-19
+all_best_mse["GATE_III"]["tke_mean"] = 5.4637034538767e-15
+all_best_mse["GATE_III"]["temperature_mean"] = 3.4386468585602514e-20
+all_best_mse["GATE_III"]["ql_mean"] = 4.2061041281769915e-14
+all_best_mse["GATE_III"]["qi_mean"] = 8.712528060056637e-17
+all_best_mse["GATE_III"]["qr_mean"] = 5.040280688212245e-13
+all_best_mse["GATE_III"]["qs_mean"] = 5.8900642113111296e-15
+all_best_mse["GATE_III"]["thetal_mean"] = 1.8476177043943748e-20
+all_best_mse["GATE_III"]["Hvar_mean"] = 5.833665028073461e-14
+all_best_mse["GATE_III"]["QTvar_mean"] = 3.666722357065889e-13
 #
 all_best_mse["life_cycle_Tan2018"] = OrderedCollections.OrderedDict()
-all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["v_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["u_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["temperature_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["thetal_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 6517.586756676236
-all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 4884.704601518459
+all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 1.1930778350948134e-12
+all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 4.856320550271249e-11
+all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 4.700330857032967e-12
+all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 3.4780115808939696e-10
+all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 3.4866062424258736e-13
+all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 1.9017271421464622e-16
+all_best_mse["life_cycle_Tan2018"]["v_mean"] = 8.440183468620281e-13
+all_best_mse["life_cycle_Tan2018"]["u_mean"] = 1.6361855585013259e-15
+all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 1.3261257800624283e-12
+all_best_mse["life_cycle_Tan2018"]["temperature_mean"] = 5.746895768249593e-16
+all_best_mse["life_cycle_Tan2018"]["thetal_mean"] = 5.810925162561764e-16
+all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 0.000215944180887789
+all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 0.0001195316954772662
 #
 all_best_mse["Nieuwstadt"] = OrderedCollections.OrderedDict()
-all_best_mse["Nieuwstadt"]["updraft_area"] = 99.95518312641288
-all_best_mse["Nieuwstadt"]["updraft_w"] = 14.234210189622093
-all_best_mse["Nieuwstadt"]["updraft_thetal"] = 117.60593544135058
-all_best_mse["Nieuwstadt"]["u_mean"] = 13.5533124537228
-all_best_mse["Nieuwstadt"]["tke_mean"] = 283.4527945078942
-all_best_mse["Nieuwstadt"]["temperature_mean"] = 1.1069166079595213e-5
-all_best_mse["Nieuwstadt"]["thetal_mean"] = 1.1163599024010392e-5
-all_best_mse["Nieuwstadt"]["Hvar_mean"] = 719.3877989748415
+all_best_mse["Nieuwstadt"]["updraft_area"] = 99.95518216478133
+all_best_mse["Nieuwstadt"]["updraft_w"] = 14.234210201962897
+all_best_mse["Nieuwstadt"]["updraft_thetal"] = 117.60593544133755
+all_best_mse["Nieuwstadt"]["u_mean"] = 13.553312453847973
+all_best_mse["Nieuwstadt"]["tke_mean"] = 283.4527946810255
+all_best_mse["Nieuwstadt"]["temperature_mean"] = 1.1069166056920466e-5
+all_best_mse["Nieuwstadt"]["thetal_mean"] = 1.116359900144678e-5
+all_best_mse["Nieuwstadt"]["Hvar_mean"] = 719.3877968155495
 #
 all_best_mse["Rico"] = OrderedCollections.OrderedDict()
-all_best_mse["Rico"]["qt_mean"] = 0.8722842742643365
-all_best_mse["Rico"]["updraft_area"] = 478.585114377104
-all_best_mse["Rico"]["updraft_w"] = 75.07963512673996
-all_best_mse["Rico"]["updraft_qt"] = 16.3639507812084
-all_best_mse["Rico"]["updraft_thetal"] = 133.8731196223532
-all_best_mse["Rico"]["v_mean"] = 0.5567836809940072
-all_best_mse["Rico"]["u_mean"] = 0.376846847076814
-all_best_mse["Rico"]["tke_mean"] = 144.5913986516131
-all_best_mse["Rico"]["temperature_mean"] = 0.0005403453265783186
-all_best_mse["Rico"]["ql_mean"] = 81.8434756745955
+all_best_mse["Rico"]["qt_mean"] = 0.8685256836011042
+all_best_mse["Rico"]["updraft_area"] = 478.56505133327573
+all_best_mse["Rico"]["updraft_w"] = 75.16164491943583
+all_best_mse["Rico"]["updraft_qt"] = 16.275399293359392
+all_best_mse["Rico"]["updraft_thetal"] = 133.87297518151618
+all_best_mse["Rico"]["v_mean"] = 0.5582957429052721
+all_best_mse["Rico"]["u_mean"] = 0.3770143595330914
+all_best_mse["Rico"]["tke_mean"] = 144.47757062599226
+all_best_mse["Rico"]["temperature_mean"] = 0.0005399337245447159
+all_best_mse["Rico"]["ql_mean"] = 81.64732868600886
 all_best_mse["Rico"]["qi_mean"] = "NA"
-all_best_mse["Rico"]["qr_mean"] = 751.7824809370619
-all_best_mse["Rico"]["thetal_mean"] = 0.0005280009571436345
-all_best_mse["Rico"]["Hvar_mean"] = 51393.014177146
-all_best_mse["Rico"]["QTvar_mean"] = 11394.994081496592
+all_best_mse["Rico"]["qr_mean"] = 751.8177494026789
+all_best_mse["Rico"]["thetal_mean"] = 0.0005275671713452303
+all_best_mse["Rico"]["Hvar_mean"] = 40691.28380130209
+all_best_mse["Rico"]["QTvar_mean"] = 9017.926388319349
 #
 all_best_mse["Soares"] = OrderedCollections.OrderedDict()
-all_best_mse["Soares"]["qt_mean"] = 0.14273283251496136
-all_best_mse["Soares"]["updraft_area"] = 97.54831237906166
-all_best_mse["Soares"]["updraft_w"] = 13.105743423539547
-all_best_mse["Soares"]["updraft_qt"] = 23.641103294226017
-all_best_mse["Soares"]["updraft_thetal"] = 65.72138701636425
-all_best_mse["Soares"]["u_mean"] = 93.89558704065772
-all_best_mse["Soares"]["tke_mean"] = 216.4518468909944
-all_best_mse["Soares"]["temperature_mean"] = 1.2381752395823234e-5
-all_best_mse["Soares"]["thetal_mean"] = 1.2064931346843553e-5
-all_best_mse["Soares"]["Hvar_mean"] = 680.9440520288947
+all_best_mse["Soares"]["qt_mean"] = 0.1427328325093554
+all_best_mse["Soares"]["updraft_area"] = 97.54831238810962
+all_best_mse["Soares"]["updraft_w"] = 13.105743423491914
+all_best_mse["Soares"]["updraft_qt"] = 23.64110329423007
+all_best_mse["Soares"]["updraft_thetal"] = 65.72138701636423
+all_best_mse["Soares"]["u_mean"] = 93.8955870406489
+all_best_mse["Soares"]["tke_mean"] = 216.4518468904579
+all_best_mse["Soares"]["temperature_mean"] = 1.2381752395697475e-5
+all_best_mse["Soares"]["thetal_mean"] = 1.2064931346715829e-5
+all_best_mse["Soares"]["Hvar_mean"] = 680.9440520288703
 #
 all_best_mse["TRMM_LBA"] = OrderedCollections.OrderedDict()
-all_best_mse["TRMM_LBA"]["qt_mean"] = 2.0998380118618782
-all_best_mse["TRMM_LBA"]["updraft_area"] = 1268.3808307624563
-all_best_mse["TRMM_LBA"]["updraft_w"] = 15532.477439550537
-all_best_mse["TRMM_LBA"]["updraft_qt"] = 285.1069426743044
-all_best_mse["TRMM_LBA"]["updraft_thetal"] = 490.0926293240085
-all_best_mse["TRMM_LBA"]["v_mean"] = 67.455540617838
-all_best_mse["TRMM_LBA"]["u_mean"] = 27.227766500480826
-all_best_mse["TRMM_LBA"]["tke_mean"] = 135664.71364592292
-all_best_mse["TRMM_LBA"]["temperature_mean"] = 0.000672844938693891
-all_best_mse["TRMM_LBA"]["ql_mean"] = 235473.98702276798
+all_best_mse["TRMM_LBA"]["qt_mean"] = 2.0998380118612485
+all_best_mse["TRMM_LBA"]["updraft_area"] = 1268.3808323837927
+all_best_mse["TRMM_LBA"]["updraft_w"] = 15532.477432102563
+all_best_mse["TRMM_LBA"]["updraft_qt"] = 285.1069426549674
+all_best_mse["TRMM_LBA"]["updraft_thetal"] = 490.09262932405204
+all_best_mse["TRMM_LBA"]["v_mean"] = 67.45554062489389
+all_best_mse["TRMM_LBA"]["u_mean"] = 27.22776649937904
+all_best_mse["TRMM_LBA"]["tke_mean"] = 135664.71363635483
+all_best_mse["TRMM_LBA"]["temperature_mean"] = 0.0006728449387071617
+all_best_mse["TRMM_LBA"]["ql_mean"] = 235473.98702879949
 all_best_mse["TRMM_LBA"]["qi_mean"] = "NA"
 all_best_mse["TRMM_LBA"]["qr_mean"] = "NA"
 all_best_mse["TRMM_LBA"]["qs_mean"] = "NA"
-all_best_mse["TRMM_LBA"]["thetal_mean"] = 0.0041402379026755035
-all_best_mse["TRMM_LBA"]["Hvar_mean"] = 430415.19984565926
-all_best_mse["TRMM_LBA"]["QTvar_mean"] = 6659.921910406574
+all_best_mse["TRMM_LBA"]["thetal_mean"] = 0.00414023790425752
+all_best_mse["TRMM_LBA"]["Hvar_mean"] = 430415.1998350132
+all_best_mse["TRMM_LBA"]["QTvar_mean"] = 6659.921911008241
 #
 all_best_mse["LES_driven_SCM"] = OrderedCollections.OrderedDict()
-all_best_mse["LES_driven_SCM"]["qt_mean"] = 3.565679133346325
-all_best_mse["LES_driven_SCM"]["v_mean"] = 1.3786051560755312
-all_best_mse["LES_driven_SCM"]["u_mean"] = 0.4764379333062847
-all_best_mse["LES_driven_SCM"]["temperature_mean"] = 0.0012830981248823014
-all_best_mse["LES_driven_SCM"]["ql_mean"] = 51738.50247731223
-all_best_mse["LES_driven_SCM"]["thetal_mean"] = 0.0015164180359003641
+all_best_mse["LES_driven_SCM"]["qt_mean"] = 3.5656791615746664
+all_best_mse["LES_driven_SCM"]["v_mean"] = 1.378605161974001
+all_best_mse["LES_driven_SCM"]["u_mean"] = 0.47643793378470234
+all_best_mse["LES_driven_SCM"]["temperature_mean"] = 0.0012830981258744414
+all_best_mse["LES_driven_SCM"]["ql_mean"] = 51738.51124851987
+all_best_mse["LES_driven_SCM"]["thetal_mean"] = 0.0015164180471682785
 #
 #################################
 #################################
diff --git a/src/EDMF_Environment.jl b/src/EDMF_Environment.jl
index 274855ea85..b541f563d5 100644
--- a/src/EDMF_Environment.jl
+++ b/src/EDMF_Environment.jl
@@ -9,11 +9,13 @@ function microphysics(
 )
 
     tendencies_pr = center_tendencies_precipitation(state)
-    p0_c = center_ref_state(state).p0
-    ρ0_c = center_ref_state(state).ρ0
     aux_en = center_aux_environment(state)
     prog_pr = center_prog_precipitation(state)
+    prog_gm = center_prog_grid_mean(state)
+    aux_gm = center_aux_grid_mean(state)
     ts_env = center_aux_environment(state).ts
+    p_c = aux_gm.p
+    ρ_c = prog_gm.ρ
     aux_en_sat = aux_en.sat
     aux_en_unsat = aux_en.unsat
 
@@ -27,7 +29,7 @@ function microphysics(
             prog_pr.q_rai[k],
             prog_pr.q_sno[k],
             aux_en.area[k],
-            ρ0_c[k],
+            ρ_c[k],
             Δt,
             ts,
         )
@@ -78,9 +80,9 @@ function quad_loop(en_thermo::SGSQuadrature, precip_model, vars, param_set, Δt:
     # θl - liquid ice potential temperature
     # _mean and ′ - subdomain mean and (co)variances
     # q_rai, q_sno - grid mean precipitation
-    UnPack.@unpack qt′qt′, qt_mean, θl′θl′, θl_mean, θl′qt′, subdomain_area, q_rai, q_sno, ρ0_c, p0_c = vars
+    UnPack.@unpack qt′qt′, qt_mean, θl′θl′, θl_mean, θl′qt′, subdomain_area, q_rai, q_sno, ρ_c, p_c = vars
 
-    FT = eltype(ρ0_c)
+    FT = eltype(ρ_c)
 
     inner_env = SA.MVector{env_len, FT}(undef)
     outer_env = SA.MVector{env_len, FT}(undef)
@@ -152,12 +154,12 @@ function quad_loop(en_thermo::SGSQuadrature, precip_model, vars, param_set, Δt:
             end
 
             # condensation
-            ts = thermo_state_pθq(param_set, p0_c, h_hat, qt_hat)
+            ts = thermo_state_pθq(param_set, p_c, h_hat, qt_hat)
             q_liq_en = TD.liquid_specific_humidity(param_set, ts)
             q_ice_en = TD.ice_specific_humidity(param_set, ts)
             T = TD.air_temperature(param_set, ts)
             # autoconversion and accretion
-            mph = precipitation_formation(param_set, precip_model, q_rai, q_sno, subdomain_area, ρ0_c, Δt, ts)
+            mph = precipitation_formation(param_set, precip_model, q_rai, q_sno, subdomain_area, ρ_c, Δt, ts)
 
             # environmental variables
             inner_env[i_ql] += q_liq_en * weights[m_h] * sqpi_inv
@@ -222,14 +224,16 @@ function microphysics(
     Δt::Real,
     param_set::APS,
 )
-    p0_c = center_ref_state(state).p0
-    ρ0_c = center_ref_state(state).ρ0
     aux_en = center_aux_environment(state)
     prog_pr = center_prog_precipitation(state)
+    prog_gm = center_prog_grid_mean(state)
+    aux_gm = center_aux_grid_mean(state)
     aux_en_unsat = aux_en.unsat
     aux_en_sat = aux_en.sat
     tendencies_pr = center_tendencies_precipitation(state)
     ts_env = center_aux_environment(state).ts
+    p_c = aux_gm.p
+    ρ_c = prog_gm.ρ
 
     #TODO - if we start using eos_smpl for the updrafts calculations
     #       we can get rid of the two categories for outer and inner quad. points
@@ -262,8 +266,8 @@ function microphysics(
                 subdomain_area = aux_en.area[k],
                 q_rai = prog_pr.q_rai[k],
                 q_sno = prog_pr.q_sno[k],
-                ρ0_c = ρ0_c[k],
-                p0_c = p0_c[k],
+                ρ_c = ρ_c[k],
+                p_c = p_c[k],
             )
             outer_env, outer_src = quad_loop(en_thermo, precip_model, vars, param_set, Δt)
 
@@ -287,7 +291,7 @@ function microphysics(
             if aux_en.cloud_fraction[k] < 1
                 aux_en_unsat.q_tot[k] = outer_env.qt_unsat / (1 - aux_en.cloud_fraction[k])
                 T_unsat = outer_env.T_unsat / (1 - aux_en.cloud_fraction[k])
-                ts_unsat = TD.PhaseEquil_pTq(param_set, p0_c[k], T_unsat, aux_en_unsat.q_tot[k])
+                ts_unsat = TD.PhaseEquil_pTq(param_set, p_c[k], T_unsat, aux_en_unsat.q_tot[k])
                 aux_en_unsat.θ_dry[k] = TD.dry_pottemp(param_set, ts_unsat)
             else
                 aux_en_unsat.q_tot[k] = 0
@@ -298,7 +302,7 @@ function microphysics(
                 aux_en_sat.T[k] = outer_env.T_sat / aux_en.cloud_fraction[k]
                 aux_en_sat.q_tot[k] = outer_env.qt_sat / aux_en.cloud_fraction[k]
                 aux_en_sat.q_vap[k] = (outer_env.qt_sat - outer_env.ql - outer_env.qi) / aux_en.cloud_fraction[k]
-                ts_sat = TD.PhaseEquil_pTq(param_set, p0_c[k], aux_en_sat.T[k], aux_en_sat.q_tot[k])
+                ts_sat = TD.PhaseEquil_pTq(param_set, p_c[k], aux_en_sat.T[k], aux_en_sat.q_tot[k])
                 aux_en_sat.θ_dry[k] = TD.dry_pottemp(param_set, ts_sat)
                 aux_en_sat.θ_liq_ice[k] = TD.liquid_ice_pottemp(param_set, ts_sat)
             else
@@ -324,7 +328,7 @@ function microphysics(
                 prog_pr.q_rai[k],
                 prog_pr.q_sno[k],
                 aux_en.area[k],
-                ρ0_c[k],
+                ρ_c[k],
                 Δt,
                 ts,
             )
diff --git a/src/EDMF_Precipitation.jl b/src/EDMF_Precipitation.jl
index 77564ff8f3..70b7488fad 100644
--- a/src/EDMF_Precipitation.jl
+++ b/src/EDMF_Precipitation.jl
@@ -18,10 +18,11 @@ function compute_precipitation_advection_tendencies(
 )
     FT = eltype(grid)
 
-    ρ0_c = center_ref_state(state).ρ0
     tendencies_pr = center_tendencies_precipitation(state)
     prog_pr = center_prog_precipitation(state)
     aux_tc = center_aux_turbconv(state)
+    prog_gm = center_prog_grid_mean(state)
+    ρ_c = prog_gm.ρ
 
     # helper to calculate the rain velocity
     # TODO: assuming w_gm = 0
@@ -38,8 +39,8 @@ function compute_precipitation_advection_tendencies(
     wvec = CC.Geometry.WVector
 
     # TODO - some positivity limiters are needed
-    @. aux_tc.qr_tendency_advection = ∇(wvec(RB(ρ0_c * q_rai * term_vel_rain))) / ρ0_c
-    @. aux_tc.qs_tendency_advection = ∇(wvec(RB(ρ0_c * q_sno * term_vel_snow))) / ρ0_c
+    @. aux_tc.qr_tendency_advection = ∇(wvec(RB(ρ_c * q_rai * term_vel_rain))) / ρ_c
+    @. aux_tc.qs_tendency_advection = ∇(wvec(RB(ρ_c * q_sno * term_vel_snow))) / ρ_c
 
     @. tendencies_pr.q_rai += aux_tc.qr_tendency_advection
     @. tendencies_pr.q_sno += aux_tc.qs_tendency_advection
@@ -67,18 +68,18 @@ function compute_precipitation_sink_tendencies(
     param_set::APS,
     Δt::Real,
 )
-    ρ0_c = center_ref_state(state).ρ0
     aux_gm = center_aux_grid_mean(state)
     aux_tc = center_aux_turbconv(state)
     prog_gm = center_prog_grid_mean(state)
     prog_pr = center_prog_precipitation(state)
+    ρ_c = prog_gm.ρ
     tendencies_pr = center_tendencies_precipitation(state)
     ts_gm = aux_gm.ts
 
     @inbounds for k in real_center_indices(grid)
         qr = prog_pr.q_rai[k]
         qs = prog_pr.q_sno[k]
-        ρ0 = ρ0_c[k]
+        ρ = ρ_c[k]
         q_tot_gm = aux_gm.q_tot[k]
         T_gm = aux_gm.T[k]
         # When we fuse loops, this should hopefully disappear
@@ -104,9 +105,9 @@ function compute_precipitation_sink_tendencies(
         # TODO - move limiters elsewhere
         # TODO - when using adaptive timestepping we are limiting the source terms
         #        with the previous timestep dt
-        S_qr_evap = -min(qr / Δt, -α_evp * CM1.evaporation_sublimation(param_set, rain_type, q, qr, ρ0, T_gm))
-        S_qs_melt = -min(qs / Δt, α_melt * CM1.snow_melt(param_set, qs, ρ0, T_gm))
-        tmp = α_dep_sub * CM1.evaporation_sublimation(param_set, snow_type, q, qs, ρ0, T_gm)
+        S_qr_evap = -min(qr / Δt, -α_evp * CM1.evaporation_sublimation(param_set, rain_type, q, qr, ρ, T_gm))
+        S_qs_melt = -min(qs / Δt, α_melt * CM1.snow_melt(param_set, qs, ρ, T_gm))
+        tmp = α_dep_sub * CM1.evaporation_sublimation(param_set, snow_type, q, qs, ρ, T_gm)
         if tmp > 0
             S_qs_sub_dep = min(qv / Δt, tmp)
         else
diff --git a/src/EDMF_Updrafts.jl b/src/EDMF_Updrafts.jl
index 434311315e..6c7e6e094f 100644
--- a/src/EDMF_Updrafts.jl
+++ b/src/EDMF_Updrafts.jl
@@ -10,19 +10,21 @@ function compute_nonequilibrium_moisture_tendencies!(
     param_set::APS,
 )
     N_up = n_updrafts(edmf)
-    p0_c = center_ref_state(state).p0
-    ρ0_c = center_ref_state(state).ρ0
+    aux_gm = center_aux_grid_mean(state)
     aux_up = center_aux_updrafts(state)
     aux_bulk = center_aux_bulk(state)
+    prog_gm = center_prog_grid_mean(state)
+    p_c = aux_gm.p
+    ρ_c = prog_gm.ρ
 
     @inbounds for i in 1:N_up
         @inbounds for k in real_center_indices(grid)
             T_up = aux_up[i].T[k]
             q_up = TD.PhasePartition(aux_up[i].q_tot[k], aux_up[i].q_liq[k], aux_up[i].q_ice[k])
-            ts_up = TD.PhaseNonEquil_pTq(param_set, p0_c[k], T_up, q_up)
+            ts_up = TD.PhaseNonEquil_pTq(param_set, p_c[k], T_up, q_up)
 
             # condensation/evaporation, deposition/sublimation
-            mph = noneq_moisture_sources(param_set, aux_up[i].area[k], ρ0_c[k], Δt, ts_up)
+            mph = noneq_moisture_sources(param_set, aux_up[i].area[k], ρ_c[k], Δt, ts_up)
             aux_up[i].ql_tendency_noneq[k] = mph.ql_tendency * aux_up[i].area[k]
             aux_up[i].qi_tendency_noneq[k] = mph.qi_tendency * aux_up[i].area[k]
         end
@@ -50,24 +52,26 @@ function compute_precipitation_formation_tendencies(
     param_set::APS,
 )
     N_up = n_updrafts(edmf)
-    p0_c = center_ref_state(state).p0
-    ρ0_c = center_ref_state(state).ρ0
+    prog_gm = center_prog_grid_mean(state)
+    aux_gm = center_aux_grid_mean(state)
     aux_up = center_aux_updrafts(state)
     aux_bulk = center_aux_bulk(state)
     prog_pr = center_prog_precipitation(state)
     tendencies_pr = center_tendencies_precipitation(state)
+    p_c = aux_gm.p
+    ρ_c = prog_gm.ρ
 
     @inbounds for i in 1:N_up
         @inbounds for k in real_center_indices(grid)
             T_up = aux_up[i].T[k]
             q_tot_up = aux_up[i].q_tot[k]
             if edmf.moisture_model isa EquilibriumMoisture
-                ts_up = TD.PhaseEquil_pTq(param_set, p0_c[k], T_up, q_tot_up)
+                ts_up = TD.PhaseEquil_pTq(param_set, p_c[k], T_up, q_tot_up)
             elseif edmf.moisture_model isa NonEquilibriumMoisture
                 q_liq_up = aux_up[i].q_liq[k]
                 q_ice_up = aux_up[i].q_ice[k]
                 q = TD.PhasePartition(q_tot_up, q_liq_up, q_ice_up)
-                ts_up = TD.PhaseNonEquil_pTq(param_set, p0_c[k], T_up, q)
+                ts_up = TD.PhaseNonEquil_pTq(param_set, p_c[k], T_up, q)
             else
                 error(
                     "Something went wrong in EDMF_Updrafts. The expected moisture model is Equilibrium or NonEquilibrium",
@@ -81,7 +85,7 @@ function compute_precipitation_formation_tendencies(
                 prog_pr.q_rai[k],
                 prog_pr.q_sno[k],
                 aux_up[i].area[k],
-                ρ0_c[k],
+                ρ_c[k],
                 Δt,
                 ts_up,
             )
diff --git a/src/EDMF_functions.jl b/src/EDMF_functions.jl
index 7a95806da4..0ebfd72a13 100755
--- a/src/EDMF_functions.jl
+++ b/src/EDMF_functions.jl
@@ -60,9 +60,9 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
     aux_up = center_aux_updrafts(state)
     aux_tc_f = face_aux_turbconv(state)
     aux_up_f = face_aux_updrafts(state)
-    ρ0_f = face_ref_state(state).ρ0
-    p0_c = center_ref_state(state).p0
-    α0_c = center_ref_state(state).α0
+    ρ_f = aux_gm_f.ρ
+    ρ_c = prog_gm.ρ
+    p_c = aux_gm.p
     kf_surf = kf_surface(grid)
     kc_surf = kc_surface(grid)
     massflux = aux_tc_f.massflux
@@ -86,9 +86,9 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
     If = CCO.InterpolateC2F(; bottom = CCO.SetValue(FT(0)), top = CCO.SetValue(FT(0)))
 
     # Compute the mass flux and associated scalar fluxes
-    @. massflux = ρ0_f * Ifae(a_en) * (w_en - w_gm)
-    @. massflux_h = ρ0_f * Ifae(a_en) * (w_en - w_gm) * (If(θ_liq_ice_en) - If(θ_liq_ice_gm))
-    @. massflux_qt = ρ0_f * Ifae(a_en) * (w_en - w_gm) * (If(q_tot_en) - If(q_tot_gm))
+    @. massflux = ρ_f * Ifae(a_en) * (w_en - w_gm)
+    @. massflux_h = ρ_f * Ifae(a_en) * (w_en - w_gm) * (If(θ_liq_ice_en) - If(θ_liq_ice_gm))
+    @. massflux_qt = ρ_f * Ifae(a_en) * (w_en - w_gm) * (If(q_tot_en) - If(q_tot_gm))
     @inbounds for i in 1:N_up
         aux_up_f_i = aux_up_f[i]
         aux_up_i = aux_up[i]
@@ -98,11 +98,11 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
         w_up = aux_up_f[i].w
         θ_liq_ice_up = aux_up_i.θ_liq_ice
         q_tot_up = aux_up_i.q_tot
-        @. aux_up_f[i].massflux = ρ0_f * Ifau(a_up) * (w_up - w_gm)
+        @. aux_up_f[i].massflux = ρ_f * Ifau(a_up) * (w_up - w_gm)
         # We know that, since W = 0 at z = 0, m = 0 also, and
         # therefore θ_liq_ice / q_tot values do not matter
-        @. massflux_h += ρ0_f * (Ifau(a_up) * (w_up - w_gm) * (If(θ_liq_ice_up) - If(θ_liq_ice_gm)))
-        @. massflux_qt += ρ0_f * (Ifau(a_up) * (w_up - w_gm) * (If(q_tot_up) - If(q_tot_gm)))
+        @. massflux_h += ρ_f * (Ifau(a_up) * (w_up - w_gm) * (If(θ_liq_ice_up) - If(θ_liq_ice_gm)))
+        @. massflux_qt += ρ_f * (Ifau(a_up) * (w_up - w_gm) * (If(q_tot_up) - If(q_tot_gm)))
     end
     massflux[kf_surf] = 0
     massflux_h[kf_surf] = 0
@@ -115,7 +115,7 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
         q_ice_en = aux_en.q_ice
         q_liq_gm = prog_gm.q_liq
         q_ice_gm = prog_gm.q_ice
-        @. massflux_en = ρ0_f * Ifae(a_en) * (w_en - w_gm)
+        @. massflux_en = ρ_f * Ifae(a_en) * (w_en - w_gm)
         @. massflux_ql = massflux_en * (If(q_liq_en) - If(q_liq_gm))
         @. massflux_qi = massflux_en * (If(q_ice_en) - If(q_ice_gm))
         @inbounds for i in 1:N_up
@@ -136,8 +136,8 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
     # Compute the  mass flux tendencies
     # Adjust the values of the grid mean variables
     # Prepare the output
-    @. massflux_tendency_h = -∇c(wvec(massflux_h)) * α0_c
-    @. massflux_tendency_qt = -∇c(wvec(massflux_qt)) * α0_c
+    @. massflux_tendency_h = -∇c(wvec(massflux_h)) / ρ_c
+    @. massflux_tendency_qt = -∇c(wvec(massflux_qt)) / ρ_c
 
     diffusive_flux_h = aux_tc_f.diffusive_flux_h
     diffusive_flux_qt = aux_tc_f.diffusive_flux_qt
@@ -158,8 +158,8 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
         massflux_tendency_ql = aux_tc.massflux_tendency_ql
         massflux_tendency_qi = aux_tc.massflux_tendency_qi
 
-        @. massflux_tendency_ql = -∇c(wvec(massflux_ql)) * α0_c
-        @. massflux_tendency_qi = -∇c(wvec(massflux_qi)) * α0_c
+        @. massflux_tendency_ql = -∇c(wvec(massflux_ql)) / ρ_c
+        @. massflux_tendency_qi = -∇c(wvec(massflux_qi)) / ρ_c
 
         diffusive_flux_ql = aux_tc_f.diffusive_flux_ql
         diffusive_flux_qi = aux_tc_f.diffusive_flux_qi
@@ -176,14 +176,15 @@ end
 
 function compute_diffusive_fluxes(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBase, param_set::APS)
     FT = eltype(grid)
-    ρ0_f = face_ref_state(state).ρ0
     aux_bulk = center_aux_bulk(state)
     aux_tc_f = face_aux_turbconv(state)
     aux_en = center_aux_environment(state)
+    aux_gm_f = face_aux_grid_mean(state)
     KM = center_aux_turbconv(state).KM
     KH = center_aux_turbconv(state).KH
     aeKM = center_aux_turbconv(state).ϕ_temporary
     aeKH = center_aux_turbconv(state).ψ_temporary
+    ρ_f = aux_gm_f.ρ
     a_en = aux_en.area
     @. aeKM = a_en * KM
     @. aeKH = a_en * KH
@@ -194,8 +195,8 @@ function compute_diffusive_fluxes(edmf::EDMFModel, grid::Grid, state::State, sur
     IfKH = CCO.InterpolateC2F(; bottom = CCO.SetValue(aeKH[kc_surf]), top = CCO.SetValue(aeKH[kc_toa]))
     IfKM = CCO.InterpolateC2F(; bottom = CCO.SetValue(aeKM[kc_surf]), top = CCO.SetValue(aeKM[kc_toa]))
 
-    @. aux_tc_f.ρ_ae_KH = IfKH(aeKH) * ρ0_f
-    @. aux_tc_f.ρ_ae_KM = IfKM(aeKM) * ρ0_f
+    @. aux_tc_f.ρ_ae_KH = IfKH(aeKH) * ρ_f
+    @. aux_tc_f.ρ_ae_KM = IfKM(aeKM) * ρ_f
 
     aeKHq_tot_bc = -surf.ρq_tot_flux / a_en[kc_surf] / aux_tc_f.ρ_ae_KH[kf_surf]
     aeKHθ_liq_ice_bc = -surf.ρθ_liq_ice_flux / a_en[kc_surf] / aux_tc_f.ρ_ae_KH[kf_surf]
@@ -265,18 +266,20 @@ function set_edmf_surface_bc(edmf::EDMFModel, grid::Grid, state::State, surf::Su
     N_up = n_updrafts(edmf)
     kc_surf = kc_surface(grid)
     kf_surf = kf_surface(grid)
-    ρ0_c = center_ref_state(state).ρ0
-    ρ0_f = face_ref_state(state).ρ0
     prog_up = center_prog_updrafts(state)
     prog_en = center_prog_environment(state)
     prog_up_f = face_prog_updrafts(state)
     aux_bulk = center_aux_bulk(state)
     aux_en = center_aux_environment(state)
+    prog_gm = center_prog_grid_mean(state)
+    aux_gm_f = face_aux_grid_mean(state)
+    ρ_c = prog_gm.ρ
+    ρ_f = aux_gm_f.ρ
     @inbounds for i in 1:N_up
         θ_surf = θ_surface_bc(surf, grid, state, edmf, i)
         q_surf = q_surface_bc(surf, grid, state, edmf, i)
         a_surf = area_surface_bc(surf, edmf, i)
-        prog_up[i].ρarea[kc_surf] = ρ0_c[kc_surf] * a_surf
+        prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_surf
         prog_up[i].ρaθ_liq_ice[kc_surf] = prog_up[i].ρarea[kc_surf] * θ_surf
         prog_up[i].ρaq_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * q_surf
         if edmf.moisture_model isa NonEquilibriumMoisture
@@ -285,7 +288,7 @@ function set_edmf_surface_bc(edmf::EDMFModel, grid::Grid, state::State, surf::Su
             prog_up[i].ρaq_liq[kc_surf] = prog_up[i].ρarea[kc_surf] * q_liq_surf
             prog_up[i].ρaq_ice[kc_surf] = prog_up[i].ρarea[kc_surf] * q_ice_surf
         end
-        prog_up_f[i].ρaw[kf_surf] = ρ0_f[kf_surf] * w_surface_bc(surf)
+        prog_up_f[i].ρaw[kf_surf] = ρ_f[kf_surf] * w_surface_bc(surf)
     end
 
     flux1 = surf.ρθ_liq_ice_flux
@@ -293,28 +296,29 @@ function set_edmf_surface_bc(edmf::EDMFModel, grid::Grid, state::State, surf::Su
     zLL = grid.zc[kc_surf].z
     ustar = surf.ustar
     oblength = surf.obukhov_length
-    α0LL = center_ref_state(state).α0[kc_surf]
+    ρLL = prog_gm.ρ[kc_surf]
     # TODO: is bulk even defined before this is called?
     ae_surf = 1 - aux_bulk.area[kc_surf] # area of environment
 
-    ρ0_ae = ρ0_c[kc_surf] * ae_surf
+    ρ_ae = ρ_c[kc_surf] * ae_surf
 
-    prog_en.ρatke[kc_surf] = ρ0_ae * get_surface_tke(param_set, surf.ustar, zLL, surf.obukhov_length)
+    prog_en.ρatke[kc_surf] = ρ_ae * get_surface_tke(param_set, surf.ustar, zLL, surf.obukhov_length)
     if edmf.thermo_covariance_model isa PrognosticThermoCovariances
-        prog_en.ρaHvar[kc_surf] = ρ0_ae * get_surface_variance(flux1 * α0LL, flux1 * α0LL, ustar, zLL, oblength)
-        prog_en.ρaQTvar[kc_surf] = ρ0_ae * get_surface_variance(flux2 * α0LL, flux2 * α0LL, ustar, zLL, oblength)
-        prog_en.ρaHQTcov[kc_surf] = ρ0_ae * get_surface_variance(flux1 * α0LL, flux2 * α0LL, ustar, zLL, oblength)
+        prog_en.ρaHvar[kc_surf] = ρ_ae * get_surface_variance(flux1 / ρLL, flux1 / ρLL, ustar, zLL, oblength)
+        prog_en.ρaQTvar[kc_surf] = ρ_ae * get_surface_variance(flux2 / ρLL, flux2 / ρLL, ustar, zLL, oblength)
+        prog_en.ρaHQTcov[kc_surf] = ρ_ae * get_surface_variance(flux1 / ρLL, flux2 / ρLL, ustar, zLL, oblength)
     end
     return nothing
 end
 
 function surface_helper(surf::SurfaceBase, grid::Grid, state::State)
     kc_surf = kc_surface(grid)
+    prog_gm = center_prog_grid_mean(state)
     zLL = grid.zc[kc_surf].z
     ustar = surf.ustar
     oblength = surf.obukhov_length
-    α0LL = center_ref_state(state).α0[kc_surf]
-    return (; ustar, zLL, oblength, α0LL)
+    ρLL = prog_gm.ρ[kc_surf]
+    return (; ustar, zLL, oblength, ρLL)
 end
 
 function a_up_boundary_conditions(surf::SurfaceBase, edmf::EDMFModel, i::Int)
@@ -349,9 +353,9 @@ function θ_surface_bc(surf::SurfaceBase{FT}, grid::Grid, state::State, edmf::ED
     surf.bflux > 0 || return FT(0)
     a_total = edmf.surface_area
     a_ = area_surface_bc(surf, edmf, i)
-    UnPack.@unpack ustar, zLL, oblength, α0LL = surface_helper(surf, grid, state)
+    UnPack.@unpack ustar, zLL, oblength, ρLL = surface_helper(surf, grid, state)
     ρθ_liq_ice_flux = surf.ρθ_liq_ice_flux
-    h_var = get_surface_variance(ρθ_liq_ice_flux * α0LL, ρθ_liq_ice_flux * α0LL, ustar, zLL, oblength)
+    h_var = get_surface_variance(ρθ_liq_ice_flux / ρLL, ρθ_liq_ice_flux / ρLL, ustar, zLL, oblength)
     surface_scalar_coeff = percentile_bounds_mean_norm(1 - a_total + (i - 1) * a_, 1 - a_total + i * a_, 1000)
     return aux_gm.θ_liq_ice[kc_surf] + surface_scalar_coeff * sqrt(h_var)
 end
@@ -362,9 +366,9 @@ function q_surface_bc(surf::SurfaceBase{FT}, grid::Grid, state::State, edmf::EDM
     surf.bflux > 0 || return aux_gm.q_tot[kc_surf]
     a_total = edmf.surface_area
     a_ = area_surface_bc(surf, edmf, i)
-    UnPack.@unpack ustar, zLL, oblength, α0LL = surface_helper(surf, grid, state)
+    UnPack.@unpack ustar, zLL, oblength, ρLL = surface_helper(surf, grid, state)
     ρq_tot_flux = surf.ρq_tot_flux
-    qt_var = get_surface_variance(ρq_tot_flux * α0LL, ρq_tot_flux * α0LL, ustar, zLL, oblength)
+    qt_var = get_surface_variance(ρq_tot_flux / ρLL, ρq_tot_flux / ρLL, ustar, zLL, oblength)
     surface_scalar_coeff = percentile_bounds_mean_norm(1 - a_total + (i - 1) * a_, 1 - a_total + i * a_, 1000)
     return aux_gm.q_tot[kc_surf] + surface_scalar_coeff * sqrt(qt_var)
 end
@@ -463,8 +467,10 @@ function compute_up_tendencies!(edmf::EDMFModel, grid::Grid, state::State, param
     prog_up = center_prog_updrafts(state)
     tendencies_up = center_tendencies_updrafts(state)
     tendencies_up_f = face_tendencies_updrafts(state)
-    ρ0_c = center_ref_state(state).ρ0
-    ρ0_f = face_ref_state(state).ρ0
+    prog_gm = center_prog_grid_mean(state)
+    aux_gm_f = face_aux_grid_mean(state)
+    ρ_c = prog_gm.ρ
+    ρ_f = aux_gm_f.ρ
     au_lim = edmf.max_area
 
     @inbounds for i in 1:N_up
@@ -500,17 +506,17 @@ function compute_up_tendencies!(edmf::EDMFModel, grid::Grid, state::State, param
         tends_ρaq_tot = tendencies_up[i].ρaq_tot
 
         @. tends_ρarea =
-            -∇c(wvec(LBF(Ic(w_up) * ρ0_c * a_up))) + (ρ0_c * a_up * Ic(w_up) * entr_turb_dyn) -
-            (ρ0_c * a_up * Ic(w_up) * detr_turb_dyn)
+            -∇c(wvec(LBF(Ic(w_up) * ρ_c * a_up))) + (ρ_c * a_up * Ic(w_up) * entr_turb_dyn) -
+            (ρ_c * a_up * Ic(w_up) * detr_turb_dyn)
 
         @. tends_ρaθ_liq_ice =
-            -∇c(wvec(LBF(Ic(w_up) * ρ0_c * a_up * θ_liq_ice_up))) +
-            (ρ0_c * a_up * Ic(w_up) * entr_turb_dyn * θ_liq_ice_en) -
-            (ρ0_c * a_up * Ic(w_up) * detr_turb_dyn * θ_liq_ice_up) + (ρ0_c * θ_liq_ice_tendency_precip_formation)
+            -∇c(wvec(LBF(Ic(w_up) * ρ_c * a_up * θ_liq_ice_up))) +
+            (ρ_c * a_up * Ic(w_up) * entr_turb_dyn * θ_liq_ice_en) -
+            (ρ_c * a_up * Ic(w_up) * detr_turb_dyn * θ_liq_ice_up) + (ρ_c * θ_liq_ice_tendency_precip_formation)
 
         @. tends_ρaq_tot =
-            -∇c(wvec(LBF(Ic(w_up) * ρ0_c * a_up * q_tot_up))) + (ρ0_c * a_up * Ic(w_up) * entr_turb_dyn * q_tot_en) -
-            (ρ0_c * a_up * Ic(w_up) * detr_turb_dyn * q_tot_up) + (ρ0_c * qt_tendency_precip_formation)
+            -∇c(wvec(LBF(Ic(w_up) * ρ_c * a_up * q_tot_up))) + (ρ_c * a_up * Ic(w_up) * entr_turb_dyn * q_tot_en) -
+            (ρ_c * a_up * Ic(w_up) * detr_turb_dyn * q_tot_up) + (ρ_c * qt_tendency_precip_formation)
 
         if edmf.moisture_model isa NonEquilibriumMoisture
 
@@ -528,16 +534,14 @@ function compute_up_tendencies!(edmf::EDMFModel, grid::Grid, state::State, param
             tends_ρaq_ice = tendencies_up[i].ρaq_ice
 
             @. tends_ρaq_liq =
-                -∇c(wvec(LBF(Ic(w_up) * ρ0_c * a_up * q_liq_up))) +
-                (ρ0_c * a_up * Ic(w_up) * entr_turb_dyn * q_liq_en) -
-                (ρ0_c * a_up * Ic(w_up) * detr_turb_dyn * q_liq_up) +
-                (ρ0_c * (ql_tendency_precip_formation + ql_tendency_noneq))
+                -∇c(wvec(LBF(Ic(w_up) * ρ_c * a_up * q_liq_up))) + (ρ_c * a_up * Ic(w_up) * entr_turb_dyn * q_liq_en) -
+                (ρ_c * a_up * Ic(w_up) * detr_turb_dyn * q_liq_up) +
+                (ρ_c * (ql_tendency_precip_formation + ql_tendency_noneq))
 
             @. tends_ρaq_ice =
-                -∇c(wvec(LBF(Ic(w_up) * ρ0_c * a_up * q_ice_up))) +
-                (ρ0_c * a_up * Ic(w_up) * entr_turb_dyn * q_ice_en) -
-                (ρ0_c * a_up * Ic(w_up) * detr_turb_dyn * q_ice_up) +
-                (ρ0_c * (qi_tendency_precip_formation + qi_tendency_noneq))
+                -∇c(wvec(LBF(Ic(w_up) * ρ_c * a_up * q_ice_up))) + (ρ_c * a_up * Ic(w_up) * entr_turb_dyn * q_ice_en) -
+                (ρ_c * a_up * Ic(w_up) * detr_turb_dyn * q_ice_up) +
+                (ρ_c * (qi_tendency_precip_formation + qi_tendency_noneq))
 
             tends_ρaq_liq[kc_surf] = 0
             tends_ρaq_ice[kc_surf] = 0
@@ -585,10 +589,10 @@ function compute_up_tendencies!(edmf::EDMFModel, grid::Grid, state::State, param
         detr_w = aux_up[i].detr_turb_dyn
         buoy = aux_up[i].buoy
 
-        @. tends_ρaw = -(∇f(wvec(LBC(Iaf(a_up) * ρ0_f * w_up * w_up))))
+        @. tends_ρaw = -(∇f(wvec(LBC(Iaf(a_up) * ρ_f * w_up * w_up))))
         @. tends_ρaw +=
-            (ρ0_f * Iaf(a_up) * w_up * (I0f(entr_w) * w_en - I0f(detr_w) * w_up)) +
-            (ρ0_f * Iaf(a_up) * I0f(buoy)) +
+            (ρ_f * Iaf(a_up) * w_up * (I0f(entr_w) * w_en - I0f(detr_w) * w_up)) +
+            (ρ_f * Iaf(a_up) * I0f(buoy)) +
             nh_pressure
         tends_ρaw[kf_surf] = 0
     end
@@ -605,11 +609,12 @@ function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::Su
 
     prog_up = center_prog_updrafts(state)
     prog_gm = center_prog_grid_mean(state)
+    aux_gm_f = face_aux_grid_mean(state)
     aux_up = center_aux_updrafts(state)
     aux_up_f = face_aux_updrafts(state)
     prog_up_f = face_prog_updrafts(state)
-    ρ0_c = center_ref_state(state).ρ0
-    ρ0_f = face_ref_state(state).ρ0
+    ρ_c = prog_gm.ρ
+    ρ_f = aux_gm_f.ρ
     a_min = edmf.minimum_area
     a_max = edmf.max_area
 
@@ -622,7 +627,7 @@ function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::Su
             @. prog_up[i].δ_nondim = max(prog_up[i].δ_nondim, 0)
         end
         @inbounds for k in real_center_indices(grid)
-            prog_up[i].ρarea[k] = min(prog_up[i].ρarea[k], ρ0_c[k] * a_max)
+            prog_up[i].ρarea[k] = min(prog_up[i].ρarea[k], ρ_c[k] * a_max)
         end
     end
     if edmf.moisture_model isa NonEquilibriumMoisture
@@ -636,7 +641,7 @@ function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::Su
         @. prog_up_f[i].ρaw = max.(prog_up_f[i].ρaw, 0)
         a_up_bcs = a_up_boundary_conditions(surf, edmf, i)
         If = CCO.InterpolateC2F(; a_up_bcs...)
-        @. prog_up_f[i].ρaw = Int(If(prog_up[i].ρarea) >= ρ0_f * a_min) * prog_up_f[i].ρaw
+        @. prog_up_f[i].ρaw = Int(If(prog_up[i].ρarea) >= ρ_f * a_min) * prog_up_f[i].ρaw
     end
 
     @inbounds for k in real_center_indices(grid)
@@ -646,7 +651,7 @@ function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::Su
             # this is needed to make sure Rico is unchanged.
             # TODO : look into it further to see why
             # a similar filtering of ρaθ_liq_ice breaks the simulation
-            if prog_up[i].ρarea[k] / ρ0_c[k] < a_min
+            if prog_up[i].ρarea[k] / ρ_c[k] < a_min
                 prog_up[i].ρaq_tot[k] = 0
             end
         end
@@ -659,7 +664,7 @@ function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::Su
                 # this is needed to make sure Rico is unchanged.
                 # TODO : look into it further to see why
                 # a similar filtering of ρaθ_liq_ice breaks the simulation
-                if prog_up[i].ρarea[k] / ρ0_c[k] < a_min
+                if prog_up[i].ρarea[k] / ρ_c[k] < a_min
                     prog_up[i].ρaq_liq[k] = 0
                     prog_up[i].ρaq_ice[k] = 0
                 end
@@ -676,7 +681,7 @@ function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::Su
         θ_surf = θ_surface_bc(surf, grid, state, edmf, i)
         q_surf = q_surface_bc(surf, grid, state, edmf, i)
         a_surf = area_surface_bc(surf, edmf, i)
-        prog_up[i].ρarea[kc_surf] = ρ0_c[kc_surf] * a_surf
+        prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_surf
         prog_up[i].ρaθ_liq_ice[kc_surf] = prog_up[i].ρarea[kc_surf] * θ_surf
         prog_up[i].ρaq_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * q_surf
     end
@@ -703,8 +708,8 @@ function compute_covariance_shear(
 ) where {covar_sym, ϕ_en_sym, ψ_en_sym}
 
     aux_tc = center_aux_turbconv(state)
-    ρ0_c = center_ref_state(state).ρ0
     prog_gm = center_prog_grid_mean(state)
+    ρ_c = prog_gm.ρ
     is_tke = covar_sym == :tke
     tke_factor = is_tke ? 0.5 : 1
     k_eddy = is_tke ? aux_tc.KM : aux_tc.KH
@@ -731,12 +736,12 @@ function compute_covariance_shear(
         @. shear =
             tke_factor *
             2 *
-            ρ0_c *
+            ρ_c *
             area_en *
             k_eddy *
             (∇c(wvec(ϕ_en)) * ∇c(wvec(ψ_en)) + (∇c(wvec(If(u))))^2 + (∇c(wvec(If(v))))^2)
     else
-        @. shear = tke_factor * 2 * ρ0_c * area_en * k_eddy * ∇c(wvec(If(ϕ_en))) * ∇c(wvec(If(ψ_en)))
+        @. shear = tke_factor * 2 * ρ_c * area_en * k_eddy * ∇c(wvec(If(ϕ_en))) * ∇c(wvec(If(ψ_en)))
     end
     return nothing
 end
@@ -784,7 +789,6 @@ function compute_covariance_entr(
 ) where {covar_sym, ϕ_sym, ψ_sym}
 
     N_up = n_updrafts(edmf)
-    ρ0_c = center_ref_state(state).ρ0
     FT = eltype(grid)
     is_tke = covar_sym == :tke
     tke_factor = is_tke ? 0.5 : 1
@@ -792,6 +796,8 @@ function compute_covariance_entr(
     aux_up_f = face_aux_updrafts(state)
     aux_gm_c = center_aux_grid_mean(state)
     prog_gm_f = face_prog_grid_mean(state)
+    prog_gm = center_prog_grid_mean(state)
+    ρ_c = prog_gm.ρ
     gm = is_tke ? prog_gm_f : aux_gm_c
     prog_up = is_tke ? aux_up_f : aux_up
     ϕ_gm = getproperty(gm, ϕ_sym)
@@ -828,16 +834,16 @@ function compute_covariance_entr(
 
         @. entr_gain +=
             Int(a_up > min_area) *
-            (tke_factor * ρ0_c * a_up * abs(Ic(w_up)) * detr_sc * (Idc(ϕ_up) - Idc(ϕ_en)) * (Idc(ψ_up) - Idc(ψ_en))) + (
+            (tke_factor * ρ_c * a_up * abs(Ic(w_up)) * detr_sc * (Idc(ϕ_up) - Idc(ϕ_en)) * (Idc(ψ_up) - Idc(ψ_en))) + (
                 tke_factor *
-                ρ0_c *
+                ρ_c *
                 a_up *
                 abs(Ic(w_up)) *
                 eps_turb *
                 ((Idc(ϕ_en) - Idc(ϕ_gm)) * (Idc(ψ_up) - Idc(ψ_en)) + (Idc(ψ_en) - Idc(ψ_gm)) * (Idc(ϕ_up) - Idc(ϕ_en)))
             )
 
-        @. detr_loss += Int(a_up > min_area) * tke_factor * ρ0_c * a_up * abs(Ic(w_up)) * (entr_sc + eps_turb) * covar
+        @. detr_loss += Int(a_up > min_area) * tke_factor * ρ_c * a_up * abs(Ic(w_up)) * (entr_sc + eps_turb) * covar
 
     end
 
@@ -854,9 +860,10 @@ function compute_covariance_dissipation(
     FT = eltype(grid)
     c_d::FT = ICP.c_d(param_set)
     aux_tc = center_aux_turbconv(state)
-    ρ0_c = center_ref_state(state).ρ0
     prog_en = center_prog_environment(state)
+    prog_gm = center_prog_grid_mean(state)
     aux_en = center_aux_environment(state)
+    ρ_c = prog_gm.ρ
     aux_en_2m = center_aux_environment_2m(state)
     aux_covar = getproperty(aux_en_2m, covar_sym)
     covar = getproperty(aux_en, covar_sym)
@@ -865,7 +872,7 @@ function compute_covariance_dissipation(
     tke_en = aux_en.tke
     mixing_length = aux_tc.mixing_length
 
-    @. dissipation = ρ0_c * area_en * covar * max(tke_en, 0)^FT(0.5) / max(mixing_length, FT(1.0e-3)) * c_d
+    @. dissipation = ρ_c * area_en * covar * max(tke_en, 0)^FT(0.5) / max(mixing_length, FT(1.0e-3)) * c_d
     return nothing
 end
 
@@ -880,9 +887,9 @@ function compute_en_tendencies!(
     N_up = n_updrafts(edmf)
     kc_surf = kc_surface(grid)
     kc_toa = kc_top_of_atmos(grid)
-    ρ0_c = center_ref_state(state).ρ0
-    ρ0_f = face_ref_state(state).ρ0
+    aux_gm_f = face_aux_grid_mean(state)
     prog_en = center_prog_environment(state)
+    prog_gm = center_prog_grid_mean(state)
     aux_en_2m = center_aux_environment_2m(state)
     tendencies_en = center_tendencies_environment(state)
     prog_covar = getproperty(tendencies_en, prog_sym)
@@ -892,6 +899,8 @@ function compute_en_tendencies!(
     aux_covar = getproperty(aux_en_2m, covar_sym)
     aux_up = center_aux_updrafts(state)
     w_en_f = face_aux_environment(state).w
+    ρ_c = prog_gm.ρ
+    ρ_f = aux_gm_f.ρ
     c_d = ICP.c_d(param_set)
     is_tke = covar_sym == :tke
     FT = eltype(grid)
@@ -940,14 +949,14 @@ function compute_en_tendencies!(
         a_up = aux_up[i].area
         # TODO: using `Int(bool) *` means that NaNs can propagate
         # into the solution. Could we somehow call `ifelse` instead?
-        @. D_env += Int(a_up > min_area) * ρ0_c * a_up * Ic(w_up) * (entr_sc + turb_entr)
+        @. D_env += Int(a_up > min_area) * ρ_c * a_up * Ic(w_up) * (entr_sc + turb_entr)
     end
 
     RB = CCO.RightBiasedC2F(; top = CCO.SetValue(FT(0)))
     @. prog_covar =
         press + buoy + shear + entr_gain + rain_src - D_env * covar -
-        (ρ0_c * area_en * c_d * sqrt(max(tke_en, 0)) / max(mixing_length, 1)) * covar -
-        ∇c(wvec(RB(ρ0_c * area_en * Ic(w_en_f) * covar))) + ∇c(ρ0_f * If(aeK) * ∇f(covar))
+        (ρ_c * area_en * c_d * sqrt(max(tke_en, 0)) / max(mixing_length, 1)) * covar -
+        ∇c(wvec(RB(ρ_c * area_en * Ic(w_en_f) * covar))) + ∇c(ρ_f * If(aeK) * ∇f(covar))
 
     return nothing
 end
@@ -963,7 +972,8 @@ function update_diagnostic_covariances!(
     N_up = n_updrafts(edmf)
     kc_surf = kc_surface(grid)
     kc_toa = kc_top_of_atmos(grid)
-    ρ0_c = center_ref_state(state).ρ0
+    prog_gm = center_prog_grid_mean(state)
+    ρ_c = prog_gm.ρ
     aux_en_2m = center_aux_environment_2m(state)
     aux_up_f = face_aux_updrafts(state)
     aux_en = center_aux_environment(state)
@@ -999,12 +1009,12 @@ function update_diagnostic_covariances!(
         a_up = aux_up[i].area
         # TODO: using `Int(bool) *` means that NaNs can propagate
         # into the solution. Could we somehow call `ifelse` instead?
-        @. D_env += Int(a_up > min_area) * ρ0_c * a_up * Ic(w_up) * (entr_sc + turb_entr)
+        @. D_env += Int(a_up > min_area) * ρ_c * a_up * Ic(w_up) * (entr_sc + turb_entr)
     end
 
     @. covar =
         (shear + entr_gain + rain_src) /
-        max(D_env + ρ0_c * area_en * c_d * sqrt(max(tke_en, 0)) / max(mixing_length, 1), covar_lim)
+        max(D_env + ρ_c * area_en * c_d * sqrt(max(tke_en, 0)) / max(mixing_length, 1), covar_lim)
     return nothing
 end
 
diff --git a/src/closures/buoyancy_gradients.jl b/src/closures/buoyancy_gradients.jl
index 1a13c5a8d6..b262071b82 100644
--- a/src/closures/buoyancy_gradients.jl
+++ b/src/closures/buoyancy_gradients.jl
@@ -20,12 +20,12 @@ function buoyancy_gradients(
     R_v = ICP.R_v(param_set)
 
     phase_part = TD.PhasePartition(0.0, 0.0, 0.0) # assuming R_d = R_m
-    Π = TD.exner_given_pressure(param_set, bg_model.p0, phase_part)
+    Π = TD.exner_given_pressure(param_set, bg_model.p, phase_part)
 
-    ∂b∂θv = g * (R_d / bg_model.alpha0 / bg_model.p0) * Π
+    ∂b∂θv = g * (R_d / bg_model.alpha0 / bg_model.p) * Π
 
     if bg_model.en_cld_frac > 0.0
-        ts_sat = thermo_state_pθq(param_set, bg_model.p0, bg_model.θ_liq_ice_sat, bg_model.qt_sat)
+        ts_sat = thermo_state_pθq(param_set, bg_model.p, bg_model.θ_liq_ice_sat, bg_model.qt_sat)
         phase_part = TD.PhasePartition(param_set, ts_sat)
         lh = TD.latent_heat_liq_ice(param_set, phase_part)
         cp_m = TD.cp_m(param_set, ts_sat)
diff --git a/src/closures/entr_detr.jl b/src/closures/entr_detr.jl
index a02752a7c0..db214f8cc8 100644
--- a/src/closures/entr_detr.jl
+++ b/src/closures/entr_detr.jl
@@ -156,10 +156,11 @@ function compute_entr_detr!(
     aux_en = center_aux_environment(state)
     aux_en_f = face_aux_environment(state)
     prog_gm_f = face_prog_grid_mean(state)
+    prog_gm = center_prog_grid_mean(state)
     aux_gm = center_aux_grid_mean(state)
     aux_tc = center_aux_turbconv(state)
-    p0_c = center_ref_state(state).p0
-    ρ0_c = center_ref_state(state).ρ0
+    p_c = aux_gm.p
+    ρ_c = prog_gm.ρ
     g::FT = ICP.grav(param_set)
     w_up_c = aux_tc.w_up_c
     w_en_c = aux_tc.w_en_c
@@ -204,7 +205,7 @@ function compute_entr_detr!(
                     a_up = aux_up[i].area[k], # updraft area fraction
                     a_en = aux_en.area[k], # environment area fraction
                     H_up = plume_scale_height[i], # plume scale height
-                    ref_H = p0_c[k] / (ρ0_c[k] * g), # reference state scale height
+                    ref_H = p_c[k] / (ρ_c[k] * g), # reference state scale height
                     RH_up = aux_up[i].RH[k], # updraft relative humidity
                     RH_en = aux_en.RH[k], # environment relative humidity
                     max_area = max_area, # maximum updraft area
@@ -263,8 +264,9 @@ function compute_entr_detr!(
     prog_gm_f = face_prog_grid_mean(state)
     aux_gm = center_aux_grid_mean(state)
     aux_tc = center_aux_turbconv(state)
-    p0_c = center_ref_state(state).p0
-    ρ0_c = center_ref_state(state).ρ0
+    prog_gm = center_prog_grid_mean(state)
+    p_c = aux_gm.p
+    ρ_c = prog_gm.ρ
     g::FT = ICP.grav(param_set)
     w_up_c = aux_tc.w_up_c
     w_en_c = aux_tc.w_en_c
@@ -312,7 +314,7 @@ function compute_entr_detr!(
                     a_up = aux_up[i].area[k], # updraft area fraction
                     a_en = aux_en.area[k], # environment area fraction
                     H_up = plume_scale_height[i], # plume scale height
-                    ref_H = p0_c[k] / (ρ0_c[k] * g), # reference state scale height
+                    ref_H = p_c[k] / (ρ_c[k] * g), # reference state scale height
                     RH_up = aux_up[i].RH[k], # updraft relative humidity
                     RH_en = aux_en.RH[k], # environment relative humidity
                     max_area = max_area, # maximum updraft area
diff --git a/src/closures/perturbation_pressure.jl b/src/closures/perturbation_pressure.jl
index 4fc4c0a5e2..88f77974cd 100644
--- a/src/closures/perturbation_pressure.jl
+++ b/src/closures/perturbation_pressure.jl
@@ -31,8 +31,9 @@ function compute_nh_pressure!(state::State, grid::Grid, edmf::EDMFModel, param_s
     ∇ = CCO.DivergenceC2F(; w_bcs...)
     aux_up = center_aux_updrafts(state)
     aux_up_f = face_aux_updrafts(state)
-    ρ0 = face_ref_state(state).ρ0
+    aux_gm_f = face_aux_grid_mean(state)
     aux_en_f = face_aux_environment(state)
+    ρ_f = aux_gm_f.ρ
     plume_scale_height = map(1:N_up) do i
         compute_plume_scale_height(grid, state, param_set, i)
     end
@@ -61,10 +62,10 @@ function compute_nh_pressure!(state::State, grid::Grid, edmf::EDMFModel, param_s
         nh_press_drag = aux_up_f[i].nh_pressure_drag
         nh_pressure = aux_up_f[i].nh_pressure
 
-        @. nh_press_buoy = Int(Ifa(a_up) > 0) * -α_b / (1 + α₂_asp_ratio²) * ρ0 * Ifa(a_up) * Ifb(b_up)
-        @. nh_press_adv = Int(Ifa(a_up) > 0) * ρ0 * Ifa(a_up) * α_a * w_up * ∇(wvec(Ifc(w_up)))
+        @. nh_press_buoy = Int(Ifa(a_up) > 0) * -α_b / (1 + α₂_asp_ratio²) * ρ_f * Ifa(a_up) * Ifb(b_up)
+        @. nh_press_adv = Int(Ifa(a_up) > 0) * ρ_f * Ifa(a_up) * α_a * w_up * ∇(wvec(Ifc(w_up)))
         # drag as w_dif and account for downdrafts
-        @. nh_press_drag = Int(Ifa(a_up) > 0) * -1 * ρ0 * Ifa(a_up) * α_d * (w_up - w_en) * abs(w_up - w_en) / H_up
+        @. nh_press_drag = Int(Ifa(a_up) > 0) * -1 * ρ_f * Ifa(a_up) * α_d * (w_up - w_en) * abs(w_up - w_en) / H_up
         @. nh_pressure = nh_press_buoy + nh_press_adv + nh_press_drag
     end
     return nothing
diff --git a/src/dycore_api.jl b/src/dycore_api.jl
index 110ca1cb97..801da4e5b3 100644
--- a/src/dycore_api.jl
+++ b/src/dycore_api.jl
@@ -38,12 +38,6 @@ tendencies(state, fl) = getproperty(state.tendencies, field_loc(fl))
 
 center_tendencies_grid_mean(state) = tendencies(state, CentField())
 
-""" Reference state fields for the host model """
-ref_state(state, fl) = aux(state, fl).ref_state
-
-face_ref_state(state) = ref_state(state, FaceField())
-center_ref_state(state) = ref_state(state, CentField())
-
 #####
 ##### TurbulenceConvection fields
 #####
diff --git a/src/microphysics_coupling.jl b/src/microphysics_coupling.jl
index e4a1aebef1..87b4d42529 100644
--- a/src/microphysics_coupling.jl
+++ b/src/microphysics_coupling.jl
@@ -2,7 +2,7 @@
 Computes the tendencies to qt and θ_liq_ice due to precipitation formation
 (autoconversion + accretion)
 """
-function noneq_moisture_sources(param_set::APS, area::FT, ρ0::FT, Δt::Real, ts) where {FT}
+function noneq_moisture_sources(param_set::APS, area::FT, ρ::FT, Δt::Real, ts) where {FT}
 
     # TODO - when using adaptive timestepping we are limiting the source terms
     #        with the previous timestep Δt
@@ -15,7 +15,7 @@ function noneq_moisture_sources(param_set::APS, area::FT, ρ0::FT, Δt::Real, ts
         q_vap = TD.vapor_specific_humidity(param_set, ts)
 
         # TODO - is that the state we want to be relaxing to?
-        ts_eq = TD.PhaseEquil_ρTq(param_set, ρ0, T, q.tot)
+        ts_eq = TD.PhaseEquil_ρTq(param_set, ρ, T, q.tot)
         q_eq = TD.PhasePartition(param_set, ts_eq)
 
         S_ql = CM1.conv_q_vap_to_q_liq_ice(param_set, liq_type, q_eq, q)
@@ -49,7 +49,7 @@ function precipitation_formation(
     qr::FT,
     qs::FT,
     area::FT,
-    ρ0::FT,
+    ρ::FT,
     Δt::Real,
     ts,
 ) where {FT}
@@ -108,21 +108,21 @@ function precipitation_formation(
             θ_liq_ice_tendency -= 1 / Π_m / c_pm * (L_v0 * S_qt_rain + L_s0 * S_qt_snow)
 
             # accretion cloud water + rain
-            S_qr = min(q.liq / Δt, CM1.accretion(param_set, liq_type, rain_type, q.liq, qr, ρ0))
+            S_qr = min(q.liq / Δt, CM1.accretion(param_set, liq_type, rain_type, q.liq, qr, ρ))
             qr_tendency += S_qr
             qt_tendency -= S_qr
             ql_tendency -= S_qr
             θ_liq_ice_tendency += S_qr / Π_m / c_pm * L_v0
 
             # accretion cloud ice + snow
-            S_qs = min(q.ice / Δt, CM1.accretion(param_set, ice_type, snow_type, q.ice, qs, ρ0))
+            S_qs = min(q.ice / Δt, CM1.accretion(param_set, ice_type, snow_type, q.ice, qs, ρ))
             qs_tendency += S_qs
             qt_tendency -= S_qs
             qi_tendency -= S_qs
             θ_liq_ice_tendency += S_qs / Π_m / c_pm * L_s0
 
             # sink of cloud water via accretion cloud water + snow
-            S_qt = -min(q.liq / Δt, CM1.accretion(param_set, liq_type, snow_type, q.liq, qs, ρ0))
+            S_qt = -min(q.liq / Δt, CM1.accretion(param_set, liq_type, snow_type, q.liq, qs, ρ))
             if T < T_fr # cloud droplets freeze to become snow)
                 qs_tendency -= S_qt
                 qt_tendency += S_qt
@@ -138,9 +138,9 @@ function precipitation_formation(
             end
 
             # sink of cloud ice via accretion cloud ice - rain
-            S_qt = -min(q.ice / Δt, CM1.accretion(param_set, ice_type, rain_type, q.ice, qr, ρ0))
+            S_qt = -min(q.ice / Δt, CM1.accretion(param_set, ice_type, rain_type, q.ice, qr, ρ))
             # sink of rain via accretion cloud ice - rain
-            S_qr = -min(qr / Δt, CM1.accretion_rain_sink(param_set, q.ice, qr, ρ0))
+            S_qr = -min(qr / Δt, CM1.accretion_rain_sink(param_set, q.ice, qr, ρ))
             qt_tendency += S_qt
             qi_tendency += S_qt
             qr_tendency += S_qr
@@ -149,9 +149,9 @@ function precipitation_formation(
 
             # accretion rain - snow
             if T < T_fr
-                S_qs = min(qr / Δt, CM1.accretion_snow_rain(param_set, snow_type, rain_type, qs, qr, ρ0))
+                S_qs = min(qr / Δt, CM1.accretion_snow_rain(param_set, snow_type, rain_type, qs, qr, ρ))
             else
-                S_qs = -min(qs / Δt, CM1.accretion_snow_rain(param_set, rain_type, snow_type, qr, qs, ρ0))
+                S_qs = -min(qs / Δt, CM1.accretion_snow_rain(param_set, rain_type, snow_type, qr, qs, ρ))
             end
             qs_tendency += S_qs
             qr_tendency -= S_qs
diff --git a/src/name_aliases.jl b/src/name_aliases.jl
index a2e3817900..83bb48885f 100644
--- a/src/name_aliases.jl
+++ b/src/name_aliases.jl
@@ -9,12 +9,10 @@ function name_aliases()
     dict = Dict(
         "zc" => ("z_half",),
         "zf" => ("z",),
-        "α0_c" => ("alpha0_half",),
-        "α0_f" => ("alpha0",),
-        "p0_c" => ("p0_half",),
-        "p0_f" => ("p0",),
-        "ρ0_c" => ("rho0_half",),
-        "ρ0_f" => ("rho0",),
+        "p_c" => ("p_half",),
+        "p_f" => ("p",),
+        "ρ_c" => ("rho0_half",),
+        "ρ_f" => ("rho0",),
         "updraft_area" => ("updraft_fraction",),
         "updraft_thetal" => ("updraft_thetali",),
         "thetal_mean" => ("thetali_mean", "theta_mean"),
diff --git a/src/types.jl b/src/types.jl
index c64d2c5f85..01719b44dd 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -120,7 +120,7 @@ Base.@kwdef struct EnvBuoyGrad{FT, EBC <: AbstractEnvBuoyGradClosure}
     "liquid ice potential temperature gradient in the saturated part"
     ∂θl∂z_sat::FT
     "reference pressure"
-    p0::FT
+    p::FT
     "cloud fraction"
     en_cld_frac::FT
     "specific volume"
@@ -163,7 +163,7 @@ Base.@kwdef struct MinDisspLen{FT}
     "turbulent Prandtl number"
     Pr::FT
     "reference pressure"
-    p0::FT
+    p::FT
     "vertical buoyancy gradient struct"
     ∇b::GradBuoy{FT}
     "env shear"
diff --git a/src/update_aux.jl b/src/update_aux.jl
index 964d1e64d8..48e498dc66 100644
--- a/src/update_aux.jl
+++ b/src/update_aux.jl
@@ -6,10 +6,6 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
     kc_surf = kc_surface(grid)
     kf_surf = kf_surface(grid)
     kc_toa = kc_top_of_atmos(grid)
-    ρ0_f = face_ref_state(state).ρ0
-    p0_c = center_ref_state(state).p0
-    ρ0_c = center_ref_state(state).ρ0
-    α0_c = center_ref_state(state).α0
     c_m = ICP.c_m(param_set)
     KM = center_aux_turbconv(state).KM
     KH = center_aux_turbconv(state).KH
@@ -21,6 +17,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
     aux_en = center_aux_environment(state)
     aux_en_f = face_aux_environment(state)
     aux_gm = center_aux_grid_mean(state)
+    aux_gm_f = face_aux_grid_mean(state)
     aux_tc_f = face_aux_turbconv(state)
     aux_tc = center_aux_turbconv(state)
     aux_bulk = center_aux_bulk(state)
@@ -28,6 +25,9 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
     aux_en_2m = center_aux_environment_2m(state)
     prog_up = center_prog_updrafts(state)
     prog_up_f = face_prog_updrafts(state)
+    ρ_f = aux_gm_f.ρ
+    p_c = aux_gm.p
+    ρ_c = prog_gm.ρ
     aux_en_unsat = aux_en.unsat
     aux_en_sat = aux_en.sat
     m_entr_detr = aux_tc.ϕ_temporary
@@ -47,7 +47,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
         #####
         @inbounds for i in 1:N_up
             if is_surface_center(grid, k)
-                if prog_up[i].ρarea[k] / ρ0_c[k] >= edmf.minimum_area
+                if prog_up[i].ρarea[k] / ρ_c[k] >= edmf.minimum_area
                     θ_surf = θ_surface_bc(surf, grid, state, edmf, i)
                     q_surf = q_surface_bc(surf, grid, state, edmf, i)
                     a_surf = area_surface_bc(surf, edmf, i)
@@ -59,10 +59,10 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
                     aux_up[i].q_tot[k] = aux_gm.q_tot[k]
                 end
             else
-                if prog_up[i].ρarea[k] / ρ0_c[k] >= edmf.minimum_area
+                if prog_up[i].ρarea[k] / ρ_c[k] >= edmf.minimum_area
                     aux_up[i].θ_liq_ice[k] = prog_up[i].ρaθ_liq_ice[k] / prog_up[i].ρarea[k]
                     aux_up[i].q_tot[k] = prog_up[i].ρaq_tot[k] / prog_up[i].ρarea[k]
-                    aux_up[i].area[k] = prog_up[i].ρarea[k] / ρ0_c[k]
+                    aux_up[i].area[k] = prog_up[i].ρarea[k] / ρ_c[k]
                 else
                     aux_up[i].θ_liq_ice[k] = aux_gm.θ_liq_ice[k]
                     aux_up[i].q_tot[k] = aux_gm.q_tot[k]
@@ -73,7 +73,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
         if edmf.moisture_model isa NonEquilibriumMoisture
             @inbounds for i in 1:N_up
                 if is_surface_center(grid, k)
-                    if prog_up[i].ρarea[k] / ρ0_c[k] >= edmf.minimum_area
+                    if prog_up[i].ρarea[k] / ρ_c[k] >= edmf.minimum_area
                         ql_surf = ql_surface_bc(surf)
                         qi_surf = qi_surface_bc(surf)
                         aux_up[i].q_liq[k] = ql_surf
@@ -83,7 +83,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
                         aux_up[i].q_ice[k] = prog_gm.q_ice[k]
                     end
                 else
-                    if prog_up[i].ρarea[k] / ρ0_c[k] >= edmf.minimum_area
+                    if prog_up[i].ρarea[k] / ρ_c[k] >= edmf.minimum_area
                         aux_up[i].q_liq[k] = prog_up[i].ρaq_liq[k] / prog_up[i].ρarea[k]
                         aux_up[i].q_ice[k] = prog_up[i].ρaq_ice[k] / prog_up[i].ρarea[k]
                     else
@@ -127,11 +127,11 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
             end
         end
         aux_en.area[k] = 1 - aux_bulk.area[k]
-        aux_en.tke[k] = prog_en.ρatke[k] / (ρ0_c[k] * aux_en.area[k])
+        aux_en.tke[k] = prog_en.ρatke[k] / (ρ_c[k] * aux_en.area[k])
         if edmf.thermo_covariance_model isa PrognosticThermoCovariances
-            aux_en.Hvar[k] = prog_en.ρaHvar[k] / (ρ0_c[k] * aux_en.area[k])
-            aux_en.QTvar[k] = prog_en.ρaQTvar[k] / (ρ0_c[k] * aux_en.area[k])
-            aux_en.HQTcov[k] = prog_en.ρaHQTcov[k] / (ρ0_c[k] * aux_en.area[k])
+            aux_en.Hvar[k] = prog_en.ρaHvar[k] / (ρ_c[k] * aux_en.area[k])
+            aux_en.QTvar[k] = prog_en.ρaQTvar[k] / (ρ_c[k] * aux_en.area[k])
+            aux_en.HQTcov[k] = prog_en.ρaHQTcov[k] / (ρ_c[k] * aux_en.area[k])
         end
 
         #####
@@ -157,7 +157,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
         else
             error("Something went wrong. The moisture_model options are equilibrium or nonequilibrium")
         end
-        ts_env[k] = thermo_state_pθq(param_set, p0_c[k], aux_en.θ_liq_ice[k], aux_en.q_tot[k], thermo_args...)
+        ts_env[k] = thermo_state_pθq(param_set, p_c[k], aux_en.θ_liq_ice[k], aux_en.q_tot[k], thermo_args...)
         ts_en = ts_env[k]
         aux_en.T[k] = TD.air_temperature(param_set, ts_en)
         aux_en.θ_virt[k] = TD.virtual_pottemp(param_set, ts_en)
@@ -165,7 +165,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
         aux_en.q_liq[k] = TD.liquid_specific_humidity(param_set, ts_en)
         aux_en.q_ice[k] = TD.ice_specific_humidity(param_set, ts_en)
         rho = TD.air_density(param_set, ts_en)
-        aux_en.buoy[k] = buoyancy_c(param_set, ρ0_c[k], rho)
+        aux_en.buoy[k] = buoyancy_c(param_set, ρ_c[k], rho)
 
         # update_sat_unsat
         if TD.has_condensate(param_set, ts_en)
@@ -189,21 +189,21 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
                 qt = aux_up[i].q_tot[k - 1]
                 h = aux_up[i].θ_liq_ice[k - 1]
                 if edmf.moisture_model isa EquilibriumMoisture
-                    ts_up = thermo_state_pθq(param_set, p0_c[k], h, qt)
+                    ts_up = thermo_state_pθq(param_set, p_c[k], h, qt)
                 elseif edmf.moisture_model isa NonEquilibriumMoisture
                     ql = aux_up[i].q_liq[k - 1]
                     qi = aux_up[i].q_ice[k - 1]
-                    ts_up = thermo_state_pθq(param_set, p0_c[k], h, qt, ql, qi)
+                    ts_up = thermo_state_pθq(param_set, p_c[k], h, qt, ql, qi)
                 else
                     error("Something went wrong. emdf.moisture_model options are equilibrium or nonequilibrium")
                 end
             else
                 if edmf.moisture_model isa EquilibriumMoisture
-                    ts_up = thermo_state_pθq(param_set, p0_c[k], aux_up[i].θ_liq_ice[k], aux_up[i].q_tot[k])
+                    ts_up = thermo_state_pθq(param_set, p_c[k], aux_up[i].θ_liq_ice[k], aux_up[i].q_tot[k])
                 elseif edmf.moisture_model isa NonEquilibriumMoisture
                     ts_up = thermo_state_pθq(
                         param_set,
-                        p0_c[k],
+                        p_c[k],
                         aux_up[i].θ_liq_ice[k],
                         aux_up[i].q_tot[k],
                         aux_up[i].q_liq[k],
@@ -217,7 +217,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
             aux_up[i].q_ice[k] = TD.ice_specific_humidity(param_set, ts_up)
             aux_up[i].T[k] = TD.air_temperature(param_set, ts_up)
             ρ = TD.air_density(param_set, ts_up)
-            aux_up[i].buoy[k] = buoyancy_c(param_set, ρ0_c[k], ρ)
+            aux_up[i].buoy[k] = buoyancy_c(param_set, ρ_c[k], ρ)
             aux_up[i].RH[k] = TD.relative_humidity(param_set, ts_up)
         end
         aux_gm.buoy[k] = (1.0 - aux_bulk.area[k]) * aux_en.buoy[k]
@@ -277,7 +277,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
         If = CCO.InterpolateC2F(; a_up_bcs...)
         a_min = edmf.minimum_area
         a_up = aux_up[i].area
-        @. aux_up_f[i].w = ifelse(If(a_up) >= a_min, max(prog_up_f[i].ρaw / (ρ0_f * If(a_up)), 0), FT(0))
+        @. aux_up_f[i].w = ifelse(If(a_up) >= a_min, max(prog_up_f[i].ρaw / (ρ_f * If(a_up)), 0), FT(0))
     end
     @inbounds for i in 1:N_up
         aux_up_f[i].w[kf_surf] = w_surface_bc(surf)
@@ -391,9 +391,9 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
                 ∂θv∂z_unsat = ∂θv∂z[k],
                 ∂qt∂z_sat = ∂qt∂z[k],
                 ∂θl∂z_sat = ∂θl∂z[k],
-                p0 = p0_c[k],
+                p = p_c[k],
                 en_cld_frac = aux_en.cloud_fraction[k],
-                alpha0 = α0_c[k],
+                alpha0 = 1 / ρ_c[k],
             )
             bg_model = EnvBuoyGrad(edmf.bg_closure; bg_kwargs...)
 
@@ -408,9 +408,9 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
                 ∂θv∂z_unsat = ∂θv∂z_unsat[k],
                 ∂qt∂z_sat = ∂qt∂z_sat[k],
                 ∂θl∂z_sat = ∂θl∂z_sat[k],
-                p0 = p0_c[k],
+                p = p_c[k],
                 en_cld_frac = aux_en.cloud_fraction[k],
-                alpha0 = α0_c[k],
+                alpha0 = 1 / ρ_c[k],
             )
             bg_model = EnvBuoyGrad(edmf.bg_closure; bg_kwargs...)
         else
@@ -429,7 +429,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
             tke_surf = aux_en.tke[kc_surf],
             ustar = surf.ustar,
             Pr = aux_tc.prandtl_nvec[k],
-            p0 = p0_c[k],
+            p = p_c[k],
             ∇b = bg,
             Shear² = Shear²[k],
             tke = aux_en.tke[k],
@@ -444,7 +444,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
         KM[k] = c_m * ml.mixing_length * sqrt(max(aux_en.tke[k], 0.0))
         KH[k] = KM[k] / aux_tc.prandtl_nvec[k]
 
-        aux_en_2m.tke.buoy[k] = -aux_en.area[k] * ρ0_c[k] * KH[k] * bg.∂b∂z
+        aux_en_2m.tke.buoy[k] = -aux_en.area[k] * ρ_c[k] * KH[k] * bg.∂b∂z
     end
 
 
@@ -476,9 +476,9 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
     # TODO defined again in compute_covariance_shear and compute_covaraince
     @inbounds for k in real_center_indices(grid)
         aux_en_2m.tke.rain_src[k] = 0
-        aux_en_2m.Hvar.rain_src[k] = ρ0_c[k] * aux_en.area[k] * 2 * aux_en.Hvar_rain_dt[k]
-        aux_en_2m.QTvar.rain_src[k] = ρ0_c[k] * aux_en.area[k] * 2 * aux_en.QTvar_rain_dt[k]
-        aux_en_2m.HQTcov.rain_src[k] = ρ0_c[k] * aux_en.area[k] * aux_en.HQTcov_rain_dt[k]
+        aux_en_2m.Hvar.rain_src[k] = ρ_c[k] * aux_en.area[k] * 2 * aux_en.Hvar_rain_dt[k]
+        aux_en_2m.QTvar.rain_src[k] = ρ_c[k] * aux_en.area[k] * 2 * aux_en.QTvar_rain_dt[k]
+        aux_en_2m.HQTcov.rain_src[k] = ρ_c[k] * aux_en.area[k] * aux_en.HQTcov_rain_dt[k]
     end
 
     get_GMV_CoVar(edmf, grid, state, Val(:tke), Val(:w), Val(:w))
@@ -496,8 +496,8 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
         prog_pr = center_prog_precipitation(state)
 
         @inbounds for k in real_center_indices(grid)
-            term_vel_rain[k] = CM1.terminal_velocity(param_set, CM1.RainType(), ρ0_c[k], prog_pr.q_rai[k])
-            term_vel_snow[k] = CM1.terminal_velocity(param_set, CM1.SnowType(), ρ0_c[k], prog_pr.q_sno[k])
+            term_vel_rain[k] = CM1.terminal_velocity(param_set, CM1.RainType(), ρ_c[k], prog_pr.q_rai[k])
+            term_vel_snow[k] = CM1.terminal_velocity(param_set, CM1.SnowType(), ρ_c[k], prog_pr.q_sno[k])
         end
     end
 
@@ -508,7 +508,8 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
         zLL = grid.zc[kc_surf].z
         ustar = surf.ustar
         oblength = surf.obukhov_length
-        α0LL = center_ref_state(state).α0[kc_surf]
+        prog_gm = center_prog_grid_mean(state)
+        ρLL = prog_gm.ρ[kc_surf]
         update_diagnostic_covariances!(edmf, grid, state, param_set, Val(:Hvar))
         update_diagnostic_covariances!(edmf, grid, state, param_set, Val(:QTvar))
         update_diagnostic_covariances!(edmf, grid, state, param_set, Val(:HQTcov))
@@ -519,9 +520,9 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
             aux_en.HQTcov[k] = min(aux_en.HQTcov[k], sqrt(aux_en.Hvar[k] * aux_en.QTvar[k]))
         end
         ae_surf = 1 - aux_bulk.area[kc_surf]
-        aux_en.Hvar[kc_surf] = ae_surf * get_surface_variance(flux1 * α0LL, flux1 * α0LL, ustar, zLL, oblength)
-        aux_en.QTvar[kc_surf] = ae_surf * get_surface_variance(flux2 * α0LL, flux2 * α0LL, ustar, zLL, oblength)
-        aux_en.HQTcov[kc_surf] = ae_surf * get_surface_variance(flux1 * α0LL, flux2 * α0LL, ustar, zLL, oblength)
+        aux_en.Hvar[kc_surf] = ae_surf * get_surface_variance(flux1 / ρLL, flux1 / ρLL, ustar, zLL, oblength)
+        aux_en.QTvar[kc_surf] = ae_surf * get_surface_variance(flux2 / ρLL, flux2 / ρLL, ustar, zLL, oblength)
+        aux_en.HQTcov[kc_surf] = ae_surf * get_surface_variance(flux1 / ρLL, flux2 / ρLL, ustar, zLL, oblength)
     end
     return nothing
 end