compute prandtl_nvec in prognostic edmfx, add mixing length to plots

post_processing/contours_and_profiles.jl

@@ -410,7 +410,11 @@ function contours_and_profiles(output_dir, ref_job_id = nothing)
     end
     if has_sgs
         if has_moisture
-            contour_variables = (contour_variables..., (:gm, :cloud_fraction))
+            contour_variables =
+                (contour_variables...,
+                (:gm, :cloud_fraction),
+                (:draft, :w_velocity),
+            )
         end
         contour_variables =
             (contour_variables..., (:draft, :area), (:env, :tke))
@@ -448,8 +452,12 @@ function contours_and_profiles(output_dir, ref_job_id = nothing)
             profile_variables =
                 (profile_variables..., ((:gm,), :cloud_fraction))
         end
-        profile_variables =
-            (profile_variables..., ((:draft,), :area), ((:env,), :tke))
+        profile_variables = (
+            profile_variables...,
+            ((:draft,), :area),
+            ((:env,), :tke),
+            ((:env,), :mixing_length),
+        )
     end
 
     ##

src/cache/edmf_precomputed_quantities.jl

@@ -226,9 +226,20 @@ function set_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
         ),
     )
 
-    @. ᶜstrain_rate_norm = $(FT(1e-4))
+    # TODO: Currently the shear production only includes vertical gradients
+    ᶠu⁰ = p.ᶠtemp_C123
+    @. ᶠu⁰ = C123(ᶠinterp(Y.c.uₕ)) + C123(ᶠu³⁰)
+    ᶜstrain_rate = p.ᶜtemp_UVWxUVW
+    compute_strain_rate!(ᶜstrain_rate, ᶠu⁰)
+    @. ᶜstrain_rate_norm = norm_sqr(ᶜstrain_rate)
+
     ᶜprandtl_nvec = p.ᶜtemp_scalar
-    @. ᶜprandtl_nvec = FT(1) / 3
+    @. ᶜprandtl_nvec = turbulent_prandtl_number(
+        params,
+        obukhov_length,
+        ᶜlinear_buoygrad,
+        ᶜstrain_rate_norm,
+    )
     ᶜtke_exch = p.ᶜtemp_scalar_2
     @. ᶜtke_exch = 0
     for j in 1:n