Skip to content

Commit

Permalink
Add cloud fraction diagnostics
Browse files Browse the repository at this point in the history
  • Loading branch information
trontrytel committed Dec 2, 2023
1 parent 4d36738 commit 21d0851
Showing 9 changed files with 107 additions and 91 deletions.
5 changes: 5 additions & 0 deletions examples/hybrid/driver.jl
Original file line number Diff line number Diff line change
@@ -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
2 changes: 1 addition & 1 deletion src/ClimaAtmos.jl
Original file line number Diff line number Diff line change
@@ -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"),
)
11 changes: 6 additions & 5 deletions src/cache/cache.jl
Original file line number Diff line number Diff line change
@@ -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,
Original file line number Diff line number Diff line change
@@ -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
9 changes: 9 additions & 0 deletions src/cache/precomputed_quantities.jl
Original file line number Diff line number Diff line change
@@ -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

52 changes: 16 additions & 36 deletions src/callbacks/callbacks.jl
Original file line number Diff line number Diff line change
@@ -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()
1 change: 1 addition & 0 deletions src/solver/type_getters.jl
Original file line number Diff line number Diff line change
@@ -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"

10 changes: 5 additions & 5 deletions src/solver/types.jl
Original file line number Diff line number Diff line change
@@ -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"
2 changes: 1 addition & 1 deletion test/coupler_compatibility.jl
Original file line number Diff line number Diff line change
@@ -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,

0 comments on commit 21d0851

Please sign in to comment.