Modify Pi groups and include in stats file output

driver/compute_diagnostics.jl

@@ -55,6 +55,12 @@ function io_dictionary_diagnostics()
         "massflux" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_diagnostics_turbconv(state).massflux),
         "rain_flux" => (; dims = ("zf", "t"), group = "profiles", field = state -> face_diagnostics_precip(state).rain_flux),
         "snow_flux" => (; dims = ("zf", "t"), group = "profiles", field = state -> face_diagnostics_precip(state).snow_flux,),
+        "pi_1" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_diagnostics_turbconv(state).Π₁),
+        "pi_2" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_diagnostics_turbconv(state).Π₂),
+        "pi_3" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_diagnostics_turbconv(state).Π₃),
+        "pi_4" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_diagnostics_turbconv(state).Π₄),
+        "pi_5" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_diagnostics_turbconv(state).Π₅),
+        "pi_6" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_diagnostics_turbconv(state).Π₆),
     )
     return io_dict
 end
@@ -273,7 +279,7 @@ function compute_diagnostics!(
     @inbounds for i in 1:N_up
         @. diag_tc.massflux += Ic(aux_up_f[i].massflux)
     end
-
+    pi_subset = TC.entrainment_Π_subset(edmf)
     @inbounds for k in TC.real_center_indices(grid)
         a_up_bulk_k = a_up_bulk[k]
         diag_tc.entr_sc[k] = 0
@@ -286,6 +292,12 @@ function compute_diagnostics!(
         diag_tc.δ_ml_nondim[k] = 0
         diag_tc.asp_ratio[k] = 0
         diag_tc.frac_turb_entr[k] = 0
+        diag_tc.Π₁[k] = 0
+        diag_tc.Π₂[k] = 0
+        diag_tc.Π₃[k] = 0
+        diag_tc.Π₄[k] = 0
+        diag_tc.Π₅[k] = 0
+        diag_tc.Π₆[k] = 0
         if a_up_bulk_k > 0.0
             @inbounds for i in 1:N_up
                 aux_up_i = aux_up[i]
@@ -299,6 +311,17 @@ function compute_diagnostics!(
                 diag_tc.δ_ml_nondim[k] += aux_up_i.area[k] * aux_up_i.δ_ml_nondim[k] / a_up_bulk_k
                 diag_tc.asp_ratio[k] += aux_up_i.area[k] * aux_up_i.asp_ratio[k] / a_up_bulk_k
                 diag_tc.frac_turb_entr[k] += aux_up_i.area[k] * aux_up_i.frac_turb_entr[k] / a_up_bulk_k
+
+                for Π_i in 1:length(pi_subset)
+                    sub_script = ""
+                    for digit in string(pi_subset[Π_i])
+                        sub_script *= Char('₀' + parse(Int, digit))
+                    end
+                    property_name = "Π" * sub_script  # Concatenate "Π" with the subscript
+                    property = getproperty(diag_tc, Symbol(property_name))
+                    property[k] = aux_up_i.Π_groups[Π_i][k]
+                end
+
             end
         end
     end

src/EDMF_functions.jl

@@ -51,6 +51,9 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
     massflux_h = aux_tc_f.massflux_h
     massflux_qt = aux_tc_f.massflux_qt
     aux_tc = center_aux_turbconv(state)
+    ∂M∂z = aux_tc.∂M∂z
+    ∂lnM∂z = aux_tc.∂lnM∂z
+    massflux_c = aux_tc.massflux
 
     wvec = CC.Geometry.WVector
     ∇c = CCO.DivergenceF2C()
@@ -64,9 +67,9 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
     q_tot_gm = aux_gm.q_tot
     q_tot_en = aux_en.q_tot
     If = CCO.InterpolateC2F(; bottom = CCO.SetValue(FT(0)), top = CCO.SetValue(FT(0)))
+    Ic = CCO.InterpolateF2C()
 
     # Compute the mass flux and associated scalar fluxes
-    @. massflux = ρ_f * ᶠinterp_a(a_en) * (w_en - toscalar(w_gm))
     @. massflux_h = ρ_f * ᶠinterp_a(a_en) * (w_en - toscalar(w_gm)) * (If(θ_liq_ice_en) - If(θ_liq_ice_gm))
     @. massflux_qt = ρ_f * ᶠinterp_a(a_en) * (w_en - toscalar(w_gm)) * (If(q_tot_en) - If(q_tot_gm))
     @inbounds for i in 1:N_up
@@ -77,10 +80,20 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
         q_tot_up = aux_up_i.q_tot
         θ_liq_ice_up = aux_up_i.θ_liq_ice
         @. aux_up_f[i].massflux = ρ_f * ᶠinterp_a(a_up) * (w_up_i - toscalar(w_gm))
+        @. massflux_c += Ic(aux_up_f[i].massflux)
         @. massflux_h += ρ_f * (ᶠinterp_a(a_up) * (w_up_i - toscalar(w_gm)) * (If(θ_liq_ice_up) - If(θ_liq_ice_gm)))
         @. massflux_qt += ρ_f * (ᶠinterp_a(a_up) * (w_up_i - toscalar(w_gm)) * (If(q_tot_up) - If(q_tot_gm)))
     end
 
+    @inbounds for i in 1:N_up
+        @. massflux += aux_up_f[i].massflux
+    end
+    @. ∂M∂z = ∇c(wvec(massflux))
+
+    @inbounds for k in real_center_indices(grid)
+        aux_tc.∂lnM∂z[k] = ∂M∂z[k] / (massflux_c[k] + eps(FT))
+    end
+
     if edmf.moisture_model isa NonEquilibriumMoisture
         massflux_en = aux_tc_f.massflux_en
         massflux_ql = aux_tc_f.massflux_ql

src/closures/entr_detr.jl

@@ -37,6 +37,7 @@ function entrainment_inv_length_scale(
     tke::FT,
     zc_i::FT,
     ref_H::FT,
+    ∂lnM∂z::FT,
     ::BuoyVelEntrDimScale,
 ) where {FT}
     Δw = get_Δw(εδ_model, w_up, w_en)
@@ -53,6 +54,7 @@ function entrainment_inv_length_scale(
     tke::FT,
     zc_i::FT,
     ref_H::FT,
+    ∂lnM∂z::FT,
     ::InvScaleHeightEntrDimScale,
 ) where {FT}
     return (1 / ref_H)
@@ -67,6 +69,7 @@ function entrainment_inv_length_scale(
     tke::FT,
     zc_i::FT,
     ref_H::FT,
+    ∂lnM∂z::FT,
     ::InvZEntrDimScale,
 ) where {FT}
     return (1 / zc_i)
@@ -81,11 +84,57 @@ function entrainment_inv_length_scale(
     tke::FT,
     zc_i::FT,
     ref_H::FT,
+    ∂lnM∂z::FT,
     ::InvMeterEntrDimScale,
 ) where {FT}
     return FT(1)
 end
 
+function entrainment_inv_length_scale(
+    εδ_model,
+    b_up::FT,
+    b_en::FT,
+    w_up::FT,
+    w_en::FT,
+    tke::FT,
+    zc_i::FT,
+    ref_H::FT,
+    ∂lnM∂z::FT,
+    ::PosMassFluxGradDimScale,
+) where {FT}
+    return max(∂lnM∂z, FT(0))
+end
+
+function entrainment_inv_length_scale(
+    εδ_model,
+    b_up::FT,
+    b_en::FT,
+    w_up::FT,
+    w_en::FT,
+    tke::FT,
+    zc_i::FT,
+    ref_H::FT,
+    ∂lnM∂z::FT,
+    ::NegMassFluxGradDimScale,
+) where {FT}
+    return abs(min(∂lnM∂z, FT(0)))
+end
+
+function entrainment_inv_length_scale(
+    εδ_model,
+    b_up::FT,
+    b_en::FT,
+    w_up::FT,
+    w_en::FT,
+    tke::FT,
+    zc_i::FT,
+    ref_H::FT,
+    ∂lnM∂z::FT,
+    ::AbsMassFluxGradDimScale,
+) where {FT}
+    return abs(∂lnM∂z)
+end
+
 """A convenience wrapper for entrainment_inv_length_scale"""
 function entrainment_inv_length_scale(εδ_model, εδ_vars, dim_scale)
     return entrainment_inv_length_scale(
@@ -97,6 +146,7 @@ function entrainment_inv_length_scale(εδ_model, εδ_vars, dim_scale)
         εδ_vars.tke_en,
         εδ_vars.zc_i,
         εδ_vars.ref_H,
+        εδ_vars.∂lnM∂z,
         dim_scale,
     )
 end
@@ -254,16 +304,26 @@ function compute_phys_entr_detr!(
                     tke_en = aux_en.tke[k], # environment tke
                     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 = 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
                     zc_i = FT(grid.zc[k].z), # vertical coordinate
+                    ∂lnM∂z = aux_tc.∂lnM∂z[k], # ln(massflux) gradient
                     Δt = Δt, # Model time step
                     ε_nondim = aux_up[i].ε_nondim[k], # nondimensional fractional dynamical entrainment
                     δ_nondim = aux_up[i].δ_nondim[k], # nondimensional fractional dynamical detrainment
+                    entr_Π_subset = entrainment_Π_subset(edmf), # indices of Pi groups to include
                 )
 
+                # store pi groups for output
+                Π = non_dimensional_groups(εδ_closure, εδ_model_vars)
+                @assert length(Π) == n_Π_groups(edmf)
+                for Π_i in 1:length(entrainment_Π_subset(edmf))
+                    aux_up[i].Π_groups[Π_i][k] = Π[Π_i]
+                end
+
                 # update fractional and turbulent entr/detr
                 if εδ_closure isa PrognosticNoisyRelaxationProcess
                     # fractional dynamical entrainment from prognostic state
@@ -392,10 +452,15 @@ function compute_ml_entr_detr!(
                     RH_en = aux_en.RH[k], # environment relative humidity
                     max_area = max_area, # maximum updraft area
                     zc_i = FT(grid.zc[k].z), # vertical coordinate
-                    wstar = surf.wstar, # convective velocity
+                    ∂lnM∂z = aux_tc.∂lnM∂z[k], # ln(massflux) gradient
                     entr_Π_subset = entrainment_Π_subset(edmf), # indices of Pi groups to include
                 )
-
+                # store pi groups for output
+                Π = non_dimensional_groups(εδ_closure, εδ_model_vars)
+                @assert length(Π) == n_Π_groups(edmf)
+                for Π_i in 1:length(entrainment_Π_subset(edmf))
+                    aux_up[i].Π_groups[Π_i][k] = Π[Π_i]
+                end
                 # update fractional and turbulent entr/detr
                 # fractional, turbulent & nondimensional entrainment
                 ε_ml_nondim, δ_ml_nondim = non_dimensional_function(εδ_closure, εδ_model_vars)
@@ -489,9 +554,9 @@ function compute_ml_entr_detr!(
                     RH_en = aux_en.RH[k], # environment relative humidity
                     max_area = max_area, # maximum updraft area
                     zc_i = FT(grid.zc[k].z), # vertical coordinate
+                    ∂lnM∂z = aux_tc.∂lnM∂z[k], # ln(massflux) gradient
                     ε_ml_nondim = aux_up[i].ε_ml_nondim[k], # nondimensional fractional dynamical entrainment
                     δ_ml_nondim = aux_up[i].δ_ml_nondim[k], # nondimensional fractional dynamical detrainment
-                    wstar = surf.wstar, # convective velocity
                     entr_Π_subset = entrainment_Π_subset(edmf), # indices of Pi groups to include
                 )
                 Π = non_dimensional_groups(εδ_model, εδ_model_vars)
@@ -522,6 +587,7 @@ function compute_ml_entr_detr!(
                 aux_en.tke[k],
                 FT(grid.zc[k].z),
                 p_c[k] / (ρ_c[k] * g),
+                aux_tc.∂lnM∂z[k],
                 edmf.entr_dim_scale,
             )
             δ_dim_scale = entrainment_inv_length_scale(
@@ -533,6 +599,7 @@ function compute_ml_entr_detr!(
                 aux_en.tke[k],
                 FT(grid.zc[k].z),
                 p_c[k] / (ρ_c[k] * g),
+                aux_tc.∂lnM∂z[k],
                 edmf.detr_dim_scale,
             )
             aux_up[i].entr_ml[k] = ε_dim_scale * aux_up[i].ε_ml_nondim[k]

src/closures/nondimensional_exchange_functions.jl

@@ -12,14 +12,13 @@ function non_dimensional_groups(εδ_model, εδ_model_vars)
     Δw = get_Δw(εδ_model, εδ_model_vars.w_up, εδ_model_vars.w_en)
     Δb = εδ_model_vars.b_up - εδ_model_vars.b_en
     Π_norm = εδ_params(εδ_model).Π_norm
-    Π₁ = (εδ_model_vars.zc_i * Δb) / (Δw^2 + εδ_model_vars.wstar^2) / Π_norm[1]
-    Π₂ =
-        (εδ_model_vars.tke_gm - εδ_model_vars.a_en * εδ_model_vars.tke_en) / (εδ_model_vars.tke_gm + eps(FT)) /
-        Π_norm[2]
+    Π₁ = (εδ_model_vars.zc_i * Δb) / (Δw^2 + eps(FT)) / Π_norm[1]
+    Π₂ = εδ_model_vars.tke_en / (Δw^2 + eps(FT)) / Π_norm[2]
     Π₃ = √(εδ_model_vars.a_up) / Π_norm[3]
     Π₄ = (εδ_model_vars.RH_up - εδ_model_vars.RH_en) / Π_norm[4]
     Π₅ = εδ_model_vars.zc_i / εδ_model_vars.H_up / Π_norm[5]
     Π₆ = εδ_model_vars.zc_i / εδ_model_vars.ref_H / Π_norm[6]
+
     Π_groups = (Π₁, Π₂, Π₃, Π₄, Π₅, Π₆)
 
     return map(i -> Π_groups[i], εδ_model_vars.entr_Π_subset)

src/diagnostics.jl

@@ -105,6 +105,8 @@ function io_dictionary_aux()
         "updraft_qt_precip" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_bulk(state).qt_tendency_precip_formation),
         "updraft_thetal_precip" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_bulk(state).θ_liq_ice_tendency_precip_formation),
 
+        "massflux_grad" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_turbconv(state).∂M∂z),
+        "ln_massflux_grad" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_turbconv(state).∂lnM∂z),
         "massflux_tendency_h" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_turbconv(state).massflux_tendency_h),
         "massflux_tendency_qt" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_turbconv(state).massflux_tendency_qt),
         "diffusive_tendency_h" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_turbconv(state).diffusive_tendency_h),

src/types.jl

@@ -132,6 +132,9 @@ struct BuoyVelEntrDimScale <: EntrDimScale end
 struct InvScaleHeightEntrDimScale <: EntrDimScale end
 struct InvZEntrDimScale <: EntrDimScale end
 struct InvMeterEntrDimScale <: EntrDimScale end
+struct PosMassFluxGradDimScale <: EntrDimScale end
+struct NegMassFluxGradDimScale <: EntrDimScale end
+struct AbsMassFluxGradDimScale <: EntrDimScale end
 
 """
     GradBuoy
@@ -683,7 +686,15 @@ function EDMFModel(::Type{FT}, namelist, precip_model, rain_formation_model) whe
         "EDMF_PrognosticTKE",
         "entr_dim_scale";
         default = "buoy_vel",
-        valid_options = ["buoy_vel", "inv_scale_height", "inv_z", "none"],
+        valid_options = [
+            "buoy_vel",
+            "inv_scale_height",
+            "inv_z",
+            "pos_massflux",
+            "neg_massflux",
+            "abs_massflux",
+            "none",
+        ],
     )
 
     pressure_model_params = PressureModelParams{FT}(;
@@ -712,6 +723,12 @@ function EDMFModel(::Type{FT}, namelist, precip_model, rain_formation_model) whe
         InvScaleHeightEntrDimScale()
     elseif entr_dim_scale == "inv_z"
         InvZEntrDimScale()
+    elseif entr_dim_scale == "pos_massflux"
+        PosMassFluxGradDimScale()
+    elseif entr_dim_scale == "neg_massflux"
+        NegMassFluxGradDimScale()
+    elseif entr_dim_scale == "abs_massflux"
+        AbsMassFluxGradDimScale()
     elseif entr_dim_scale == "none"
         InvMeterEntrDimScale()
     else
@@ -724,7 +741,15 @@ function EDMFModel(::Type{FT}, namelist, precip_model, rain_formation_model) whe
         "EDMF_PrognosticTKE",
         "detr_dim_scale";
         default = "buoy_vel",
-        valid_options = ["buoy_vel", "inv_scale_height", "inv_z", "none"],
+        valid_options = [
+            "buoy_vel",
+            "inv_scale_height",
+            "inv_z",
+            "pos_massflux",
+            "neg_massflux",
+            "abs_massflux",
+            "none",
+        ],
     )
 
     detr_dim_scale = if detr_dim_scale == "buoy_vel"
@@ -733,6 +758,12 @@ function EDMFModel(::Type{FT}, namelist, precip_model, rain_formation_model) whe
         InvScaleHeightEntrDimScale()
     elseif detr_dim_scale == "inv_z"
         InvZEntrDimScale()
+    elseif detr_dim_scale == "pos_massflux"
+        PosMassFluxGradDimScale()
+    elseif detr_dim_scale == "neg_massflux"
+        NegMassFluxGradDimScale()
+    elseif detr_dim_scale == "abs_massflux"
+        AbsMassFluxGradDimScale()
     elseif detr_dim_scale == "none"
         InvMeterEntrDimScale()
     else

src/variables.jl

@@ -144,6 +144,7 @@ cent_aux_vars_edmf(::Type{FT}, local_geometry, edmf) where {FT} = (;
         KH = FT(0),
         KQ = FT(0),
         mixing_length = FT(0),
+        massflux = FT(0),
         massflux_tendency_h = FT(0),
         massflux_tendency_qt = FT(0),
         diffusive_tendency_h = FT(0),
@@ -158,6 +159,8 @@ cent_aux_vars_edmf(::Type{FT}, local_geometry, edmf) where {FT} = (;
         w_up_c = FT(0),
         w_en_c = FT(0),
         Shear² = FT(0),
+        ∂M∂z = FT(0),
+        ∂lnM∂z = FT(0),
         ∂θv∂z = FT(0),
         ∂qt∂z = FT(0),
         ∂θl∂z = FT(0),
@@ -227,6 +230,13 @@ cent_diagnostic_vars_edmf(FT, local_geometry, edmf) = (;
         δ_ml_nondim = FT(0),
         massflux = FT(0),
         frac_turb_entr = FT(0),
+        Π_groups = ntuple(i -> FT(0), n_Π_groups(edmf)),
+        Π₁ = FT(0),
+        Π₂ = FT(0),
+        Π₃ = FT(0),
+        Π₄ = FT(0),
+        Π₅ = FT(0),
+        Π₆ = FT(0),
     )
 )