modified: src/prognostic_equations/vertical_diffusion_boundary_layer.jl

src/prognostic_equations/vertical_diffusion_boundary_layer.jl

@@ -23,18 +23,19 @@ function compute_boundary_layer_height(Ri_c::FT) where {FT}
     return nothing
 end
 function compute_bulk_richardson_number(θᵥ::FT, θₐ::FT, zₐ::FT, norm_va, grav::FT, z::FT) where {FT}
-    norm²_va = norm_va^2
+    norm²_va = norm_va^2 + FT(10)
     return (grav*z)*(θᵥ-θₐ)/(θₐ*norm²_va)
 end
-function compute_eddy_diffusivity_coefficient(z::FT, h::FT, norm_va, Ri_a::FT) where {FT}
+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 (?)
-    if Ri_a < FT(0) # Unstable BL
-        C_E = κ^2 * (log(zₐ/z₀))^(-2)
-        K_b(z) = κ*norm_va*sqrt(C_E)*z
-    else # Stable or Neutral BL
+
+
+    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
@@ -94,7 +95,7 @@ function vertical_diffusion_boundary_layer_tendency!(
     interior_uₕ = Fields.level(Y.c.uₕ, 1)
     ᶜΔz_surface = Fields.Δz_field(interior_uₕ)
 
-    boundary_layer_height = similar(zₐ)
+    boundary_layer_height = FT(0)
     Ri_c = FT(1)
     ### Detect 𝒽, boundary layer height per column.
     for il in 1:Spaces.nlevels(axes(ᶜz))
@@ -105,15 +106,16 @@ function vertical_diffusion_boundary_layer_tendency!(
             norm(interior_uₕ[colidx]),
             grav,
             Spaces.level(ᶜz[colidx], il)[],
-        ) > Ri_c
-            @. boundary_layer_height[colidx] = Spaces.level(ᶜz[colidx], il)[] - zₐ[colidx] - ᶜΔz_surface[colidx]/2
+        ) < Ri_c
+             boundary_layer_height = Spaces.level(ᶜz[colidx], il)[]
         end
     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] * compute_eddy_diffusivity_coefficient(
             ᶜz[colidx],
-            boundary_layer_height[colidx], 
+            zₐ[colidx],
+            boundary_layer_height, 
             norm(interior_uₕ[colidx]),
             Ri_a[colidx]
            ))