Add cloud fraction diagnostics

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"

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:

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"

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

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

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

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

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,

src/cache/cloud_fraction.jl

@@ -0,0 +1,234 @@
+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 smagorinsky_length_scale(c_smag, N_eff, dz, Pr, ϵ_st)
+    FT = eltype(c_smag)
+    return N_eff > FT(0) ?
+           c_smag *
+           dz *
+           max(0, 1 - N_eff^2 / Pr / 2 / max(ϵ_st, eps(FT)))^(1 / 4) :
+           c_smag * dz
+end
+
+function compute_gm_mixing_length!(ᶜmixing_length, Y, p)
+    (; params) = p
+    thermo_params = CAP.thermodynamics_params(params)
+
+    FT = eltype(params)
+    ᶜdz = Fields.Δz_field(axes(Y.c))
+    ᶜlg = Fields.local_geometry_field(Y.c)
+    (; ᶜts, ᶜp, ᶠu³) = 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)
+
+    @. ᶜmixing_length = smagorinsky_length_scale(
+        CAP.c_smag(params),
+        sqrt(max(ᶜlinear_buoygrad, 0)),   #N_eff
+        ᶜdz,
+        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(
+        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

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(