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

@@ -169,8 +169,14 @@ amd_les:
 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`]"
+smagorinsky:
+  help: "Smagorinsky diffusive closure in the horizontal and vertical [`false` (default), `true`]"
+  value: false
+smagorinsky_horizontal:
+  help: "Smagorinsky diffusive closure in the horizontal [`false` (default), `true`]"
+  value: false
+smagorinsky_vertical:
+  help: "Smagorinsky diffusive closure in the vertical [`false` (default), `true`]"
   value: false
 bubble:
   help: "Enable bubble correction for more accurate surface areas"

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_horizontal isa SmagorinskyLilly ||
+           atmos.smagorinsky_vertical 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),
+                ᶜD_h = similar(Y.c, FT), ᶜD_v = similar(Y.c, FT),
+                ᶜstrain_rate_norm_h = similar(Y.c, FT),
+                ᶜstrain_rate_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,7 +612,8 @@ 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
+    if p.atmos.smagorinsky_horizontal isa SmagorinskyLilly ||
+       p.atmos.smagorinsky_vertical isa SmagorinskyLilly
         set_smagorinsky_lilly_precomputed_quantities!(Y, p)
     end
 

src/diagnostics/core_diagnostics.jl

@@ -343,6 +343,65 @@ add_diagnostic_variable!(
     end,
 )
 
+###
+# Smagorinski Lilly diffusivity
+###
+add_diagnostic_variable!(
+    short_name = "Dh_smag",
+    long_name = "Horizontal smagorinski diffusivity",
+    units = "m^2 s^-1",
+    compute! = (out, _, cache, _) -> begin
+        isnothing(out) ? copy(cache.precomputed.ᶜD_h) : (out .= cache.precomputed.ᶜD_h)
+    end,
+)
+add_diagnostic_variable!(
+    short_name = "Dv_smag",
+    long_name = "Horizontal smagorinski diffusivity",
+    units = "m^2 s^-1",
+    compute! = (out, _, cache, _) -> begin
+        isnothing(out) ? copy(cache.precomputed.ᶜD_v) : (out .= cache.precomputed.ᶜD_v)
+    end,
+)
+add_diagnostic_variable!(
+    short_name = "dtmaxh_smag",
+    long_name = "Max horizontal diffusion dt",
+    units = "s",
+    compute! = (out, state, cache, _) -> begin
+        Δh = Spaces.node_horizontal_length_scale(Spaces.horizontal_space(axes(state.f)))
+        dtmax = @. lazy(Δh^2 / cache.precomputed.ᶜD_h)
+        isnothing(out) ? copy(dtmax) : (out .= dtmax)
+    end,
+)
+add_diagnostic_variable!(
+    short_name = "dtmaxv_smag",
+    long_name = "Max diffusion dt",
+    units = "s",
+    compute! = (out, state, cache, _) -> begin
+        ᶜΔ_z = Fields.Δz_field(state.c)
+        dtmax = @. lazy(ᶜΔ_z^2 / cache.precomputed.ᶜD_v)
+        isnothing(out) ? copy(dtmax) : (out .= dtmax)
+    end,
+)
+import ..projected_strain_rate_norm
+add_diagnostic_variable!(
+    short_name = "strainrateh_smag",
+    long_name = "Strain rate horizontal",
+    units = "s",
+    compute! = (out, state, cache, _) -> begin
+        ᶜS_norm_h = projected_strain_rate_norm(cache.precomputed.ᶜS, Geometry.UVAxis())
+        isnothing(out) ? copy(ᶜS_norm_h) : (out .= ᶜS_norm_h)
+    end,
+)
+add_diagnostic_variable!(
+    short_name = "strainratev_smag",
+    long_name = "Strain rate vertical",
+    units = "s",
+    compute! = (out, state, cache, _) -> begin
+        ᶜS_norm_v = projected_strain_rate_norm(cache.precomputed.ᶜS, Geometry.WAxis())
+        isnothing(out) ? copy(ᶜS_norm_v) : (out .= ᶜS_norm_v)
+    end,
+)
+
 ###
 # Relative humidity (3d)
 ###