add a simple edmf test case

config/default_configs/default_edmf_config.yml

@@ -11,6 +11,9 @@ gs_tendency:
 edmf_coriolis:
   help: "EDMF coriolis [`nothing` (default), `Bomex`,`LifeCycleTan2018`,`Rico`,`ARM_SGP`,`DYCOMS_RF01`,`DYCOMS_RF02`,`GABLS`]"
   value: ~
+edmfx_velocity_relaxation:
+  help: "If set to true, it switches on the relaxation of negative velocity in EDMFX.  [`true`, `false` (default)]"
+  value: false
 edmfx_nh_pressure:
   help: "If set to true, it switches on EDMFX pressure drag closure.  [`true`, `false` (default)]"
   value: false

config/model_configs/prognostic_edmfx_bomex_box.yml

@@ -23,7 +23,8 @@ y_elem: 2
 z_elem: 60
 z_stretch: false
 perturb_initstate: false
-dt: "40secs"
+dt: "20secs"
+dt: "20secs"
 t_end: "6hours"
 dt_save_state_to_disk: "30mins"
 toml: [toml/prognostic_edmfx_bomex_box.toml]

config/model_configs/prognostic_edmfx_simpleplume_column.yml

@@ -9,6 +9,7 @@ edmfx_detr_model: "Generalized"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
+edmfx_velocity_relaxation: true
 prognostic_tke: true
 moist: "equil"
 config: "column"
@@ -17,7 +18,7 @@ z_elem: 80
 z_stretch: false
 perturb_initstate: false
 dt: "1secs"
-t_end: "6hours"
+t_end: "12hours"
 dt_save_to_disk: "10mins"
 toml: [toml/prognostic_edmfx_simpleplume.toml]
 netcdf_output_at_levels: true

src/cache/precomputed_quantities.jl

@@ -77,7 +77,7 @@ function precomputed_quantities(Y, atmos)
             ᶜ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}),
-            ᶜKᵥʲs = similar(Y.c, NTuple{n, FT}),
+            ᶠKᵥʲs = similar(Y.f, 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}),

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -62,22 +62,24 @@ function set_prognostic_edmf_precomputed_quantities_draft_and_bc!(Y, p, ᶠuₕ
 
     (; ᶜΦ,) = p.core
     (; ᶜspecific, ᶜp, ᶜh_tot, ᶜK) = p.precomputed
-    (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜKᵥʲs, ᶜtsʲs, ᶜρʲs) = p.precomputed
+    (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶠKᵥʲs, ᶜtsʲs, ᶜρʲs) = p.precomputed
     (; ustar, obukhov_length, buoyancy_flux) = p.precomputed.sfc_conditions
+    ᶜinterp_lb = Operators.LeftBiasedF2C()
+    ᶜinterp_rb = Operators.RightBiasedF2C()
 
     for j in 1:n
         ᶜuʲ = ᶜuʲs.:($j)
         ᶠu³ʲ = ᶠu³ʲs.:($j)
         ᶜKʲ = ᶜKʲs.:($j)
-        ᶜKᵥʲ = ᶜKᵥʲs.:($j)
+        ᶠKᵥʲ = ᶠKᵥʲs.:($j)
         ᶠu₃ʲ = Y.f.sgsʲs.:($j).u₃
         ᶜtsʲ = ᶜtsʲs.:($j)
         ᶜρʲ = ᶜρʲs.:($j)
         ᶜmseʲ = Y.c.sgsʲs.:($j).mse
         ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot
 
         set_velocity_quantities!(ᶜuʲ, ᶠu³ʲ, ᶜKʲ, ᶠu₃ʲ, Y.c.uₕ, ᶠuₕ³)
-        @. ᶜKᵥʲ = ᶜinterp(adjoint(CT3(ᶠu₃ʲ)) * ᶠu₃ʲ) / 2
+        @. ᶠKᵥʲ = (adjoint(CT3(ᶠu₃ʲ)) * ᶠu₃ʲ) / 2
         @. ᶜtsʲ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmseʲ - ᶜΦ, ᶜq_totʲ)
         @. ᶜρʲ = TD.air_density(thermo_params, ᶜtsʲ)
 

src/prognostic_equations/advection.jl

@@ -76,7 +76,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     (; ᶜf³, ᶠf¹², ᶜΦ) = p.core
     (; ᶜu, ᶠu³, ᶜK) = p.precomputed
     (; edmfx_upwinding) = n > 0 || advect_tke ? p.atmos.numerics : all_nothing
-    (; ᶜuʲs, ᶜKʲs, ᶜKᵥʲs) = n > 0 ? p.precomputed : all_nothing
+    (; ᶜuʲs, ᶜKʲs, ᶠKᵥʲs) = n > 0 ? p.precomputed : all_nothing
     (; ᶠu³⁰) = advect_tke ? p.precomputed : all_nothing
     (; energy_upwinding, tracer_upwinding) = p.atmos.numerics
     (; ᶜspecific) = p.precomputed
@@ -157,7 +157,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
             for j in 1:n
                 @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] -=
                     ᶠω¹²ʲs.:($$j)[colidx] × ᶠinterp(CT12(ᶜuʲs.:($$j)[colidx])) +
-                    ᶠgradᵥ(ᶜKʲs.:($$j)[colidx] - ᶜKᵥʲs.:($$j)[colidx])
+                    ᶠgradᵥ(ᶜKʲs.:($$j)[colidx] - ᶜinterp(ᶠKᵥʲs.:($$j)[colidx]))
             end
         else
             # deep atmosphere
@@ -174,7 +174,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
                 @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] -=
                     (ᶠf¹²[colidx] + ᶠω¹²ʲs.:($$j)[colidx]) ×
                     ᶠinterp(CT12(ᶜuʲs.:($$j)[colidx])) +
-                    ᶠgradᵥ(ᶜKʲs.:($$j)[colidx] - ᶜKᵥʲs.:($$j)[colidx])
+                    ᶠgradᵥ(ᶜKʲs.:($$j)[colidx] - ᶜinterp(ᶠKᵥʲs.:($$j)[colidx]))
             end
         end
 
@@ -210,17 +210,28 @@ function edmfx_sgs_vertical_advection_tendency!(
     (; dt) = p
     ᶜJ = Fields.local_geometry_field(Y.c).J
     (; edmfx_upwinding) = p.atmos.numerics
-    (; ᶠu³ʲs, ᶜKᵥʲs, ᶜρʲs) = p.precomputed
+    (; ᶠu³ʲs, ᶠKᵥʲs, ᶜρʲs) = p.precomputed
     (; ᶠgradᵥ_ᶜΦ) = p.core
 
     ᶠz = Fields.coordinate_field(Y.f).z
     ᶜa_scalar = p.scratch.ᶜtemp_scalar
+    ᶜu₃ʲ = p.scratch.ᶜtemp_C3
+    ᶜKᵥʲ = p.scratch.ᶜtemp_scalar_2
+    ᶜinterp_lb = Operators.LeftBiasedF2C()
+    ᶜinterp_rb = Operators.RightBiasedF2C()
     for j in 1:n
+        # TODO: Add a biased GradientF2F operator in ClimaCore
+        @. ᶜu₃ʲ[colidx] = ᶜinterp(Y.f.sgsʲs.:($$j).u₃[colidx])
+        @. ᶜKᵥʲ[colidx] = ifelse(
+            ᶜu₃ʲ[colidx].components.data.:1 > 0,
+            ᶜinterp_lb(ᶠKᵥʲs.:($$j)[colidx]),
+            ᶜinterp_rb(ᶠKᵥʲs.:($$j)[colidx]),
+        )
         # For the updraft u_3 equation, we assume the grid-mean to be hydrostatic
         # and calcuate the buoyancy term relative to the grid-mean density.
         @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] -=
             (ᶠinterp(ᶜρʲs.:($$j)[colidx] - Y.c.ρ[colidx]) * ᶠgradᵥ_ᶜΦ[colidx]) /
-            ᶠinterp(ᶜρʲs.:($$j)[colidx]) + ᶠgradᵥ(ᶜKᵥʲs.:($$j)[colidx])
+            ᶠinterp(ᶜρʲs.:($$j)[colidx]) + ᶠgradᵥ(ᶜKᵥʲ[colidx])
 
         # buoyancy term in mse equation
         @. Yₜ.c.sgsʲs.:($$j).mse[colidx] +=

src/prognostic_equations/edmfx_closures.jl

@@ -310,3 +310,25 @@ function turbulent_prandtl_number(
     end
     return prandtl_nvec
 end
+
+edmfx_velocity_relaxation_tendency!(Yₜ, Y, p, t, colidx, turbconv_model) =
+    nothing
+function edmfx_velocity_relaxation_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    colidx,
+    turbconv_model::PrognosticEDMFX,
+)
+
+    n = n_mass_flux_subdomains(turbconv_model)
+    (; dt) = p
+
+    if p.atmos.edmfx_velocity_relaxation
+        for j in 1:n
+            @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] -=
+                C3(min(Y.f.sgsʲs.:($$j).u₃[colidx].components.data.:1, 0)) / dt
+        end
+    end
+end

src/prognostic_equations/remaining_tendency.jl

@@ -61,6 +61,14 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
             p.atmos.turbconv_model,
         )
         edmfx_nh_pressure_tendency!(Yₜ, Y, p, t, colidx, p.atmos.turbconv_model)
+        edmfx_velocity_relaxation_tendency!(
+            Yₜ,
+            Y,
+            p,
+            t,
+            colidx,
+            p.atmos.turbconv_model,
+        )
         edmfx_tke_tendency!(Yₜ, Y, p, t, colidx, p.atmos.turbconv_model)
         edmfx_precipitation_tendency!(
             Yₜ,

src/solver/type_getters.jl

@@ -41,6 +41,9 @@ function get_atmos(config::AtmosConfig, params)
     edmfx_nh_pressure = parsed_args["edmfx_nh_pressure"]
     @assert edmfx_nh_pressure in (false, true)
 
+    edmfx_velocity_relaxation = parsed_args["edmfx_velocity_relaxation"]
+    @assert edmfx_velocity_relaxation in (false, true)
+
     implicit_diffusion = parsed_args["implicit_diffusion"]
     @assert implicit_diffusion in (true, false)
 
@@ -62,6 +65,7 @@ function get_atmos(config::AtmosConfig, params)
         edmfx_sgs_mass_flux,
         edmfx_sgs_diffusive_flux,
         edmfx_nh_pressure,
+        edmfx_velocity_relaxation,
         precip_model,
         cloud_model,
         forcing_type,

src/solver/types.jl

@@ -333,6 +333,7 @@ Base.@kwdef struct AtmosModel{
     ESMF,
     ESDF,
     ENP,
+    EVR,
     TCM,
     NOGW,
     OGW,
@@ -363,6 +364,7 @@ Base.@kwdef struct AtmosModel{
     edmfx_sgs_mass_flux::ESMF = nothing
     edmfx_sgs_diffusive_flux::ESDF = nothing
     edmfx_nh_pressure::ENP = nothing
+    edmfx_velocity_relaxation::EVR = nothing
     turbconv_model::TCM = nothing
     non_orographic_gravity_wave::NOGW = nothing
     orographic_gravity_wave::OGW = nothing