Make some diagnostic vars ClimaCore fields

integration_tests/utils/main.jl

@@ -146,12 +146,26 @@ face_aux_vars(FT, n_up) = (; aux_vars_ref_state(FT)..., face_aux_vars_gm(FT)...,
 
 # Center only
 cent_diagnostic_vars_gm(FT) = ()
-cent_diagnostic_vars_edmf(FT, n_up) = ()
+cent_diagnostic_vars_edmf(FT, n_up) = (
+    asp_ratio = FT(0),
+    entr_sc = FT(0),
+    detr_sc = FT(0),
+    massflux = FT(0),
+    frac_turb_entr = FT(0),
+    horiz_K_eddy = FT(0),
+    sorting_function = FT(0),
+    b_mix = FT(0),
+)
 cent_diagnostic_vars(FT, n_up) = (; cent_diagnostic_vars_gm(FT)..., cent_diagnostic_vars_edmf(FT, n_up)...)
 
 # Face only
 face_diagnostic_vars_gm(FT) = ()
-face_diagnostic_vars_edmf(FT, n_up) = ()
+face_diagnostic_vars_edmf(FT, n_up) = (
+    nh_pressure = FT(0),
+    nh_pressure_adv = FT(0),
+    nh_pressure_drag = FT(0),
+    nh_pressure_b = FT(0),
+)
 face_diagnostic_vars(FT, n_up) = (; face_diagnostic_vars_gm(FT)..., face_diagnostic_vars_edmf(FT, n_up)...)
 
 #####

src/diagnostics.jl

@@ -112,6 +112,17 @@ function io_dictionary_aux(state)
 
         "qr_mean" => (; dims = ("zc", "t"), group = "profiles", field = center_prog_precipitation(state).qr),
 
+        "nh_pressure" => (; dims = ("zf", "t"), group = "profiles", field = face_diagnostics_turbconv(state).nh_pressure),
+        "nh_pressure_adv" => (; dims = ("zf", "t"), group = "profiles", field = face_diagnostics_turbconv(state).nh_pressure_adv),
+        "nh_pressure_drag" => (; dims = ("zf", "t"), group = "profiles", field = face_diagnostics_turbconv(state).nh_pressure_drag),
+        "nh_pressure_b" => (; dims = ("zf", "t"), group = "profiles", field = face_diagnostics_turbconv(state).nh_pressure_b),
+        "turbulent_entrainment" => (; dims = ("zc", "t"), group = "profiles", field = center_diagnostics_turbconv(state).frac_turb_entr),
+        "horiz_K_eddy" => (; dims = ("zc", "t"), group = "profiles", field = center_diagnostics_turbconv(state).horiz_K_eddy),
+        "entrainment_sc" => (; dims = ("zc", "t"), group = "profiles", field = center_diagnostics_turbconv(state).entr_sc),
+        "detrainment_sc" => (; dims = ("zc", "t"), group = "profiles", field = center_diagnostics_turbconv(state).detr_sc),
+        "asp_ratio" => (; dims = ("zc", "t"), group = "profiles", field = center_diagnostics_turbconv(state).asp_ratio),
+        "massflux" => (; dims = ("zc", "t"), group = "profiles", field = center_diagnostics_turbconv(state).massflux),
+
         "updraft_cloud_fraction" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_turbconv(state).bulk.cloud_fraction),
 
     )
@@ -154,6 +165,7 @@ function compute_diagnostics!(edmf, gm, grid, state, Case, TS)
     aux_tc_f = face_aux_turbconv(state)
     aux_gm_f = face_aux_grid_mean(state)
     prog_pr = center_prog_precipitation(state)
+    a_up_bulk = center_aux_turbconv(state).bulk.area
     kc_toa = kc_top_of_atmos(grid)
     gm.cloud_base = grid.zc[kc_toa]
     gm.cloud_top = 0.0
@@ -165,6 +177,8 @@ function compute_diagnostics!(edmf, gm, grid, state, Case, TS)
     en_thermo = edmf.EnvThermo
     up_thermo = edmf.UpdThermo
     n_updrafts = up.n_updrafts
+    diag_tc = center_diagnostics_turbconv(state)
+    diag_tc_f = face_diagnostics_turbconv(state)
 
     @inbounds for k in real_center_indices(grid)
         gm.lwp += ρ0_c[k] * aux_gm.q_liq[k] * grid.Δz
@@ -256,5 +270,41 @@ function compute_diagnostics!(edmf, gm, grid, state, Case, TS)
         end
     end
 
+    @inbounds for k in real_center_indices(grid)
+        if a_up_bulk[k] > 0.0
+            @inbounds for i in 1:(edmf.n_updrafts)
+                diag_tc.massflux[k] += interpf2c(edmf.m, grid, k, i)
+                diag_tc.entr_sc[k] += aux_up[i].area[k] * edmf.entr_sc[i, k] / a_up_bulk[k]
+                diag_tc.detr_sc[k] += aux_up[i].area[k] * edmf.detr_sc[i, k] / a_up_bulk[k]
+                diag_tc.asp_ratio[k] += aux_up[i].area[k] * edmf.asp_ratio[i, k] / a_up_bulk[k]
+                diag_tc.frac_turb_entr[k] += aux_up[i].area[k] * edmf.frac_turb_entr[i, k] / a_up_bulk[k]
+                diag_tc.horiz_K_eddy[k] += aux_up[i].area[k] * edmf.horiz_K_eddy[i, k] / a_up_bulk[k]
+                diag_tc.sorting_function[k] += aux_up[i].area[k] * edmf.sorting_function[i, k] / a_up_bulk[k]
+                diag_tc.b_mix[k] += aux_up[i].area[k] * edmf.b_mix[i, k] / a_up_bulk[k]
+            end
+        end
+    end
+
+    @inbounds for k in real_face_indices(grid)
+        a_up_bulk_f =
+            interpc2f(a_up_bulk, grid, k; bottom = SetValue(sum(edmf.area_surface_bc)), top = SetZeroGradient())
+        if a_up_bulk_f > 0.0
+            @inbounds for i in 1:(edmf.n_updrafts)
+                a_up_f = interpc2f(
+                    aux_up[i].area,
+                    grid,
+                    k;
+                    bottom = SetValue(edmf.area_surface_bc[i]),
+                    top = SetZeroGradient(),
+                )
+                diag_tc.nh_pressure[k] += a_up_f * edmf.nh_pressure[i, k] / a_up_bulk_f
+                diag_tc.nh_pressure_b[k] += a_up_f * edmf.nh_pressure_b[i, k] / a_up_bulk_f
+                diag_tc.nh_pressure_adv[k] += a_up_f * edmf.nh_pressure_adv[i, k] / a_up_bulk_f
+                diag_tc.nh_pressure_drag[k] += a_up_f * edmf.nh_pressure_drag[i, k] / a_up_bulk_f
+            end
+        end
+    end
+
+
     return
 end

src/dycore_api.jl

@@ -36,6 +36,12 @@ tendencies(state, fl) = getproperty(state.tendencies, field_loc(fl))
 
 center_tendencies_grid_mean(state) = tendencies(state, CentField())
 
+""" Purely diagnostic fields for the host model """
+diagnostics(state, fl) = getproperty(state.diagnostics, field_loc(fl))
+
+center_diagnostics_grid_mean(state) = diagnostics(state, CentField())
+center_diagnostics_turbconv(state) = diagnostics(state, CentField()).turbconv
+
 """ Reference state fields for the host model """
 ref_state(state, fl) = aux(state, fl).ref_state
 

src/io.jl

@@ -111,23 +111,12 @@ function initialize_io(edmf::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats)
     initialize_io(edmf.EnvVar, Stats)
     initialize_io(edmf.Precip, Stats)
 
-    add_profile(Stats, "entrainment_sc")
-    add_profile(Stats, "detrainment_sc")
-    add_profile(Stats, "nh_pressure")
-    add_profile(Stats, "nh_pressure_adv")
-    add_profile(Stats, "nh_pressure_drag")
-    add_profile(Stats, "nh_pressure_b")
-    add_profile(Stats, "asp_ratio")
-    add_profile(Stats, "horiz_K_eddy")
     add_profile(Stats, "sorting_function")
-    add_profile(Stats, "b_mix")
     add_ts(Stats, "rd")
-    add_profile(Stats, "turbulent_entrainment")
     add_profile(Stats, "turbulent_entrainment_full")
     add_profile(Stats, "turbulent_entrainment_W")
     add_profile(Stats, "turbulent_entrainment_H")
     add_profile(Stats, "turbulent_entrainment_QT")
-    add_profile(Stats, "massflux")
     add_profile(Stats, "massflux_h")
     add_profile(Stats, "massflux_qt")
     add_profile(Stats, "massflux_tendency_h")
@@ -151,75 +140,10 @@ function initialize_io(edmf::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats)
 end
 
 function io(edmf::EDMF_PrognosticTKE, grid, state, Stats::NetCDFIO_Stats, TS::TimeStepping, param_set)
-
-    mean_nh_pressure = face_field(grid)
-    mean_nh_pressure_adv = face_field(grid)
-    mean_nh_pressure_drag = face_field(grid)
-    mean_nh_pressure_b = face_field(grid)
-
-    mean_asp_ratio = center_field(grid)
-    mean_entr_sc = center_field(grid)
-    mean_detr_sc = center_field(grid)
-    massflux = center_field(grid)
-    mean_frac_turb_entr = center_field(grid)
-    mean_horiz_K_eddy = center_field(grid)
-    mean_sorting_function = center_field(grid)
-    mean_b_mix = center_field(grid)
-
     io(edmf.UpdVar, grid, state, Stats)
     io(edmf.EnvVar, grid, state, Stats)
     io(edmf.Precip, grid, state, Stats)
-
-    aux_up = center_aux_updrafts(state)
-    a_up_bulk = center_aux_turbconv(state).bulk.area
-
     write_ts(Stats, "rd", StatsBase.mean(edmf.pressure_plume_spacing))
-
-    @inbounds for k in real_center_indices(grid)
-        if a_up_bulk[k] > 0.0
-            @inbounds for i in 1:(edmf.n_updrafts)
-                massflux[k] += interpf2c(edmf.m, grid, k, i)
-                mean_entr_sc[k] += aux_up[i].area[k] * edmf.entr_sc[i, k] / a_up_bulk[k]
-                mean_detr_sc[k] += aux_up[i].area[k] * edmf.detr_sc[i, k] / a_up_bulk[k]
-                mean_asp_ratio[k] += aux_up[i].area[k] * edmf.asp_ratio[i, k] / a_up_bulk[k]
-                mean_frac_turb_entr[k] += aux_up[i].area[k] * edmf.frac_turb_entr[i, k] / a_up_bulk[k]
-                mean_horiz_K_eddy[k] += aux_up[i].area[k] * edmf.horiz_K_eddy[i, k] / a_up_bulk[k]
-                mean_sorting_function[k] += aux_up[i].area[k] * edmf.sorting_function[i, k] / a_up_bulk[k]
-                mean_b_mix[k] += aux_up[i].area[k] * edmf.b_mix[i, k] / a_up_bulk[k]
-            end
-        end
-    end
-
-    @inbounds for k in real_face_indices(grid)
-        a_up_bulk_f =
-            interpc2f(a_up_bulk, grid, k; bottom = SetValue(sum(edmf.area_surface_bc)), top = SetZeroGradient())
-        if a_up_bulk_f > 0.0
-            @inbounds for i in 1:(edmf.n_updrafts)
-                a_up_f = interpc2f(
-                    aux_up[i].area,
-                    grid,
-                    k;
-                    bottom = SetValue(edmf.area_surface_bc[i]),
-                    top = SetZeroGradient(),
-                )
-                mean_nh_pressure[k] += a_up_f * edmf.nh_pressure[i, k] / a_up_bulk_f
-                mean_nh_pressure_b[k] += a_up_f * edmf.nh_pressure_b[i, k] / a_up_bulk_f
-                mean_nh_pressure_adv[k] += a_up_f * edmf.nh_pressure_adv[i, k] / a_up_bulk_f
-                mean_nh_pressure_drag[k] += a_up_f * edmf.nh_pressure_drag[i, k] / a_up_bulk_f
-            end
-        end
-    end
-
-    write_profile(Stats, "turbulent_entrainment", mean_frac_turb_entr)
-    write_profile(Stats, "horiz_K_eddy", mean_horiz_K_eddy)
-    write_profile(Stats, "entrainment_sc", mean_entr_sc)
-    write_profile(Stats, "detrainment_sc", mean_detr_sc)
-    write_profile(Stats, "nh_pressure", mean_nh_pressure)
-    write_profile(Stats, "nh_pressure_adv", mean_nh_pressure_adv)
-    write_profile(Stats, "nh_pressure_drag", mean_nh_pressure_drag)
-    write_profile(Stats, "nh_pressure_b", mean_nh_pressure_b)
-    write_profile(Stats, "asp_ratio", mean_asp_ratio)
-    write_profile(Stats, "massflux", massflux)
     write_profile(Stats, "massflux_h", edmf.massflux_h)
     write_profile(Stats, "massflux_qt", edmf.massflux_qt)
     write_profile(Stats, "massflux_tendency_h", edmf.massflux_tendency_h)