Make more variables ClimaCore fields

integration_tests/utils/main.jl

@@ -154,6 +154,10 @@ cent_aux_vars_edmf(FT, n_up) = (;
         θ_virt = FT(0),
         ∇ρ_ae_K∇ϕ = FT(0),
         ρaew_en_ϕ = FT(0),
+        massflux_tendency_h = FT(0),
+        massflux_tendency_qt = FT(0),
+        diffusive_tendency_h = FT(0),
+        diffusive_tendency_qt = FT(0),
     ),
 )
 cent_aux_vars(FT, n_up) = (; aux_vars_ref_state(FT)..., cent_aux_vars_gm(FT)..., cent_aux_vars_edmf(FT, n_up)...)
@@ -178,6 +182,12 @@ face_aux_vars_edmf(FT, n_up) = (;
         ρ_ae_K∇ϕ = FT(0),
         en = (; w = FT(0)),
         up = ntuple(i -> face_aux_vars_up(FT), n_up),
+        massflux_h = FT(0),
+        massflux_qt = FT(0),
+        diffusive_flux_h = FT(0),
+        diffusive_flux_qt = FT(0),
+        diffusive_flux_u = FT(0),
+        diffusive_flux_v = FT(0),
     ),
 )
 face_aux_vars(FT, n_up) = (; aux_vars_ref_state(FT)..., face_aux_vars_gm(FT)..., face_aux_vars_edmf(FT, n_up)...)

src/Turbulence_PrognosticTKE.jl

@@ -25,6 +25,8 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm,
     aux_en_f = face_aux_environment(state)
     aux_up = center_aux_updrafts(state)
     aux_bulk = center_aux_bulk(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
@@ -100,9 +102,9 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm,
             aux_tc.θ_liq_ice_tendency_precip_sinks[k]
     end
 
-    aux_up_f = face_aux_updrafts(state)
-    edmf.massflux_h .= 0.0
-    edmf.massflux_qt .= 0.0
+    # TODO: we shouldn't need to call parent here
+    parent(aux_tc_f.massflux_h) .= 0
+    parent(aux_tc_f.massflux_qt) .= 0
     # Compute the mass flux and associated scalar fluxes
     @inbounds for i in 1:(up.n_updrafts)
         aux_up_f[i].massflux[kf_surf] = 0.0
@@ -115,8 +117,8 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm,
     end
 
     @inbounds for k in real_face_indices(grid)
-        edmf.massflux_h[k] = 0.0
-        edmf.massflux_qt[k] = 0.0
+        aux_tc_f.massflux_h[k] = 0.0
+        aux_tc_f.massflux_qt[k] = 0.0
         # We know that, since W = 0 at z = 0, m = 0 also, and
         # therefore θ_liq_ice / q_tot values do not matter
         m_bcs = (; bottom = SetValue(0), top = SetValue(0))
@@ -125,16 +127,16 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm,
         @inbounds for i in 1:(up.n_updrafts)
             h_up_f = interpc2f(aux_up[i].θ_liq_ice, grid, k; m_bcs...)
             qt_up_f = interpc2f(aux_up[i].q_tot, grid, k; m_bcs...)
-            edmf.massflux_h[k] += aux_up_f[i].massflux[k] * (h_up_f - h_en_f)
-            edmf.massflux_qt[k] += aux_up_f[i].massflux[k] * (qt_up_f - qt_en_f)
+            aux_tc_f.massflux_h[k] += aux_up_f[i].massflux[k] * (h_up_f - h_en_f)
+            aux_tc_f.massflux_qt[k] += aux_up_f[i].massflux[k] * (qt_up_f - qt_en_f)
         end
     end
 
     # Compute the  mass flux tendencies
     # Adjust the values of the grid mean variables
     @inbounds for k in real_center_indices(grid)
-        mf_tend_h_dual = dual_faces(edmf.massflux_h, grid, k)
-        mf_tend_qt_dual = dual_faces(edmf.massflux_qt, grid, k)
+        mf_tend_h_dual = dual_faces(aux_tc_f.massflux_h, grid, k)
+        mf_tend_qt_dual = dual_faces(aux_tc_f.massflux_qt, grid, k)
 
         ∇mf_tend_h = ∇f2c(mf_tend_h_dual, grid, k)
         ∇mf_tend_qt = ∇f2c(mf_tend_qt_dual, grid, k)
@@ -143,8 +145,8 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm,
         mf_tend_qt = -∇mf_tend_qt * α0_c[k]
 
         # Prepare the output
-        edmf.massflux_tendency_h[k] = mf_tend_h
-        edmf.massflux_tendency_qt[k] = mf_tend_qt
+        aux_tc.massflux_tendency_h[k] = mf_tend_h
+        aux_tc.massflux_tendency_qt[k] = mf_tend_qt
         tendencies_gm.θ_liq_ice[k] += mf_tend_h
         tendencies_gm.q_tot[k] += mf_tend_qt
     end
@@ -154,19 +156,19 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm,
     aeKMu_bc = Case.Sur.rho_uflux / aux_en.area[kc_surf]
     aeKMv_bc = Case.Sur.rho_vflux / aux_en.area[kc_surf]
     @inbounds for k in real_center_indices(grid)
-        aeKH_q_tot_cut = dual_faces(edmf.diffusive_flux_qt, grid, k)
+        aeKH_q_tot_cut = dual_faces(aux_tc_f.diffusive_flux_qt, grid, k)
         ∇aeKH_q_tot = ∇f2c(aeKH_q_tot_cut, grid, k; bottom = SetValue(aeKHq_tot_bc), top = SetValue(0))
         tendencies_gm.q_tot[k] += -α0_c[k] * aux_en.area[k] * ∇aeKH_q_tot
 
-        aeKH_θ_liq_ice_cut = dual_faces(edmf.diffusive_flux_h, grid, k)
+        aeKH_θ_liq_ice_cut = dual_faces(aux_tc_f.diffusive_flux_h, grid, k)
         ∇aeKH_θ_liq_ice = ∇f2c(aeKH_θ_liq_ice_cut, grid, k; bottom = SetValue(aeKHθ_liq_ice_bc), top = SetValue(0))
         tendencies_gm.θ_liq_ice[k] += -α0_c[k] * aux_en.area[k] * ∇aeKH_θ_liq_ice
 
-        aeKM_u_cut = dual_faces(edmf.diffusive_flux_u, grid, k)
+        aeKM_u_cut = dual_faces(aux_tc_f.diffusive_flux_u, grid, k)
         ∇aeKM_u = ∇f2c(aeKM_u_cut, grid, k; bottom = SetValue(aeKMu_bc), top = SetValue(0))
         tendencies_gm.u[k] += -α0_c[k] * aux_en.area[k] * ∇aeKM_u
 
-        aeKM_v_cut = dual_faces(edmf.diffusive_flux_v, grid, k)
+        aeKM_v_cut = dual_faces(aux_tc_f.diffusive_flux_v, grid, k)
         ∇aeKM_v = ∇f2c(aeKM_v_cut, grid, k; bottom = SetValue(aeKMv_bc), top = SetValue(0))
         tendencies_gm.v[k] += -α0_c[k] * aux_en.area[k] * ∇aeKM_v
     end
@@ -210,19 +212,19 @@ function compute_diffusive_fluxes(
     @inbounds for k in real_face_indices(grid)
         q_dual = dual_centers(aux_en.q_tot, grid, k)
         ∇q_tot_f = ∇c2f(q_dual, grid, k; bottom = SetGradient(aeKHq_tot_bc), top = SetGradient(0))
-        edmf.diffusive_flux_qt[k] = -aux_tc_f.ρ_ae_KH[k] * ∇q_tot_f
+        aux_tc_f.diffusive_flux_qt[k] = -aux_tc_f.ρ_ae_KH[k] * ∇q_tot_f
 
         θ_liq_ice_dual = dual_centers(aux_en.θ_liq_ice, grid, k)
         ∇θ_liq_ice_f = ∇c2f(θ_liq_ice_dual, grid, k; bottom = SetGradient(aeKHθ_liq_ice_bc), top = SetGradient(0))
-        edmf.diffusive_flux_h[k] = -aux_tc_f.ρ_ae_KH[k] * ∇θ_liq_ice_f
+        aux_tc_f.diffusive_flux_h[k] = -aux_tc_f.ρ_ae_KH[k] * ∇θ_liq_ice_f
 
         u_dual = dual_centers(prog_gm.u, grid, k)
         ∇u_f = ∇c2f(u_dual, grid, k; bottom = SetGradient(aeKMu_bc), top = SetGradient(0))
-        edmf.diffusive_flux_u[k] = -aux_tc_f.ρ_ae_KM[k] * ∇u_f
+        aux_tc_f.diffusive_flux_u[k] = -aux_tc_f.ρ_ae_KM[k] * ∇u_f
 
         v_dual = dual_centers(prog_gm.v, grid, k)
         ∇v_f = ∇c2f(v_dual, grid, k; bottom = SetGradient(aeKMv_bc), top = SetGradient(0))
-        edmf.diffusive_flux_v[k] = -aux_tc_f.ρ_ae_KM[k] * ∇v_f
+        aux_tc_f.diffusive_flux_v[k] = -aux_tc_f.ρ_ae_KM[k] * ∇v_f
     end
     return
 end

src/diagnostics.jl

@@ -131,6 +131,14 @@ function io_dictionary_aux(state)
         "updraft_qt_precip" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_bulk(state).qt_tendency_precip_formation),
         "updraft_thetal_precip" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_bulk(state).θ_liq_ice_tendency_precip_formation),
 
+        "massflux_tendency_h" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_turbconv(state).massflux_tendency_h),
+        "massflux_tendency_qt" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_turbconv(state).massflux_tendency_qt),
+        "diffusive_tendency_h" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_turbconv(state).diffusive_tendency_h),
+        "diffusive_tendency_qt" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_turbconv(state).diffusive_tendency_qt),
+
+        "total_flux_h" => (; dims = ("zf", "t"), group = "profiles", field = face_aux_turbconv(state).diffusive_flux_h .+ face_aux_turbconv(state).massflux_h),
+        "total_flux_qt" => (; dims = ("zf", "t"), group = "profiles", field = face_aux_turbconv(state).diffusive_flux_qt .+ face_aux_turbconv(state).massflux_qt),
+
     )
     return io_dict
 end

src/io.jl

@@ -114,22 +114,6 @@ function initialize_io(edmf::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats)
     initialize_io(edmf.Precip, Stats)
 
     add_ts(Stats, "rd")
-    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_h")
-    add_profile(Stats, "massflux_qt")
-    add_profile(Stats, "massflux_tendency_h")
-    add_profile(Stats, "massflux_tendency_qt")
-    add_profile(Stats, "diffusive_flux_h")
-    add_profile(Stats, "diffusive_flux_u")
-    add_profile(Stats, "diffusive_flux_v")
-    add_profile(Stats, "diffusive_flux_qt")
-    add_profile(Stats, "diffusive_tendency_h")
-    add_profile(Stats, "diffusive_tendency_qt")
-    add_profile(Stats, "total_flux_h")
-    add_profile(Stats, "total_flux_qt")
     # Diff mixing lengths: Ignacio
     add_profile(Stats, "ed_length_scheme")
     add_profile(Stats, "mixing_length_ratio")
@@ -142,19 +126,6 @@ function io(edmf::EDMF_PrognosticTKE, grid, state, Stats::NetCDFIO_Stats, TS::Ti
     io(edmf.EnvVar, grid, state, Stats)
     io(edmf.Precip, grid, state, Stats)
     write_ts(Stats, "rd", StatsBase.mean(edmf.pressure_plume_spacing))
-    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)
-    write_profile(Stats, "massflux_tendency_qt", edmf.massflux_tendency_qt)
-    write_profile(Stats, "diffusive_flux_h", edmf.diffusive_flux_h)
-    write_profile(Stats, "diffusive_flux_qt", edmf.diffusive_flux_qt)
-    write_profile(Stats, "diffusive_flux_u", edmf.diffusive_flux_u)
-    write_profile(Stats, "diffusive_flux_v", edmf.diffusive_flux_v)
-    write_profile(Stats, "diffusive_tendency_h", edmf.diffusive_tendency_h)
-    write_profile(Stats, "diffusive_tendency_qt", edmf.diffusive_tendency_qt)
-    write_profile(Stats, "total_flux_h", edmf.massflux_h .+ edmf.diffusive_flux_h)
-    write_profile(Stats, "total_flux_qt", edmf.massflux_qt .+ edmf.diffusive_flux_qt)
-
     #Different mixing lengths : Ignacio
     write_profile(Stats, "ed_length_scheme", edmf.mls)
     write_profile(Stats, "mixing_length_ratio", edmf.ml_ratio)

src/types.jl

@@ -443,17 +443,6 @@ mutable struct EDMF_PrognosticTKE{A1}
     h_surface_bc::A1
     qt_surface_bc::A1
     pressure_plume_spacing::A1
-    massflux_tendency_h::A1
-    massflux_tendency_qt::A1
-    diffusive_tendency_h::A1
-    diffusive_tendency_qt::A1
-    massflux_h::A1
-    massflux_qt::A1
-    diffusive_flux_h::A1
-    diffusive_flux_qt::A1
-    diffusive_flux_u::A1
-    diffusive_flux_v::A1
-    massflux_tke::A1 # remove
     prandtl_nvec::A1
     prandtl_number::Float64 # remove
     mls::A1
@@ -527,23 +516,6 @@ mutable struct EDMF_PrognosticTKE{A1}
         qt_surface_bc = zeros(n_updrafts)
         pressure_plume_spacing = zeros(n_updrafts)
 
-        # Mass flux tendencies of mean scalars (for output)
-        massflux_tendency_h = center_field(grid)
-        massflux_tendency_qt = center_field(grid)
-
-        # (Eddy) diffusive tendencies of mean scalars (for output)
-        diffusive_tendency_h = center_field(grid)
-        diffusive_tendency_qt = center_field(grid)
-
-        # Vertical fluxes for output
-        massflux_h = face_field(grid)
-        massflux_qt = face_field(grid)
-        diffusive_flux_h = face_field(grid)
-        diffusive_flux_qt = face_field(grid)
-        diffusive_flux_u = face_field(grid)
-        diffusive_flux_v = face_field(grid)
-        massflux_tke = center_field(grid)
-
         # Initialize SDE parameters
         dt = parse_namelist(namelist, "time_stepping", "dt"; default = 1.0)
         closure = parse_namelist(
@@ -592,7 +564,7 @@ mutable struct EDMF_PrognosticTKE{A1}
         wstar = 0
         entr_surface_bc = 0
         detr_surface_bc = 0
-        A1 = typeof(massflux_tendency_h)
+        A1 = typeof(area_surface_bc)
         return new{A1}(
             Ri_bulk_crit,
             zi,
@@ -616,17 +588,6 @@ mutable struct EDMF_PrognosticTKE{A1}
             h_surface_bc,
             qt_surface_bc,
             pressure_plume_spacing,
-            massflux_tendency_h,
-            massflux_tendency_qt,
-            diffusive_tendency_h,
-            diffusive_tendency_qt,
-            massflux_h,
-            massflux_qt,
-            diffusive_flux_h,
-            diffusive_flux_qt,
-            diffusive_flux_u,
-            diffusive_flux_v,
-            massflux_tke,
             prandtl_nvec,
             prandtl_number,
             mls,