remove ref_state, move ref variables to grid mean

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

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

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)

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,27 +311,27 @@ 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
 
     # TODO: Move rho_cloud_liq to CLIMAParameters
     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