From 313f2b079b41bd5de4462182130085a98e46f848 Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Thu, 30 Nov 2023 13:42:21 -0800 Subject: [PATCH 1/4] Add cloud fraction diagnostics --- config/default_configs/default_config.yml | 3 + .../single_column_precipitation_test.yml | 1 + examples/hybrid/driver.jl | 5 + src/ClimaAtmos.jl | 2 +- src/cache/cache.jl | 11 +- src/cache/cloud_fraction.jl | 225 ++++++++++++++++++ .../diagnostic_edmf_precomputed_quantities.jl | 4 +- src/cache/precomputed_quantities.jl | 5 + .../prognostic_edmf_precomputed_quantities.jl | 2 + src/callbacks/callbacks.jl | 55 ++--- src/diagnostics/core_diagnostics.jl | 15 ++ .../buoyancy_gradients.jl | 24 +- src/prognostic_equations/cloud_fraction.jl | 113 --------- src/prognostic_equations/edmfx_closures.jl | 13 - src/solver/type_getters.jl | 4 + src/solver/types.jl | 14 +- test/coupler_compatibility.jl | 2 +- 17 files changed, 323 insertions(+), 175 deletions(-) create mode 100644 src/cache/cloud_fraction.jl delete mode 100644 src/prognostic_equations/cloud_fraction.jl diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml index 90e2df29e6..ea9473c786 100644 --- a/config/default_configs/default_config.yml +++ b/config/default_configs/default_config.yml @@ -93,6 +93,9 @@ idealized_clouds: idealized_insolation: help: "Use idealized insolation in radiation model [`false`, `true` (default)]" value: true +dt_cloud_fraction: + help: "Time between calling cloud fraction update" + value: "3hours" config: diff --git a/config/model_configs/single_column_precipitation_test.yml b/config/model_configs/single_column_precipitation_test.yml index c81c92da08..9c31ee109f 100644 --- a/config/model_configs/single_column_precipitation_test.yml +++ b/config/model_configs/single_column_precipitation_test.yml @@ -7,6 +7,7 @@ z_stretch: false dt: "10secs" t_end: "1500secs" dt_save_to_disk: "500secs" +dt_cloud_fraction: "60secs" moist: "nonequil" precip_model: "1M" precip_upwinding: "first_order" diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl index 11e8e31108..14a96e0dda 100644 --- a/examples/hybrid/driver.jl +++ b/examples/hybrid/driver.jl @@ -175,5 +175,10 @@ if config.parsed_args["check_precipitation"] @test all( ClimaCore.isapprox(Yₜ_ρ[colidx], Yₜ_ρqₜ[colidx], rtol = eps(FT)), ) + + # cloud fraction diagnostics + @assert !any(isnan, sol.prob.p.precomputed.ᶜcloud_fraction[colidx]) + @test minimum(sol.prob.p.precomputed.ᶜcloud_fraction[colidx]) >= FT(0) + @test maximum(sol.prob.p.precomputed.ᶜcloud_fraction[colidx]) <= FT(1) end end diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl index dfdf69930a..51b7e197e8 100644 --- a/src/ClimaAtmos.jl +++ b/src/ClimaAtmos.jl @@ -25,6 +25,7 @@ include(joinpath("cache", "prognostic_edmf_precomputed_quantities.jl")) include(joinpath("cache", "diagnostic_edmf_precomputed_quantities.jl")) include(joinpath("cache", "precipitation_precomputed_quantities.jl")) include(joinpath("cache", "precomputed_quantities.jl")) +include(joinpath("cache", "cloud_fraction.jl")) include(joinpath("initial_conditions", "InitialConditions.jl")) include( @@ -80,7 +81,6 @@ include(joinpath("prognostic_equations", "edmfx_entr_detr.jl")) include(joinpath("prognostic_equations", "edmfx_tke.jl")) include(joinpath("prognostic_equations", "edmfx_sgs_flux.jl")) include(joinpath("prognostic_equations", "edmfx_boundary_condition.jl")) -include(joinpath("prognostic_equations", "cloud_fraction.jl")) include( joinpath("parameterized_tendencies", "microphysics", "precipitation.jl"), ) diff --git a/src/cache/cache.jl b/src/cache/cache.jl index 354e1b2813..4e7cd73b0d 100644 --- a/src/cache/cache.jl +++ b/src/cache/cache.jl @@ -7,7 +7,7 @@ struct AtmosCache{ COR, SFC, GHOST, - ENV, + SGQ, PREC, SCRA, HYPE, @@ -49,7 +49,8 @@ struct AtmosCache{ """Center and face ghost buffers used by DSS""" ghost_buffer::GHOST - env_thermo_quad::ENV + """Struct with sub-grid sampling quadrature""" + SG_quad::SGQ """Quantities that are updated with set_precomputed_quantities!""" precomputed::PREC @@ -153,11 +154,11 @@ function build_cache(Y, atmos, params, surface_setup, dt, start_date) sfc_setup = surface_setup(params) scratch = temporary_quantities(Y, atmos) - env_thermo_quad = SGSQuadrature(FT) + SG_quad = SGSQuadrature(FT) precomputed = precomputed_quantities(Y, atmos) precomputing_arguments = - (; atmos, core, params, sfc_setup, precomputed, scratch, dt) + (; atmos, core, params, sfc_setup, precomputed, scratch, dt, SG_quad) # Coupler compatibility isnothing(precomputing_arguments.sfc_setup) && @@ -190,7 +191,7 @@ function build_cache(Y, atmos, params, surface_setup, dt, start_date) core, sfc_setup, ghost_buffer, - env_thermo_quad, + SG_quad, precomputed, scratch, hyperdiff, diff --git a/src/cache/cloud_fraction.jl b/src/cache/cloud_fraction.jl new file mode 100644 index 0000000000..79d9e83123 --- /dev/null +++ b/src/cache/cloud_fraction.jl @@ -0,0 +1,225 @@ +import StaticArrays as SA +import ClimaCore.RecursiveApply: rzero, ⊞, ⊠ + +# TODO: write a test with scalars that are linear with z +""" + Diagnose horizontal covariances based on vertical gradients + (i.e. taking turbulence production as the only term) +""" +function covariance_from_grad(coeff, mixing_length, ∇Φ, ∇Ψ) + return 2 * coeff * mixing_length^2 * dot(∇Φ, ∇Ψ) +end + +""" + Compute the Smagorinsky length scale from + - c_smag coefficient + - N_eff - buoyancy frequency = sqrt(max(ᶜlinear_buoygrad, 0)) + - dz - vertical grid scale + - Pr - Prandtl number + - ϵ_st - strain rate norm +""" +function compute_smagorinsky_length_scale(c_smag, N_eff, dz, Pr, ϵ_st) + FT = eltype(c_smag) + #return N_eff > FT(0) && N_eff < sqrt(2 * Pr * ϵ_st) ? + return N_eff > FT(0) ? + c_smag * dz * max(0, (1 - N_eff^2 / Pr / 2 / ϵ_st))^(1 / 4) : + c_smag * dz +end + +""" + Compute the grid scale cloud fraction based on sub-grid scale properties +""" +function set_cloud_fraction!(Y, p, ::DryModel) + @. p.precomputed.ᶜcloud_fraction = 0 +end +function set_cloud_fraction!(Y, p, ::Union{EquilMoistModel, NonEquilMoistModel}) + (; SG_quad, params) = p + + FT = eltype(params) + thermo_params = CAP.thermodynamics_params(params) + ᶜdz = Fields.Δz_field(axes(Y.c)) + ᶜlg = Fields.local_geometry_field(Y.c) + (; ᶜts, ᶜp, ᶠu³, ᶜcloud_fraction) = p.precomputed + (; obukhov_length) = p.precomputed.sfc_conditions + + ᶜlinear_buoygrad = p.scratch.ᶜtemp_scalar + @. ᶜlinear_buoygrad = buoyancy_gradients( + params, + p.atmos.moisture_model, + EnvBuoyGrad( + BuoyGradMean(), + TD.air_temperature(thermo_params, ᶜts), # t_sat + TD.vapor_specific_humidity(thermo_params, ᶜts), # qv_sat + TD.total_specific_humidity(thermo_params, ᶜts), # qt_sat + TD.liquid_specific_humidity(thermo_params, ᶜts), # q_liq + TD.ice_specific_humidity(thermo_params, ᶜts), # q_ice + TD.dry_pottemp(thermo_params, ᶜts), # θ_sat + TD.liquid_ice_pottemp(thermo_params, ᶜts), # θ_liq_ice_sat + projected_vector_data( + C3, + ᶜgradᵥ(ᶠinterp(TD.virtual_pottemp(thermo_params, ᶜts))), + ᶜlg, + ), # ∂θv∂z_unsat + projected_vector_data( + C3, + ᶜgradᵥ(ᶠinterp(TD.total_specific_humidity(thermo_params, ᶜts))), + ᶜlg, + ), # ∂qt∂z_sat + projected_vector_data( + C3, + ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts))), + ᶜlg, + ), # ∂θl∂z_sat + ᶜp, # p + ifelse(TD.has_condensate(thermo_params, ᶜts), 1, 0),# en_cld_frac + Y.c.ρ, # ρ + ), + ) + + ᶠu = p.scratch.ᶠtemp_C123 + @. ᶠu = C123(ᶠinterp(Y.c.uₕ)) + C123(ᶠu³) + + ᶜstrain_rate = p.scratch.ᶜtemp_UVWxUVW + compute_strain_rate_center!(ᶜstrain_rate, ᶠu) + + ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar + @. ᶜprandtl_nvec = turbulent_prandtl_number( + params, + obukhov_length, + ᶜlinear_buoygrad, + norm_sqr(ᶜstrain_rate), + ) + + ᶜl_smag = p.scratch.ᶜtemp_scalar_2 + @. ᶜl_smag = compute_smagorinsky_length_scale( + CAP.c_smag(params), + sqrt(max(ᶜlinear_buoygrad, 0)), #N_eff + ᶜdz, + ᶜprandtl_nvec, + norm_sqr(ᶜstrain_rate), + ) + + coeff = FT(2.1) # TODO - move to parameters + @. ᶜcloud_fraction = quad_loop( + SG_quad, + ᶜp, + TD.total_specific_humidity(thermo_params, ᶜts), + TD.liquid_ice_pottemp(thermo_params, ᶜts), + Geometry.WVector( + ᶜgradᵥ(ᶠinterp(TD.total_specific_humidity(thermo_params, ᶜts))), + ), + Geometry.WVector( + ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts))), + ), + coeff, + ᶜl_smag, # replace with mixing_length when using EDMF SGS + thermo_params, + ) +end + +""" + function quad_loop(SG_quad, p_c, q_mean, θ_mean, ᶜ∇q, ᶜ∇θ, + coeff, ᶜlength_scale, thermo_params) + +where: + - SG_quad is a struct containing information about quadrature type and order + - p_c is the atmospheric pressure + - q_mean, θ_mean is the grid mean q_tot and liquid ice potential temperature + - ᶜ∇q, ᶜ∇θ are the gradients of q_tot and liquid ice potential temperature + - coeff - a free parameter (to be moved into params) + - ᶜlength_scale - mixing length for simulations with EDMF and Smagorinsky + length scale for simulations without EDMF + - thermo params - thermodynamics parameters + +The function imposes additional limits on the quadrature points and +returns cloud fraction computed as a sum over quadrature points. +""" +function quad_loop( + SG_quad::SGSQuadrature, + p_c, + q_mean, + θ_mean, + ᶜ∇q, + ᶜ∇θ, + coeff, + ᶜlength_scale, + thermo_params, +) + # Returns the physical values based on quadrature sampling points + # and limited covarainces + function get_x_hat(χ1, χ2) + + @assert SG_quad.quadrature_type isa GaussianQuad + FT = eltype(χ1) + + q′q′ = covariance_from_grad(coeff, ᶜlength_scale, ᶜ∇q, ᶜ∇q) + θ′θ′ = covariance_from_grad(coeff, ᶜlength_scale, ᶜ∇θ, ᶜ∇θ) + θ′q′ = covariance_from_grad(coeff, ᶜlength_scale, ᶜ∇θ, ᶜ∇q) + + # Epsilon defined per typical variable fluctuation + eps_q = eps(FT) * max(eps(FT), q_mean) + eps_θ = eps(FT) + + # limit σ_q to prevent negative q_tot_hat + σ_q_lim = -q_mean / (sqrt(FT(2)) * SG_quad.a[1]) + σ_q = min(sqrt(q′q′), σ_q_lim) + # Do we also have to try to limit θ in the same way as q?? + σ_θ = sqrt(θ′θ′) + + # Enforce Cauchy-Schwarz inequality, numerically stable compute + _corr = (θ′q′ / max(σ_q, eps_q)) + corr = max(min(_corr / max(σ_θ, eps_θ), FT(1)), FT(-1)) + + # Conditionals + σ_c = sqrt(max(1 - corr * corr, 0)) * σ_θ + + μ_c = θ_mean + sqrt(FT(2)) * corr * σ_θ * χ1 + θ_hat = μ_c + sqrt(FT(2)) * σ_c * χ2 + q_hat = q_mean + sqrt(FT(2)) * σ_q * χ1 + # The σ_q_lim limits q_tot_hat to be close to zero + # for the negative sampling points. However due to numerical erros + # we sometimes still get small negative numers here + return (θ_hat, max(FT(0), q_hat)) + end + + function f(x1_hat, x2_hat) + FT = eltype(x1_hat) + @assert(x1_hat >= FT(0)) + @assert(x2_hat >= FT(0)) + ts = thermo_state(thermo_params; p = p_c, θ = x1_hat, q_tot = x2_hat) + hc = TD.has_condensate(thermo_params, ts) + return (; + cf = hc ? FT(1) : FT(0), # cloud fraction + q_tot_sat = hc ? x2_hat : FT(0), # cloudy/dry for buoyancy in TKE + ) + end + + return quad(f, get_x_hat, SG_quad).cf +end + +""" + Compute f(θ, q) as a sum over inner and outer quadrature points + that approximate the sub-grid scale variability of θ and q. + + θ - liquid ice potential temperature + q - total water specific humidity +""" +function quad(f::F, get_x_hat::F1, quad) where {F <: Function, F1 <: Function} + χ = quad.a + weights = quad.w + quad_order = quadrature_order(quad) + FT = eltype(χ) + # zero outer quadrature points + T = typeof(f(get_x_hat(χ[1], χ[1])...)) + outer_env = rzero(T) + @inbounds for m_q in 1:quad_order + # zero inner quadrature points + inner_env = rzero(T) + for m_h in 1:quad_order + x_hat = get_x_hat(χ[m_q], χ[m_h]) + inner_env = inner_env ⊞ f(x_hat...) ⊠ weights[m_h] ⊠ FT(1 / sqrt(π)) + end + outer_env = outer_env ⊞ inner_env ⊠ weights[m_q] ⊠ FT(1 / sqrt(π)) + end + return outer_env +end diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl index 06edfee53c..6ed7bcc1f4 100644 --- a/src/cache/diagnostic_edmf_precomputed_quantities.jl +++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl @@ -449,7 +449,7 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t) nh_pressure³ʲ_data_prev_halflevel ) - # get u³ʲ to calculate divergence term for detrainment, + # get u³ʲ to calculate divergence term for detrainment, # u³ʲ will be clipped later after we get area fraction minimum_value = FT(1e-6) @. u³ʲ_halflevel = ifelse( @@ -717,6 +717,8 @@ function set_diagnostic_edmf_precomputed_quantities_env_closures!(Y, p, t) TD.air_temperature(thermo_params, ᶜts), # t_sat TD.vapor_specific_humidity(thermo_params, ᶜts), # qv_sat q_tot, # qt_sat + TD.liquid_specific_humidity(thermo_params, ᶜts), + TD.ice_specific_humidity(thermo_params, ᶜts), TD.dry_pottemp(thermo_params, ᶜts), # θ_sat TD.liquid_ice_pottemp(thermo_params, ᶜts), # θ_liq_ice_sat projected_vector_data( diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index 6f647b9013..759b94fd03 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -52,6 +52,7 @@ function precomputed_quantities(Y, atmos) ᶜh_tot = similar(Y.c, FT), sfc_conditions = Fields.Field(SCT, Spaces.level(axes(Y.f), half)), ) + cloud_diagnostics = (; ᶜcloud_fraction = similar(Y.c, FT),) advective_sgs_quantities = atmos.turbconv_model isa PrognosticEDMFX ? (; @@ -124,6 +125,7 @@ function precomputed_quantities(Y, atmos) diagnostic_sgs_quantities..., vert_diff_quantities..., precipitation_quantities..., + cloud_diagnostics..., ) end @@ -374,6 +376,9 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) set_precipitation_precomputed_quantities!(Y, p, t) end + # TODO + #set_cloud_fraction!(Y, p, moisture_model) + return nothing end diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl index 989f162395..2878395888 100644 --- a/src/cache/prognostic_edmf_precomputed_quantities.jl +++ b/src/cache/prognostic_edmf_precomputed_quantities.jl @@ -240,6 +240,8 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t) TD.air_temperature(thermo_params, ᶜts⁰), # t_sat TD.vapor_specific_humidity(thermo_params, ᶜts⁰), # qv_sat ᶜq_tot⁰, # qt_sat + TD.liquid_specific_humidity(thermo_params, ᶜts⁰), + TD.ice_specific_humidity(thermo_params, ᶜts⁰), TD.dry_pottemp(thermo_params, ᶜts⁰), # θ_sat TD.liquid_ice_pottemp(thermo_params, ᶜts⁰), # θ_liq_ice_sat projected_vector_data( diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index f2cc7c0b2b..5b055c5c3d 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -37,6 +37,12 @@ function flux_accumulation!(integrator) return nothing end +function cloud_fraction_model_callback!(integrator) + Y = integrator.u + p = integrator.p + set_cloud_fraction!(Y, p, p.atmos.moisture_model) +end + NVTX.@annotate function rrtmgp_model_callback!(integrator) Y = integrator.u p = integrator.p @@ -174,7 +180,7 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator) end function common_diagnostics(p, ᶜu, ᶜts) - (; env_thermo_quad, params) = p + (; params) = p thermo_params = CAP.thermodynamics_params(params) ᶜρ = TD.air_density.(thermo_params, ᶜts) diagnostics = (; @@ -186,6 +192,7 @@ function common_diagnostics(p, ᶜu, ᶜts) specific_enthalpy = TD.specific_enthalpy.(thermo_params, ᶜts), buoyancy = CAP.grav(p.params) .* (p.core.ᶜρ_ref .- ᶜρ) ./ ᶜρ, density = TD.air_density.(thermo_params, ᶜts), + cloud_fraction_gm = p.precomputed.ᶜcloud_fraction, ) if !(p.atmos.moisture_model isa DryModel) diagnostics = (; @@ -195,12 +202,6 @@ function common_diagnostics(p, ᶜu, ᶜts) q_ice = TD.ice_specific_humidity.(thermo_params, ᶜts), q_tot = TD.total_specific_humidity.(thermo_params, ᶜts), relative_humidity = TD.relative_humidity.(thermo_params, ᶜts), - cloud_fraction_gm = get_cloud_fraction.( - thermo_params, - env_thermo_quad, - p.precomputed.ᶜp, - ᶜts, - ), ) end return diagnostics @@ -211,13 +212,15 @@ end add_prefix(diagnostics::NamedTuple{names}, prefix) where {names} = NamedTuple{Symbol.(prefix, names)}(values(diagnostics)) -cloud_fraction(thermo_params, ts, area::FT) where {FT} = +# TODO - temporary crutch untill we move this to the new diagnostics +# and start using quadratures in EDMFX +draft_cloud_fraction(thermo_params, ts, area::FT) where {FT} = TD.has_condensate(thermo_params, ts) && area > 1e-3 ? FT(1) : FT(0) NVTX.@annotate function compute_diagnostics(integrator) (; t, u, p) = integrator Y = u - (; params, env_thermo_quad) = p + (; params) = p FT = eltype(params) thermo_params = CAP.thermodynamics_params(params) @@ -270,58 +273,40 @@ NVTX.@annotate function compute_diagnostics(integrator) env_diagnostics = (; common_diagnostics(p, ᶜu⁰, ᶜts⁰)..., area = ᶜa⁰, - cloud_fraction = get_cloud_fraction.( - thermo_params, - env_thermo_quad, - ᶜp, - ᶜts⁰, - ), + cloud_fraction = draft_cloud_fraction.(thermo_params, ᶜts⁰, ᶜa⁰), tke = ᶜtke⁰, mixing_length = ᶜmixing_length, ) draft_diagnostics = (; common_diagnostics(p, ᶜu⁺, ᶜts⁺)..., area = ᶜa⁺, - cloud_fraction = cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺), + cloud_fraction = draft_cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺), ) turbulence_convection_diagnostic = (; add_prefix(env_diagnostics, :env_)..., add_prefix(draft_diagnostics, :draft_)..., cloud_fraction = ᶜa⁰ .* - get_cloud_fraction.( - thermo_params, - env_thermo_quad, - ᶜp, - ᶜts⁰, - ) .+ ᶜa⁺ .* cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺), + draft_cloud_fraction.(thermo_params, ᶜts⁰, ᶜa⁰) .+ ᶜa⁺ .* + draft_cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺), ) elseif p.atmos.turbconv_model isa DiagnosticEDMFX (; ᶜtke⁰, ᶜmixing_length) = p.precomputed (; ᶜu⁺, ᶜts⁺, ᶜa⁺) = output_diagnostic_sgs_quantities(Y, p, t) env_diagnostics = (; - cloud_fraction = get_cloud_fraction.( - thermo_params, - env_thermo_quad, - ᶜp, - ᶜts, - ), + cloud_fraction = draft_cloud_fraction.(thermo_params, ᶜts, FT(1)), tke = ᶜtke⁰, mixing_length = ᶜmixing_length, ) draft_diagnostics = (; common_diagnostics(p, ᶜu⁺, ᶜts⁺)..., area = ᶜa⁺, - cloud_fraction = cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺), + cloud_fraction = draft_cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺), ) turbulence_convection_diagnostic = (; add_prefix(env_diagnostics, :env_)..., add_prefix(draft_diagnostics, :draft_)..., - cloud_fraction = get_cloud_fraction.( - thermo_params, - env_thermo_quad, - ᶜp, - ᶜts, - ) .+ ᶜa⁺ .* cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺), + cloud_fraction = draft_cloud_fraction.(thermo_params, ᶜts, FT(1)) .+ ᶜa⁺ .* + draft_cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺), ) else turbulence_convection_diagnostic = NamedTuple() diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl index 007b74da34..bdf7b9ce60 100644 --- a/src/diagnostics/core_diagnostics.jl +++ b/src/diagnostics/core_diagnostics.jl @@ -210,6 +210,21 @@ add_diagnostic_variable!( end, ) +### +# Cloud fraction (3d) +### +add_diagnostic_variable!( + short_name = "cl", + long_name = "Cloud fraction", + units = "%", + compute! = (out, state, cache, time) -> begin + if isnothing(out) + return copy(cache.precomputed.ᶜcloud_fraction) .* 100 + else + out .= cache.precomputed.ᶜcloud_fraction .* 100 + end + end, +) ### # Relative humidity (3d) diff --git a/src/prognostic_equations/buoyancy_gradients.jl b/src/prognostic_equations/buoyancy_gradients.jl index b8c4b196f4..e99cf44229 100644 --- a/src/prognostic_equations/buoyancy_gradients.jl +++ b/src/prognostic_equations/buoyancy_gradients.jl @@ -38,12 +38,34 @@ function buoyancy_gradients( bg_model.qt_sat, ) elseif moisture_model isa NonEquilMoistModel + TD.PhaseNonEquil_pθq( + thermo_params, + bg_model.p, + bg_model.θ_liq_ice_sat, + TD.PhasePartition( + bg_model.qt_sat, + bg_model.ql_sat, + bg_model.qi_sat, + ), + ) + else error("Unsupported moisture model") end - phase_part = TD.PhasePartition(thermo_params, ts_sat) + phase_part = + if moisture_model isa DryModel || moisture_model isa EquilMoistModel + TD.PhasePartition(thermo_params, ts_sat) + elseif moisture_model isa NonEquilMoistModel + TD.PhasePartition( + bg_model.qt_sat, + bg_model.ql_sat, + bg_model.qi_sat, + ) + end + lh = TD.latent_heat_liq_ice(thermo_params, phase_part) cp_m = TD.cp_m(thermo_params, ts_sat) + # TODO - double check if this is assuming liquid only? ∂b∂θl_sat = ( ∂b∂θv * ( 1 + diff --git a/src/prognostic_equations/cloud_fraction.jl b/src/prognostic_equations/cloud_fraction.jl deleted file mode 100644 index 6ed6b91aa2..0000000000 --- a/src/prognostic_equations/cloud_fraction.jl +++ /dev/null @@ -1,113 +0,0 @@ -import StaticArrays as SA - -function get_cloud_fraction(thermo_params, env_thermo_quad, ᶜp, ᶜts) - q_tot = TD.total_specific_humidity(thermo_params, ᶜts) - FT = eltype(thermo_params) - θ_liq_ice = TD.liquid_ice_pottemp(thermo_params, ᶜts) - qt′qt′ = (FT(0.05) * q_tot)^2 - θl′θl′ = FT(5) - θl′qt′ = FT(0) - return compute_cloud_fraction( - env_thermo_quad, - thermo_params, - ᶜp, - q_tot, - θ_liq_ice, - qt′qt′, - θl′θl′, - θl′qt′, - ) -end - -function compute_cloud_fraction( - env_thermo_quad, - thermo_params, - ᶜp::FT, - q_tot, - θ_liq_ice, - qt′qt′, - θl′θl′, - θl′qt′, -) where {FT} - vars = (; - p_c = ᶜp, - qt_mean = q_tot, - θl_mean = θ_liq_ice, - qt′qt′, - θl′θl′, - θl′qt′, - ) - return quad_loop(env_thermo_quad, vars, thermo_params)::FT -end - -function quad_loop(env_thermo_quad::SGSQuadrature, vars, thermo_params) - - # qt - total water specific humidity - # θl - liquid ice potential temperature - # _mean and ′ - subdomain mean and (co)variances - (; qt′qt′, qt_mean, θl′θl′, θl_mean, θl′qt′, p_c) = vars - - FT = eltype(qt_mean) - - sqrt2 = FT(sqrt(2)) - - # Epsilon defined per typical variable fluctuation - eps_q = FT(eps(FT)) * max(FT(eps(FT)), qt_mean) - eps_θ = FT(eps(FT)) - - # limit σ_q to prevent negative q_tot_hat - σ_q_lim = -qt_mean / (sqrt2 * env_thermo_quad.a[1]) # TODO: is this correct? - σ_q::FT = min(sqrt(qt′qt′), σ_q_lim) - σ_θ::FT = sqrt(θl′θl′) - - # Enforce Cauchy-Schwarz inequality, numerically stable compute - _corr::FT = (θl′qt′ / max(σ_q, eps_q)) - corr::FT = max(min(_corr / max(σ_θ, eps_θ), 1), -1) - - # Conditionals - σ_c = sqrt(max(1 - corr * corr, 0)) * σ_θ - - function get_x_hat(χ::Tuple{<:Real, <:Real}) - μ_c = θl_mean + sqrt2 * corr * σ_θ * χ[1] - θ_hat = μ_c + sqrt2 * σ_c * χ[2] - q_tot_hat = qt_mean + sqrt2 * σ_q * χ[1] - return (θ_hat, q_tot_hat) - end - - # cloudy/dry categories for buoyancy in TKE - f_q_tot_sat(x_hat::Tuple{<:Real, <:Real}, hc) = - hc ? x_hat[2] : eltype(x_hat)(0) - - get_ts(x_hat::Tuple{<:Real, <:Real}) = - thermo_state(thermo_params; p = p_c, θ = x_hat[1], q_tot = x_hat[2]) - f_cf(x_hat::Tuple{<:Real, <:Real}, hc) = - hc ? eltype(x_hat)(1) : eltype(x_hat)(0) - function f(x_hat::Tuple{<:Real, <:Real}) - ts = get_ts(x_hat) - hc = TD.has_condensate(thermo_params, ts) - return (; cf = f_cf(x_hat, hc), q_tot_sat = f_q_tot_sat(x_hat, hc)) - end - - return quad(f, get_x_hat, env_thermo_quad).cf -end - -import ClimaCore.RecursiveApply: rzero, ⊞, ⊠ -function quad(f, get_x_hat::F, quad_type) where {F <: Function} - χ = quad_type.a - weights = quad_type.w - quad_order = quadrature_order(quad_type) - FT = eltype(χ) - # zero outer quadrature points - T = typeof(f(get_x_hat((χ[1], χ[1])))) - outer_env = rzero(T) - @inbounds for m_q in 1:quad_order - # zero inner quadrature points - inner_env = rzero(T) - for m_h in 1:quad_order - x_hat = get_x_hat((χ[m_q], χ[m_h])) - inner_env = inner_env ⊞ f(x_hat) ⊠ weights[m_h] ⊠ FT(1 / sqrt(π)) - end - outer_env = outer_env ⊞ inner_env ⊠ weights[m_q] ⊠ FT(1 / sqrt(π)) - end - return outer_env -end diff --git a/src/prognostic_equations/edmfx_closures.jl b/src/prognostic_equations/edmfx_closures.jl index 7e478bb45f..09ad08a40f 100644 --- a/src/prognostic_equations/edmfx_closures.jl +++ b/src/prognostic_equations/edmfx_closures.jl @@ -248,19 +248,6 @@ function mixing_length( l_N = l_z end - # compute l_smag - the Smagorinsky length scale. - # TODO: This should be added to ClimaParameters - c_smag = CAP.c_smag(params) - N_eff = sqrt(max(ᶜlinear_buoygrad, 0)) - if N_eff > 0.0 - l_smag = - c_smag * - ᶜdz * - max(0, 1 - N_eff^2 / ᶜPr / (2 * ᶜstrain_rate_norm))^(1 / 4) - else - l_smag = c_smag * ᶜdz - end - # add limiters l = SA.SVector( (l_N < eps(FT) || l_N > l_z) ? l_z : l_N, diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index aac1f9866b..7cf2aff14e 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -512,6 +512,10 @@ function get_callbacks(parsed_args, sim_info, atmos, params, comms_ctx) (callbacks..., call_every_dt(rrtmgp_model_callback!, dt_rad)) end + dt_cf = FT(time_to_seconds(parsed_args["dt_cloud_fraction"])) + callbacks = + (callbacks..., call_every_dt(cloud_fraction_model_callback!, dt_cf)) + return callbacks end diff --git a/src/solver/types.jl b/src/solver/types.jl index 5933a18516..e1b8cc8afc 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -121,6 +121,10 @@ Base.@kwdef struct EnvBuoyGrad{FT, EBC <: AbstractEnvBuoyGradClosure} qv_sat::FT "total specific humidity in the saturated part" qt_sat::FT + "liquid specific humidity in the saturated part" + ql_sat::FT + "ice specific humidity in the saturated part" + qi_sat::FT "potential temperature in the saturated part" θ_sat::FT "liquid ice potential temperature in the saturated part" @@ -190,9 +194,9 @@ abstract type AbstractQuadratureType end struct LogNormalQuad <: AbstractQuadratureType end struct GaussianQuad <: AbstractQuadratureType end -abstract type AbstractEnvThermo end -struct SGSMean <: AbstractEnvThermo end -struct SGSQuadrature{N, QT, A, W} <: AbstractEnvThermo +abstract type AbstractSGSamplingType end +struct SGSMean <: AbstractSGSamplingType end +struct SGSQuadrature{N, QT, A, W} <: AbstractSGSamplingType quadrature_type::QT a::A w::W @@ -218,7 +222,7 @@ struct SGSQuadrature{N, QT, A, W} <: AbstractEnvThermo end end quadrature_order(::SGSQuadrature{N}) where {N} = N -quad_type(::SGSQuadrature{N}) where {N} = N +quad_type(::SGSQuadrature{N}) where {N} = N #TODO - this seems wrong? abstract type AbstractSurfaceThermoState end struct GCMSurfaceThermoState <: AbstractSurfaceThermoState end @@ -232,7 +236,7 @@ Base.broadcastable(x::PrognosticEDMFX) = tuple(x) Base.broadcastable(x::DiagnosticEDMFX) = tuple(x) Base.broadcastable(x::AbstractEntrainmentModel) = tuple(x) Base.broadcastable(x::AbstractDetrainmentModel) = tuple(x) -Base.broadcastable(x::AbstractEnvThermo) = tuple(x) +Base.broadcastable(x::AbstractSGSamplingType) = tuple(x) Base.@kwdef struct RadiationDYCOMS_RF01{FT} "Large-scale divergence" diff --git a/test/coupler_compatibility.jl b/test/coupler_compatibility.jl index 4cd9f977b5..9918a62b8d 100644 --- a/test/coupler_compatibility.jl +++ b/test/coupler_compatibility.jl @@ -65,7 +65,7 @@ const T2 = 290 p.core, sfc_setup, p.ghost_buffer, - p.env_thermo_quad, + p.SG_quad, p.precomputed, p.scratch, p.hyperdiff, From 837e2285f13d11a562bfd4471e28ce127ecc0db4 Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Wed, 13 Dec 2023 16:28:17 -0800 Subject: [PATCH 2/4] Add mixing length to diagnostics --- ...uarez_rhoe_equilmoist_topography_dcmip.yml | 15 +++++ src/cache/cloud_fraction.jl | 59 +++++++++++-------- src/diagnostics/Diagnostics.jl | 1 + src/diagnostics/core_diagnostics.jl | 21 ++++++- src/diagnostics/netcdf_writer.jl | 1 + 5 files changed, 71 insertions(+), 26 deletions(-) diff --git a/config/model_configs/sphere_held_suarez_rhoe_equilmoist_topography_dcmip.yml b/config/model_configs/sphere_held_suarez_rhoe_equilmoist_topography_dcmip.yml index d41ed96e0b..793d9dfb9a 100644 --- a/config/model_configs/sphere_held_suarez_rhoe_equilmoist_topography_dcmip.yml +++ b/config/model_configs/sphere_held_suarez_rhoe_equilmoist_topography_dcmip.yml @@ -11,3 +11,18 @@ topography: "DCMIP200" job_id: "sphere_held_suarez_rhoe_equilmoist_topography_dcmip" moist: "equil" netcdf_interpolate_z_over_msl: true +diagnostics: + - short_name: cl + period: 3hours + - short_name: cli + period: 3hours + - short_name: clw + period: 3hours + - short_name: ta + period: 3hours + - short_name: hus + period: 3hours + - short_name: hur + period: 3hours + - short_name: mixlgm + period: 3hours diff --git a/src/cache/cloud_fraction.jl b/src/cache/cloud_fraction.jl index 79d9e83123..3cb3fa2004 100644 --- a/src/cache/cloud_fraction.jl +++ b/src/cache/cloud_fraction.jl @@ -18,28 +18,23 @@ end - Pr - Prandtl number - ϵ_st - strain rate norm """ -function compute_smagorinsky_length_scale(c_smag, N_eff, dz, Pr, ϵ_st) +function smagorinsky_length_scale(c_smag, N_eff, dz, Pr, ϵ_st) FT = eltype(c_smag) - #return N_eff > FT(0) && N_eff < sqrt(2 * Pr * ϵ_st) ? return N_eff > FT(0) ? - c_smag * dz * max(0, (1 - N_eff^2 / Pr / 2 / ϵ_st))^(1 / 4) : + c_smag * + dz * + max(0, 1 - N_eff^2 / Pr / 2 / max(ϵ_st, eps(FT)))^(1 / 4) : c_smag * dz end -""" - Compute the grid scale cloud fraction based on sub-grid scale properties -""" -function set_cloud_fraction!(Y, p, ::DryModel) - @. p.precomputed.ᶜcloud_fraction = 0 -end -function set_cloud_fraction!(Y, p, ::Union{EquilMoistModel, NonEquilMoistModel}) - (; SG_quad, params) = p +function compute_gm_mixing_length!(ᶜmixing_length, Y, p) + (; params) = p + thermo_params = CAP.thermodynamics_params(params) FT = eltype(params) - thermo_params = CAP.thermodynamics_params(params) ᶜdz = Fields.Δz_field(axes(Y.c)) ᶜlg = Fields.local_geometry_field(Y.c) - (; ᶜts, ᶜp, ᶠu³, ᶜcloud_fraction) = p.precomputed + (; ᶜts, ᶜp, ᶠu³) = p.precomputed (; obukhov_length) = p.precomputed.sfc_conditions ᶜlinear_buoygrad = p.scratch.ᶜtemp_scalar @@ -78,26 +73,40 @@ function set_cloud_fraction!(Y, p, ::Union{EquilMoistModel, NonEquilMoistModel}) ᶠu = p.scratch.ᶠtemp_C123 @. ᶠu = C123(ᶠinterp(Y.c.uₕ)) + C123(ᶠu³) - ᶜstrain_rate = p.scratch.ᶜtemp_UVWxUVW compute_strain_rate_center!(ᶜstrain_rate, ᶠu) - ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar - @. ᶜprandtl_nvec = turbulent_prandtl_number( - params, - obukhov_length, - ᶜlinear_buoygrad, - norm_sqr(ᶜstrain_rate), - ) - - ᶜl_smag = p.scratch.ᶜtemp_scalar_2 - @. ᶜl_smag = compute_smagorinsky_length_scale( + @. ᶜmixing_length = smagorinsky_length_scale( CAP.c_smag(params), sqrt(max(ᶜlinear_buoygrad, 0)), #N_eff ᶜdz, - ᶜprandtl_nvec, + turbulent_prandtl_number( + params, + obukhov_length, + ᶜlinear_buoygrad, + norm_sqr(ᶜstrain_rate), + ), norm_sqr(ᶜstrain_rate), ) +end + +""" + Compute the grid scale cloud fraction based on sub-grid scale properties +""" +function set_cloud_fraction!(Y, p, ::DryModel) + @. p.precomputed.ᶜcloud_fraction = 0 +end +function set_cloud_fraction!(Y, p, ::Union{EquilMoistModel, NonEquilMoistModel}) + (; SG_quad, params) = p + + FT = eltype(params) + thermo_params = CAP.thermodynamics_params(params) + (; ᶜts, ᶜp, ᶜcloud_fraction) = p.precomputed + + # temp_scalar is already in use inside the compute_gm_mixing_length + # this is not a great pattern, but I don't have a more elegant idea + ᶜl_smag = p.scratch.ᶜtemp_scalar_2 + compute_gm_mixing_length!(ᶜl_smag, Y, p) coeff = FT(2.1) # TODO - move to parameters @. ᶜcloud_fraction = quad_loop( diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl index 562778626f..51d3669608 100644 --- a/src/diagnostics/Diagnostics.jl +++ b/src/diagnostics/Diagnostics.jl @@ -35,6 +35,7 @@ import ..DiagnosticEDMFX # functions used to calculate diagnostics import ..draft_area +import ..compute_gm_mixing_length! # We need the abbreviations for symbols like curl, grad, and so on include(joinpath("..", "utils", "abbreviations.jl")) diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl index bdf7b9ce60..1add87444a 100644 --- a/src/diagnostics/core_diagnostics.jl +++ b/src/diagnostics/core_diagnostics.jl @@ -4,7 +4,7 @@ # # In addition to the metadata (names, comments, ...), the most important step in adding a # new DiagnosticVariable is defining its compute! function. `compute!` has to take four -# arguments: (out, state, cache, time), and as to write the diagnostic in place into the +# arguments: (out, state, cache, time), and has to write the diagnostic in place into the # `out` variable. # # Often, it is possible to compute certain diagnostics only for specific models (e.g., @@ -226,6 +226,25 @@ add_diagnostic_variable!( end, ) +### +# Grid mean mixing length (3d) +### +function compute_mixlgm!(out, state, cache, time) + if isnothing(out) + out = similar(state.c.ρ) + end + compute_gm_mixing_length!(out, state, cache) +end + +add_diagnostic_variable!( + short_name = "mixlgm", + long_name = "Grid mean mixing length", + standard_name = "grid_mean_mixing_length", + units = "m", + comments = "Smagorinsky length scale, to be used as an approximation of mixing length in simulations without EDMF SGS model", + compute! = compute_mixlgm!, +) + ### # Relative humidity (3d) ### diff --git a/src/diagnostics/netcdf_writer.jl b/src/diagnostics/netcdf_writer.jl index 6c9592ca99..3b22375df4 100644 --- a/src/diagnostics/netcdf_writer.jl +++ b/src/diagnostics/netcdf_writer.jl @@ -597,6 +597,7 @@ function write_field!( deflatelevel = writer.compression_level, ) v.attrib["long_name"] = diagnostic.output_long_name + v.attrib["short_name"] = var.short_name v.attrib["units"] = var.units v.attrib["comments"] = var.comments From 45319a496ebb907266713f7d7f7649c34eaf4954 Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Wed, 13 Dec 2023 17:40:16 -0800 Subject: [PATCH 3/4] Bump allocation limits in flame test --- perf/flame.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/perf/flame.jl b/perf/flame.jl index d62e064b28..f75924c7d9 100644 --- a/perf/flame.jl +++ b/perf/flame.jl @@ -36,17 +36,17 @@ ProfileCanvas.html_file(joinpath(output_dir, "flame.html"), results) ##### allocs_limit = Dict() -allocs_limit["flame_perf_target"] = 145_184 -allocs_limit["flame_perf_target_tracers"] = 177_440 +allocs_limit["flame_perf_target"] = 147_520 +allocs_limit["flame_perf_target_tracers"] = 179_776 allocs_limit["flame_perf_target_edmfx"] = 7_005_552 allocs_limit["flame_perf_diagnostics"] = 25_356_928 -allocs_limit["flame_perf_target_diagnostic_edmfx"] = 1_309_104 +allocs_limit["flame_perf_target_diagnostic_edmfx"] = 1_309_968 allocs_limit["flame_sphere_baroclinic_wave_rhoe_equilmoist_expvdiff"] = 4_018_252_656 allocs_limit["flame_perf_target_threaded"] = 1_276_864 allocs_limit["flame_perf_target_callbacks"] = 37_277_112 -allocs_limit["flame_perf_gw"] = 3_226_427_872 -allocs_limit["flame_perf_target_prognostic_edmfx_aquaplanet"] = 1_251_728 +allocs_limit["flame_perf_gw"] = 3_226_428_736 +allocs_limit["flame_perf_target_prognostic_edmfx_aquaplanet"] = 1_257_712 # Ideally, we would like to track all the allocations, but this becomes too # expensive there is too many of them. Here, we set the default sample rate to From 72b09d4d7cde9b4700e00bbdfb8c68d97b2a0ec5 Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Wed, 13 Dec 2023 16:45:48 -0800 Subject: [PATCH 4/4] Bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index c144ed44f6..0ccdcfeb9e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ClimaAtmos" uuid = "b2c96348-7fb7-4fe0-8da9-78d88439e717" authors = ["Climate Modeling Alliance"] -version = "0.17.2" +version = "0.18.0" [deps] ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"