From 8c7d298f4032394756c8e63444f7a35c98f04b57 Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Thu, 30 Nov 2023 13:42:21 -0800 Subject: [PATCH] wip on cloud fraction diagnostics [ci skip] --- ...uarez_rhoe_equilmoist_topography_dcmip.yml | 1 + src/ClimaAtmos.jl | 2 +- src/cache/precomputed_quantities.jl | 11 ++++ src/solver/model_getters.jl | 11 ++++ src/solver/type_getters.jl | 2 + src/solver/types.jl | 5 ++ .../cloud_fraction.jl | 60 ++++++++++--------- 7 files changed, 62 insertions(+), 30 deletions(-) rename src/{prognostic_equations => utils}/cloud_fraction.jl (65%) 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 9f0cc86e451..c99f303fd6e 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 @@ -7,6 +7,7 @@ vert_diff: "true" use_reference_state: false forcing: "held_suarez" precip_model: "0M" +cloud_diagnostics: "diag_covariance" topography: "DCMIP200" job_id: "sphere_held_suarez_rhoe_equilmoist_topography_dcmip" moist: "equil" diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl index 2e62dfefcee..4b52932873e 100644 --- a/src/ClimaAtmos.jl +++ b/src/ClimaAtmos.jl @@ -13,6 +13,7 @@ include(joinpath("utils", "utilities.jl")) include(joinpath("utils", "debug_utils.jl")) include(joinpath("topography", "topography.jl")) include(joinpath("utils", "variable_manipulations.jl")) +include(joinpath("utils", "cloud_fraction.jl")) include( joinpath("parameterized_tendencies", "radiation", "radiation_utilities.jl"), @@ -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/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index 876bac8bfc8..737f10a7dba 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -112,11 +112,22 @@ function precomputed_quantities(Y, atmos) precipitation_quantities = atmos.precip_model isa Microphysics1Moment ? (; ᶜwᵣ = similar(Y.c, FT), ᶜwₛ = similar(Y.c, FT)) : (;) + cloud_diagnostics_quantities = + atmos.cloud_diagnostics isa DiagnosticCovariance ? + (; + ᶜqₜ′qₜ′ = similar(Y.c, FT), + ᶜθ′θ′ = similar(Y.c, FT), + ᶜθ′qₜ′ = similar(Y.c, FT), + cloud_fraction = similar(Y.c, FT) + ᶜqₜ_grad = (Y.f, Geometry.WVector{FT}), + ᶜθ_grad = (Y.f, Geometry.WVector{FT}), + ) : (;) return (; gs_quantities..., advective_sgs_quantities..., diagnostic_sgs_quantities..., precipitation_quantities..., + cloud_diagnostics_quantities..., ) end diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl index ea76d26efef..50dc3cbf64f 100644 --- a/src/solver/model_getters.jl +++ b/src/solver/model_getters.jl @@ -221,6 +221,17 @@ function get_precipitation_model(parsed_args) end end +function get_cloud_diagnostics_type(parsed_args) + cloud_diag = parsed_args["cloud_diagnostics"] + return if cloud_diag == nothing + nothing + elseif cloud_diag = "diag_covariance" + DiagnosticCovariance() + else + error("Invalid cloud diagnostics type $(cloud_diagnostics)") + end +end + function get_forcing_type(parsed_args) forcing = parsed_args["forcing"] @assert forcing in (nothing, "held_suarez") diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index 5c2e19265d4..273207e34a6 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -17,6 +17,7 @@ function get_atmos(config::AtmosConfig, params) FT = eltype(config) moisture_model = get_moisture_model(parsed_args) precip_model = get_precipitation_model(parsed_args) + cloud_diagnostics_type = get_cloud_diagnostics_type(parsed_args) radiation_mode = get_radiation_mode(parsed_args, FT) forcing_type = get_forcing_type(parsed_args) @@ -58,6 +59,7 @@ function get_atmos(config::AtmosConfig, params) edmfx_sgs_diffusive_flux, edmfx_nh_pressure, precip_model, + cloud_diagnostics, forcing_type, turbconv_model = get_turbconv_model(FT, parsed_args, turbconv_params), non_orographic_gravity_wave = get_non_orographic_gravity_wave_model( diff --git a/src/solver/types.jl b/src/solver/types.jl index 5933a18516a..8d6606c42cf 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -11,6 +11,9 @@ struct NoPrecipitation <: AbstractPrecipitationModel end struct Microphysics0Moment <: AbstractPrecipitationModel end struct Microphysics1Moment <: AbstractPrecipitationModel end +abstract type AbstractCloudDiagnosticsType end +struct DiagnosticCovariance <: AbstractCloudDiagnosticsType end + abstract type AbstractModelConfig end struct SingleColumnModel <: AbstractModelConfig end struct SphericalModel <: AbstractModelConfig end @@ -325,6 +328,7 @@ Base.@kwdef struct AtmosModel{ PEM, MM, PM, + CD, F, S, RM, @@ -352,6 +356,7 @@ Base.@kwdef struct AtmosModel{ perf_mode::PEM = nothing moisture_model::MM = nothing precip_model::PM = nothing + cloud_diagnostics::CD = nothing forcing_type::F = nothing subsidence::S = nothing radiation_mode::RM = nothing diff --git a/src/prognostic_equations/cloud_fraction.jl b/src/utils/cloud_fraction.jl similarity index 65% rename from src/prognostic_equations/cloud_fraction.jl rename to src/utils/cloud_fraction.jl index 6ed6b91aa29..64212293529 100644 --- a/src/prognostic_equations/cloud_fraction.jl +++ b/src/utils/cloud_fraction.jl @@ -1,34 +1,36 @@ import StaticArrays as SA +import ClimaCore.RecursiveApply: rzero, ⊞, ⊠ -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′, - ) +# TODO: write a test with scalars that are linear with z +function get_covariance(coeff, ᶜmixing_length, ᶜscalar1_grad, ᶜscalar2_grad) + return 2 * coeff * ᶜmixing_length^2 * dot(ᶜscalar1_grad, ᶜscalar2_grad) end -function compute_cloud_fraction( - env_thermo_quad, - thermo_params, - ᶜp::FT, - q_tot, - θ_liq_ice, - qt′qt′, - θl′θl′, - θl′qt′, -) where {FT} +function get_cloud_fraction(Y, p, colidx, ::SGSQuadrature) + (; params) = p + FT = eltype(params) + thermo_params = CAP.thermodynamics_params(params) + ᶜdz = Fields.Δz_field(axes(Y.c)) + (; ᶜts, ᶜp) = p.precomputed + + (; qt′qt′, θl′θl′, θl′qt′) = p.precomputed.cloud_diagnostics_quantities + (; cloud_fraction) = p.precomputed.cloud_diagnostics_quantities + (; ᶜqₜ_grad, ᶜθ_grad) = p.precomputed.cloud_diagnostics_quantities + + ᶜ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) + + @. ᶜqₜ_grad = Geometry.WVector(ᶜgradᵥ(ᶠinterp(ᶜqₜ))) + @. ᶜθ_grad = Geometry.WVector(ᶜgradᵥ(ᶠinterp(ᶜθ))) + + coeff = FT(2.1) # TODO - move to parameters + # TODO - replace dz with mixinglength when using EDMF SGS + @. qt′qt′ = get_covariance(coeff, ᶜdz, ᶜqₜ_grad, ᶜqₜ_grad) + @. θl′θl′ = get_covariance(coeff, ᶜdz, ᶜθ_grad, ᶜθ_grad) + @. θl′qt′ = get_covariance(coeff, ᶜdz, ᶜθ_grad, ᶜqₜ_grad) + vars = (; p_c = ᶜp, qt_mean = q_tot, @@ -37,7 +39,8 @@ function compute_cloud_fraction( θl′θl′, θl′qt′, ) - return quad_loop(env_thermo_quad, vars, thermo_params)::FT + + @. cloud_fraction = quad_loop(::SGSQuadrature, vars, thermo_params)::FT end function quad_loop(env_thermo_quad::SGSQuadrature, vars, thermo_params) @@ -91,7 +94,6 @@ function quad_loop(env_thermo_quad::SGSQuadrature, vars, thermo_params) 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