From 5f6bf1a2f12bbf69d753554c79882d766178b054 Mon Sep 17 00:00:00 2001 From: akshaysridhar Date: Thu, 16 Nov 2023 11:30:16 -0800 Subject: [PATCH] Add Frierson(2006) diffusion profile modified: src/prognostic_equations/vertical_diffusion_boundary_layer.jl modified: src/utils/common_spaces.jl +Formatter --- .buildkite/pipeline.yml | 1 + .../vertical_diffusion_boundary_layer.jl | 173 ++++++++++++++---- src/utils/common_spaces.jl | 6 +- 3 files changed, 140 insertions(+), 40 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 7974b5f12bc..9cbf0d47c8b 100644 --- a/.buildkite/pipeline.yml +++ b/.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()'" diff --git a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl b/src/prognostic_equations/vertical_diffusion_boundary_layer.jl index a949eb96173..c6090101db7 100644 --- a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl +++ b/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( diff --git a/src/utils/common_spaces.jl b/src/utils/common_spaces.jl index 9d1b18c8796..1b25821e2cc 100644 --- a/src/utils/common_spaces.jl +++ b/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