Merge pull request #3245 from CliMA/cc/update_entr_pi_groups

Add entrainment parameter vectors and turbulent entrainment

regression_tests/ref_counter.jl

@@ -1,6 +1,9 @@
-178
+179
 
 #=
+179:
+- Update Pi entr groups/ add parameter vectors
+
 178:
 - Added aerosol to one of the regression tests
 

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -242,6 +242,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
         ᶜρʲs,
         ᶜentrʲs,
         ᶜdetrʲs,
+        ᶜturb_entrʲs,
         ᶠnh_pressure³ʲs,
     ) = p.precomputed
     (; ᶠu³⁰, ᶜK⁰, ᶜtke⁰) = p.precomputed
@@ -324,6 +325,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
             ᶜq_totʲ = ᶜq_totʲs.:($j)
             ᶜentrʲ = ᶜentrʲs.:($j)
             ᶜdetrʲ = ᶜdetrʲs.:($j)
+            ᶜturb_entrʲ = ᶜturb_entrʲs.:($j)
             ᶠnh_pressure³ʲ = ᶠnh_pressure³ʲs.:($j)
 
             if precip_model isa Union{Microphysics0Moment, Microphysics1Moment}
@@ -356,6 +358,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
             tsʲ_prev_level = Fields.field_values(Fields.level(ᶜtsʲ, i - 1))
             entrʲ_prev_level = Fields.field_values(Fields.level(ᶜentrʲ, i - 1))
             detrʲ_prev_level = Fields.field_values(Fields.level(ᶜdetrʲ, i - 1))
+            turb_entrʲ_prev_level =
+                Fields.field_values(Fields.level(ᶜturb_entrʲ, i - 1))
             nh_pressure³ʲ_prev_halflevel =
                 Fields.field_values(Fields.level(ᶠnh_pressure³ʲ, i - 1 - half))
             scale_height =
@@ -421,6 +425,21 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                 dt,
             )
 
+            @. turb_entrʲ_prev_level = turbulent_entrainment(
+                params,
+                draft_area(ρaʲ_prev_level, ρʲ_prev_level),
+            )
+
+            @. turb_entrʲ_prev_level = limit_turb_entrainment(
+                entrʲ_prev_level,
+                turb_entrʲ_prev_level,
+                get_physical_w(
+                    u³ʲ_prev_halflevel,
+                    local_geometry_prev_halflevel,
+                ),
+                dz_prev_level,
+            )
+
             # TODO: use updraft top instead of scale height
             if p.atmos.edmfx_model.nh_pressure isa Val{true}
                 @. nh_pressure³ʲ_prev_halflevel = ᶠupdraft_nh_pressure(
@@ -508,8 +527,10 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                     local_geometry_prev_level.J *
                     2 *
                     (
-                        entrʲ_prev_level * u³⁰_data_prev_halflevel -
-                        entrʲ_prev_level * u³ʲ_data_prev_halflevel
+                        (entrʲ_prev_level + turb_entrʲ_prev_level) *
+                        u³⁰_data_prev_halflevel -
+                        (entrʲ_prev_level + turb_entrʲ_prev_level) *
+                        u³ʲ_data_prev_halflevel
                     )
                 )
 
@@ -639,8 +660,10 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                     local_geometry_prev_level.J *
                     ρaʲ_prev_level *
                     (
-                        entrʲ_prev_level * (h_tot_prev_level - K_prev_level) -
-                        detrʲ_prev_level * mseʲ_prev_level
+                        (entrʲ_prev_level + turb_entrʲ_prev_level) *
+                        (h_tot_prev_level - K_prev_level) -
+                        (detrʲ_prev_level + turb_entrʲ_prev_level) *
+                        mseʲ_prev_level
                     )
                 )
             if precip_model isa Microphysics0Moment
@@ -686,8 +709,10 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                     local_geometry_prev_level.J *
                     ρaʲ_prev_level *
                     (
-                        entrʲ_prev_level * q_tot_prev_level -
-                        detrʲ_prev_level * q_totʲ_prev_level
+                        (entrʲ_prev_level + turb_entrʲ_prev_level) *
+                        q_tot_prev_level -
+                        (detrʲ_prev_level + turb_entrʲ_prev_level) *
+                        q_totʲ_prev_level
                     )
                 )
             if precip_model isa Union{Microphysics0Moment, Microphysics1Moment}
@@ -791,7 +816,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_top_bc!(
     t,
 )
     n = n_mass_flux_subdomains(p.atmos.turbconv_model)
-    (; ᶜentrʲs, ᶜdetrʲs) = p.precomputed
+    (; ᶜentrʲs, ᶜdetrʲs, ᶜturb_entrʲs) = p.precomputed
     (; ᶠu³⁰, ᶠu³ʲs, ᶜuʲs, ᶠnh_pressure³ʲs) = p.precomputed
     (; precip_model) = p.atmos
 
@@ -806,6 +831,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_top_bc!(
         ᶠnh_pressure³ʲ = ᶠnh_pressure³ʲs.:($j)
         ᶜentrʲ = ᶜentrʲs.:($j)
         ᶜdetrʲ = ᶜdetrʲs.:($j)
+        ᶜturb_entrʲ = ᶜturb_entrʲs.:($j)
 
         u³ʲ_halflevel = Fields.field_values(Fields.level(ᶠu³ʲ, i_top + half))
         @. u³ʲ_halflevel = CT3(0)
@@ -818,8 +844,10 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_top_bc!(
 
         entrʲ_level = Fields.field_values(Fields.level(ᶜentrʲ, i_top))
         detrʲ_level = Fields.field_values(Fields.level(ᶜdetrʲ, i_top))
+        turb_entrʲ_level = Fields.field_values(Fields.level(ᶜturb_entrʲ, i_top))
         fill!(entrʲ_level, RecursiveApply.rzero(eltype(entrʲ_level)))
         fill!(detrʲ_level, RecursiveApply.rzero(eltype(detrʲ_level)))
+        fill!(turb_entrʲ_level, RecursiveApply.rzero(eltype(turb_entrʲ_level)))
         @. ᶜuʲ = C123(Y.c.uₕ) + ᶜinterp(C123(ᶠu³ʲ))
 
         if precip_model isa Union{Microphysics0Moment, Microphysics1Moment}

src/cache/precomputed_quantities.jl

@@ -101,6 +101,7 @@ function precomputed_quantities(Y, atmos)
             ᶜρʲs = similar(Y.c, NTuple{n, FT}),
             ᶜentrʲs = similar(Y.c, NTuple{n, FT}),
             ᶜdetrʲs = similar(Y.c, NTuple{n, FT}),
+            ᶜturb_entrʲs = similar(Y.c, NTuple{n, FT}),
             ᶠnh_pressure₃ʲs = similar(Y.f, NTuple{n, C3{FT}}),
             ᶜgradᵥ_θ_virt⁰ = Fields.Field(C3{FT}, cspace),
             ᶜgradᵥ_q_tot⁰ = Fields.Field(C3{FT}, cspace),
@@ -126,6 +127,7 @@ function precomputed_quantities(Y, atmos)
             ᶜq_totʲs = similar(Y.c, NTuple{n, FT}),
             ᶜentrʲs = similar(Y.c, NTuple{n, FT}),
             ᶜdetrʲs = similar(Y.c, NTuple{n, FT}),
+            ᶜturb_entrʲs = similar(Y.c, NTuple{n, FT}),
             ᶠnh_pressure³ʲs = similar(Y.f, NTuple{n, CT3{FT}}),
             ᶠu³⁰ = similar(Y.f, CT3{FT}),
             ᶜu⁰ = similar(Y.c, C123{FT}),

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -204,7 +204,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_closures!(
         ᶜK_h,
         ρatke_flux,
     ) = p.precomputed
-    (; ᶜuʲs, ᶜtsʲs, ᶠu³ʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs) = p.precomputed
+    (; ᶜuʲs, ᶜtsʲs, ᶠu³ʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs, ᶜturb_entrʲs) = p.precomputed
     (; ustar, obukhov_length) = p.precomputed.sfc_conditions
 
     ᶜz = Fields.coordinate_field(Y.c).z
@@ -232,11 +232,21 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_closures!(
             max(ᶜtke⁰, 0),
             p.atmos.edmfx_model.entr_model,
         )
+
         @. ᶜentrʲs.:($$j) = limit_entrainment(
             ᶜentrʲs.:($$j),
             draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
             dt,
         )
+
+        @. ᶜturb_entrʲs.:($$j) = turbulent_entrainment(
+            params,
+            draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
+        )
+
+        @. ᶜturb_entrʲs.:($$j) =
+            limit_turb_entrainment(ᶜentrʲs.:($$j), ᶜturb_entrʲs.:($$j), dt)
+
         @. ᶜvert_div = ᶜdivᵥ(ᶠinterp(ᶜρʲs.:($$j)) * ᶠu³ʲs.:($$j)) / ᶜρʲs.:($$j)
         @. ᶜmassflux_vert_div =
             ᶜdivᵥ(ᶠinterp(Y.c.sgsʲs.:($$j).ρa) * ᶠu³ʲs.:($$j))
@@ -260,6 +270,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_closures!(
             ᶜtke⁰,
             p.atmos.edmfx_model.detr_model,
         )
+
         @. ᶜdetrʲs.:($$j) = limit_detrainment(
             ᶜdetrʲs.:($$j),
             draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),

src/diagnostics/edmfx_diagnostics.jl

@@ -442,6 +442,33 @@ add_diagnostic_variable!(
     compute! = compute_entr!,
 )
 
+compute_turbentr!(out, state, cache, time) =
+    compute_turbentr!(out, state, cache, time, cache.atmos.turbconv_model)
+compute_turbentr!(_, _, _, _, turbconv_model::T) where {T} =
+    error_diagnostic_variable("turbentr", turbconv_model)
+
+function compute_turbentr!(
+    out,
+    state,
+    cache,
+    time,
+    turbconv_model::Union{PrognosticEDMFX, DiagnosticEDMFX},
+)
+    if isnothing(out)
+        return copy(cache.precomputed.ᶜturb_entrʲs.:1)
+    else
+        out .= cache.precomputed.ᶜturb_entrʲs.:1
+    end
+end
+
+add_diagnostic_variable!(
+    short_name = "turbentr",
+    long_name = "Turbulent entrainment rate",
+    units = "s^-1",
+    comments = "Turbulent entrainment rate of the first updraft",
+    compute! = compute_turbentr!,
+)
+
 ###
 # Detrainment (3d)
 ###

src/parameters/Parameters.jl

@@ -19,7 +19,7 @@ Base.broadcastable(param_set::ACAP) = tuple(param_set)
 Base.broadcastable(param_set::ATCP) = tuple(param_set)
 Base.broadcastable(param_set::ASTP) = tuple(param_set)
 
-Base.@kwdef struct TurbulenceConvectionParameters{FT} <: ATCP
+Base.@kwdef struct TurbulenceConvectionParameters{FT, VFT1, VFT2} <: ATCP
     surface_area::FT
     max_area::FT
     min_area::FT
@@ -42,6 +42,8 @@ Base.@kwdef struct TurbulenceConvectionParameters{FT} <: ATCP
     detr_coeff::FT
     detr_buoy_coeff::FT
     detr_vertdiv_coeff::FT
+    entr_param_vec::VFT1
+    turb_entr_param_vec::VFT2
     detr_massflux_vertdiv_coeff::FT
     min_area_limiter_scale::FT
     min_area_limiter_power::FT

src/parameters/create_parameters.jl

@@ -5,6 +5,12 @@ import SurfaceFluxes.UniversalFunctions as UF
 import Insolation.Parameters.InsolationParameters
 import Thermodynamics.Parameters.ThermodynamicsParameters
 import CloudMicrophysics as CM
+import StaticArrays as SA
+
+
+to_svec(x::AbstractArray) = SA.SVector{length(x)}(x)
+to_svec(x) = x
+to_svec(x::NamedTuple) = map(x -> to_svec(x), x)
 
 function TurbulenceConvectionParameters(toml_dict::CP.AbstractTOMLDict)
     name_map = (;
@@ -20,6 +26,8 @@ function TurbulenceConvectionParameters(toml_dict::CP.AbstractTOMLDict)
         :mixing_length_Ri_crit => :Ri_crit,
         :detr_coeff => :detr_coeff,
         :EDMF_surface_area => :surface_area,
+        :entr_param_vec => :entr_param_vec,
+        :turb_entr_param_vec => :turb_entr_param_vec,
         :minimum_updraft_top => :min_updraft_top,
         :mixing_length_eddy_viscosity_coefficient => :tke_ed_coeff,
         :mixing_length_smin_ub => :smin_ub,
@@ -39,7 +47,10 @@ function TurbulenceConvectionParameters(toml_dict::CP.AbstractTOMLDict)
     )
     parameters = CP.get_parameter_values(toml_dict, name_map, "ClimaAtmos")
     FT = CP.float_type(toml_dict)
-    CAP.TurbulenceConvectionParameters{FT}(; parameters...)
+    parameters = to_svec(parameters)
+    VFT1 = typeof(parameters.entr_param_vec)
+    VFT2 = typeof(parameters.turb_entr_param_vec)
+    CAP.TurbulenceConvectionParameters{FT, VFT1, VFT2}(; parameters...)
 end
 
 function SurfaceTemperatureParameters(toml_dict::CP.AbstractTOMLDict)