modified: ../src/parameterized_tendencies/les_sgs_models/smagorinsky…

…_lilly.jl

src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl

@@ -16,7 +16,7 @@ function smagorinsky_lilly_cache(Y, sl::SmagorinskyLilly)
     (; Cs) = sl
     FT = eltype(Y)
     h_space = Spaces.horizontal_space(axes(Y.c))
-    Δ_filter = Fields.Δz_field(axes(Y.f))
+    Δ_filter = Spaces.node_horizontal_length_scale(h_space)
     ᶜtemp_scalar_3 = zero(Fields.Field(Float32, axes(Y.c.ρ)))
     v_t = ᶜtemp_scalar_3
     return (; v_t, Δ_filter)
@@ -27,9 +27,11 @@ vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::Nothing) = nothing
 
 function compute_strain_rate_f!(ϵ::Fields.Field, u::Fields.Field)
     @assert eltype(u) <: C123
+    ∇ᵥuvw_boundary =
+        Geometry.outer(Geometry.WVector(0), Geometry.UVWVector(0, 0, 0))
     ᶠgradᵥ = Operators.GradientC2F(
        bottom = Operators.SetValue(Geometry.UVWVector(0,0,0)),
-       top = Operators.SetValue(Geometry.UVWVector(0,0,0)),
+       top = Operators.SetGradient(∇ᵥuvw_boundary),
     )
     axis_uvw = Geometry.UVWAxis()
     @. ϵ =
@@ -78,8 +80,7 @@ function horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl::SmagorinskyLi
     CA.compute_strain_rate_center!(ᶜϵ, Geometry.Covariant123Vector.(ᶠu))
     compute_strain_rate_f!(ᶠϵ, Geometry.Covariant123Vector.(localu))
     @. ᶜS = (ᶜS + adjoint(ᶜS)) + ᶜϵ
-    @. ᶠS = ᶠinterp(ᶜS)
-
+    @. ᶠS = ᶠinterp(ᶜS + adjoint(ᶜS)) + ᶠϵ
     thermo_params = CAP.thermodynamics_params(p.params)
 
     ### Interp Checks
@@ -93,7 +94,6 @@ function horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl::SmagorinskyLi
     ∇θ = @. ᶠinterp(ᶜ∇θ)
     ###
 
-    θ_v = @. TD.virtual_pottemp(thermo_params, ᶜts)
     N² = @. grav / ᶠinterp(θ_v) * Geometry.WVector(∇θ).components.data.:1
     @. ᶠfb = (max(FT(0), 
                   1 - 3*(N²) / (CA.norm_sqr(ᶠS) + eps(FT))))^(1/2)
@@ -102,15 +102,15 @@ function horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl::SmagorinskyLi
     ᶠρ = @. ᶠwinterp(ᶜJ, Y.c.ρ)
     ᶠv_t = @. (Cs * Δ_filter)^2 * sqrt(2 * CA.norm_sqr(ᶠS)) * ᶠfb
     ᶜv_t = @. ᶜinterp(ᶠv_t)
-    uc2f = Operators.InterpolateC2F(;top = Operators.Extrapolate(), 
-                                    bottom = Operators.Extrapolate())
-    ᶠv_t = @. uc2f(ᶜv_t)
     ᶜD = @. FT(3) * ᶜv_t
 
     @. v_t = ᶜv_t
 
-    @. Yₜ.c.uₕ += C12(wdivₕ(Y.c.ρ * ᶜv_t * ᶜS)) / Y.c.ρ
-    @. Yₜ.f.u₃ += C3(wdivₕ(ᶠρ * ᶠv_t * ᶠS)) / ᶠρ
+    ᶜτ = @. FT(2) * ᶜv_t * ᶜS
+    ᶠτ = @. FT(2) * ᶠv_t * ᶠS
+
+    @. Yₜ.c.uₕ += C12(wdivₕ(Y.c.ρ * ᶜτ)) / Y.c.ρ
+    @. Yₜ.f.u₃ += C3(wdivₕ(ᶠρ * ᶠτ)) / ᶠρ
 
     # energy adjustment
     (; ᶜspecific) = p.precomputed
@@ -127,7 +127,6 @@ function horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl::SmagorinskyLi
         @. ᶜρχₜ += wdivₕ(Y.c.ρ * ᶜD * gradₕ(ᶜχ)) 
         @. Yₜ.c.ρ += wdivₕ(Y.c.ρ * ᶜD * gradₕ(ᶜχ)) 
     end
-
 end
 
 function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl::SmagorinskyLilly) 
@@ -163,7 +162,7 @@ function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl::SmagorinskyLill
     CA.compute_strain_rate_center!(ᶜϵ, Geometry.Covariant123Vector.(ᶠu))
     compute_strain_rate_f!(ᶠϵ, Geometry.Covariant123Vector.(localu))
     @. ᶜS = (ᶜS + adjoint(ᶜS)) + ᶜϵ
-    @. ᶠS = ᶠinterp(ᶜS)
+    @. ᶠS = ᶠinterp(ᶜS + adjoint(ᶜS)) + ᶠϵ
 
     thermo_params = CAP.thermodynamics_params(p.params)
     θ_v = @. TD.virtual_pottemp(thermo_params, ᶜts)
@@ -184,9 +183,6 @@ function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl::SmagorinskyLill
     ᶠρ = @. ᶠwinterp(ᶜJ, Y.c.ρ)
     ᶠv_t = @. (Cs * Δ_filter)^2 * sqrt(2 * CA.norm_sqr(ᶠS)) * ᶠfb
     ᶜv_t = @. ᶜinterp(ᶠv_t)
-    uc2f = Operators.InterpolateC2F(;top = Operators.Extrapolate(), 
-                                    bottom = Operators.Extrapolate())
-    ᶠv_t = @. uc2f(ᶜv_t)
     ᶜD = @. FT(3) * ᶜv_t
 
     @. Yₜ.c.uₕ -= C12(