From 37c98006c1c022213063c6c8b08d0bffcb1a69c2 Mon Sep 17 00:00:00 2001
From: Haakon Ludvig Langeland Ervik
 <45243236+haakon-e@users.noreply.github.com>
Date: Wed, 4 Dec 2024 15:01:57 -0800
Subject: [PATCH] Use ClimaParams c_smag and Pr_t for Smagorinsky LES closure

---
 config/default_configs/default_config.yml              |  6 ------
 config/model_configs/box_density_current_test.yml      |  2 --
 config/model_configs/les_isdac_box.yml                 |  2 --
 .../les_sgs_models/smagorinsky_lilly.jl                |  5 +++--
 src/solver/model_getters.jl                            | 10 ++--------
 src/solver/type_getters.jl                             |  6 +-----
 src/solver/types.jl                                    |  5 +----
 .../les_closures/smagorinsky_lilly.jl                  |  8 +-------
 toml/box_density_current_test.toml                     |  7 +++++++
 toml/les_isdac.toml                                    |  4 ++++
 10 files changed, 19 insertions(+), 36 deletions(-)
 create mode 100644 toml/box_density_current_test.toml

diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml
index 53b82d41ec..42b5b8f254 100644
--- a/config/default_configs/default_config.yml
+++ b/config/default_configs/default_config.yml
@@ -148,12 +148,6 @@ hyperdiff:
 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 312b536640..8d8f55e40c 100644
--- a/config/model_configs/box_density_current_test.yml
+++ b/config/model_configs/box_density_current_test.yml
@@ -2,7 +2,6 @@ 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
@@ -14,7 +13,6 @@ z_elem: 45
 z_stretch: false
 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:
diff --git a/config/model_configs/les_isdac_box.yml b/config/model_configs/les_isdac_box.yml
index 5a6f4c9a59..2164302c88 100644
--- a/config/model_configs/les_isdac_box.yml
+++ b/config/model_configs/les_isdac_box.yml
@@ -15,8 +15,6 @@ approximate_linear_solve_iters: 2
 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
diff --git a/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl b/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl
index fb63815d12..0e3130fe0a 100644
--- a/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl
+++ b/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl
@@ -28,7 +28,8 @@ These quantities are computed for both cell centers and faces, with prefixes `
 function set_smagorinsky_lilly_precomputed_quantities!(Y, p)
 
     (; atmos, precomputed, scratch, params) = p
-    (; Cs, Pr_t) = atmos.smagorinsky_lilly
+    c_smag = CAP.c_smag(params)
+    Pr_t = CAP.Prandtl_number_0(CAP.turbconv_params(params))
     (; ᶜu, ᶠu³, ᶜts, ᶜτ_smag, ᶠτ_smag, ᶜD_smag, ᶠD_smag) = precomputed
     FT = eltype(Y)
     grav = CAP.grav(params)
@@ -77,7 +78,7 @@ function set_smagorinsky_lilly_precomputed_quantities!(Y, p)
     ᶜΔ = @. ᶜtemp_scalar = ∛(Δ_xy * ᶜΔ_z) * ᶜfb
 
     # Smagorinsky-Lilly eddy viscosity
-    ᶜνₜ = @. ᶜtemp_scalar = Cs^2 * ᶜΔ^2 * ᶜS_norm
+    ᶜνₜ = @. ᶜtemp_scalar = c_smag^2 * ᶜΔ^2 * ᶜS_norm
     ᶠνₜ = @. ᶠtemp_scalar = ᶠinterp(ᶜνₜ)
 
     # Subgrid-scale momentum flux tensor, `τ = -2 νₜ ∘ S`
diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl
index 9b84033f32..bf199a0a84 100644
--- a/src/solver/model_getters.jl
+++ b/src/solver/model_getters.jl
@@ -156,16 +156,10 @@ 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}
+function get_smagorinsky_lilly_model(parsed_args)
     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
+    return is_model_active ? SmagorinskyLilly() : nothing
 end
 
 function get_rayleigh_sponge_model(parsed_args, params, ::Type{FT}) where {FT}
diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl
index 8513af2734..5610d9fc05 100644
--- a/src/solver/type_getters.jl
+++ b/src/solver/type_getters.jl
@@ -85,11 +85,7 @@ 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,
-        ),
+        smagorinsky_lilly = get_smagorinsky_lilly_model(parsed_args),
         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 9693216804..11a04ff430 100644
--- a/src/solver/types.jl
+++ b/src/solver/types.jl
@@ -127,10 +127,7 @@ Base.@kwdef struct ViscousSponge{FT} <: AbstractSponge
 end
 
 abstract type AbstractEddyViscosityModel end
-Base.@kwdef struct SmagorinskyLilly{FT} <: AbstractEddyViscosityModel
-    Cs::FT = 0.2
-    Pr_t::FT = 1 / 3
-end
+struct SmagorinskyLilly <: AbstractEddyViscosityModel end
 
 Base.@kwdef struct RayleighSponge{FT} <: AbstractSponge
     zd::FT
diff --git a/test/parameterized_tendencies/les_closures/smagorinsky_lilly.jl b/test/parameterized_tendencies/les_closures/smagorinsky_lilly.jl
index 434267f760..1afa3f508c 100644
--- a/test/parameterized_tendencies/les_closures/smagorinsky_lilly.jl
+++ b/test/parameterized_tendencies/les_closures/smagorinsky_lilly.jl
@@ -57,13 +57,7 @@ import Thermodynamics as TD
     Y.c.uₕ .= Geometry.Covariant12Vector(vel)
     Y.f.u₃ .= Geometry.Covariant3Vector(w)
 
-    horizontal_smagorinsky_lilly_tendency(
-        Yₜ,
-        Y,
-        p,
-        t,
-        SmagorinskyLilly(FT(0.2)),
-    )
+    horizontal_smagorinsky_lilly_tendency(Yₜ, Y, p, t, SmagorinskyLilly())
 
     ### Component test begins here
 end
diff --git a/toml/box_density_current_test.toml b/toml/box_density_current_test.toml
new file mode 100644
index 0000000000..e9e27d1e03
--- /dev/null
+++ b/toml/box_density_current_test.toml
@@ -0,0 +1,7 @@
+[c_smag]
+value = 0.25
+type = "float"
+
+[mixing_length_Prandtl_number_0]
+value = 0.33333
+type = "float"
diff --git a/toml/les_isdac.toml b/toml/les_isdac.toml
index 2a70aea536..ff3461477f 100644
--- a/toml/les_isdac.toml
+++ b/toml/les_isdac.toml
@@ -1,2 +1,6 @@
 [zd_rayleigh]
 value = 2000.0
+
+[mixing_length_Prandtl_number_0]
+value = 0.33333
+type = "float"