diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl index 7e289d00416..82875ca5fb2 100644 --- a/examples/hybrid/driver.jl +++ b/examples/hybrid/driver.jl @@ -174,5 +174,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 2e62dfefcee..00b99f99d8f 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 b5b84cef6c1..1863aa488f8 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 @@ -146,11 +147,11 @@ function build_cache(Y, atmos, params, surface_setup, simulation) 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, simulation) + (; atmos, core, params, sfc_setup, precomputed, scratch, simulation, SG_quad) # Coupler compatibility isnothing(precomputing_arguments.sfc_setup) && @@ -183,7 +184,7 @@ function build_cache(Y, atmos, params, surface_setup, simulation) core, sfc_setup, ghost_buffer, - env_thermo_quad, + SG_quad, precomputed, scratch, hyperdiff, diff --git a/src/prognostic_equations/cloud_fraction.jl b/src/cache/cloud_fraction.jl similarity index 52% rename from src/prognostic_equations/cloud_fraction.jl rename to src/cache/cloud_fraction.jl index 6ed6b91aa29..3cd75dfd36b 100644 --- a/src/prognostic_equations/cloud_fraction.jl +++ b/src/cache/cloud_fraction.jl @@ -1,52 +1,72 @@ import StaticArrays as SA +import ClimaCore.RecursiveApply: rzero, ⊞, ⊠ + +# TODO: write a test with scalars that are linear with z +function covariance_from_grad(coeff, ᶜmixing_length, ᶜ∇Φ, ᶜ∇Ψ) + return 2 * coeff * ᶜmixing_length^2 * dot(ᶜ∇Φ, ᶜ∇Ψ) +end + +function set_cloud_fraction!(Y, p, ::DryModel) + FT = eltype(Y) + (; cloud_fraction) = p.precomputed + @. cloud_fraction = FT(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)) + (; ᶜts, ᶜp) = p.precomputed + + (; ᶜqₜ′qₜ′, ᶜθ′θ′, ᶜθ′qₜ′) = p.precomputed + (; cloud_fraction) = p.precomputed + + ᶜqₜ = p.scratch.ᶜtemp_scalar + ᶜθ = p.scratch.ᶜtemp_scalar_2 + @. ᶜqₜ = TD.total_specific_humidity(thermo_params, ᶜts) + @. ᶜθ = TD.liquid_ice_pottemp(thermo_params, ᶜts) -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, + coeff = FT(2.1) # TODO - move to parameters + + # TODO - replace dz with mixinglength when using EDMF SGS + @. cloud_fraction = quad_loop( + SG_quad, ᶜp, - q_tot, - θ_liq_ice, - qt′qt′, - θl′θl′, - θl′qt′, + ᶜqₜ, + ᶜθ, + covariance_from_grad( + coeff, + ᶜdz, + Geometry.WVector(ᶜgradᵥ(ᶠinterp(ᶜqₜ))), + Geometry.WVector(ᶜgradᵥ(ᶠinterp(ᶜqₜ))), + ), + covariance_from_grad( + coeff, + ᶜdz, + Geometry.WVector(ᶜgradᵥ(ᶠinterp(ᶜθ))), + Geometry.WVector(ᶜgradᵥ(ᶠinterp(ᶜθ))), + ), + covariance_from_grad( + coeff, + ᶜdz, + Geometry.WVector(ᶜgradᵥ(ᶠinterp(ᶜθ))), + Geometry.WVector(ᶜgradᵥ(ᶠinterp(ᶜqₜ))) + ), + thermo_params ) end -function compute_cloud_fraction( - env_thermo_quad, - thermo_params, - ᶜp::FT, - q_tot, - θ_liq_ice, +function quad_loop( + SG_quad::SGSQuadrature, + p_c, + qt_mean, + θl_mean, 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 - + thermo_params +) FT = eltype(qt_mean) sqrt2 = FT(sqrt(2)) @@ -56,7 +76,7 @@ function quad_loop(env_thermo_quad::SGSQuadrature, vars, thermo_params) 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_lim = -qt_mean / (sqrt2 * SG_quad.a[1]) # TODO: is this correct? σ_q::FT = min(sqrt(qt′qt′), σ_q_lim) σ_θ::FT = sqrt(θl′θl′) @@ -88,10 +108,10 @@ function quad_loop(env_thermo_quad::SGSQuadrature, vars, thermo_params) 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 + return quad(f, get_x_hat, SG_quad).cf end -import ClimaCore.RecursiveApply: rzero, ⊞, ⊠ +# TODO f::F1 <: Functon function quad(f, get_x_hat::F, quad_type) where {F <: Function} χ = quad_type.a weights = quad_type.w diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index 2331af0afa8..1c302a244ff 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -54,6 +54,12 @@ function precomputed_quantities(Y, atmos) ᶜh_tot = similar(Y.c, FT), sfc_conditions = Fields.Field(SCT, Spaces.level(axes(Y.f), half)), ) + cloud_diagnostics = (; + ᶜqₜ′qₜ′ = similar(Y.c, FT), + ᶜθ′θ′ = similar(Y.c, FT), + ᶜθ′qₜ′ = similar(Y.c, FT), + cloud_fraction = similar(Y.c, FT), + ) advective_sgs_quantities = atmos.turbconv_model isa PrognosticEDMFX ? (; @@ -125,6 +131,7 @@ function precomputed_quantities(Y, atmos) diagnostic_sgs_quantities..., vert_diff_quantities..., precipitation_quantities..., + cloud_diagnostics..., ) end @@ -375,6 +382,8 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) set_precipitation_precomputed_quantities!(Y, p, t) end + set_cloud_fraction!(Y, p, moisture_model) + return nothing end diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 02e9078115a..5473639d0de 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -174,7 +174,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 +186,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 +196,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 +206,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 +267,41 @@ 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⁺), + cloud_fraction = + ᶜ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/solver/type_getters.jl b/src/solver/type_getters.jl index 5c2e19265d4..cf5ce261ed4 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -751,6 +751,7 @@ function get_integrator(config::AtmosConfig) s = @timed_str begin p = build_cache(Y, atmos, params, surface_setup, simulation) + @info(propertynames(p)) end @info "Allocating cache (p): $s" diff --git a/src/solver/types.jl b/src/solver/types.jl index 5933a18516a..4ae032d3498 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -190,9 +190,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 +218,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 +232,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 02c2b0adb31..3c43bd94158 100644 --- a/test/coupler_compatibility.jl +++ b/test/coupler_compatibility.jl @@ -64,7 +64,7 @@ const T2 = 290 p.core, sfc_setup, p.ghost_buffer, - p.env_thermo_quad, + p.SG_quad, p.precomputed, p.scratch, p.hyperdiff,