wip [skip ci]

config/model_configs/advective_edmfx_bomex_box.yml

@@ -0,0 +1,30 @@
+job_id: "advectiveedmfx_bomex_box"
+initial_condition: "Bomex"
+subsidence: "Bomex"
+edmf_coriolis: "Bomex"
+ls_adv: "Bomex"
+surface_setup: "Bomex"
+turbconv: "advective_edmfx"
+edmfx_upwinding: first_order
+edmfx_entr_model: "ConstantCoefficient"
+edmfx_detr_model: "ConstantCoefficient"
+edmfx_sgs_mass_flux: true
+edmfx_sgs_diffusive_flux: true
+edmfx_nh_pressure: false
+prognostic_tke: false
+moist: "equil"
+config: "box"
+hyperdiff: "true"
+kappa_4: 1.0e12
+x_max: 1e5
+y_max: 1e5
+z_max: 3e3
+x_elem: 2
+y_elem: 2
+z_elem: 60
+z_stretch: false
+perturb_initstate: false
+dt: "1secs"
+t_end: "10800secs"
+dt_save_to_disk: "10secs"
+toml: [toml/edmfx_box.toml]

src/ClimaAtmos.jl

@@ -26,6 +26,7 @@ include(joinpath("TurbulenceConvection_deprecated", "TurbulenceConvection.jl"))
 import .TurbulenceConvection as TC
 
 include(joinpath("cache", "edmf_precomputed_quantities.jl"))
+include(joinpath("cache", "advective_edmf_precomputed_quantities.jl"))
 include(joinpath("cache", "diagnostic_edmf_precomputed_quantities.jl"))
 include(joinpath("cache", "precomputed_quantities.jl"))
 

src/cache/advective_edmf_precomputed_quantities.jl

@@ -23,44 +23,65 @@ function set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
     thermo_args = (thermo_params, energy_form, moisture_model)
 
     (; ᶜspecific, ᶜp, ᶜΦ, ᶜh_tot, ᶜρ_ref) = p
-    (; ᶜspecific⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜh_tot⁰) = p
+    (; ᶜtke⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜh_tot⁰, ᶜq_tot⁰) = p
     (; ᶜmixing_length, ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜK_u, ᶜK_h) = p
     (;
-        ᶜspecificʲs,
         ᶜuʲs,
         ᶠu³ʲs,
         ᶜKʲs,
         ᶜtsʲs,
         ᶜρʲs,
-        ᶜh_totʲs,
         ᶜentrʲs,
         ᶜdetrʲs,
     ) = p
     (; ustar, obukhov_length, buoyancy_flux) = p.sfc_conditions
 
-    @. ᶜspecific⁰ = specific_full_sgs⁰(Y.c, turbconv_model)
     @. ᶜρa⁰ = ρa⁰(Y.c)
+    @. ᶜtke⁰ = divide_by_ρa(Y.c.sgs⁰.ρatke, ᶜρa⁰, 0, Y.c.ρ, turbconv_model)
+    # TODO: use all updrafts
+    @. ᶜh_tot⁰ = divide_by_ρa(
+        Y.c.ρ * ᶜh_tot - (Y.c.sgsʲs.:1).ρa * (Y.c.sgsʲs.:1).h_tot,
+        ᶜρa⁰,
+        Y.c.ρ * ᶜh_tot,
+        Y.c.ρ,
+        turbconv_model,
+    )
+    @. ᶜq_tot⁰ = divide_by_ρa(
+        Y.c.ρq_tot - (Y.c.sgsʲs.:1).ρa * (Y.c.sgsʲs.:1).q_tot,
+        ᶜρa⁰,
+        Y.c.ρq_tot,
+        Y.c.ρ,
+        turbconv_model,
+    )
     set_sgs_ᶠu₃!(u₃⁰, ᶠu₃⁰, Y, turbconv_model)
     set_velocity_quantities!(ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶠu₃⁰, Y.c.uₕ, ᶠuₕ³)
-    @. ᶜK⁰ += ᶜspecific⁰.tke
-    @. ᶜts⁰ = ts_sgs(thermo_args..., ᶜspecific⁰, ᶜK⁰, ᶜΦ, ᶜp)
+    @. ᶜK⁰ += ᶜtke⁰
+    @. ᶜts⁰ = TD.PhaseEquil_phq(
+        thermo_params,
+        ᶜp,
+        ᶜh_tot⁰ - ᶜK⁰ - ᶜΦ,
+        ᶜq_tot⁰,
+    )
     @. ᶜρ⁰ = TD.air_density(thermo_params, ᶜts⁰)
-    @. ᶜh_tot⁰ =
-        TD.total_specific_enthalpy(thermo_params, ᶜts⁰, ᶜspecific⁰.e_tot)
-    @. ᶜspecificʲs = specific_sgsʲs(Y.c, turbconv_model)
+
     for j in 1:n
         ᶜuʲ = ᶜuʲs.:($j)
         ᶠu³ʲ = ᶠu³ʲs.:($j)
         ᶜKʲ = ᶜKʲs.:($j)
         ᶠu₃ʲ = Y.f.sgsʲs.:($j).u₃
-        ᶜspecificʲ = ᶜspecificʲs.:($j)
         ᶜtsʲ = ᶜtsʲs.:($j)
         ᶜρʲ = ᶜρʲs.:($j)
-        ᶜh_totʲ = Y.c.sgsʲs.h_tot.:($j)
-        ᶜq_totʲ = Y.c.sgsʲs.q_tot.:($j)
+        ᶜh_totʲ = Y.c.sgsʲs.:($j).h_tot
+        ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot
 
         set_velocity_quantities!(ᶜuʲ, ᶠu³ʲ, ᶜKʲ, ᶠu₃ʲ, Y.c.uₕ, ᶠuₕ³)
-        @. ᶜtsʲ = ts_sgs(thermo_args..., ᶜspecificʲ, ᶜKʲ, ᶜΦ, ᶜp)
+        # FIXME: get ts
+        @. ᶜtsʲ = TD.PhaseEquil_phq(
+            thermo_params,
+            ᶜp,
+            ᶜh_totʲ - ᶜKʲ - ᶜΦ,
+            ᶜq_totʲ,
+        )
         @. ᶜρʲ = TD.air_density(thermo_params, ᶜtsʲ)
 
         # When ᶜe_intʲ = ᶜe_int and ᶜq_totʲ = ᶜq_tot, we still observe that
@@ -96,8 +117,6 @@ function set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
 
         # Based on boundary conditions for updrafts we overwrite
         # the first interior point for EDMFX ᶜh_totʲ...
-        @. ᶜh_totʲ = 
-            TD.total_specific_enthalpy(thermo_params, ᶜtsʲ, ᶜspecificʲ.e_tot)
         ᶜh_tot_int_val = Fields.field_values(Fields.level(ᶜh_tot, 1))
         ᶜh_totʲ_int_val = Fields.field_values(Fields.level(ᶜh_totʲ, 1))
         @. ᶜh_totʲ_int_val = sgs_scalar_first_interior_bc(
@@ -237,11 +256,11 @@ function set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
                 1 / 2 *
                 (
                     get_physical_w(ᶜuʲs.:($$j), ᶜlg) - get_physical_w(ᶜu⁰, ᶜlg)
-                )^2 - ᶜspecific⁰.tke
+                )^2 - ᶜtke⁰
             )
     end
 
-    sfc_tke = Fields.level(ᶜspecific⁰.tke, 1)
+    sfc_tke = Fields.level(ᶜtke⁰, 1)
     @. ᶜmixing_length = mixing_length(
         p.params,
         ustar,
@@ -250,7 +269,7 @@ function set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
         ᶜdz,
         sfc_tke,
         ᶜlinear_buoygrad,
-        ᶜspecific⁰.tke,
+        ᶜtke⁰,
         obukhov_length,
         ᶜstrain_rate_norm,
         ᶜprandtl_nvec,
@@ -259,7 +278,7 @@ function set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
 
     turbconv_params = CAP.turbconv_params(params)
     c_m = TCP.tke_ed_coeff(turbconv_params)
-    @. ᶜK_u = c_m * ᶜmixing_length * sqrt(max(ᶜspecific⁰.tke, 0))
+    @. ᶜK_u = c_m * ᶜmixing_length * sqrt(max(ᶜtke⁰, 0))
     # TODO: add Prantdl number
     @. ᶜK_h = ᶜK_u / ᶜprandtl_nvec
 

src/cache/precomputed_quantities.jl

@@ -95,34 +95,28 @@ function precomputed_quantities(Y, atmos)
     advective_sgs_quantities =
         atmos.turbconv_model isa AdvectiveEDMFX ?
         (;
-            ᶜspecific⁰ = specific_full_sgs⁰.(Y.c, atmos.turbconv_model),
+            ᶜtke⁰ = similar(Y.c, FT),
             ᶜρa⁰ = similar(Y.c, FT),
             ᶠu₃⁰ = similar(Y.f, C3{FT}),
             ᶜu⁰ = similar(Y.c, C123{FT}),
             ᶠu³⁰ = similar(Y.f, CT3{FT}),
             ᶜK⁰ = similar(Y.c, FT),
+            ᶜh_tot⁰ = similar(Y.c, FT),
+            ᶜq_tot⁰ = similar(Y.c, FT),
             ᶜts⁰ = similar(Y.c, TST),
             ᶜρ⁰ = similar(Y.c, FT),
             ᶜlinear_buoygrad = similar(Y.c, FT),
             ᶜstrain_rate_norm = similar(Y.c, FT),
             ᶜmixing_length = similar(Y.c, FT),
             ᶜK_u = similar(Y.c, FT),
             ᶜK_h = similar(Y.c, FT),
-            ᶜspecificʲs = specific_sgsʲs.(Y.c, atmos.turbconv_model),
             ᶜuʲs = similar(Y.c, NTuple{n, C123{FT}}),
             ᶠu³ʲs = similar(Y.f, NTuple{n, CT3{FT}}),
             ᶜKʲs = similar(Y.c, NTuple{n, FT}),
             ᶜtsʲs = similar(Y.c, NTuple{n, TST}),
             ᶜρʲs = similar(Y.c, NTuple{n, FT}),
             ᶜentrʲs = similar(Y.c, NTuple{n, FT}),
             ᶜdetrʲs = similar(Y.c, NTuple{n, FT}),
-            (
-                atmos.energy_form isa TotalEnergy ?
-                (;
-                    ᶜh_totʲs = similar(Y.c, NTuple{n, FT}),
-                    ᶜh_tot⁰ = similar(Y.c, FT),
-                ) : (;)
-            )...,
         ) : (;)
     diagnostic_sgs_quantities =
         atmos.turbconv_model isa DiagnosticEDMFX ?
@@ -389,6 +383,27 @@ function output_sgs_quantities(Y, p, t)
     return (; ᶜspecific⁺, ᶠu₃⁺, ᶜu⁺, ᶠu³⁺, ᶜK⁺, ᶜts⁺, ᶜa⁺, ᶜa⁰)
 end
 
+"""
+    output_advective_sgs_quantities(Y, p, t)
+
+Sets `ᶜu⁺`, `ᶠu³⁺`, `ᶜts⁺` and `ᶜa⁺` to be the same as the
+values of the first updraft.
+"""
+function output_advective_sgs_quantities(Y, p, t)
+    (; turbconv_model) = p.atmos
+    thermo_params = CAP.thermodynamics_params(p.params)
+    (; ᶜp, ᶜρa⁰, ᶜρ⁰, ᶜΦ, ᶜtsʲs) = p
+    ᶠuₕ³ = p.ᶠtemp_CT3
+    set_ᶠuₕ³!(ᶠuₕ³, Y)
+    (ᶠu₃⁺, ᶜu⁺, ᶠu³⁺, ᶜK⁺) = similar.((p.ᶠu₃⁰, p.ᶜu⁰, p.ᶠu³⁰, p.ᶜK⁰))
+    set_sgs_ᶠu₃!(u₃⁺, ᶠu₃⁺, Y, turbconv_model)
+    set_velocity_quantities!(ᶜu⁺, ᶠu³⁺, ᶜK⁺, ᶠu₃⁺, Y.c.uₕ, ᶠuₕ³)
+    ᶜts⁺ = ᶜtsʲs.:1
+    ᶜa⁺ = @. draft_area(ρa⁺(Y.c), TD.air_density(thermo_params, ᶜts⁺))
+    ᶜa⁰ = @. draft_area(ᶜρa⁰, ᶜρ⁰)
+    return (; ᶠu₃⁺, ᶜu⁺, ᶠu³⁺, ᶜK⁺, ᶜts⁺, ᶜa⁺, ᶜa⁰)
+end
+
 """
     output_diagnostic_sgs_quantities(Y, p, t)
 

src/prognostic_equations/advection.jl

@@ -101,13 +101,14 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     ᶜJ = Fields.local_geometry_field(Y.c).J
     (; ᶜu, ᶠu³, ᶜK, ᶜf) = p
     (; edmfx_upwinding) = n > 0 || advect_tke ? p : all_nothing
-    (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜρʲs, ᶜspecificʲs) = n > 0 ? p : all_nothing
+    (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜρʲs) = n > 0 ? p : all_nothing
+    (; ᶜspecificʲs) = turbconv_model isa EDMFX ? p : all_nothing
     (; ᶜp, ᶜp_ref, ᶜρ_ref, ᶠgradᵥ_ᶜΦ) = n > 0 ? p : all_nothing
-    (; ᶜh_totʲs) = n > 0 && is_total_energy ? p : all_nothing
+    (; ᶜh_totʲs) = turbconv_model isa EDMFX && is_total_energy ? p : all_nothing
     (; ᶠu³⁰) = advect_tke ? p : all_nothing
     ᶜρa⁰ = advect_tke ? (n > 0 ? p.ᶜρa⁰ : Y.c.ρ) : nothing
     ᶜρ⁰ = advect_tke ? (n > 0 ? p.ᶜρ⁰ : Y.c.ρ) : nothing
-    ᶜtke⁰ = advect_tke ? (n > 0 ? p.ᶜspecific⁰.tke : p.ᶜtke⁰) : nothing
+    ᶜtke⁰ = advect_tke ? (turbconv_model isa EDMFX ? p.ᶜspecific⁰.tke : p.ᶜtke⁰) : nothing
     ᶜa_scalar = p.ᶜtemp_scalar
     ᶜω³ = p.ᶜtemp_CT3
     ᶠω¹² = p.ᶠtemp_CT12

src/prognostic_equations/edmfx_entr_detr.jl

@@ -318,8 +318,8 @@ end
 function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::AdvectiveEDMFX)
 
     n = n_mass_flux_subdomains(turbconv_model)
-    (; ᶜspecificʲs, ᶜh_totʲs, ᶜentrʲs, ᶜdetrʲs) = p
-    (; ᶜspecific⁰, ᶜh_tot⁰, ᶠu₃⁰) = p
+    (; ᶜentrʲs, ᶜdetrʲs) = p
+    (; ᶜq_tot⁰, ᶜh_tot⁰, ᶠu₃⁰) = p
 
     for j in 1:n
 
@@ -329,11 +329,11 @@ function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::Advect
 
         @. Yₜ.c.sgsʲs.:($$j).h_tot[colidx] +=
             ᶜentrʲs.:($$j)[colidx] *
-            (ᶜh_tot⁰[colidx] - ᶜh_totʲs.:($$j)[colidx])
+            (ᶜh_tot⁰[colidx] - Y.c.sgsʲs.:($$j).h_tot[colidx])
 
         @. Yₜ.c.sgsʲs.:($$j).q_tot[colidx] +=
             ᶜentrʲs.:($$j)[colidx] *
-            (ᶜq_tot⁰[colidx] - ᶜq_totʲs.:($$j)[colidx])
+            (ᶜq_tot⁰[colidx] - Y.c.sgsʲs.:($$j).q_tot[colidx])
 
         @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] +=
             ᶠinterp(ᶜentrʲs.:($$j)[colidx]) *

src/solver/model_getters.jl

@@ -341,7 +341,7 @@ function get_turbconv_model(
     turbconv_params,
 )
     turbconv = parsed_args["turbconv"]
-    @assert turbconv in (nothing, "edmf", "edmfx", "diagnostic_edmfx")
+    @assert turbconv in (nothing, "edmf", "edmfx", "advective_edmfx", "diagnostic_edmfx")
 
     return if turbconv == "edmf"
         TC.EDMFModel(