From aeb97fc1a68aafc886083e084314194a22ca3ef3 Mon Sep 17 00:00:00 2001 From: Akshay Sridhar Date: Thu, 25 Jul 2024 09:03:01 -0700 Subject: [PATCH] Use C_E_min in surface exchange coeff and min diffusivity in surface layer --- src/cache/precomputed_quantities.jl | 61 +++++++++++++++++------------ src/solver/model_getters.jl | 2 +- src/solver/types.jl | 4 +- 3 files changed, 40 insertions(+), 27 deletions(-) diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index 5aac993507..d0b9b0b560 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -336,7 +336,7 @@ function eddy_diffusivity_coefficient( κ::FT, ) where {FT} # Equations (17), (18) - if z < f_b * h + if z <= f_b * h K_b = compute_surface_layer_diffusivity(z, z₀, κ, C_E, Ri, Ri_a, Ri_c, uₐ) return K_b @@ -360,14 +360,13 @@ end function compute_boundary_layer_height!( h_boundary_layer, - f_b::FT, dz, Ri_local, Ri_c::FT, Ri_a, ) where {FT} nlevels = Spaces.nlevels(Spaces.axes(Ri_local)) - for level in 1:nlevels + for level in 1:(nlevels - 1) h_boundary_layer .= ifelse.( Fields.Field( @@ -375,7 +374,7 @@ function compute_boundary_layer_height!( axes(h_boundary_layer), ) .< Ri_c, Fields.Field( - Fields.field_values(Fields.level(dz, level)), + Fields.field_values(Fields.level(dz, level + 1)), axes(h_boundary_layer), ), h_boundary_layer, @@ -390,17 +389,24 @@ function compute_bulk_richardson_number( grav, z::FT, ) where {FT} - # TODO Gustiness from params - return (grav * z) * (θ_v - θ_v_a) / (θ_v_a * (norm_ua)^2 + FT(10)) + # TODO Gustiness from ClimaParams + return (grav * z) * (θ_v - θ_v_a) / (θ_v_a * (max((norm_ua)^2, FT(10)))) end -function compute_exchange_coefficient(Ri_a, Ri_c, zₐ, z₀, κ::FT) where {FT} +function compute_exchange_coefficient( + Ri_a, + Ri_c, + zₐ, + z₀, + κ::FT, + C_E_min::FT, +) where {FT} # Equations (12), (13), (14) - if Ri_a < FT(0) + 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) + return FT(C_E_min) end end @@ -415,17 +421,19 @@ function compute_surface_layer_diffusivity( norm_uₐ, ) where {FT} # Equations (19), (20) - if Ri_a < FT(0) - return κ * norm_uₐ * sqrt(C_E) * z + if Ri_a <= FT(0) + return max(κ * norm_uₐ * sqrt(C_E) * z, FT(1)) else - return κ * - norm_uₐ * - sqrt(C_E) * - z * - (1 + Ri / Ri_c * (log(z / z₀) / (1 - Ri / Ri_c)))^(-1) + return max( + κ * + norm_uₐ * + sqrt(C_E) * + z * + (1 + Ri / Ri_c * (log(z / z₀) / (1 - Ri / Ri_c)))^(-1), + FT(1), + ) end end -### """ set_precomputed_quantities!(Y, p, t) @@ -551,9 +559,9 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) z₀ = FT(1e-5) Ri_c = FT(1.0) f_b = FT(0.1) + C_E_min = p.atmos.vert_diff.C_E # Prepare scratch vars - ᶠρK_E = p.scratch.ᶠtemp_scalar θ_v = p.scratch.ᶜtemp_scalar Ri = p.scratch.ᶜtemp_scalar_2 dz_local = p.scratch.ᶜtemp_scalar_3 @@ -571,15 +579,13 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) @. θ_v_sfc = TD.virtual_pottemp(thermo_params, ᶠts_sfc) θ_v_a = Fields.level(θ_v, 1) - ## Compute boundary layer height - - ## TODO: Cache elevation field? z_local .= Fields.field_values(Fields.coordinate_field(Y.c).z) z_sfc .= Fields.field_values( Fields.level(Fields.coordinate_field(Y.f).z, half), ) @. z_local = z_local - z_sfc dz_local .= Fields.Field(z_local, axes(Y.c)) + zₐ = Fields.level(dz_local, 1) ᶜθ_v_sfc .= Fields.Field(Fields.field_values(θ_v_sfc), axes(interior_uₕ)) @@ -599,10 +605,9 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) ) #### Detect 𝒽, boundary layer height per column - h_boundary_layer = f_b .* Fields.level(ᶜz, Spaces.nlevels(axes(Y.c))) + h_boundary_layer = ᶜΔz_surface ./ 2 .+ FT(1000) compute_boundary_layer_height!( h_boundary_layer, - f_b, dz_local, Ri, Ri_c, @@ -610,8 +615,14 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) ) ## Exchange coefficients - @. C_E = - compute_exchange_coefficient(Ri_a, Ri_c, ᶜΔz_surface ./ 2, z₀, κ) + @. C_E = compute_exchange_coefficient( + Ri_a, + Ri_c, + ᶜΔz_surface ./ 2, + z₀, + κ, + C_E_min, + ) @. ᶜK_h = eddy_diffusivity_coefficient( dz_local, z₀, diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl index 284cafc14a..93093f14d4 100644 --- a/src/solver/model_getters.jl +++ b/src/solver/model_getters.jl @@ -108,7 +108,7 @@ function get_vertical_diffusion_model( elseif vert_diff_name in ("true", true, "VerticalDiffusion") VerticalDiffusion{diffuse_momentum, FT}(; C_E = params.C_E) elseif vert_diff_name in ("FriersonDiffusion",) - FriersonDiffusion{diffuse_momentum, FT}() + FriersonDiffusion{diffuse_momentum, FT}(; C_E = params.C_E) else error("Uncaught diffusion model `$vert_diff_name`.") end diff --git a/src/solver/types.jl b/src/solver/types.jl index 43a83e1260..f54400042e 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -59,7 +59,9 @@ Base.@kwdef struct VerticalDiffusion{DM, FT} <: AbstractVerticalDiffusion C_E::FT end diffuse_momentum(::VerticalDiffusion{DM}) where {DM} = DM -struct FriersonDiffusion{DM, FT} <: AbstractVerticalDiffusion end +Base.@kwdef struct FriersonDiffusion{DM, FT} <: AbstractVerticalDiffusion + C_E::FT +end diffuse_momentum(::FriersonDiffusion{DM}) where {DM} = DM diffuse_momentum(::Nothing) = false