-
Notifications
You must be signed in to change notification settings - Fork 19
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
1 parent
4d36738
commit 84a1573
Showing
9 changed files
with
197 additions
and
161 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,152 @@ | ||
| import StaticArrays as SA | ||
| import ClimaCore.RecursiveApply: rzero, ⊞, ⊠ | ||
|
|
||
| # TODO: write a test with scalars that are linear with z | ||
| """ | ||
| Diagnose horizntal 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 grid scale cloud fraction based on sub-grid scale properties | ||
| """ | ||
| 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, 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) | ||
|
|
||
| coeff = FT(2.1) # TODO - move to parameters | ||
|
|
||
| @. cloud_fraction = quad_loop( | ||
| SG_quad, | ||
| ᶜp, | ||
| ᶜqₜ, | ||
| ᶜθ, | ||
| covariance_from_grad( | ||
| coeff, | ||
| ᶜdz, # TODO - replace dz with mixinglength when using EDMF SGS | ||
| 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 quad_loop(SG_quad, p_c, qt_mean, θl_mean, qt′qt′, θl′θl′, θl′qt′, thermo_params) | ||
| where: | ||
| - SG_quad is a struct containing information about quadrature type and order | ||
| - p_c is the atmospheic pressure | ||
| - qt_mean, θl_mean is the grid mean q_tot and liquid ice potential temperature | ||
| - qt′qt′, θl′θl′, θl′qt′ are the (co)variances of q_tot and liquid ice potential temperature | ||
| - thermo params are the thermodynamics parameters | ||
| The function trnasforms and imposes additional limits on the quadrature points. | ||
| It returns cloud fraction computed as a sum over quadrature points. | ||
| """ | ||
| function quad_loop( | ||
| SG_quad::SGSQuadrature, | ||
| p_c, | ||
| qt_mean, | ||
| θl_mean, | ||
| qt′qt′, | ||
| θl′θl′, | ||
| θl′qt′, | ||
| thermo_params, | ||
| ) | ||
| FT = eltype(qt_mean) | ||
|
|
||
| sqrt2 = FT(sqrt(2)) | ||
|
|
||
| # Epsilon defined per typical variable fluctuation | ||
| eps_q = FT(eps(FT)) * max(FT(eps(FT)), qt_mean) | ||
| eps_θ = FT(eps(FT)) | ||
|
|
||
| # limit σ_q to prevent negative q_tot_hat | ||
| σ_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′) | ||
|
|
||
| # Enforce Cauchy-Schwarz inequality, numerically stable compute | ||
| _corr::FT = (θl′qt′ / max(σ_q, eps_q)) | ||
| corr::FT = max(min(_corr / max(σ_θ, eps_θ), 1), -1) | ||
|
|
||
| # Conditionals | ||
| σ_c = sqrt(max(1 - corr * corr, 0)) * σ_θ | ||
|
|
||
| function get_x_hat(χ::Tuple{<:Real, <:Real}) | ||
| μ_c = θl_mean + sqrt2 * corr * σ_θ * χ[1] | ||
| θ_hat = μ_c + sqrt2 * σ_c * χ[2] | ||
| q_tot_hat = qt_mean + sqrt2 * σ_q * χ[1] | ||
| return (θ_hat, q_tot_hat) | ||
| end | ||
|
|
||
| # cloudy/dry categories for buoyancy in TKE | ||
| f_q_tot_sat(x_hat::Tuple{<:Real, <:Real}, hc) = | ||
| hc ? x_hat[2] : eltype(x_hat)(0) | ||
|
|
||
| get_ts(x_hat::Tuple{<:Real, <:Real}) = | ||
| thermo_state(thermo_params; p = p_c, θ = x_hat[1], q_tot = x_hat[2]) | ||
| f_cf(x_hat::Tuple{<:Real, <:Real}, hc) = | ||
| hc ? eltype(x_hat)(1) : eltype(x_hat)(0) | ||
| function f(x_hat::Tuple{<:Real, <:Real}) | ||
| ts = get_ts(x_hat) | ||
| hc = TD.has_condensate(thermo_params, ts) | ||
| return (; cf = f_cf(x_hat, hc), q_tot_sat = f_q_tot_sat(x_hat, hc)) | ||
| end | ||
|
|
||
| return quad(f, get_x_hat, SG_quad).cf | ||
| end | ||
|
|
||
| """ | ||
| Compute f(A, B) as a sum over inner and outer quadrature points | ||
| that approximate the sub-grid scale variability of A and B | ||
| """ | ||
| 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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.