Add minimum exchange coeff + exponentially decaying diffusion to Frierson functions

src/cache/precomputed_quantities.jl

@@ -336,7 +336,7 @@ function eddy_diffusivity_coefficient(
     κ::FT,
 ) where {FT}
     # Equations (17), (18)
-    if z < f_b * h
+    if z <= f_b * h
         K_b =
             compute_surface_layer_diffusivity(z, z₀, κ, C_E, Ri, Ri_a, Ri_c, uₐ)
         return K_b
@@ -360,22 +360,21 @@ end
 
 function compute_boundary_layer_height!(
     h_boundary_layer,
-    f_b::FT,
     dz,
     Ri_local,
     Ri_c::FT,
     Ri_a,
 ) where {FT}
     nlevels = Spaces.nlevels(Spaces.axes(Ri_local))
-    for level in 1:nlevels
+    for level in 1:(nlevels - 1)
         h_boundary_layer .=
             ifelse.(
                 Fields.Field(
                     Fields.field_values(Fields.level(Ri_local, level)),
                     axes(h_boundary_layer),
                 ) .< Ri_c,
                 Fields.Field(
-                    Fields.field_values(Fields.level(dz, level)),
+                    Fields.field_values(Fields.level(dz, level + 1)),
                     axes(h_boundary_layer),
                 ),
                 h_boundary_layer,
@@ -390,17 +389,24 @@ function compute_bulk_richardson_number(
     grav,
     z::FT,
 ) where {FT}
-    # TODO Gustiness from params
-    return (grav * z) * (θ_v - θ_v_a) / (θ_v_a * (norm_ua)^2 + FT(10))
+    # TODO Gustiness from ClimaParams
+    return (grav * z) * (θ_v - θ_v_a) / (θ_v_a * (max((norm_ua)^2, FT(10))))
 end
-function compute_exchange_coefficient(Ri_a, Ri_c, zₐ, z₀, κ::FT) where {FT}
+function compute_exchange_coefficient(
+    Ri_a,
+    Ri_c,
+    zₐ,
+    z₀,
+    κ::FT,
+    C_E_min::FT,
+) where {FT}
     # Equations (12), (13), (14)
-    if Ri_a < FT(0)
+    if Ri_a <= FT(0)
         return κ^2 * (log(zₐ / z₀))^(-2)
     elseif FT(0) < Ri_a < Ri_c
         return κ^2 * (log(zₐ / z₀))^(-2) * (1 - Ri_a / Ri_c)^2
     else
-        return FT(0)
+        return FT(C_E_min)
     end
 end
 
@@ -415,17 +421,19 @@ function compute_surface_layer_diffusivity(
     norm_uₐ,
 ) where {FT}
     # Equations (19), (20)
-    if Ri_a < FT(0)
-        return κ * norm_uₐ * sqrt(C_E) * z
+    if Ri_a <= FT(0)
+        return max(κ * norm_uₐ * sqrt(C_E) * z, FT(1))
     else
-        return κ *
-               norm_uₐ *
-               sqrt(C_E) *
-               z *
-               (1 + Ri / Ri_c * (log(z / z₀) / (1 - Ri / Ri_c)))^(-1)
+        return max(
+            κ *
+            norm_uₐ *
+            sqrt(C_E) *
+            z *
+            (1 + Ri / Ri_c * (log(z / z₀) / (1 - Ri / Ri_c)))^(-1),
+            FT(1),
+        )
     end
 end
-###
 
 """
     set_precomputed_quantities!(Y, p, t)
@@ -551,9 +559,9 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
         z₀ = FT(1e-5)
         Ri_c = FT(1.0)
         f_b = FT(0.1)
+        C_E_min = p.atmos.vert_diff.C_E
 
         # Prepare scratch vars
-        ᶠρK_E = p.scratch.ᶠtemp_scalar
         θ_v = p.scratch.ᶜtemp_scalar
         Ri = p.scratch.ᶜtemp_scalar_2
         dz_local = p.scratch.ᶜtemp_scalar_3
@@ -571,15 +579,13 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
         @. θ_v_sfc = TD.virtual_pottemp(thermo_params, ᶠts_sfc)
         θ_v_a = Fields.level(θ_v, 1)
 
-        ## Compute boundary layer height
-
-        ## TODO: Cache elevation field?
         z_local .= Fields.field_values(Fields.coordinate_field(Y.c).z)
         z_sfc .= Fields.field_values(
             Fields.level(Fields.coordinate_field(Y.f).z, half),
         )
         @. z_local = z_local - z_sfc
         dz_local .= Fields.Field(z_local, axes(Y.c))
+        zₐ = Fields.level(dz_local, 1)
         ᶜθ_v_sfc .=
             Fields.Field(Fields.field_values(θ_v_sfc), axes(interior_uₕ))
 
@@ -599,19 +605,24 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
         )
 
         #### Detect 𝒽, boundary layer height per column
-        h_boundary_layer = f_b .* Fields.level(ᶜz, Spaces.nlevels(axes(Y.c)))
+        h_boundary_layer = ᶜΔz_surface ./ 2 .+ FT(1000)
         compute_boundary_layer_height!(
             h_boundary_layer,
-            f_b,
             dz_local,
             Ri,
             Ri_c,
             Ri_a,
         )
 
         ## Exchange coefficients
-        @. C_E =
-            compute_exchange_coefficient(Ri_a, Ri_c, ᶜΔz_surface ./ 2, z₀, κ)
+        @. C_E = compute_exchange_coefficient(
+            Ri_a,
+            Ri_c,
+            ᶜΔz_surface ./ 2,
+            z₀,
+            κ,
+            C_E_min,
+        )
         @. ᶜK_h = eddy_diffusivity_coefficient(
             dz_local,
             z₀,

src/solver/model_getters.jl

@@ -108,7 +108,7 @@ function get_vertical_diffusion_model(
     elseif vert_diff_name in ("true", true, "VerticalDiffusion")
         VerticalDiffusion{diffuse_momentum, FT}(; C_E = params.C_E)
     elseif vert_diff_name in ("FriersonDiffusion",)
-        FriersonDiffusion{diffuse_momentum, FT}()
+        FriersonDiffusion{diffuse_momentum, FT}(; C_E = params.C_E)
     else
         error("Uncaught diffusion model `$vert_diff_name`.")
     end

src/solver/types.jl

@@ -59,7 +59,9 @@ Base.@kwdef struct VerticalDiffusion{DM, FT} <: AbstractVerticalDiffusion
     C_E::FT
 end
 diffuse_momentum(::VerticalDiffusion{DM}) where {DM} = DM
-struct FriersonDiffusion{DM, FT} <: AbstractVerticalDiffusion end
+Base.@kwdef struct FriersonDiffusion{DM, FT} <: AbstractVerticalDiffusion
+    C_E::FT
+end
 diffuse_momentum(::FriersonDiffusion{DM}) where {DM} = DM
 diffuse_momentum(::Nothing) = false