From e2c5939bcab0d16b814e63a1e0706ef77656854a Mon Sep 17 00:00:00 2001 From: Zhaoyi Shen <11598433+szy21@users.noreply.github.com> Date: Thu, 11 Apr 2024 16:50:03 -0700 Subject: [PATCH] make edmf boundary condition dependent on surface area fraction --- Project.toml | 2 ++ docs/Manifest.toml | 2 +- examples/Manifest.toml | 2 +- perf/Manifest.toml | 2 +- regression_tests/ref_counter.jl | 9 ++++++-- .../diagnostic_edmf_precomputed_quantities.jl | 2 ++ .../prognostic_edmf_precomputed_quantities.jl | 12 ++++++++--- .../edmfx_boundary_condition.jl | 21 ++++++++++++++++--- 8 files changed, 41 insertions(+), 11 deletions(-) diff --git a/Project.toml b/Project.toml index 5e5aec717a..b1dc799ba7 100644 --- a/Project.toml +++ b/Project.toml @@ -20,6 +20,7 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" Insolation = "e98cc03f-d57e-4e3c-b70c-8d51efe9e0d8" @@ -61,6 +62,7 @@ Dates = "1" Dierckx = "0.5" DiffEqBase = "6" DiffEqCallbacks = "2" +Distributions = "0.25" DocStringExtensions = "0.8, 0.9" FastGaussQuadrature = "0.4, 0.5, 1" Insolation = "0.9.2" diff --git a/docs/Manifest.toml b/docs/Manifest.toml index 1dc5db8db8..48c3d46560 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -319,7 +319,7 @@ weakdeps = ["SparseArrays"] ChainRulesCoreSparseArraysExt = "SparseArrays" [[deps.ClimaAtmos]] -deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "ClimaComms", "ClimaCore", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Colors", "Dates", "Dierckx", "DiffEqBase", "DiffEqCallbacks", "DocStringExtensions", "FastGaussQuadrature", "Insolation", "Interpolations", "IntervalSets", "Krylov", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "Pkg", "Printf", "RRTMGP", "Random", "RootSolvers", "SciMLBase", "StaticArrays", "Statistics", "StatsBase", "SurfaceFluxes", "Thermodynamics", "YAML"] +deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "ClimaComms", "ClimaCore", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Colors", "Dates", "Dierckx", "DiffEqBase", "DiffEqCallbacks", "Distributions", "DocStringExtensions", "FastGaussQuadrature", "Insolation", "Interpolations", "IntervalSets", "Krylov", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "Pkg", "Printf", "RRTMGP", "Random", "RootSolvers", "SciMLBase", "StaticArrays", "Statistics", "StatsBase", "SurfaceFluxes", "Thermodynamics", "YAML"] path = ".." uuid = "b2c96348-7fb7-4fe0-8da9-78d88439e717" version = "0.23.0" diff --git a/examples/Manifest.toml b/examples/Manifest.toml index 7381aab9b7..caee870569 100644 --- a/examples/Manifest.toml +++ b/examples/Manifest.toml @@ -307,7 +307,7 @@ weakdeps = ["CairoMakie"] CairoMakieExt = "CairoMakie" [[deps.ClimaAtmos]] -deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "ClimaComms", "ClimaCore", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Colors", "Dates", "Dierckx", "DiffEqBase", "DiffEqCallbacks", "DocStringExtensions", "FastGaussQuadrature", "Insolation", "Interpolations", "IntervalSets", "Krylov", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "Pkg", "Printf", "RRTMGP", "Random", "RootSolvers", "SciMLBase", "StaticArrays", "Statistics", "StatsBase", "SurfaceFluxes", "Thermodynamics", "YAML"] +deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "ClimaComms", "ClimaCore", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Colors", "Dates", "Dierckx", "DiffEqBase", "DiffEqCallbacks", "Distributions", "DocStringExtensions", "FastGaussQuadrature", "Insolation", "Interpolations", "IntervalSets", "Krylov", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "Pkg", "Printf", "RRTMGP", "Random", "RootSolvers", "SciMLBase", "StaticArrays", "Statistics", "StatsBase", "SurfaceFluxes", "Thermodynamics", "YAML"] path = ".." uuid = "b2c96348-7fb7-4fe0-8da9-78d88439e717" version = "0.23.0" diff --git a/perf/Manifest.toml b/perf/Manifest.toml index 96bcbd232d..a4657a5245 100644 --- a/perf/Manifest.toml +++ b/perf/Manifest.toml @@ -322,7 +322,7 @@ version = "0.5.0" GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6" [[deps.ClimaAtmos]] -deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "ClimaComms", "ClimaCore", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Colors", "Dates", "Dierckx", "DiffEqBase", "DiffEqCallbacks", "DocStringExtensions", "FastGaussQuadrature", "Insolation", "Interpolations", "IntervalSets", "Krylov", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "Pkg", "Printf", "RRTMGP", "Random", "RootSolvers", "SciMLBase", "StaticArrays", "Statistics", "StatsBase", "SurfaceFluxes", "Thermodynamics", "YAML"] +deps = ["Adapt", "ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "ClimaComms", "ClimaCore", "ClimaParams", "ClimaTimeSteppers", "ClimaUtilities", "CloudMicrophysics", "Colors", "Dates", "Dierckx", "DiffEqBase", "DiffEqCallbacks", "Distributions", "DocStringExtensions", "FastGaussQuadrature", "Insolation", "Interpolations", "IntervalSets", "Krylov", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "Pkg", "Printf", "RRTMGP", "Random", "RootSolvers", "SciMLBase", "StaticArrays", "Statistics", "StatsBase", "SurfaceFluxes", "Thermodynamics", "YAML"] path = ".." uuid = "b2c96348-7fb7-4fe0-8da9-78d88439e717" version = "0.23.0" diff --git a/regression_tests/ref_counter.jl b/regression_tests/ref_counter.jl index 7a9aee41ba..3bfd645718 100644 --- a/regression_tests/ref_counter.jl +++ b/regression_tests/ref_counter.jl @@ -1,8 +1,13 @@ -158 +159 + +# 159: +# - Changed the boundary condition of edmf updraft properties +# to be dependent on the surface area # 158: # - Switched back the precipitation threshold defintion in the -# 0-moment scheme to specific humidity +# 0-moment scheme to specific humidity + # 157: # - For the grid mean precipitation tendency in the 0-moment scheme: # - added limiting by q_tot/dt diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl index c5f3971ccd..3b4ac876b7 100644 --- a/src/cache/diagnostic_edmf_precomputed_quantities.jl +++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl @@ -155,6 +155,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!( @. mseʲ_int_level = sgs_scalar_first_interior_bc( z_int_level - z_sfc_halflevel, ρ_int_level, + turbconv_model.a_int, h_tot_int_level - K_int_level, buoyancy_flux_sfc_halflevel, ρ_flux_h_tot_sfc_halflevel, @@ -165,6 +166,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!( @. q_totʲ_int_level = sgs_scalar_first_interior_bc( z_int_level - z_sfc_halflevel, ρ_int_level, + turbconv_model.a_int, q_tot_int_level, buoyancy_flux_sfc_halflevel, ρ_flux_q_tot_sfc_halflevel, diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl index 8542205849..0988224d46 100644 --- a/src/cache/prognostic_edmf_precomputed_quantities.jl +++ b/src/cache/prognostic_edmf_precomputed_quantities.jl @@ -70,6 +70,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft_and_bc! (; params) = p thermo_params = CAP.thermodynamics_params(params) + turbconv_params = CAP.turbconv_params(params) (; ᶜΦ,) = p.core (; ᶜspecific, ᶜp, ᶜh_tot, ᶜK) = p.precomputed @@ -96,7 +97,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft_and_bc! # We need field_values everywhere because we are mixing # information from surface and first interior inside the - # sgs_h/q_tot_first_interior_bc call. + # sgs_scalar_first_interior_bc call. ᶜz_int_val = Fields.field_values(Fields.level(Fields.coordinate_field(Y.c).z, 1)) z_sfc_val = Fields.field_values( @@ -116,13 +117,18 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft_and_bc! ) # Based on boundary conditions for updrafts we overwrite - # the first interior point for EDMFX ᶜh_totʲ... + # the first interior point for EDMFX ᶜmseʲ... + ᶜaʲ_int_val = p.scratch.temp_data_level + # TODO: replace this with the actual surface area fraction when + # using prognostic surface area + @. ᶜaʲ_int_val = FT(turbconv_params.surface_area) ᶜh_tot_int_val = Fields.field_values(Fields.level(ᶜh_tot, 1)) ᶜK_int_val = Fields.field_values(Fields.level(ᶜK, 1)) ᶜmseʲ_int_val = Fields.field_values(Fields.level(ᶜmseʲ, 1)) @. ᶜmseʲ_int_val = sgs_scalar_first_interior_bc( ᶜz_int_val - z_sfc_val, ᶜρ_int_val, + ᶜaʲ_int_val, ᶜh_tot_int_val - ᶜK_int_val, buoyancy_flux_val, ρ_flux_h_tot_val, @@ -137,6 +143,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft_and_bc! @. ᶜq_totʲ_int_val = sgs_scalar_first_interior_bc( ᶜz_int_val - z_sfc_val, ᶜρ_int_val, + ᶜaʲ_int_val, ᶜq_tot_int_val, buoyancy_flux_val, ρ_flux_q_tot_val, @@ -158,7 +165,6 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft_and_bc! sgsʲs_ρa_int_val = Fields.field_values(Fields.level(Y.c.sgsʲs.:($j).ρa, 1)) - turbconv_params = CAP.turbconv_params(params) @. sgsʲs_ρ_int_val = TD.air_density(thermo_params, ᶜtsʲ_int_val) @. sgsʲs_ρa_int_val = $(FT(turbconv_params.surface_area)) * diff --git a/src/prognostic_equations/edmfx_boundary_condition.jl b/src/prognostic_equations/edmfx_boundary_condition.jl index a54c7d2b20..b10c6f4749 100644 --- a/src/prognostic_equations/edmfx_boundary_condition.jl +++ b/src/prognostic_equations/edmfx_boundary_condition.jl @@ -2,9 +2,12 @@ ##### EDMFX SGS boundary condition ##### +import Distributions + function sgs_scalar_first_interior_bc( ᶜz_int::FT, ᶜρ_int::FT, + ᶜaʲ_int::FT, ᶜscalar_int::FT, sfc_buoyancy_flux, sfc_ρ_flux_scalar, @@ -13,8 +16,6 @@ function sgs_scalar_first_interior_bc( sfc_local_geometry, ) where {FT} sfc_buoyancy_flux > 0 || return ᶜscalar_int - # a_total = edmf.surface_area - # a_ = area_surface_bc(surface_conditions, edmf) scalar_var = get_first_interior_variance( sfc_ρ_flux_scalar / ᶜρ_int, ustar, @@ -26,7 +27,8 @@ function sgs_scalar_first_interior_bc( # 1 - a_total + (i - 1) * a_, # 1 - a_total + i * a_, # ) - surface_scalar_coeff = FT(1.75) + # TODO: This assumes that there is only one updraft + surface_scalar_coeff = percentile_bounds_mean_norm(1 - ᶜaʲ_int, FT(1)) return ᶜscalar_int + surface_scalar_coeff * sqrt(scalar_var) end @@ -46,3 +48,16 @@ function get_first_interior_variance( return 4 * Geometry._norm_sqr(c_star, local_geometry) end end + +function percentile_bounds_mean_norm( + low_percentile::FT, + high_percentile::FT, +) where {FT <: Real} + gauss_int(x) = -exp(-x * x / 2) / sqrt(2 * pi) + xp_low = Distributions.quantile(Distributions.Normal(), low_percentile) + xp_high = Distributions.quantile(Distributions.Normal(), high_percentile) + return (gauss_int(xp_high) - gauss_int(xp_low)) / ( + Distributions.cdf(Distributions.Normal(), xp_high) - + Distributions.cdf(Distributions.Normal(), xp_low) + ) +end