modified: src/prognostic_equations/vertical_diffusion_boundary_layer.jl

src/prognostic_equations/vertical_diffusion_boundary_layer.jl

@@ -17,36 +17,31 @@ 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(θᵥ, θₐ, zₐ, norm_va, grav, z) 
+    norm²_va = norm_va^2 + eltype(θᵥ)(10)
+    return (grav * z) * (θᵥ-θₐ)/(θₐ*norm²_va)
+end
+function compute_eddy_diffusivity_coefficient(z::FT, zₐ::FT, h::FT, norm_va, Ri_a::FT) where {FT}
+    Ri_c = FT(1) # Critical Richardson number Table 1 Frierson et al (2006)
+    f_b = FT(0.1) # Surface layer fraction Table 1 Frierson et al (2006)
+    κ = FT(0.4) #TODO This is available in CLIMAParameters
+    z₀ = FT(3.21e-5) #TODO This is taken from Frierson et al (2006) but should be in sfc_conditions (?)
 
-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)
+
+    C_E = κ^2 * (log(zₐ/z₀))^(-2)
+    K_b(z) = κ*norm_va*sqrt(C_E)*z
+    if Ri_a >= FT(0) # Stable or Neutral BL
+        C_E = Ri_a > Ri_c ? FT(0) : κ^2 * (log(zₐ/z₀))^(-2)*(1-Ri_a/Ri_c)^2
+        K_b(z) = κ*norm_va*sqrt(C_E)*z*(1 + (Ri_a * log(z / z₀))/(Ri_c * (1-Ri_a/Ri_c)))^(-1)
+    end
+    if f_b * h < z < h
+        K_E = K_b(f_b*h)*(z/(f_b*h))*(1-(z-f_b*h)/(1-f_b)/h)^2
     else
-        error("Unsupported thermo option")
+        K_E = K_b(z)
     end
 end
 
@@ -77,23 +72,42 @@ 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
+    (; params) = p
+    thermo_params = CAP.thermodynamics_params(params)
+    (; ᶜts, ᶜp, ᶜspecific, sfc_conditions) = p.precomputed # assume ᶜts and ᶜp have been updated
     (; C_E) = p.atmos.vert_diff
 
     ᶠgradᵥ = Operators.GradientC2F() # apply BCs to ᶜdivᵥ, which wraps ᶠgradᵥ
+
+    ᶠρK_E = p.scratch.ᶠtemp_scalar
+    Ri_a = θᵥ = 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)
+    @. θᵥ = TD.virtual_pottemp(thermo_params, ᶜts)
+    θₐ = Fields.level(θᵥ,1)
+    θᵥ_sfc = @. TD.virtual_pottemp(thermo_params, ᶜts_sfc)
     interior_uₕ = Fields.level(Y.c.uₕ, 1)
-    ᶠp = ᶠρK_E = p.scratch.ᶠtemp_scalar
-    @. ᶠp[colidx] = ᶠinterp(ᶜp[colidx])
     ᶜΔz_surface = Fields.Δz_field(interior_uₕ)
+
+    boundary_layer_height = zₐ
+    Ri_c = FT(1)
+    ### Detect 𝒽, boundary layer height per column.
+    for il in 1:Spaces.nlevels(axes(ᶜz))
+    end
+    @. Ri_a[colidx] = compute_bulk_richardson_number(θₐ[colidx], θᵥ_sfc[colidx], zₐ[colidx], norm(interior_uₕ[colidx]), grav, ᶜz[colidx])
     @. ᶠρK_E[colidx] =
-        ᶠinterp(Y.c.ρ[colidx]) * eddy_diffusivity_coefficient(
-            C_E,
+        ᶠinterp(Y.c.ρ[colidx] * compute_eddy_diffusivity_coefficient(
+            ᶜz[colidx],
+            zₐ[colidx],
+            boundary_layer_height[colidx], 
             norm(interior_uₕ[colidx]),
-            ᶜΔz_surface[colidx] / 2,
-            ᶠp[colidx],
-        )
+            Ri_a[colidx]
+           ))
 
     if diffuse_momentum(p.atmos.vert_diff)
         ᶜdivᵥ_uₕ = Operators.DivergenceF2C(