rft: split horizontal and vertical smagorinsky

- add diffusivity diagnostics
- implicit vertical smagorinsky
- les ci_plots edits
- use vertical strain for lilly correction
- full strain rate method

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`, `UV`, `W`, `UV_W_decoupled`] (!! experimental !!)"
+  value: ~
 bubble:
   help: "Enable bubble correction for more accurate surface areas"
   value: false

post_processing/ci_plots.jl

@@ -1310,17 +1310,10 @@ 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", "dtmaxh_smag", "strainrateh_smag",  # smag horizontal
+        "Dv_smag", "dtmaxv_smag", "strainratev_smag",  # smag vertical
+        "edt"  # DecayWithHeight vertical diffusivity
     ]
     short_names = short_names ∩ collect(keys(simdirs[1].vars))
 
@@ -1333,10 +1326,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

@@ -245,14 +245,17 @@ function precomputed_quantities(Y, atmos)
             precipitation_sgs_quantities...,
             diagnostic_precipitation_sgs_quantities...,
         ) : (;)
-    smagorinsky_lilly_quantities =
-        if atmos.smagorinsky_lilly isa SmagorinskyLilly
+    smagorinsky_quantities =
+        if atmos.smagorinsky 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')),
+                ᶜ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),
+                ᶜS_norm_h = similar(Y.c, FT),
+                ᶜS_norm_v = similar(Y.c, FT),
             )
         else
             (;)
@@ -281,8 +284,9 @@ function precomputed_quantities(Y, atmos)
         precipitation_quantities...,
         surface_precip_fluxes...,
         cloud_diagnostics_tuple,
-        smagorinsky_lilly_quantities...,
-        amd_les_quantities...)
+        smagorinsky_quantities...,
+        amd_les_quantities...,
+    )
 end
 
 # Interpolates the third contravariant component of Y.c.uₕ to cell faces.
@@ -608,9 +612,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 = "strainrateh_smag",
+    long_name = "Strain rate horizontal 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 = "strainratev_smag",
+    long_name = "Strain rate vertical 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)
 ###