rft: smagorinsky-lilly for future h-v splitting

Author: haakon-e

config/default_configs/default_config.yml

@@ -170,8 +170,8 @@ c_amd:
   help: "Model coefficient for AMD-LES closure (TODO: Move to parameters.toml)"
   value: 0.29
 smagorinsky_lilly:
-  help: "Smagorinsky-Lilly diffusive closure [`false` (default), `true`]"
-  value: false
+  help: "Smagorinsky-Lilly diffusive closure [`nothing` (default), `UVW`]"
+  value: ~
 bubble:
   help: "Enable bubble correction for more accurate surface areas"
   value: false

config/model_configs/box_density_current_test.yml

@@ -1,7 +1,7 @@
 reference_job_id: "les_box"
 initial_condition: "DryDensityCurrentProfile"
 config: "box"
-smagorinsky_lilly: true
+smagorinsky_lilly: "UVW"
 discrete_hydrostatic_balance: true
 hyperdiff: "false"
 x_max: 51200.0

config/model_configs/les_isdac_box.yml

@@ -14,7 +14,7 @@ implicit_diffusion: false
 approximate_linear_solve_iters: 2
 hyperdiff: "false"
 apply_limiter: false
-smagorinsky_lilly: true
+smagorinsky_lilly: "UVW"
 # time- and spatial discretization
 x_elem: 10
 x_max: 3.2e3

post_processing/ci_plots.jl

@@ -1310,17 +1310,11 @@ function make_plots(
 
     reduction = "inst"
     short_names = [
-        "wa",
-        "ua",
-        "va",
-        "ta",
-        "thetaa",
-        "ha",
-        "hus",
-        "hur",
-        "cl",
-        "clw",
-        "cli",
+        "wa", "ua", "va", "ta", "thetaa", "ha",
+        "hus", "hur", "cl", "clw", "cli", "ke",
+        "Dh_smag", "strainrateh_smag",  # smag horizontal
+        "Dv_smag", "strainratev_smag",  # smag vertical
+        "edt"  # DecayWithHeight vertical diffusivity
     ]
     short_names = short_names ∩ collect(keys(simdirs[1].vars))
 
@@ -1333,10 +1327,7 @@ function make_plots(
         window_end = last(var.dims["time"])
         window_start = window_end - 2hours
         var_window = ClimaAnalysis.window(
-            var,
-            "time";
-            left = window_start,
-            right = window_end,
+            var, "time"; left = window_start, right = window_end,
         )
         var_reduced = horizontal_average(average_time(var_window))
         return var_reduced

src/cache/precomputed_quantities.jl

@@ -249,10 +249,12 @@ function precomputed_quantities(Y, atmos)
         if atmos.smagorinsky_lilly isa SmagorinskyLilly
             uvw_vec = UVW(FT(0), FT(0), FT(0))
             (;
-                ᶜτ_smag = similar(Y.c, typeof(uvw_vec * uvw_vec')),
-                ᶠτ_smag = similar(Y.f, typeof(uvw_vec * uvw_vec')),
-                ᶜD_smag = similar(Y.c, FT),
-                ᶠD_smag = similar(Y.f, FT),
+                ᶜS = similar(Y.c, typeof(uvw_vec * uvw_vec')),
+                ᶠS = similar(Y.f, typeof(uvw_vec * uvw_vec')),
+                ᶜS_norm_h = similar(Y.c, FT), ᶜS_norm_v = similar(Y.c, FT),
+                ᶜL_h = similar(Y.c, FT), ᶜL_v = similar(Y.c, FT),
+                ᶜνₜ_h = similar(Y.c, FT), ᶜνₜ_v = similar(Y.c, FT),
+                ᶜD_h = similar(Y.c, FT), ᶜD_v = similar(Y.c, FT),
             )
         else
             (;)
@@ -608,9 +610,7 @@ NVTX.@annotate function set_explicit_precomputed_quantities_part2!(Y, p, t)
         set_cloud_fraction!(Y, p, moisture_model, cloud_model)
     end
 
-    if p.atmos.smagorinsky_lilly isa SmagorinskyLilly
-        set_smagorinsky_lilly_precomputed_quantities!(Y, p)
-    end
+    set_smagorinsky_lilly_precomputed_quantities!(Y, p, p.atmos.smagorinsky)
 
     if p.atmos.amd_les isa AnisotropicMinimumDissipation
         set_amd_precomputed_quantities!(Y, p)

src/diagnostics/core_diagnostics.jl

@@ -343,6 +343,46 @@ add_diagnostic_variable!(
     end,
 )
 
+###
+# Smagorinsky Lilly diffusivity
+###
+add_diagnostic_variable!(
+    short_name = "Dh_smag",
+    long_name = "Horizontal smagorinsky diffusivity",
+    units = "m^2 s^-1",
+    compute! = (out, _, cache, _) -> begin
+        (; ᶜD_h) = cache.precomputed
+        isnothing(out) ? copy(ᶜD_h) : (out .= ᶜD_h)
+    end,
+)
+add_diagnostic_variable!(
+    short_name = "Dv_smag",
+    long_name = "Horizontal smagorinsky diffusivity",
+    units = "m^2 s^-1",
+    compute! = (out, _, cache, _) -> begin
+        (; ᶜD_v) = cache.precomputed
+        isnothing(out) ? copy(ᶜD_v) : (out .= ᶜD_v)
+    end,
+)
+add_diagnostic_variable!(
+    short_name = "strainh_smag",
+    long_name = "Horizontal strain rate magnitude (for Smagorinsky)",
+    units = "s",
+    compute! = (out, state, cache, _) -> begin
+        (; ᶜS_norm_h) = cache.precomputed
+        isnothing(out) ? copy(ᶜS_norm_h) : (out .= ᶜS_norm_h)
+    end,
+)
+add_diagnostic_variable!(
+    short_name = "strainv_smag",
+    long_name = "Vertical strain rate magnitude (for Smagorinsky)",
+    units = "s",
+    compute! = (out, state, cache, _) -> begin
+        (; ᶜS_norm_v) = cache.precomputed
+        isnothing(out) ? copy(ᶜS_norm_v) : (out .= ᶜS_norm_v)
+    end,
+)
+
 ###
 # Relative humidity (3d)
 ###