From 69ac28e8aa990fd3ebc3b79a6a6f4d6aae3a9619 Mon Sep 17 00:00:00 2001
From: Akshay Sridhar <asridhar@clima.gps.caltech.edu>
Date: Mon, 3 Jun 2024 10:43:00 -0700
Subject: [PATCH] implement smagorinsky lilly closure

---
 .gitignore                                    |   2 +
 config/default_configs/default_config.yml     |   9 +
 .../box_density_current_test.yml              |  28 +-
 config/model_configs/les_isdac_box.yml        |   6 +-
 examples/Manifest.toml                        | 294 ++++++++++--------
 examples/Project.toml                         |   2 +-
 post_processing/ci_plots.jl                   |  98 +++++-
 src/ClimaAtmos.jl                             |   7 +
 src/cache/precomputed_quantities.jl           |  18 ++
 src/cache/temporary_quantities.jl             |  15 +-
 .../les_sgs_models/smagorinsky_lilly.jl       | 174 +++++++++++
 .../radiation/radiation.jl                    |   4 +-
 .../remaining_tendency.jl                     |   4 +
 src/solver/model_getters.jl                   |  12 +
 src/solver/type_getters.jl                    |   5 +
 src/solver/types.jl                           |   8 +
 .../gravity_wave/gw_plotutils.jl              |   9 +-
 .../les_closures/smagorinsky_lilly.jl         |  69 ++++
 toml/{isdac_box.toml => les_isdac.toml}       |   0
 19 files changed, 596 insertions(+), 168 deletions(-)
 create mode 100644 src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl
 create mode 100644 test/parameterized_tendencies/les_closures/smagorinsky_lilly.jl
 rename toml/{isdac_box.toml => les_isdac.toml} (100%)

diff --git a/.gitignore b/.gitignore
index 4cf65a220c..ed7ca00073 100644
--- a/.gitignore
+++ b/.gitignore
@@ -32,6 +32,8 @@ docs/src/config.md
 test/Manifest*.toml
 calibration/test/Manifest*.toml
 
+Artifacts.toml
+
 # ignore vscode artifacts
 *.vscode
 *.DS
diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml
index 708cdf0801..53b82d41ec 100644
--- a/config/default_configs/default_config.yml
+++ b/config/default_configs/default_config.yml
@@ -145,6 +145,15 @@ vert_diff:
 hyperdiff:
   help: "Hyperdiffusion [`ClimaHyperdiffusion` (or `true`; default), `none` (or `false`)]"
   value: "CAM_SE"
+smagorinsky_lilly:
+  help: "Smagorinsky-Lilly diffusive closure [`false` (default), `true`]"
+  value: false
+c_smag: 
+  help: "Smagorinsky coefficient"
+  value: 0.2
+prandtl_turbulent_neutral:
+  help: "Turbulent Prandtl number for neutral stratification"
+  value: 0.333
 bubble:
   help: "Enable bubble correction for more accurate surface areas"
   value: true
diff --git a/config/model_configs/box_density_current_test.yml b/config/model_configs/box_density_current_test.yml
index 0c3e346d21..312b536640 100644
--- a/config/model_configs/box_density_current_test.yml
+++ b/config/model_configs/box_density_current_test.yml
@@ -1,16 +1,22 @@
-dt_save_state_to_disk: "10secs"
+reference_job_id: "les_box"
 initial_condition: "DryDensityCurrentProfile"
+config: "box"
+smagorinsky_lilly: true
+c_smag: 0.25
+discrete_hydrostatic_balance: true
+hyperdiff: "false"
 x_max: 51200.0
+y_max: 51200.0
+z_max: 6400.0
+x_elem: 45
+y_elem: 45
 z_elem: 45
-dt: "0.1secs"
-t_end: "10.0secs"
-discrete_hydrostatic_balance: true
-y_max: 6400.0
-y_elem: 15
 z_stretch: false
-x_elem: 45
-config: "box"
-z_max: 6400.0
+dt: "0.25secs"
+t_end: "20secs"
+prandtl_turbulent_neutral: 0.333
+dt_save_state_to_disk: "5secs"
+netcdf_interpolation_num_points: [40, 40, 80]
 diagnostics:
-  - short_name: thetaa
-    period: 10secs
+  - short_name: [wa, ua, va, ta, thetaa, ha]
+    period: 5secs 
diff --git a/config/model_configs/les_isdac_box.yml b/config/model_configs/les_isdac_box.yml
index 06bdbf55a9..5a6f4c9a59 100644
--- a/config/model_configs/les_isdac_box.yml
+++ b/config/model_configs/les_isdac_box.yml
@@ -1,4 +1,5 @@
 # ISDAC config
+reference_job_id: "les_box"  # for plotting
 initial_condition: "ISDAC" 
 subsidence: "ISDAC" 
 surface_setup: "ISDAC"
@@ -15,6 +16,7 @@ hyperdiff: "false"
 apply_limiter: false
 smagorinsky_lilly: true
 c_smag: 0.20
+prandtl_turbulent_neutral: 0.333
 # time- and spatial discretization
 x_elem: 10
 x_max: 3.2e3
@@ -24,11 +26,13 @@ z_elem: 15
 z_max: 2.5e3
 z_stretch: false
 rayleigh_sponge: true
-toml: [toml/isdac_box.toml]  # sponge height
+toml: [toml/les_isdac.toml]  # sponge height
 ode_algo: "SSPKnoth"
 dt: "0.05secs" 
 t_end: "2mins" 
 dt_cloud_fraction: "10mins"
+restart_file: "/groups/esm/hervik/climaatmos-les-artifacts/les_isdac_day0.0.hdf5"
+# diagnostics
 dt_save_state_to_disk: "1mins"
 netcdf_interpolation_num_points: [30, 30, 150]
 diagnostics:
diff --git a/examples/Manifest.toml b/examples/Manifest.toml
index 909a799dcd..ba93fa277f 100644
--- a/examples/Manifest.toml
+++ b/examples/Manifest.toml
@@ -1,13 +1,13 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.6"
+julia_version = "1.10.5"
 manifest_format = "2.0"
-project_hash = "b0cd3fa0e3153549529fa43102bd0bdae16ff23d"
+project_hash = "13bc74579bbfd1b4bab2e5d6cddd05f9a01dba06"
 
 [[deps.ADTypes]]
-git-tree-sha1 = "30bb95a372787af850addf28ac937f1be7b79173"
+git-tree-sha1 = "72af59f5b8f09faee36b4ec48e014a79210f2f4f"
 uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
-version = "1.10.0"
+version = "1.11.0"
 weakdeps = ["ChainRulesCore", "ConstructionBase", "EnzymeCore"]
 
     [deps.ADTypes.extensions]
@@ -79,9 +79,9 @@ version = "1.1.3"
 
 [[deps.Animations]]
 deps = ["Colors"]
-git-tree-sha1 = "e81c509d2c8e49592413bfb0bb3b08150056c79d"
+git-tree-sha1 = "e092fa223bf66a3c41f9c022bd074d916dc303e7"
 uuid = "27a7e980-b3e6-11e9-2bcd-0b925532e340"
-version = "0.4.1"
+version = "0.4.2"
 
 [[deps.ArgParse]]
 deps = ["Logging", "TextWrap"]
@@ -95,9 +95,9 @@ version = "1.1.1"
 
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra"]
-git-tree-sha1 = "d60a1922358aa203019b7857a2c8c37329b8736c"
+git-tree-sha1 = "d5140b60b87473df18cf4fe66382b7c3596df047"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.17.0"
+version = "7.17.1"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
@@ -201,9 +201,9 @@ version = "0.1.6"
 
 [[deps.BlockArrays]]
 deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra"]
-git-tree-sha1 = "d434647f798823bcae510aee0bc0401927f64391"
+git-tree-sha1 = "62551a412126a8c979febd1a05a986928719639e"
 uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-version = "1.1.1"
+version = "1.2.0"
 weakdeps = ["BandedMatrices"]
 
     [deps.BlockArrays.extensions]
@@ -254,10 +254,10 @@ uuid = "4e9b3aee-d8a1-5a3d-ad8b-7d824db253f0"
 version = "1.0.1+0"
 
 [[deps.CUDA]]
-deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "Crayons", "DataFrames", "ExprTools", "GPUArrays", "GPUCompiler", "KernelAbstractions", "LLVM", "LLVMLoopInfo", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "NVTX", "Preferences", "PrettyTables", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "StaticArrays", "Statistics"]
-git-tree-sha1 = "cdbdca28f19c2c7fcf34ffb48bfdaca404dcd18a"
+deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "Crayons", "DataFrames", "ExprTools", "GPUArrays", "GPUCompiler", "KernelAbstractions", "LLVM", "LLVMLoopInfo", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "NVTX", "Preferences", "PrettyTables", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "StaticArrays", "Statistics", "demumble_jll"]
+git-tree-sha1 = "e0725a467822697171af4dae15cec10b4fc19053"
 uuid = "052768ef-5323-5732-b1bb-66c8b64840ba"
-version = "5.5.0"
+version = "5.5.2"
 weakdeps = ["ChainRulesCore", "EnzymeCore", "SpecialFunctions"]
 
     [deps.CUDA.extensions]
@@ -267,9 +267,9 @@ weakdeps = ["ChainRulesCore", "EnzymeCore", "SpecialFunctions"]
 
 [[deps.CUDA_Driver_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "ccd1e54610c222fadfd4737dac66bff786f63656"
+git-tree-sha1 = "14996d716a2eaaeccfc8d7bc854dd87fde720ac1"
 uuid = "4ee394cb-3365-5eb0-8335-949819d2adfc"
-version = "0.10.3+0"
+version = "0.10.4+0"
 
 [[deps.CUDA_Runtime_Discovery]]
 deps = ["Libdl"]
@@ -279,9 +279,9 @@ version = "0.3.5"
 
 [[deps.CUDA_Runtime_jll]]
 deps = ["Artifacts", "CUDA_Driver_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"]
-git-tree-sha1 = "ed8a056a88f5b852df94046060f6770a57334728"
+git-tree-sha1 = "17f1536c600133f7c4113bae0a2d98dbf27c7ebc"
 uuid = "76a88914-d11a-5bdc-97e0-2f5a05c973a2"
-version = "0.15.4+0"
+version = "0.15.5+0"
 
 [[deps.Cairo]]
 deps = ["Cairo_jll", "Colors", "Glib_jll", "Graphics", "Libdl", "Pango_jll"]
@@ -290,10 +290,10 @@ uuid = "159f3aea-2a34-519c-b102-8c37f9878175"
 version = "1.1.1"
 
 [[deps.CairoMakie]]
-deps = ["CRC32c", "Cairo", "Colors", "FileIO", "FreeType", "GeometryBasics", "LinearAlgebra", "Makie", "PrecompileTools"]
-git-tree-sha1 = "d69c7593fe9d7d617973adcbe4762028c6899b2c"
+deps = ["CRC32c", "Cairo", "Cairo_jll", "Colors", "FileIO", "FreeType", "GeometryBasics", "LinearAlgebra", "Makie", "PrecompileTools"]
+git-tree-sha1 = "c3161fbfe99d9d7ee121cf2017d49966b152857c"
 uuid = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
-version = "0.11.11"
+version = "0.12.16"
 
 [[deps.Cairo_jll]]
 deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
@@ -312,14 +312,14 @@ weakdeps = ["SparseArrays"]
     ChainRulesCoreSparseArraysExt = "SparseArrays"
 
 [[deps.ClimaAnalysis]]
-deps = ["Dates", "Interpolations", "NCDatasets", "NaNStatistics", "OrderedCollections", "Reexport", "Statistics"]
-git-tree-sha1 = "e13d742cd5a5ad287cf72ab22cb0927780cfe596"
+deps = ["Artifacts", "Dates", "Interpolations", "NCDatasets", "NaNStatistics", "OrderedCollections", "Reexport", "Statistics", "Unitful"]
+git-tree-sha1 = "b67f8f5f754fde6132bae1a2add06e51cf26227b"
 uuid = "29b5916a-a76c-4e73-9657-3c8fd22e65e6"
-version = "0.5.7"
+version = "0.5.12"
 
     [deps.ClimaAnalysis.extensions]
-    GeoMakieExt = "GeoMakie"
-    MakieExt = "Makie"
+    ClimaAnalysisGeoMakieExt = "GeoMakie"
+    ClimaAnalysisMakieExt = "Makie"
 
     [deps.ClimaAnalysis.weakdeps]
     GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"
@@ -417,15 +417,14 @@ version = "0.7.38"
     StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [[deps.ClimaUtilities]]
-deps = ["Artifacts", "Dates"]
-git-tree-sha1 = "f99f60059af1de49a24d16049dfa8d65e3d3907d"
+deps = ["Artifacts", "ClimaComms", "Dates"]
+git-tree-sha1 = "ae6e4d3223ecb29d3f3d9c5f82b622bae12f3bff"
 uuid = "b3f4f4ca-9299-4f7f-bd9b-81e1242a7513"
-version = "0.1.18"
-weakdeps = ["Adapt", "CUDA", "ClimaComms", "ClimaCore", "ClimaCoreTempestRemap", "Interpolations", "NCDatasets"]
+version = "0.1.19"
+weakdeps = ["Adapt", "CUDA", "ClimaCore", "ClimaCoreTempestRemap", "Interpolations", "NCDatasets"]
 
     [deps.ClimaUtilities.extensions]
-    ClimaUtilitiesClimaCommsCUDAExt = ["ClimaComms", "CUDA"]
-    ClimaUtilitiesClimaCommsExt = "ClimaComms"
+    ClimaUtilitiesCUDAExt = "CUDA"
     ClimaUtilitiesClimaCoreExt = "ClimaCore"
     ClimaUtilitiesClimaCoreInterpolationsExt = ["ClimaCore", "Interpolations"]
     ClimaUtilitiesClimaCoreNCDatasetsExt = ["ClimaCore", "NCDatasets"]
@@ -440,9 +439,9 @@ version = "0.1.13"
 
 [[deps.CloudMicrophysics]]
 deps = ["ClimaParams", "DocStringExtensions", "ForwardDiff", "HCubature", "LazyArtifacts", "QuadGK", "RootSolvers", "SpecialFunctions", "Thermodynamics"]
-git-tree-sha1 = "8f93d0c730dd4f41fade18a1954abf1eb5ac69ce"
+git-tree-sha1 = "e3b2ae212b68aea23d11c03300abb1268a56af87"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
-version = "0.22.4"
+version = "0.22.5"
 
     [deps.CloudMicrophysics.extensions]
     EmulatorModelsExt = ["DataFrames", "MLJ"]
@@ -519,6 +518,11 @@ git-tree-sha1 = "cda2cfaebb4be89c9084adaca7dd7333369715c5"
 uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
 version = "0.3.1"
 
+[[deps.CommonWorldInvalidations]]
+git-tree-sha1 = "ae52d1c52048455e85a387fbee9be553ec2b68d0"
+uuid = "f70d9fcc-98c5-4d4a-abd7-e4cdeebd8ca8"
+version = "1.0.0"
+
 [[deps.Compat]]
 deps = ["TOML", "UUIDs"]
 git-tree-sha1 = "8ae8d32e09f0dcf42a36b90d4e17f5dd2e4c4215"
@@ -608,10 +612,10 @@ deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
 [[deps.DelaunayTriangulation]]
-deps = ["AdaptivePredicates", "EnumX", "ExactPredicates", "PrecompileTools", "Random"]
-git-tree-sha1 = "89df54fbe66e5872d91d8c2cd3a375f660c3fd64"
+deps = ["AdaptivePredicates", "EnumX", "ExactPredicates", "Random"]
+git-tree-sha1 = "e1371a23fd9816080c828d0ce04373857fe73d33"
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
-version = "1.6.1"
+version = "1.6.3"
 
 [[deps.DelimitedFiles]]
 deps = ["Mmap"]
@@ -620,10 +624,10 @@ uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 version = "1.9.1"
 
 [[deps.DiffEqBase]]
-deps = ["ArrayInterface", "ConcreteStructs", "DataStructures", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PreallocationTools", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Static", "StaticArraysCore", "Statistics", "Tricks", "TruncatedStacktraces"]
-git-tree-sha1 = "d1e8a4642e28b0945bde6e2e1ac569b9e0abd728"
+deps = ["ArrayInterface", "ConcreteStructs", "DataStructures", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "FastPower", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PreallocationTools", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "Setfield", "Static", "StaticArraysCore", "Statistics", "TruncatedStacktraces"]
+git-tree-sha1 = "b7dbeaa770bad0980ddddf606de814cff2acb3bc"
 uuid = "2b5f629d-d688-5b77-993f-72d75c75574e"
-version = "6.151.5"
+version = "6.160.0"
 
     [deps.DiffEqBase.extensions]
     DiffEqBaseCUDAExt = "CUDA"
@@ -635,6 +639,7 @@ version = "6.151.5"
     DiffEqBaseMeasurementsExt = "Measurements"
     DiffEqBaseMonteCarloMeasurementsExt = "MonteCarloMeasurements"
     DiffEqBaseReverseDiffExt = "ReverseDiff"
+    DiffEqBaseSparseArraysExt = "SparseArrays"
     DiffEqBaseTrackerExt = "Tracker"
     DiffEqBaseUnitfulExt = "Unitful"
 
@@ -648,6 +653,7 @@ version = "6.151.5"
     Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
     MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
     Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
     Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
@@ -712,9 +718,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 version = "1.0.4"
 
 [[deps.EnzymeCore]]
-git-tree-sha1 = "8f205a601760f4798a10f138c3940f0451d95188"
+git-tree-sha1 = "2c366bfe21936e8f44b607eb86a1f4c110f0d841"
 uuid = "f151be2c-9106-41f4-ab19-57ee4f262869"
-version = "0.7.8"
+version = "0.8.7"
 weakdeps = ["Adapt"]
 
     [deps.EnzymeCore.extensions]
@@ -788,6 +794,27 @@ git-tree-sha1 = "fd923962364b645f3719855c88f7074413a6ad92"
 uuid = "442a2c76-b920-505d-bb47-c5924d526838"
 version = "1.0.2"
 
+[[deps.FastPower]]
+git-tree-sha1 = "58c3431137131577a7c379d00fea00be524338fb"
+uuid = "a4df4552-cc26-4903-aec0-212e50a0e84b"
+version = "1.1.1"
+
+    [deps.FastPower.extensions]
+    FastPowerEnzymeExt = "Enzyme"
+    FastPowerForwardDiffExt = "ForwardDiff"
+    FastPowerMeasurementsExt = "Measurements"
+    FastPowerMonteCarloMeasurementsExt = "MonteCarloMeasurements"
+    FastPowerReverseDiffExt = "ReverseDiff"
+    FastPowerTrackerExt = "Tracker"
+
+    [deps.FastPower.weakdeps]
+    Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
+    ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+    Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
+    MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
+    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+    Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
+
 [[deps.FastRounding]]
 deps = ["ErrorfreeArithmetic", "LinearAlgebra"]
 git-tree-sha1 = "6344aa18f654196be82e62816935225b3b9abe44"
@@ -796,9 +823,15 @@ version = "0.3.1"
 
 [[deps.FileIO]]
 deps = ["Pkg", "Requires", "UUIDs"]
-git-tree-sha1 = "91e0e5c68d02bcdaae76d3c8ceb4361e8f28d2e9"
+git-tree-sha1 = "2dd20384bf8c6d411b5c7370865b1e9b26cb2ea3"
 uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.16.5"
+version = "1.16.6"
+
+    [deps.FileIO.extensions]
+    HTTPExt = "HTTP"
+
+    [deps.FileIO.weakdeps]
+    HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 
 [[deps.FilePaths]]
 deps = ["FilePathsBase", "MacroTools", "Reexport", "Requires"]
@@ -867,15 +900,15 @@ version = "4.1.1"
 
 [[deps.FreeType2_jll]]
 deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Zlib_jll"]
-git-tree-sha1 = "5c1d8ae0efc6c2e7b1fc502cbe25def8f661b7bc"
+git-tree-sha1 = "fa8e19f44de37e225aa0f1695bc223b05ed51fb4"
 uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
-version = "2.13.2+0"
+version = "2.13.3+0"
 
 [[deps.FreeTypeAbstraction]]
 deps = ["ColorVectorSpace", "Colors", "FreeType", "GeometryBasics"]
-git-tree-sha1 = "77e2b094e61d939f9626181ab23d0b76e78f9fd3"
+git-tree-sha1 = "d52e255138ac21be31fa633200b65e4e71d26802"
 uuid = "663a7486-cb36-511b-a19d-713bb74d65c9"
-version = "0.10.5"
+version = "0.10.6"
 
 [[deps.FriBidi_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -898,11 +931,6 @@ version = "0.1.3"
 deps = ["Random"]
 uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
 
-[[deps.GMP_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"
-version = "6.2.1+6"
-
 [[deps.GPUArrays]]
 deps = ["Adapt", "GPUArraysCore", "LLVM", "LinearAlgebra", "Printf", "Random", "Reexport", "Serialization", "Statistics"]
 git-tree-sha1 = "62ee71528cca49be797076a76bdc654a170a523e"
@@ -963,21 +991,15 @@ version = "0.1.0"
 
 [[deps.Glib_jll]]
 deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"]
-git-tree-sha1 = "674ff0db93fffcd11a3573986e550d66cd4fd71f"
+git-tree-sha1 = "b36c7e110080ae48fdef61b0c31e6b17ada23b33"
 uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
-version = "2.80.5+0"
+version = "2.82.2+0"
 
 [[deps.Glob]]
 git-tree-sha1 = "97285bbd5230dd766e9ef6749b80fc617126d496"
 uuid = "c27321d9-0574-5035-807b-f59d2c89b15c"
 version = "1.3.1"
 
-[[deps.GnuTLS_jll]]
-deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Nettle_jll", "P11Kit_jll", "Zlib_jll"]
-git-tree-sha1 = "383db7d3f900f4c1f47a8a04115b053c095e48d3"
-uuid = "0951126a-58fd-58f1-b5b3-b08c7c4a876d"
-version = "3.8.4+0"
-
 [[deps.Graphics]]
 deps = ["Colors", "LinearAlgebra", "NaNMath"]
 git-tree-sha1 = "a641238db938fff9b2f60d08ed9030387daf428c"
@@ -986,15 +1008,15 @@ version = "1.1.3"
 
 [[deps.Graphite2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011"
+git-tree-sha1 = "01979f9b37367603e2848ea225918a3b3861b606"
 uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
-version = "1.3.14+0"
+version = "1.3.14+1"
 
 [[deps.GridLayoutBase]]
 deps = ["GeometryBasics", "InteractiveUtils", "Observables"]
-git-tree-sha1 = "6f93a83ca11346771a93bbde2bdad2f65b61498f"
+git-tree-sha1 = "dc6bed05c15523624909b3953686c5f5ffa10adc"
 uuid = "3955a311-db13-416c-9275-1d80ed98e5e9"
-version = "0.10.2"
+version = "0.11.1"
 
 [[deps.Grisu]]
 git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
@@ -1130,13 +1152,11 @@ deps = ["Adapt", "AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArr
 git-tree-sha1 = "88a101217d7cb38a7b481ccd50d21876e1d1b0e0"
 uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 version = "0.15.1"
+weakdeps = ["Unitful"]
 
     [deps.Interpolations.extensions]
     InterpolationsUnitfulExt = "Unitful"
 
-    [deps.Interpolations.weakdeps]
-    Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
-
 [[deps.IntervalArithmetic]]
 deps = ["CRlibm", "FastRounding", "LinearAlgebra", "Markdown", "Random", "RecipesBase", "RoundingEmulator", "SetRounding", "StaticArrays"]
 git-tree-sha1 = "5ab7744289be503d76a944784bac3f2df7b809af"
@@ -1206,9 +1226,9 @@ version = "0.9.12"
 
 [[deps.JLD2]]
 deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "PrecompileTools", "Requires", "TranscodingStreams"]
-git-tree-sha1 = "ce5737c0d4490b0e0040b5dc77fbb6a351ddf188"
+git-tree-sha1 = "f1a1c1037af2a4541ea186b26b0c0e7eeaad232b"
 uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-version = "0.5.8"
+version = "0.5.10"
 
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
@@ -1236,9 +1256,9 @@ version = "3.0.4+0"
 
 [[deps.JuliaInterpreter]]
 deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
-git-tree-sha1 = "fc8504eca188aaae4345649ca6105806bc584b70"
+git-tree-sha1 = "10da5154188682e5c0726823c2b5125957ec3778"
 uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
-version = "0.9.37"
+version = "0.9.38"
 
 [[deps.JuliaNVTXCallbacks_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -1413,9 +1433,9 @@ version = "1.17.0+1"
 
 [[deps.Libmount_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "0c4f9c4f1a50d8f35048fa0532dabbadf702f81e"
+git-tree-sha1 = "84eef7acd508ee5b3e956a2ae51b05024181dee0"
 uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
-version = "2.40.1+0"
+version = "2.40.2+0"
 
 [[deps.Libtiff_jll]]
 deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"]
@@ -1425,9 +1445,9 @@ version = "4.7.0+0"
 
 [[deps.Libuuid_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "5ee6203157c120d79034c748a2acba45b82b8807"
+git-tree-sha1 = "edbf5309f9ddf1cab25afc344b1e8150b7c832f9"
 uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
-version = "2.40.1+0"
+version = "2.40.2+0"
 
 [[deps.LightXML]]
 deps = ["Libdl", "XML2_jll"]
@@ -1482,9 +1502,9 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
 [[deps.LoweredCodeUtils]]
 deps = ["JuliaInterpreter"]
-git-tree-sha1 = "260dc274c1bc2cb839e758588c63d9c8b5e639d1"
+git-tree-sha1 = "688d6d9e098109051ae33d126fcfc88c4ce4a021"
 uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
-version = "3.0.5"
+version = "3.1.0"
 
 [[deps.Lz4_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1542,16 +1562,16 @@ uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 version = "0.5.13"
 
 [[deps.Makie]]
-deps = ["Animations", "Base64", "CRC32c", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG_jll", "FileIO", "FilePaths", "FixedPointNumbers", "Format", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageIO", "InteractiveUtils", "IntervalSets", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "MakieCore", "Markdown", "MathTeXEngine", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Scratch", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun"]
-git-tree-sha1 = "4d49c9ee830eec99d3e8de2425ff433ece7cc1bc"
+deps = ["Animations", "Base64", "CRC32c", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "Dates", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG_jll", "FileIO", "FilePaths", "FixedPointNumbers", "Format", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageBase", "ImageIO", "InteractiveUtils", "Interpolations", "IntervalSets", "InverseFunctions", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "MakieCore", "Markdown", "MathTeXEngine", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Scratch", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun", "Unitful"]
+git-tree-sha1 = "5e4e0e027642293da251bf35dac408d692ccba8b"
 uuid = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
-version = "0.20.10"
+version = "0.21.16"
 
 [[deps.MakieCore]]
-deps = ["Observables", "REPL"]
-git-tree-sha1 = "248b7a4be0f92b497f7a331aed02c1e9a878f46b"
+deps = ["ColorTypes", "GeometryBasics", "IntervalSets", "Observables"]
+git-tree-sha1 = "ae4dbe0fcf1594ed98594e5f4ee685295a2a6f74"
 uuid = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
-version = "0.7.3"
+version = "0.8.10"
 
 [[deps.ManualMemory]]
 git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd"
@@ -1569,9 +1589,9 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
 [[deps.MathTeXEngine]]
 deps = ["AbstractTrees", "Automa", "DataStructures", "FreeTypeAbstraction", "GeometryBasics", "LaTeXStrings", "REPL", "RelocatableFolders", "UnicodeFun"]
-git-tree-sha1 = "96ca8a313eb6437db5ffe946c457a401bbb8ce1d"
+git-tree-sha1 = "f45c8916e8385976e1ccd055c9874560c257ab13"
 uuid = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53"
-version = "0.5.7"
+version = "0.6.2"
 
 [[deps.MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -1580,9 +1600,9 @@ version = "2.28.2+1"
 
 [[deps.MicrosoftMPI_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "f12a29c4400ba812841c6ace3f4efbb6dbb3ba01"
+git-tree-sha1 = "bc95bf4149bf535c09602e3acdf950d9b4376227"
 uuid = "9237b28f-5490-5468-be7b-bb81f5f5e6cf"
-version = "10.1.4+2"
+version = "10.1.4+3"
 
 [[deps.Missings]]
 deps = ["DataAPI"]
@@ -1625,9 +1645,9 @@ version = "0.14.6"
 
 [[deps.NVTX]]
 deps = ["Colors", "JuliaNVTXCallbacks_jll", "Libdl", "NVTX_jll"]
-git-tree-sha1 = "53046f0483375e3ed78e49190f1154fa0a4083a1"
+git-tree-sha1 = "6a6f8bfaa91bb2e40ff562ab9f30dc827741daef"
 uuid = "5da4648a-3479-48b8-97b9-01cb529c0a1f"
-version = "0.3.4"
+version = "0.3.5"
 
 [[deps.NVTX_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -1643,9 +1663,9 @@ version = "1.0.2"
 
 [[deps.NaNStatistics]]
 deps = ["PrecompileTools", "Static", "StaticArrayInterface"]
-git-tree-sha1 = "d86d1c878cdab73ec06da5eb348967992879507f"
+git-tree-sha1 = "f2a685149d4e7d62865babd606092f8e89a6d8f8"
 uuid = "b946abbf-3ea7-4610-9019-9858bfdeaf2d"
-version = "0.6.42"
+version = "0.6.43"
 
 [[deps.NetCDF_jll]]
 deps = ["Artifacts", "Blosc_jll", "Bzip2_jll", "HDF5_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "OpenMPI_jll", "XML2_jll", "Zlib_jll", "Zstd_jll", "libzip_jll"]
@@ -1659,12 +1679,6 @@ git-tree-sha1 = "d92b107dbb887293622df7697a2223f9f8176fcd"
 uuid = "f09324ee-3d7c-5217-9330-fc30815ba969"
 version = "1.1.1"
 
-[[deps.Nettle_jll]]
-deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "eca63e3847dad608cfa6a3329b95ef674c7160b4"
-uuid = "4c82536e-c426-54e4-b420-14f461c4ed8b"
-version = "3.7.2+0"
-
 [[deps.NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 version = "1.2.0"
@@ -1725,9 +1739,9 @@ version = "0.8.1+2"
 
 [[deps.OpenMPI_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML", "Zlib_jll"]
-git-tree-sha1 = "bfce6d523861a6c562721b262c0d1aaeead2647f"
+git-tree-sha1 = "2dace87e14256edb1dd0724ab7ba831c779b96bd"
 uuid = "fe0851c0-eecd-5654-98d4-656369965a5c"
-version = "5.0.5+0"
+version = "5.0.6+0"
 
 [[deps.OpenSSL_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1748,15 +1762,9 @@ uuid = "91d4177d-7536-5919-b921-800302f37372"
 version = "1.3.3+0"
 
 [[deps.OrderedCollections]]
-git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
+git-tree-sha1 = "12f1439c4f986bb868acda6ea33ebc78e19b95ad"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.6.3"
-
-[[deps.P11Kit_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "2cd396108e178f3ae8dedbd8e938a18726ab2fbf"
-uuid = "c2071276-7c44-58a7-b746-946036e04d0a"
-version = "0.24.1+0"
+version = "1.7.0"
 
 [[deps.PCRE2_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -1783,9 +1791,9 @@ weakdeps = ["Requires", "TOML"]
 
 [[deps.Packing]]
 deps = ["GeometryBasics"]
-git-tree-sha1 = "ec3edfe723df33528e085e632414499f26650501"
+git-tree-sha1 = "bc5bf2ea3d5351edf285a06b0016788a121ce92c"
 uuid = "19eb6ba3-879d-56ad-ad62-d5c202156566"
-version = "0.5.0"
+version = "0.5.1"
 
 [[deps.PaddedViews]]
 deps = ["OffsetArrays"]
@@ -1996,9 +2004,9 @@ version = "1.3.4"
 
 [[deps.RecursiveArrayTools]]
 deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "IteratorInterfaceExtensions", "LinearAlgebra", "RecipesBase", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"]
-git-tree-sha1 = "6f4dca5fd8e97087a76b7ab8384d1c3086ace0b7"
+git-tree-sha1 = "32f824db4e5bab64e25a12b22483a30a6b813d08"
 uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
-version = "3.27.3"
+version = "3.27.4"
 
     [deps.RecursiveArrayTools.extensions]
     RecursiveArrayToolsFastBroadcastExt = "FastBroadcast"
@@ -2007,6 +2015,7 @@ version = "3.27.3"
     RecursiveArrayToolsMonteCarloMeasurementsExt = "MonteCarloMeasurements"
     RecursiveArrayToolsReverseDiffExt = ["ReverseDiff", "Zygote"]
     RecursiveArrayToolsSparseArraysExt = ["SparseArrays"]
+    RecursiveArrayToolsStructArraysExt = "StructArrays"
     RecursiveArrayToolsTrackerExt = "Tracker"
     RecursiveArrayToolsZygoteExt = "Zygote"
 
@@ -2017,6 +2026,7 @@ version = "3.27.3"
     MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
     SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+    StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
     Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
     Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
@@ -2083,9 +2093,9 @@ version = "0.1.0"
 
 [[deps.SciMLBase]]
 deps = ["ADTypes", "Accessors", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "Expronicon", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface"]
-git-tree-sha1 = "cacc7bc54bab8749b1fc1032c4911fe80cffb959"
+git-tree-sha1 = "6f156f48d3603e8023e59917547b141218f2656d"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-version = "2.61.0"
+version = "2.64.0"
 
     [deps.SciMLBase.extensions]
     SciMLBaseChainRulesCoreExt = "ChainRulesCore"
@@ -2119,9 +2129,9 @@ weakdeps = ["SparseArrays", "StaticArraysCore"]
 
 [[deps.SciMLStructures]]
 deps = ["ArrayInterface"]
-git-tree-sha1 = "25514a6f200219cd1073e4ff23a6324e4a7efe64"
+git-tree-sha1 = "0444a37a25fab98adbd90baa806ee492a3af133a"
 uuid = "53ae85a6-f571-4167-b2af-e1d143709226"
-version = "1.5.0"
+version = "1.6.1"
 
 [[deps.Scratch]]
 deps = ["Dates"]
@@ -2220,10 +2230,10 @@ uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15"
 version = "0.1.1"
 
 [[deps.Static]]
-deps = ["IfElse"]
-git-tree-sha1 = "d2fdac9ff3906e27f7a618d47b676941baa6c80c"
+deps = ["CommonWorldInvalidations", "IfElse", "PrecompileTools"]
+git-tree-sha1 = "87d51a3ee9a4b0d2fe054bdd3fc2436258db2603"
 uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
-version = "0.8.10"
+version = "1.1.1"
 
 [[deps.StaticArrayInterface]]
 deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Static"]
@@ -2423,11 +2433,6 @@ git-tree-sha1 = "0c45878dcfdcfa8480052b6ab162cdd138781742"
 uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
 version = "0.11.3"
 
-[[deps.Tricks]]
-git-tree-sha1 = "7822b97e99a1672bfb1b49b668a6d46d58d8cbcb"
-uuid = "410a4b4d-49e4-4fbc-ab6d-cb71b17b3775"
-version = "0.1.9"
-
 [[deps.TriplotBase]]
 git-tree-sha1 = "4d4ed7f294cda19382ff7de4c137d24d16adc89b"
 uuid = "981d1d27-644d-49a2-9326-4793e63143c3"
@@ -2457,6 +2462,17 @@ git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf"
 uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1"
 version = "0.4.1"
 
+[[deps.Unitful]]
+deps = ["Dates", "LinearAlgebra", "Random"]
+git-tree-sha1 = "d95fe458f26209c66a187b1114df96fd70839efd"
+uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
+version = "1.21.0"
+weakdeps = ["ConstructionBase", "InverseFunctions"]
+
+    [deps.Unitful.extensions]
+    ConstructionBaseUnitfulExt = "ConstructionBase"
+    InverseFunctionsUnitfulExt = "InverseFunctions"
+
 [[deps.Unrolled]]
 deps = ["MacroTools"]
 git-tree-sha1 = "6cc9d682755680e0f0be87c56392b7651efc2c7b"
@@ -2520,9 +2536,9 @@ version = "2.13.5+0"
 
 [[deps.XSLT_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "XML2_jll", "Zlib_jll"]
-git-tree-sha1 = "a54ee957f4c86b526460a720dbc882fa5edcbefc"
+git-tree-sha1 = "7d1671acbe47ac88e981868a078bd6b4e27c5191"
 uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
-version = "1.1.41+0"
+version = "1.1.42+0"
 
 [[deps.XZ_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -2538,15 +2554,15 @@ version = "1.8.6+0"
 
 [[deps.Xorg_libXau_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "6035850dcc70518ca32f012e46015b9beeda49d8"
+git-tree-sha1 = "2b0e27d52ec9d8d483e2ca0b72b3cb1a8df5c27a"
 uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
-version = "1.0.11+0"
+version = "1.0.11+1"
 
 [[deps.Xorg_libXdmcp_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "34d526d318358a859d7de23da945578e8e8727b7"
+git-tree-sha1 = "02054ee01980c90297412e4c809c8694d7323af3"
 uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
-version = "1.1.4+0"
+version = "1.1.4+1"
 
 [[deps.Xorg_libXext_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"]
@@ -2562,9 +2578,9 @@ version = "0.9.11+0"
 
 [[deps.Xorg_libpthread_stubs_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "8fdda4c692503d44d04a0603d9ac0982054635f9"
+git-tree-sha1 = "fee57a273563e273f0f53275101cd41a8153517a"
 uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
-version = "0.1.1+0"
+version = "0.1.1+1"
 
 [[deps.Xorg_libxcb_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
@@ -2574,9 +2590,9 @@ version = "1.17.0+0"
 
 [[deps.Xorg_xtrans_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "e92a1a012a10506618f10b7047e478403a046c77"
+git-tree-sha1 = "b9ead2d2bdb27330545eb14234a2e300da61232e"
 uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
-version = "1.5.0+0"
+version = "1.5.0+1"
 
 [[deps.YAML]]
 deps = ["Base64", "Dates", "Printf", "StringEncodings"]
@@ -2595,6 +2611,12 @@ git-tree-sha1 = "555d1076590a6cc2fdee2ef1469451f872d8b41b"
 uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
 version = "1.5.6+1"
 
+[[deps.demumble_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "6498e3581023f8e530f34760d18f75a69e3a4ea8"
+uuid = "1e29f10c-031c-5a83-9565-69cddfc27673"
+version = "1.3.0+0"
+
 [[deps.isoband_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "51b5eeb3f98367157a7a12a1fb0aa5328946c03c"
@@ -2655,10 +2677,10 @@ uuid = "c5f90fcd-3b7e-5836-afba-fc50a0988cb2"
 version = "1.4.0+0"
 
 [[deps.libzip_jll]]
-deps = ["Artifacts", "Bzip2_jll", "GnuTLS_jll", "JLLWrappers", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"]
-git-tree-sha1 = "668ac0297e6bd8f4d53dfdcd3ace71f2e00f4a35"
+deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "OpenSSL_jll", "XZ_jll", "Zlib_jll", "Zstd_jll"]
+git-tree-sha1 = "e797fa066eba69f4c0585ffbd81bc780b5118ce2"
 uuid = "337d8026-41b4-5cde-a456-74a10e5b31d1"
-version = "1.11.1+0"
+version = "1.11.2+0"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]
diff --git a/examples/Project.toml b/examples/Project.toml
index 9f60d9354c..af7c481e38 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -61,7 +61,7 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
 
 [compat]
-CairoMakie = "0.10, 0.11"
+CairoMakie = "0.10, 0.11, 0.12"
 ClimaAnalysis = "0.5.7"
 ClimaCoreMakie = "0.3, 0.4"
 ClimaCorePlots = "0.2"
diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl
index df2c26072e..fe8940b844 100644
--- a/post_processing/ci_plots.jl
+++ b/post_processing/ci_plots.jl
@@ -42,7 +42,7 @@ import CairoMakie
 import CairoMakie.Makie
 import ClimaAnalysis
 import ClimaAnalysis: Visualize as viz
-import ClimaAnalysis: SimDir, slice, read_var, average_xy
+import ClimaAnalysis: SimDir, slice, read_var, average_xy, window, average_time
 import ClimaAnalysis.Utils: kwargs as ca_kwargs
 
 import ClimaCoreSpectra: power_spectrum_2d
@@ -577,6 +577,7 @@ end
 
 function make_plots(
     ::Val{:box_density_current_test},
+    ::Val{:box_rising_thermal_test},
     output_paths::Vector{<:AbstractString},
 )
     simdirs = SimDir.(output_paths)
@@ -1065,6 +1066,100 @@ function make_plots(::Aquaplanet1MPlots, output_paths::Vector{<:AbstractString})
     )
 end
 
+LESBoxPlots = Union{Val{:les_box}}
+
+"""
+    plot_les_vert_profile!(grid_loc, var_group)
+
+Helper function for `make_plots_generic`. Takes a list of variables and plots
+them on the same axis.
+"""
+function plot_les_vert_profile!(grid_loc, var_group)
+    z = var_group[1].dims["z"]
+    units = var_group[1].attributes["units"]
+    ax = CairoMakie.Axis(
+        grid_loc[1, 1],
+        ylabel = "z [$(var_group[1].dim_attributes["z"]["units"])]",
+        xlabel = "$(short_name(var_group[1])) [$units]",
+        title = parse_var_attributes(var_group[1]),
+    )
+
+    for var in var_group
+        CairoMakie.lines!(ax, var.data, z, label = short_name(var))
+    end
+    length(var_group) > 1 && Makie.axislegend(ax)
+end
+
+function make_plots(
+    sim_type::Union{LESBoxPlots},
+    output_paths::Vector{<:AbstractString},
+)
+    simdirs = SimDir.(output_paths)
+
+    reduction = "inst"
+    short_names = [
+        "wa",
+        "ua",
+        "va",
+        "ta",
+        "thetaa",
+        "ha",
+        "hus",
+        "hur",
+        "cl",
+        "clw",
+        "cli",
+    ]
+    short_names = short_names ∩ collect(keys(simdirs[1].vars))
+
+    # Window average from instantaneous snapshots?
+    function horizontal_average(var)
+        return average_xy(var)
+    end
+    function windowed_reduction(var)
+        hours = 3600.0
+        window_end = last(var.dims["time"])
+        window_start = window_end - 2hours
+        var_window = ClimaAnalysis.window(
+            var,
+            "time";
+            left = window_start,
+            right = window_end,
+        )
+        var_reduced = horizontal_average(average_time(var_window))
+        return var_reduced
+    end
+
+    var_groups_xyt_reduced =
+        map_comparison(simdirs, short_names) do simdir, short_name
+            return [get(simdir; short_name, reduction) |> windowed_reduction]
+        end
+
+    var_groups_xy_reduced =
+        map_comparison(simdirs, short_names) do simdir, short_name
+            return [get(simdir; short_name, reduction) |> horizontal_average]
+        end
+
+    tmp_file = make_plots_generic(
+        output_paths,
+        var_groups_xyt_reduced,
+        output_name = "tmp";
+        plot_fn = plot_les_vert_profile!,
+        MAX_NUM_COLS = 2,
+        MAX_NUM_ROWS = 4,
+    )
+
+    make_plots_generic(
+        output_paths,
+        vcat(var_groups_xy_reduced...),
+        plot_fn = plot_parsed_attribute_title!,
+        summary_files = [tmp_file],
+        MAX_NUM_COLS = 2,
+        MAX_NUM_ROWS = 4,
+    )
+end
+
+
 EDMFBoxPlots = Union{
     Val{:diagnostic_edmfx_gabls_box},
     Val{:diagnostic_edmfx_bomex_box},
@@ -1097,7 +1192,6 @@ EDMFBoxPlotsWithPrecip = Union{
     Val{:diagnostic_edmfx_trmm_stretched_box},
 }
 
-
 """
     plot_edmf_vert_profile!(grid_loc, var_group)
 
diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl
index 13a447ba66..416d29f62b 100644
--- a/src/ClimaAtmos.jl
+++ b/src/ClimaAtmos.jl
@@ -111,6 +111,13 @@ include(
 )
 include(joinpath("parameterized_tendencies", "sponge", "rayleigh_sponge.jl"))
 include(joinpath("parameterized_tendencies", "sponge", "viscous_sponge.jl"))
+include(
+    joinpath(
+        "parameterized_tendencies",
+        "les_sgs_models",
+        "smagorinsky_lilly.jl",
+    ),
+)
 include(joinpath("prognostic_equations", "advection.jl"))
 
 include(joinpath("cache", "temporary_quantities.jl"))
diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
index 5d41912fd4..6acf4fa2ad 100644
--- a/src/cache/precomputed_quantities.jl
+++ b/src/cache/precomputed_quantities.jl
@@ -156,6 +156,19 @@ function precomputed_quantities(Y, atmos)
             ᶜqᵣ = similar(Y.c, FT),
             ᶜqₛ = similar(Y.c, FT),
         ) : (;)
+    smagorinsky_lilly_quantities =
+        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),
+            )
+        else
+            (;)
+        end
+
     return (;
         gs_quantities...,
         sgs_quantities...,
@@ -164,6 +177,7 @@ function precomputed_quantities(Y, atmos)
         vert_diff_quantities...,
         precipitation_quantities...,
         cloud_diagnostics_tuple,
+        smagorinsky_lilly_quantities...,
     )
 end
 
@@ -633,6 +647,10 @@ NVTX.@annotate function set_precomputed_quantities!(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
+
     return nothing
 end
 
diff --git a/src/cache/temporary_quantities.jl b/src/cache/temporary_quantities.jl
index 4482e86147..cee88bc51d 100644
--- a/src/cache/temporary_quantities.jl
+++ b/src/cache/temporary_quantities.jl
@@ -10,6 +10,7 @@ function temporary_quantities(Y, atmos)
 
     FT = Spaces.undertype(center_space)
     CTh = CTh_vector_type(Y.c)
+    uvw_vec = UVW(FT(0), FT(0), FT(0))
     return (;
         ᶠtemp_scalar = Fields.Field(FT, face_space), # ᶠp, ᶠρK_E
         ᶜtemp_scalar = Fields.Field(FT, center_space), # ᶜ1
@@ -45,14 +46,12 @@ function temporary_quantities(Y, atmos)
             face_space,
         ), # ᶠω¹²ʲs
         ᶠtemp_C123 = Fields.Field(C123{FT}, face_space), # χ₁₂₃
-        ᶜtemp_UVWxUVW = Fields.Field(
-            typeof(UVW(FT(0), FT(0), FT(0)) * UVW(FT(0), FT(0), FT(0))'),
-            center_space,
-        ), # ᶜstrain_rate
-        ᶠtemp_UVWxUVW = Fields.Field(
-            typeof(UVW(FT(0), FT(0), FT(0)) * UVW(FT(0), FT(0), FT(0))'),
-            face_space,
-        ), # ᶠstrain_rate
+        ᶜtemp_UVW = Fields.Field(typeof(uvw_vec), center_space), # UVW(ᶜu)
+        ᶠtemp_UVW = Fields.Field(typeof(uvw_vec), face_space), # UVW(ᶠu³)
+        ᶜtemp_UVWxUVW = Fields.Field(typeof(uvw_vec * uvw_vec'), center_space), # ᶜstrain_rate
+        ᶠtemp_UVWxUVW = Fields.Field(typeof(uvw_vec * uvw_vec'), face_space), # ᶠstrain_rate
+        ᶜtemp_strain = Fields.Field(typeof(uvw_vec * uvw_vec'), center_space), # ᶜstrain_rate
+        ᶠtemp_strain = Fields.Field(typeof(uvw_vec * uvw_vec'), face_space), # ᶠstrain_rate
         # TODO: Remove this hack
         sfc_temp_C3 = Fields.Field(C3{FT}, Spaces.level(face_space, half)), # ρ_flux_χ
         # Implicit solver cache:
diff --git a/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl b/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl
new file mode 100644
index 0000000000..fb63815d12
--- /dev/null
+++ b/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl
@@ -0,0 +1,174 @@
+#####
+##### Smagorinsky Lilly Diffusion
+#####
+
+import ClimaCore.Fields as Fields
+import ClimaCore.Operators as Operators
+import ClimaCore: Geometry
+
+"""
+    set_smagorinsky_lilly_precomputed_quantities!(Y, p)
+
+Compute the Smagorinsky-Lilly diffusivity tensors, `ᶜτ_smag`, `ᶠτ_smag`, `ᶜD_smag`, and `ᶠD_smag`. 
+Store in the precomputed quantities `p.precomputed`.
+
+The subgrid-scale momentum flux tensor is defined by `τ = -2 νₜ ∘ S`, where `νₜ` is the Smagorinsky-Lilly eddy viscosity 
+and `S` is the strain rate tensor. 
+
+The turbulent diffusivity is defined as `D = νₜ / Pr_t`, where `Pr_t` is the turbulent Prandtl number for neutral 
+stratification. 
+
+These quantities are computed for both cell centers and faces, with prefixes `ᶜ` and `ᶠ`, respectively.
+
+
+# Arguments
+- `Y`: The model state.
+- `p`: The model parameters, e.g. `AtmosCache`.
+"""
+function set_smagorinsky_lilly_precomputed_quantities!(Y, p)
+
+    (; atmos, precomputed, scratch, params) = p
+    (; Cs, Pr_t) = atmos.smagorinsky_lilly
+    (; ᶜu, ᶠu³, ᶜts, ᶜτ_smag, ᶠτ_smag, ᶜD_smag, ᶠD_smag) = precomputed
+    FT = eltype(Y)
+    grav = CAP.grav(params)
+    thermo_params = CAP.thermodynamics_params(params)
+    (; ᶜtemp_UVWxUVW, ᶠtemp_UVWxUVW, ᶜtemp_strain, ᶠtemp_strain) = scratch
+    (; ᶜtemp_scalar, ᶜtemp_scalar_2, ᶠtemp_scalar, ᶜtemp_UVW, ᶠtemp_UVW) =
+        scratch
+
+    ∇ᵥuvw_boundary = Geometry.outer(Geometry.WVector(0), UVW(0, 0, 0))
+    ᶠgradᵥ_uvw = Operators.GradientC2F(
+        bottom = Operators.SetGradient(∇ᵥuvw_boundary),
+        top = Operators.SetGradient(∇ᵥuvw_boundary),
+    )
+    axis_uvw = (Geometry.UVWAxis(),)
+
+    # Compute UVW velocities
+    ᶜu_uvw = @. ᶜtemp_UVW = UVW(ᶜu)
+    ᶠu_uvw = @. ᶠtemp_UVW = UVW(ᶠinterp(Y.c.uₕ)) + UVW(ᶠu³)
+
+    # Gradients
+    ## cell centers
+    ∇ᶜu_uvw = @. ᶜtemp_UVWxUVW = Geometry.project(axis_uvw, ᶜgradᵥ(ᶠu_uvw))  # vertical component
+    @. ∇ᶜu_uvw += Geometry.project(axis_uvw, gradₕ(ᶜu_uvw))  # horizontal component
+    ## cell faces
+    ∇ᶠu_uvw = @. ᶠtemp_UVWxUVW = Geometry.project(axis_uvw, ᶠgradᵥ_uvw(ᶜu_uvw))  # vertical component
+    @. ∇ᶠu_uvw += Geometry.project(axis_uvw, gradₕ(ᶠu_uvw))  # horizontal component
+
+    # Strain rate tensor
+    ᶜS = @. ᶜtemp_strain = (∇ᶜu_uvw + adjoint(∇ᶜu_uvw)) / 2
+    ᶠS = @. ᶠtemp_strain = (∇ᶠu_uvw + adjoint(∇ᶠu_uvw)) / 2
+
+    # Stratification correction
+    ᶜθ_v = @. ᶜtemp_scalar = TD.virtual_pottemp(thermo_params, ᶜts)
+    ᶜ∇ᵥθ = @. ᶜtemp_scalar_2 =
+        Geometry.WVector(ᶜgradᵥ(ᶠinterp(ᶜθ_v))).components.data.:1
+    ᶜN² = @. ᶜtemp_scalar = grav / ᶜθ_v * ᶜ∇ᵥθ
+    ᶜS_norm = @. ᶜtemp_scalar_2 = √(2 * CA.norm_sqr(ᶜS))
+
+    ᶜRi = @. ᶜtemp_scalar = ᶜN² / (ᶜS_norm^2 + eps(FT))  # Ri = N² / |S|²
+    ᶜfb = @. ᶜtemp_scalar = ifelse(ᶜRi ≤ 0, 1, max(0, 1 - ᶜRi / Pr_t)^(1 / 4))
+
+    # filter scale
+    h_space = Spaces.horizontal_space(axes(Y.c))
+    Δ_xy = Spaces.node_horizontal_length_scale(h_space)^2 # Δ_x * Δ_y
+    ᶜΔ_z = Fields.Δz_field(Y.c)
+    ᶜΔ = @. ᶜtemp_scalar = ∛(Δ_xy * ᶜΔ_z) * ᶜfb
+
+    # Smagorinsky-Lilly eddy viscosity
+    ᶜνₜ = @. ᶜtemp_scalar = Cs^2 * ᶜΔ^2 * ᶜS_norm
+    ᶠνₜ = @. ᶠtemp_scalar = ᶠinterp(ᶜνₜ)
+
+    # Subgrid-scale momentum flux tensor, `τ = -2 νₜ ∘ S`
+    @. ᶜτ_smag = -2 * ᶜνₜ * ᶜS
+    @. ᶠτ_smag = -2 * ᶠνₜ * ᶠS
+
+    # Turbulent diffusivity
+    @. ᶜD_smag = ᶜνₜ / Pr_t
+    @. ᶠD_smag = ᶠνₜ / Pr_t
+
+    nothing
+end
+
+horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::Nothing) = nothing
+vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::Nothing) = nothing
+
+function horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::SmagorinskyLilly)
+    (; ᶜτ_smag, ᶠτ_smag, ᶜD_smag, ᶜspecific, ᶜh_tot) = p.precomputed
+
+    ## Momentum tendencies
+    ᶠρ = @. p.scratch.ᶠtemp_scalar = ᶠinterp(Y.c.ρ)
+    @. Yₜ.c.uₕ -= C12(wdivₕ(Y.c.ρ * ᶜτ_smag) / Y.c.ρ)
+    @. Yₜ.f.u₃ -= C3(wdivₕ(ᶠρ * ᶠτ_smag) / ᶠρ)
+
+    ## Total energy tendency
+    @. Yₜ.c.ρe_tot += wdivₕ(Y.c.ρ * ᶜD_smag * gradₕ(ᶜh_tot))
+
+    ## Tracer diffusion and associated mass changes
+    for (ᶜρχₜ, ᶜχ, χ_name) in CA.matching_subfields(Yₜ.c, ᶜspecific)
+        χ_name == :e_tot && continue
+        ᶜρχₜ_diffusion =
+            @. p.scratch.ᶜtemp_scalar = wdivₕ(Y.c.ρ * ᶜD_smag * gradₕ(ᶜχ))
+        @. ᶜρχₜ += ᶜρχₜ_diffusion
+        # Rain and snow does not affect the mass
+        if χ_name ∉ (:q_rai, :q_sno)
+            @. Yₜ.c.ρ += ᶜρχₜ_diffusion
+        end
+    end
+
+end
+
+function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::SmagorinskyLilly)
+    FT = eltype(Y)
+    (; sfc_temp_C3, ᶠtemp_scalar, ᶜtemp_scalar) = p.scratch
+    (; ᶜτ_smag, ᶠτ_smag, ᶠD_smag, ᶜspecific, ᶜh_tot, sfc_conditions) =
+        p.precomputed
+    (; ρ_flux_uₕ, ρ_flux_h_tot) = sfc_conditions
+
+    # Define operators
+    ᶠgradᵥ = Operators.GradientC2F() # apply BCs to ᶜdivᵥ, which wraps ᶠgradᵥ
+    ᶜdivᵥ_uₕ = Operators.DivergenceF2C(
+        top = Operators.SetValue(C3(FT(0)) ⊗ C12(FT(0), FT(0))),
+        bottom = Operators.SetValue(ρ_flux_uₕ),
+    )
+    ᶠdivᵥ = Operators.DivergenceC2F(
+        bottom = Operators.SetDivergence(FT(0)),
+        top = Operators.SetDivergence(FT(0)),
+    )
+    top = Operators.SetValue(C3(FT(0)))
+    ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(;
+        top,
+        bottom = Operators.SetValue(ρ_flux_h_tot),
+    )
+
+    # Apply to tendencies
+    ## Horizontal momentum tendency
+    ᶠρ = @. ᶠtemp_scalar = ᶠinterp(Y.c.ρ)
+    @. Yₜ.c.uₕ -= C12(ᶜdivᵥ(ᶠρ * ᶠτ_smag) / Y.c.ρ)
+    ## Apply boundary condition for momentum flux
+    @. Yₜ.c.uₕ -= ᶜdivᵥ_uₕ(-(FT(0) * ᶠgradᵥ(Y.c.uₕ))) / Y.c.ρ
+    ## Vertical momentum tendency
+    @. Yₜ.f.u₃ -= C3(ᶠdivᵥ(Y.c.ρ * ᶜτ_smag) / ᶠρ)
+
+    ## Total energy tendency
+    @. Yₜ.c.ρe_tot -= ᶜdivᵥ_ρe_tot(-(ᶠρ * ᶠD_smag * ᶠgradᵥ(ᶜh_tot)))
+
+    ## Tracer diffusion and associated mass changes
+    sfc_zero = @. sfc_temp_C3 = C3(FT(0))
+    for (ᶜρχₜ, ᶜχ, χ_name) in CA.matching_subfields(Yₜ.c, ᶜspecific)
+        χ_name == :e_tot && continue
+
+        bottom = Operators.SetValue(
+            χ_name == :q_tot ? sfc_conditions.ρ_flux_q_tot : sfc_zero,
+        )
+        ᶜdivᵥ_ρχ = Operators.DivergenceF2C(; top, bottom)
+
+        ᶜ∇ᵥρD∇χₜ = @. ᶜtemp_scalar = ᶜdivᵥ_ρχ(-(ᶠρ * ᶠD_smag * ᶠgradᵥ(ᶜχ)))
+        @. ᶜρχₜ -= ᶜ∇ᵥρD∇χₜ
+        # Rain and snow does not affect the mass
+        if χ_name ∉ (:q_rai, :q_sno)
+            @. Yₜ.c.ρ -= ᶜ∇ᵥρD∇χₜ
+        end
+    end
+end
diff --git a/src/parameterized_tendencies/radiation/radiation.jl b/src/parameterized_tendencies/radiation/radiation.jl
index 019f593077..8804243bc0 100644
--- a/src/parameterized_tendencies/radiation/radiation.jl
+++ b/src/parameterized_tendencies/radiation/radiation.jl
@@ -365,9 +365,11 @@ function radiation_tendency!(Yₜ, Y, p, t, radiation_mode::RadiationDYCOMS)
 
     # Find the values of (z, ρ, q_tot) at the q_tot = 0.008 isoline, i.e., at
     # the level whose value of q_tot is closest to 0.008.
+    q_tot_isoline = FT(0.008)
     Operators.column_reduce!(
         (nt1, nt2) ->
-            abs(nt1.q_tot - FT(0.008)) < abs(nt2.q_tot - FT(0.008)) ? nt1 : nt2,
+            abs(nt1.q_tot - q_tot_isoline) < abs(nt2.q_tot - q_tot_isoline) ?
+            nt1 : nt2,
         isoline_z_ρ_q,
         Base.broadcasted(NT ∘ tuple, ᶜz, Y.c.ρ, ᶜspecific.q_tot),
     )
diff --git a/src/prognostic_equations/remaining_tendency.jl b/src/prognostic_equations/remaining_tendency.jl
index d56ece72f8..a31344f1d9 100644
--- a/src/prognostic_equations/remaining_tendency.jl
+++ b/src/prognostic_equations/remaining_tendency.jl
@@ -86,6 +86,10 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
     # NOTE: All ρa tendencies should be applied before calling this function
     pressure_work_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
 
+    sl = p.atmos.smagorinsky_lilly
+    horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl)
+    vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl)
+
     # NOTE: This will zero out all momentum tendencies in the edmfx advection test
     # please DO NOT add additional velocity tendencies after this function
     zero_velocity_tendency!(Yₜ, Y, p, t)
diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl
index e2402bc76e..9b84033f32 100644
--- a/src/solver/model_getters.jl
+++ b/src/solver/model_getters.jl
@@ -156,6 +156,18 @@ function get_viscous_sponge_model(parsed_args, params, ::Type{FT}) where {FT}
     end
 end
 
+function get_smagorinsky_lilly_model(parsed_args, params, ::Type{FT}) where {FT}
+    is_model_active = parsed_args["smagorinsky_lilly"]
+    Cs = parsed_args["c_smag"]
+    Pr_t = parsed_args["prandtl_turbulent_neutral"]  # Turbulent Prandtl number for neutral stratification
+    @assert is_model_active in (true, false)
+    return if is_model_active == true
+        SmagorinskyLilly{FT}(; Cs, Pr_t)
+    else
+        nothing
+    end
+end
+
 function get_rayleigh_sponge_model(parsed_args, params, ::Type{FT}) where {FT}
     rs_name = parsed_args["rayleigh_sponge"]
     return if rs_name in ("false", false)
diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl
index 498133a29e..8513af2734 100644
--- a/src/solver/type_getters.jl
+++ b/src/solver/type_getters.jl
@@ -85,6 +85,11 @@ function get_atmos(config::AtmosConfig, params)
         diff_mode = implicit_diffusion ? Implicit() : Explicit(),
         sgs_adv_mode = implicit_sgs_advection ? Implicit() : Explicit(),
         viscous_sponge = get_viscous_sponge_model(parsed_args, params, FT),
+        smagorinsky_lilly = get_smagorinsky_lilly_model(
+            parsed_args,
+            params,
+            FT,
+        ),
         rayleigh_sponge = get_rayleigh_sponge_model(parsed_args, params, FT),
         sfc_temperature = get_sfc_temperature_form(parsed_args),
         insolation = get_insolation_form(parsed_args),
diff --git a/src/solver/types.jl b/src/solver/types.jl
index b566d30419..9693216804 100644
--- a/src/solver/types.jl
+++ b/src/solver/types.jl
@@ -126,6 +126,12 @@ Base.@kwdef struct ViscousSponge{FT} <: AbstractSponge
     κ₂::FT
 end
 
+abstract type AbstractEddyViscosityModel end
+Base.@kwdef struct SmagorinskyLilly{FT} <: AbstractEddyViscosityModel
+    Cs::FT = 0.2
+    Pr_t::FT = 1 / 3
+end
+
 Base.@kwdef struct RayleighSponge{FT} <: AbstractSponge
     zd::FT
     α_uₕ::FT
@@ -436,6 +442,7 @@ Base.@kwdef struct AtmosModel{
     DM,
     SAM,
     VS,
+    SL,
     RS,
     ST,
     IN,
@@ -469,6 +476,7 @@ Base.@kwdef struct AtmosModel{
     diff_mode::DM = nothing
     sgs_adv_mode::SAM = nothing
     viscous_sponge::VS = nothing
+    smagorinsky_lilly::SL = nothing
     rayleigh_sponge::RS = nothing
     sfc_temperature::ST = nothing
     insolation::IN = nothing
diff --git a/test/parameterized_tendencies/gravity_wave/gw_plotutils.jl b/test/parameterized_tendencies/gravity_wave/gw_plotutils.jl
index b7a9b40c62..8ef787dae9 100644
--- a/test/parameterized_tendencies/gravity_wave/gw_plotutils.jl
+++ b/test/parameterized_tendencies/gravity_wave/gw_plotutils.jl
@@ -40,14 +40,7 @@ function create_plot!(
     if Z == nothing
         generic_axis = fig[p_loc[1], p_loc[2]] = GridLayout()
         axis = Axis(generic_axis[1, 1]; title, xlabel, ylabel, xscale, yscale)
-        CairoMakie.lines!(
-            X,
-            Y;
-            title,
-            linewidth,
-            label = label[1],
-            linestyle = :solid,
-        )
+        CairoMakie.lines!(X, Y; linewidth, label = label[1], linestyle = :solid)
     else
         generic_axis = fig[p_loc[1], p_loc[2]] = GridLayout() # Generic Axis Layout
         Axis(generic_axis[1, 1]; title, xlabel, ylabel, yscale, yreversed) # Plot specific attributes
diff --git a/test/parameterized_tendencies/les_closures/smagorinsky_lilly.jl b/test/parameterized_tendencies/les_closures/smagorinsky_lilly.jl
new file mode 100644
index 0000000000..434267f760
--- /dev/null
+++ b/test/parameterized_tendencies/les_closures/smagorinsky_lilly.jl
@@ -0,0 +1,69 @@
+
+using Test
+import ClimaAtmos as CA
+import ClimaAtmos.Parameters as CAP
+import ClimaCore as CC
+import Thermodynamics as TD
+
+### Common Objects ###
+@testset begin
+    "Smagorinsky Lilly function"
+    ### Boilerplate default integrator objects
+    config = CA.AtmosConfig(
+        Dict(
+            "initial_condition" => "DryDensityCurrentProfile",
+            "moist" => "dry",
+            "precip_model" => "0M",
+            "config" => "box",
+            "x_max" => π,
+            "y_max" => π,
+            "z_max" => 1.0,
+            "z_stretch" => false,
+            "output_default_diagnostics" => false,
+        ),
+    )
+    parsed_args = config.parsed_args
+    simulation = CA.get_simulation(config)
+    (; integrator) = simulation
+    Y = integrator.u
+    Yₜ = similar(Y)
+    p = integrator.p
+    params = p.params
+    cm_params = CAP.microphysics_params(params)
+    thermo_params = CAP.thermodynamics_params(params)
+
+    FT = eltype(Y)
+    ᶜYₜ = Y .* FT(0)
+    c_xyz = CC.Fields.coordinate_field(Y.c)
+    f_xyz = CC.Fields.coordinate_field(Y.f)
+    x = c_xyz.x
+    y = c_xyz.y
+    z = f_xyz.z
+
+    u = @. sin(x) * cos(y)
+    v = @. sin(y) * cos(x)
+    w = @. Geometry.WVector(z ./ maximum(z))
+
+    u_x = @. cos(x) * cos(y)
+    v_x = @. -sin(y) * sin(x)
+    u_y = @. -sin(x) * sin(y)
+    v_y = @. cos(y) * cos(x)
+    w_x = zeros(axes(Y.c))
+    w_y = zeros(axes(Y.c))
+
+    vel = @. CC.Geometry.UVVector(u, v)
+
+    Y.c.ρ .= FT(1)
+    Y.c.uₕ .= Geometry.Covariant12Vector(vel)
+    Y.f.u₃ .= Geometry.Covariant3Vector(w)
+
+    horizontal_smagorinsky_lilly_tendency(
+        Yₜ,
+        Y,
+        p,
+        t,
+        SmagorinskyLilly(FT(0.2)),
+    )
+
+    ### Component test begins here
+end
diff --git a/toml/isdac_box.toml b/toml/les_isdac.toml
similarity index 100%
rename from toml/isdac_box.toml
rename to toml/les_isdac.toml