Add Frierson(2006) diffusion profile

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

.buildkite/pipeline.yml

@@ -42,6 +42,7 @@ steps:
       - "julia --project=test -e 'using Pkg; Pkg.status()'"
 
       - echo "--- Instantiate examples"
+      - "julia --project=examples -e 'using Pkg; Pkg.develop(path = \".\")'"
       - "julia --project=examples -e 'using Pkg; Pkg.instantiate(;verbose=true)'"
       - "julia --project=examples -e 'using Pkg; Pkg.precompile()'"
       - "julia --project=examples -e 'using Pkg; Pkg.status()'"

src/prognostic_equations/vertical_diffusion_boundary_layer.jl

@@ -17,37 +17,89 @@ 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)
@@ -77,23 +129,70 @@ 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ₕ)
-    @. ᶠρK_E[colidx] =
-        ᶠinterp(Y.c.ρ[colidx]) * eddy_diffusivity_coefficient(
-            C_E,
+
+    ### 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] * compute_eddy_diffusivity_coefficient(
+            ᶜz[colidx],
+            z₀,
+            f_b,
+            h_boundary_layer,
             norm(interior_uₕ[colidx]),
-            ᶜΔz_surface[colidx] / 2,
-            ᶠp[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