wip on cloud fraction diagnostics [ci skip]

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"

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"),
 )

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
 

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")

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(

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

src/prognostic_equations/cloud_fraction.jl → 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