Merge pull request #2898 from CliMA/zs/edmf_scalar_bc

make edmf boundary condition dependent on surface area fraction

Author: szy21

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"

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"

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"

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"

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

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,

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

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