wip

CliMA · Oct 10, 2023 · d061ef0 · d061ef0
1 parent f8853c4
commit d061ef0
Show file tree

Hide file tree

Showing 6 changed files with 598 additions and 495 deletions.
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
diff --git a/config/model_configs/advective_edmfx_bomex_box.yml b/config/model_configs/advective_edmfx_bomex_box.yml
@@ -1,4 +1,4 @@
-job_id: "advectiveedmfx_bomex_box"
+job_id: "advective_edmfx_bomex_box"
 initial_condition: "Bomex"
 subsidence: "Bomex"
 edmf_coriolis: "Bomex"

diff --git a/config/model_configs/advective_edmfx_gabls_box.yml b/config/model_configs/advective_edmfx_gabls_box.yml
@@ -0,0 +1,30 @@
+job_id: "advective_edmfx_gabls_box"
+initial_condition: GABLS
+edmf_coriolis: GABLS
+surface_setup: GABLS
+turbconv: "advective_edmfx"
+edmfx_upwinding: "first_order"
+edmfx_entr_model: "ConstantCoefficient"
+edmfx_detr_model: "ConstantCoefficient"
+edmfx_sgs_mass_flux: false
+edmfx_sgs_diffusive_flux: false
+edmfx_nh_pressure: false
+prognostic_tke: false
+moist: "equil"
+config: "box"
+hyperdiff: "true"
+vert_diff: "true"
+kappa_4: 1e12
+x_max: 1e5
+y_max: 1e5
+z_max: 400
+x_elem: 2
+y_elem: 2
+z_elem: 8
+z_stretch: false
+dt: "5secs"
+t_end: "9hours"
+dt_save_to_disk: "10mins"
+perturb_initstate: false
+FLOAT_TYPE: "Float64"
+toml: [toml/edmfx_box.toml]
diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl
@@ -305,14 +305,16 @@ NVTX.@annotate function compute_diagnostics(integrator)
             ) .+ ᶜa⁺ .* cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺),
         )
     elseif p.atmos.turbconv_model isa AdvectiveEDMFX
-        (; ᶜtke⁰, ᶜmixing_length) = p
-        (; ᶜu⁺, ᶜts⁺, ᶜa⁺) = output_advective_sgs_quantities(Y, p, t)
+        (; ᶜtke⁰, ᶜu⁰, ᶜts⁰, ᶜmixing_length) = p
+        (; ᶜu⁺, ᶜts⁺, ᶜa⁺, ᶜa⁰) = output_advective_sgs_quantities(Y, p, t)
         env_diagnostics = (;
+            common_diagnostics(p, ᶜu⁰, ᶜts⁰)...,
+            area = ᶜa⁰,
             cloud_fraction = get_cloud_fraction.(
                 thermo_params,
                 env_thermo_quad,
                 ᶜp,
-                ᶜts,
+                ᶜts⁰,
             ),
             tke = ᶜtke⁰,
             mixing_length = ᶜmixing_length,
@@ -325,11 +327,12 @@ NVTX.@annotate function compute_diagnostics(integrator)
         turbulence_convection_diagnostic = (;
             add_prefix(env_diagnostics, :env_)...,
             add_prefix(draft_diagnostics, :draft_)...,
-            cloud_fraction = get_cloud_fraction.(
+            cloud_fraction = ᶜa⁰ .*
+                             get_cloud_fraction.(
                 thermo_params,
                 env_thermo_quad,
                 ᶜp,
-                ᶜts,
+                ᶜts⁰,
             ) .+ ᶜa⁺ .* cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺),
         )
     elseif p.atmos.turbconv_model isa DiagnosticEDMFX

diff --git a/src/prognostic_equations/edmfx_sgs_flux.jl b/src/prognostic_equations/edmfx_sgs_flux.jl
@@ -230,15 +230,67 @@ function edmfx_sgs_diffusive_flux_tendency!(
     return nothing
 end
 
-# FIXME: add SGS diffusive flux
-edmfx_sgs_diffusive_flux_tendency!(
+function edmfx_sgs_diffusive_flux_tendency!(
     Yₜ,
     Y,
     p,
     t,
     colidx,
     turbconv_model::AdvectiveEDMFX,
-) = nothing
+)
+
+    FT = Spaces.undertype(axes(Y.c))
+    (; sfc_conditions) = p
+    (; ᶜρa⁰, ᶠu³⁰, ᶜu⁰, ᶜh_tot⁰, ᶜq_tot⁰) = p
+    (; ᶜK_u, ᶜK_h) = p
+    ᶠgradᵥ = Operators.GradientC2F()
+
+    if p.atmos.edmfx_sgs_diffusive_flux
+        # energy
+        ᶠρaK_h = p.ᶠtemp_scalar
+        @. ᶠρaK_h[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * ᶠinterp(ᶜK_h[colidx])
+
+        ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
+            top = Operators.SetValue(C3(FT(0))),
+            bottom = Operators.SetValue(sfc_conditions.ρ_flux_h_tot[colidx]),
+        )
+        @. Yₜ.c.ρe_tot[colidx] -=
+            ᶜdivᵥ_ρe_tot(-(ᶠρaK_h[colidx] * ᶠgradᵥ(ᶜh_tot⁰[colidx])))
+
+        if !(p.atmos.moisture_model isa DryModel)
+            # specific humidity
+            ᶜdivᵥ_ρq_tot = Operators.DivergenceF2C(
+                top = Operators.SetValue(C3(FT(0))),
+                bottom = Operators.SetValue(
+                    sfc_conditions.ρ_flux_q_tot[colidx],
+                ),
+            )
+            @. Yₜ.c.ρq_tot[colidx] -= ᶜdivᵥ_ρq_tot(
+                -(ᶠρaK_h[colidx] * ᶠgradᵥ(ᶜq_tot⁰[colidx])),
+            )
+        end
+
+        # momentum
+        ᶠρaK_u = p.ᶠtemp_scalar
+        @. ᶠρaK_u[colidx] = ᶠinterp(ᶜρa⁰[colidx]) * ᶠinterp(ᶜK_u[colidx])
+        ᶠstrain_rate = p.ᶠtemp_UVWxUVW
+        compute_strain_rate_face!(ᶠstrain_rate[colidx], ᶜu⁰[colidx])
+        @. Yₜ.c.uₕ[colidx] -= C12(
+            ᶜdivᵥ(-(2 * ᶠρaK_u[colidx] * ᶠstrain_rate[colidx])) / Y.c.ρ[colidx],
+        )
+        # apply boundary condition for momentum flux
+        ᶜdivᵥ_uₕ = Operators.DivergenceF2C(
+            top = Operators.SetValue(C3(FT(0)) ⊗ C12(FT(0), FT(0))),
+            bottom = Operators.SetValue(sfc_conditions.ρ_flux_uₕ[colidx]),
+        )
+        @. Yₜ.c.uₕ[colidx] -=
+            ᶜdivᵥ_uₕ(-(FT(0) * ᶠgradᵥ(Y.c.uₕ[colidx]))) / Y.c.ρ[colidx]
+    end
+
+    # TODO: Add tracer flux
+
+    return nothing
+end
 
 function edmfx_sgs_diffusive_flux_tendency!(
     Yₜ,

diff --git a/toml/edmfx_box.toml b/toml/edmfx_box.toml
@@ -16,5 +16,5 @@ description = "Minimum area fraction per updraft. Parameter not described in the
 
 [entr_coeff]
 alias = "entr_coeff"
-value = 1.0e-1
+value = 1.0
 type = "float"