Merge pull request #4097 from CliMA/sa/fix_edmf_diffusion

Fix edmf diffusion;

Author: sajjadazimi

src/prognostic_equations/edmfx_sgs_flux.jl

@@ -390,7 +390,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
     thermo_params = CAP.thermodynamics_params(params)
     turbconv_params = CAP.turbconv_params(params)
     c_d = CAP.tke_diss_coeff(turbconv_params)
-    (; ᶜK, ᶜu⁰, ᶜK⁰, ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜts) = p.precomputed
+    (; ᶜK, ᶜu⁰, ᶜK⁰, ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜρʲs, ᶜts, ᶜts⁰) = p.precomputed
     (; ρatke_flux) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F()
     ᶜρa⁰ = @. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))
@@ -475,6 +475,8 @@ function edmfx_sgs_diffusive_flux_tendency!(
                 top = Operators.SetValue(C3(FT(0))),
                 bottom = Operators.SetValue(C3(FT(0))),
             )
+            ᶜρʲ = ᶜρʲs.:(1)
+            ᶜρ⁰ = @. lazy(TD.air_density(thermo_params, ᶜts⁰))
 
             microphysics_tracers = (
                 (@name(c.ρq_liq), @name(c.sgsʲs.:(1).q_liq), @name(q_liq), FT(1)),
@@ -495,15 +497,19 @@ function edmfx_sgs_diffusive_flux_tendency!(
                 ᶜχ⁰ = ᶜspecific_env_value(χ_name, Y, p)
 
                 # add updraft diffusion
-                @. ᶜρχₜ_diffusion = ᶜdivᵥ_ρq(
-                    -(ᶠinterp(Y.c.sgsʲs.:(1).ρa) * ᶠinterp(ᶜK_h) * α * ᶠgradᵥ(ᶜχʲ)),
-                )
+                @. ᶜρχₜ_diffusion =
+                    draft_area(Y.c.sgsʲs.:(1).ρa, ᶜρʲ) *
+                    ᶜdivᵥ_ρq(
+                        -(ᶠinterp(ᶜρʲ) * ᶠinterp(ᶜK_h) * α * ᶠgradᵥ(ᶜχʲ)),
+                    )
                 @. ᶜρχₜ -= ᶜρχₜ_diffusion
 
                 # add environment diffusion
-                @. ᶜρχₜ_diffusion = ᶜdivᵥ_ρq(
-                    -(ᶠinterp(ᶜρa⁰) * ᶠinterp(ᶜK_h) * α * ᶠgradᵥ(ᶜχ⁰)),
-                )
+                @. ᶜρχₜ_diffusion =
+                    draft_area(ᶜρa⁰, ᶜρ⁰) *
+                    ᶜdivᵥ_ρq(
+                        -(ᶠinterp(ᶜρ⁰) * ᶠinterp(ᶜK_h) * α * ᶠgradᵥ(ᶜχ⁰)),
+                    )
                 @. ᶜρχₜ -= ᶜρχₜ_diffusion
             end
         end

src/prognostic_equations/mass_flux_closures.jl

@@ -129,9 +129,8 @@ function edmfx_vertical_diffusion_tendency!(
 )
     if p.atmos.edmfx_model.vertical_diffusion isa Val{true}
         (; params) = p
-        (; ᶜts, ᶜK, ᶜρʲs) = p.precomputed
+        (; ᶜρʲs) = p.precomputed
         FT = eltype(p.params)
-        thermo_params = CAP.thermodynamics_params(params)
         turbconv_params = CAP.turbconv_params(params)
         n = n_mass_flux_subdomains(turbconv_model)
         ᶜdivᵥ_mse = Operators.DivergenceF2C(
@@ -142,15 +141,6 @@ function edmfx_vertical_diffusion_tendency!(
             top = Operators.SetValue(C3(0)),
             bottom = Operators.SetValue(C3(0)),
         )
-        ᶜinv_ρ̂ = (@. lazy(
-            specific(
-                FT(1),
-                Y.c.sgsʲs.:(1).ρa,
-                FT(0),
-                ᶜρʲs.:(1),
-                turbconv_model,
-            ),
-        ))
 
         (; ᶜlinear_buoygrad, ᶜstrain_rate_norm) = p.precomputed
         ᶜtke⁰ = @. lazy(specific(Y.c.sgs⁰.ρatke, Y.c.ρ))
@@ -165,7 +155,6 @@ function edmfx_vertical_diffusion_tendency!(
 
         for j in 1:n
             ᶜρʲ = ᶜρʲs.:($j)
-            ᶜρaʲ = Y.c.sgsʲs.:($j).ρa
             ᶜmseʲ = Y.c.sgsʲs.:($j).mse
             ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot
             # Note: For this and other diffusive tendencies, we should use ρaʲ instead of ρʲ,
@@ -181,7 +170,7 @@ function edmfx_vertical_diffusion_tendency!(
             p.atmos.microphysics_model isa Microphysics2Moment
         )
             α_precip = CAP.α_vert_diff_tracer(params)
-            ᶜρaʲ = Y.c.sgsʲs.:(1).ρa
+            ᶜρʲ = ᶜρʲs.:(1)
             ᶜdivᵥ_q = Operators.DivergenceF2C(
                 top = Operators.SetValue(C3(FT(0))),
                 bottom = Operators.SetValue(C3(FT(0))),
@@ -205,8 +194,7 @@ function edmfx_vertical_diffusion_tendency!(
                 ᶜχʲ = MatrixFields.get_field(Y, χʲ_name)
                 ᶜχʲₜ = MatrixFields.get_field(Yₜ, χʲ_name)
 
-                @. ᶜχʲₜ -=
-                    ᶜinv_ρ̂ * ᶜdivᵥ_q(-(ᶠinterp(ᶜρaʲ) * ᶠinterp(ᶜK_h) * α * ᶠgradᵥ(ᶜχʲ)))
+                @. ᶜχʲₜ -= ᶜdivᵥ_q(-(ᶠinterp(ᶜρʲ) * ᶠinterp(ᶜK_h) * α * ᶠgradᵥ(ᶜχʲ))) / ᶜρʲ
             end
         end
     end