change N_eff logic in l_smag

src/cache/cloud_fraction.jl

@@ -10,18 +10,18 @@ function covariance_from_grad(coeff, mixing_length, ∇Φ, ∇Ψ)
     return 2 * coeff * mixing_length^2 * dot(∇Φ, ∇Ψ)
 end
 
-# c_smag = CAP.c_smag(params)
-# N_eff = sqrt(max(ᶜlinear_buoygrad, 0))
-function compute_smagorinsky_length_scale(
-    c_smag,
-    N_eff,
-    dz,
-    Pr,
-    strain_rate_norm,
-)
+"""
+    Compute the Smagorinsky length scale from
+    - c_smag coefficient
+    - N_eff - buoyancy frequency = sqrt(max(ᶜlinear_buoygrad, 0))
+    - dz - verticla grid scale
+    - Pr - Prandtl number
+    - ϵ_st - strain rate norm
+"""
+function compute_smagorinsky_length_scale(c_smag, N_eff, dz, Pr, ϵ_st)
     FT = eltype(c_smag)
-    return N_eff > FT(0) ?
-       c_smag * dz * max(0, 1 - N_eff^2 / Pr / (2 * strain_rate_norm))^(1 / 4) :
+    return N_eff > FT(0) && N_eff < sqrt(2 * Pr * ϵ_st) ?
+       c_smag * dz * (1 - N_eff^2 / Pr / 2 / ϵ_st)^(1 / 4) :
        c_smag * dz
 end
 
@@ -90,18 +90,15 @@ function set_cloud_fraction!(Y, p, ::Union{EquilMoistModel, NonEquilMoistModel})
     )
 
     ᶜl_smag = p.scratch.ᶜtemp_scalar_2
-    # N_eff = sqrt(max(ᶜlinear_buoygrad, 0))
     @. ᶜl_smag = compute_smagorinsky_length_scale(
         CAP.c_smag(params),
         sqrt(max(ᶜlinear_buoygrad, 0)),   #N_eff
         ᶜdz,
         ᶜprandtl_nvec,
         norm_sqr(ᶜstrain_rate),
     )
-    @info(maximum(ᶜl_smag))
 
     coeff = FT(2.1) # TODO - move to parameters
-    # TODO - replace l_smag with mixing_length when using EDMF SGS
     @. ᶜcloud_fraction = quad_loop(
         SG_quad,
         ᶜp,
@@ -114,7 +111,7 @@ function set_cloud_fraction!(Y, p, ::Union{EquilMoistModel, NonEquilMoistModel})
             TD.liquid_ice_pottemp(thermo_params, ᶜts)
         ))),
         coeff,
-        ᶜl_smag,
+        ᶜl_smag, # replace with mixing_length when using EDMF SGS
         thermo_params,
     )
 end