Changes to be committed:

modified:   src/prognostic_equations/vertical_diffusion_boundary_layer.jl
modified:   src/utils/common_spaces.jl

src/prognostic_equations/vertical_diffusion_boundary_layer.jl

@@ -17,38 +17,45 @@ import ClimaCore.Operators as Operators
 # 1) turn the liquid_theta into theta version
 # 2) have a total energy version (primary goal)
 
-function eddy_diffusivity_coefficient(C_E::FT, norm_v_a, z_a, p) where {FT}
-    p_pbl = FT(85000)
-    p_strato = FT(10000)
-    K_E = C_E * norm_v_a * z_a
-    return p > p_pbl ? K_E : K_E * exp(-((p_pbl - p) / p_strato)^2)
+function compute_boundary_layer_height(Ri_c::FT) where {FT}
+    # Needs a column search given colidx
+    return nothing
+end
+function compute_bulk_richardson_number(θ_v::FT, θ_a, norm_ua, grav, z) where {FT}
+    norm²_ua = norm_ua^2# SurfaceFluxes uses FT(1) gustiness #TODO Add to ClimaParameters?
+    return (grav * z) * (θ_v-θ_a)/(θ_a*norm²_ua + eps(FT))
+end
+function compute_diffusion_coefficient(Ri_a::FT, Ri_c::FT, zₐ::FT, z₀::FT, κ::FT) where {FT}
+    # Equations (12), (13), (14)
+    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)
+    end
 end
 
-function surface_thermo_state(
-    ::GCMSurfaceThermoState,
-    thermo_params,
-    T_sfc,
-    ts_int,
-    t,
-)
-    ρ_sfc =
-        TD.air_density(thermo_params, ts_int) *
-        (
-            T_sfc / TD.air_temperature(thermo_params, ts_int)
-        )^(
-            TD.cv_m(thermo_params, ts_int) /
-            TD.gas_constant_air(thermo_params, ts_int)
-        )
-    q_sfc =
-        TD.q_vap_saturation_generic(thermo_params, T_sfc, ρ_sfc, TD.Liquid())
-    if ts_int isa TD.PhaseDry
-        return TD.PhaseDry_ρT(thermo_params, ρ_sfc, T_sfc)
-    elseif ts_int isa TD.PhaseEquil
-        return TD.PhaseEquil_ρTq(thermo_params, ρ_sfc, T_sfc, q_sfc)
+function compute_surface_layer_diffusivity(z::FT, z₀::FT, κ::FT, C::FT, Ri::FT, Ri_c::FT, Ri_a::FT, norm_uₐ) where {FT}
+    # Equations (19), (20)
+    if Ri_a < FT(0)
+        return κ * norm_uₐ * sqrt(C) * z
     else
-        error("Unsupported thermo option")
+        return κ * norm_uₐ * sqrt(C) * z * (1 + Ri/Ri_c * (log(z/z₀)/(1 - Ri/Ri_c)))^(-1)
     end
 end
+function compute_eddy_diffusivity_coefficient(z::FT, z₀, f_b::FT, h::FT, uₐ, C::FT, Ri::FT, Ri_a::FT, Ri_c::FT, κ::FT) where {FT}
+    # Equations (17), (18)
+    K = FT(0)
+    if z < f_b*h
+        K_b = compute_surface_layer_diffusivity(z, z₀, κ,C, Ri, Ri_c, Ri_a, uₐ)
+        K = K_b
+    elseif f_b*h < z < h
+        K_b = compute_surface_layer_diffusivity(f_b*h, z₀, κ, C, Ri, Ri_c, Ri_a, uₐ)
+        K =  K_b * (z/f_b/h)*(1-(z-f_b*h)/(1-f_b)/h)^2
+    end
+    return K
+end
 
 function vertical_diffusion_boundary_layer_tendency!(Yₜ, Y, p, t)
     Fields.bycolumn(axes(Y.c.uₕ)) do colidx
@@ -77,23 +84,46 @@ function vertical_diffusion_boundary_layer_tendency!(
 )
     ᶜρ = Y.c.ρ
     FT = Spaces.undertype(axes(ᶜρ))
-    (; ᶜp, ᶜspecific, sfc_conditions) = p.precomputed # assume ᶜts and ᶜp have been updated
-    (; C_E) = p.atmos.vert_diff
+    (; params) = p
+    thermo_params = CAP.thermodynamics_params(params)
+    (; ᶜts, ᶜp, ᶜspecific, sfc_conditions) = p.precomputed # assume ᶜts and ᶜp have been updated
+
+    z₀ = FT(3.21e-5)
+    f_b = FT(0.1)
+    κ = FT(0.4)
+    Ri_c = FT(1)
 
     ᶠgradᵥ = Operators.GradientC2F() # apply BCs to ᶜdivᵥ, which wraps ᶠgradᵥ
+
+    ᶠρK_E = p.scratch.ᶠtemp_scalar
+    θ_v = p.scratch.ᶜtemp_scalar
+    grav = FT(9.81)
 
+    # Compute boundary layer height
     FT = eltype(Y)
+    ᶜts_sfc = sfc_conditions.ts
+    ᶜz = Fields.coordinate_field(Y.c).z
+    zₐ = Fields.level(ᶜz, 1)
     interior_uₕ = Fields.level(Y.c.uₕ, 1)
-    ᶠp = ᶠρK_E = p.scratch.ᶠtemp_scalar
-    @. ᶠp[colidx] = ᶠinterp(ᶜp[colidx])
+    @. θ_v[colidx] = TD.virtual_pottemp(thermo_params, ᶜts[colidx])
+    θ_a = Fields.level(θ_v, 1)
+    θ_v_sfc = @. TD.virtual_pottemp(thermo_params, ᶜts_sfc[colidx])
     ᶜΔz_surface = Fields.Δz_field(interior_uₕ)
+
+    ### Detect 𝒽, boundary layer height per column.
+    Ri = @. compute_bulk_richardson_number(θ_v[colidx], θ_a[colidx], norm(Y.c.uₕ[colidx]), grav, ᶜz[colidx])
+    Ri_a = @. compute_bulk_richardson_number(θ_a[colidx], θ_v_sfc, norm(interior_uₕ[colidx]), grav, zₐ[colidx])
+
+    h_boundary_layer = FT(0)
+    for level in 1:Spaces.nlevels(axes(ᶜz))
+        if Spaces.level(Ri[colidx], level)[] < Ri_c
+            h_boundary_layer = Spaces.level(ᶜz[colidx], level)[]
+        end
+    end
+
+    C = @. compute_diffusion_coefficient(Ri_a, Ri_c, zₐ[colidx], z₀, κ)
     @. ᶠρK_E[colidx] =
-        ᶠinterp(Y.c.ρ[colidx]) * eddy_diffusivity_coefficient(
-            C_E,
-            norm(interior_uₕ[colidx]),
-            ᶜΔz_surface[colidx] / 2,
-            ᶠp[colidx],
-        )
+    ᶠinterp(Y.c.ρ[colidx] * compute_eddy_diffusivity_coefficient(ᶜz[colidx], z₀, f_b, h_boundary_layer, norm(interior_uₕ[colidx]), C[colidx], Ri[colidx], Ri_a, Ri_c, κ))
 
     if diffuse_momentum(p.atmos.vert_diff)
         ᶜdivᵥ_uₕ = Operators.DivergenceF2C(

src/utils/common_spaces.jl

@@ -75,12 +75,12 @@ end
 
 function make_hybrid_spaces(
     h_space,
-    z_max,
+    z_max::FT,
     z_elem,
     z_stretch;
     surface_warp = nothing,
     topo_smoothing = false,
-)
+) where {FT}
     z_domain = Domains.IntervalDomain(
         Geometry.ZPoint(zero(z_max)),
         Geometry.ZPoint(z_max);
@@ -103,7 +103,7 @@ function make_hybrid_spaces(
         face_space = Spaces.ExtrudedFiniteDifferenceSpace(
             h_space,
             z_face_space,
-            Hypsography.LinearAdaption(z_surface),
+            Hypsography.SLEVEAdaption(z_surface, FT(0.75), FT(10)),
         )
         center_space = Spaces.CenterExtrudedFiniteDifferenceSpace(face_space)
     end