implement bcs for pedmf; flux at the bottom boundary

Author: sajjadazimi

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -129,132 +129,19 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_bottom_bc!(
 
     FT = Spaces.undertype(axes(Y.c))
     n = n_mass_flux_subdomains(turbconv_model)
-    thermo_params = CAP.thermodynamics_params(p.params)
     turbconv_params = CAP.turbconv_params(p.params)
+    ᶜaʲ_int_val = p.scratch.temp_data_level
 
-    (; ᶜΦ,) = p.core
-    (; ᶜp, ᶜK, ᶜtsʲs, ᶜρʲs, ᶜts) = p.precomputed
-    (; ustar, obukhov_length, buoyancy_flux) = p.precomputed.sfc_conditions
+    (; ᶜρʲs) = p.precomputed
 
     for j in 1:n
-        ᶜtsʲ = ᶜtsʲs.:($j)
-        ᶜmseʲ = Y.c.sgsʲs.:($j).mse
-        ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot
-        if moisture_model isa NonEquilMoistModel && (
-            microphysics_model isa Microphysics1Moment ||
-            microphysics_model isa Microphysics2Moment
-        )
-            ᶜq_liqʲ = Y.c.sgsʲs.:($j).q_liq
-            ᶜq_iceʲ = Y.c.sgsʲs.:($j).q_ice
-            ᶜq_raiʲ = Y.c.sgsʲs.:($j).q_rai
-            ᶜq_snoʲ = Y.c.sgsʲs.:($j).q_sno
-        end
-
-        # We need field_values everywhere because we are mixing
-        # information from surface and first interior inside the
-        # 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(
-            Fields.level(Fields.coordinate_field(Y.f).z, Fields.half),
-        )
-        ᶜρ_int_val = Fields.field_values(Fields.level(Y.c.ρ, 1))
-        ᶜp_int_val = Fields.field_values(Fields.level(ᶜp, 1))
-        (; ρ_flux_h_tot, ρ_flux_q_tot, ustar, obukhov_length) =
-            p.precomputed.sfc_conditions
-
-        buoyancy_flux_val = Fields.field_values(buoyancy_flux)
-        ρ_flux_h_tot_val = Fields.field_values(ρ_flux_h_tot)
-        ρ_flux_q_tot_val = Fields.field_values(ρ_flux_q_tot)
-
-        ustar_val = Fields.field_values(ustar)
-        obukhov_length_val = Fields.field_values(obukhov_length)
-        sfc_local_geometry_val = Fields.field_values(
-            Fields.local_geometry_field(Fields.level(Y.f, Fields.half)),
-        )
-
-        # Based on boundary conditions for updrafts we overwrite
-        # 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 = @. lazy(
-            TD.total_specific_enthalpy(
-                thermo_params,
-                ᶜts,
-                specific(Y.c.ρe_tot, Y.c.ρ),
-            ),
-        )
-        ᶜ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,
-            ustar_val,
-            obukhov_length_val,
-            sfc_local_geometry_val,
-        )
-
-        # ... and the first interior point for EDMFX ᶜq_totʲ.
-
-        ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
-        ᶜq_tot_int_val = Fields.field_values(Fields.level(ᶜq_tot, 1))
-        ᶜq_totʲ_int_val = Fields.field_values(Fields.level(ᶜq_totʲ, 1))
-        @. ᶜ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,
-            ustar_val,
-            obukhov_length_val,
-            sfc_local_geometry_val,
-        )
-
-        # Then overwrite the prognostic variables at first inetrior point.
-        ᶜΦ_int_val = Fields.field_values(Fields.level(ᶜΦ, 1))
-        ᶜtsʲ_int_val = Fields.field_values(Fields.level(ᶜtsʲ, 1))
-        if moisture_model isa NonEquilMoistModel && (
-            microphysics_model isa Microphysics1Moment ||
-            microphysics_model isa Microphysics2Moment
-        )
-            ᶜq_liqʲ_int_val = Fields.field_values(Fields.level(ᶜq_liqʲ, 1))
-            ᶜq_iceʲ_int_val = Fields.field_values(Fields.level(ᶜq_iceʲ, 1))
-            ᶜq_raiʲ_int_val = Fields.field_values(Fields.level(ᶜq_raiʲ, 1))
-            ᶜq_snoʲ_int_val = Fields.field_values(Fields.level(ᶜq_snoʲ, 1))
-            @. ᶜtsʲ_int_val = TD.PhaseNonEquil_phq(
-                thermo_params,
-                ᶜp_int_val,
-                ᶜmseʲ_int_val - ᶜΦ_int_val,
-                TD.PhasePartition(
-                    ᶜq_totʲ_int_val,
-                    ᶜq_liqʲ_int_val + ᶜq_raiʲ_int_val,
-                    ᶜq_iceʲ_int_val + ᶜq_snoʲ_int_val,
-                ),
-            )
-        else
-            @. ᶜtsʲ_int_val = TD.PhaseEquil_phq(
-                thermo_params,
-                ᶜp_int_val,
-                ᶜmseʲ_int_val - ᶜΦ_int_val,
-                ᶜq_totʲ_int_val,
-            )
-        end
         sgsʲs_ρ_int_val = Fields.field_values(Fields.level(ᶜρʲs.:($j), 1))
         sgsʲs_ρa_int_val =
             Fields.field_values(Fields.level(Y.c.sgsʲs.:($j).ρa, 1))
+        @. ᶜaʲ_int_val = draft_area(sgsʲs_ρa_int_val, sgsʲs_ρ_int_val)
 
-        @. sgsʲs_ρ_int_val = TD.air_density(thermo_params, ᶜtsʲ_int_val)
         @. sgsʲs_ρa_int_val =
-            $(FT(turbconv_params.surface_area)) *
-            TD.air_density(thermo_params, ᶜtsʲ_int_val)
+            $(FT(turbconv_params.surface_area)) * sgsʲs_ρ_int_val
     end
     return nothing
 end

src/prognostic_equations/edmfx_boundary_condition.jl

@@ -2,6 +2,74 @@
 ##### EDMFX SGS boundary condition
 #####
 
+"""
+    sgs_scalar_flux_bc(
+        z_sfc, ρ_sfc, χ_sfc, ᶜχ_int, ᶜaʲ_int, ᶜχʲ_int,
+        ρ_flux_χ, ustar, obukhov_length, sfc_local_geometry,
+    )
+
+Computes the surface scalar flux for an EDMF updraft subdomain based on the
+grid-mean scalar flux and an estimated subgrid (SGS) scalar fluctuation at the
+surface.
+
+This function takes the resolved (grid-mean) surface flux of a scalar `χ`
+(e.g., `mse` or `q_tot`) and returns the corresponding flux for an EDMF updraft.
+The updraft scalar value at the surface is approximated as 
+χʲ_sfc = χ_sfc + C * σ,
+where `χ_sfc` is the grid-mean surface scalar, `σ` is the estimated SGS scalar
+variance at the surface (computed from Monin–Obukhov similarity), and `C` is a
+coefficient based on the updraft fractional area. The resulting updraft flux is
+obtained by scaling the grid-mean flux according to the ratio of surface–interior
+scalar contrasts in the updraft and grid-mean fields.
+
+# Arguments
+- `z_sfc`: Surface height of the scalar exchange level.
+- `ρ_sfc`: Air density at the surface.
+- `χ_sfc`: Grid-mean scalar value at the surface.
+- `ᶜχ_int`: Grid-mean interior scalar value at the first model level.
+- `ᶜaʲ_int`: Updraft fractional area at the first model level.
+- `ᶜχʲ_int`: Updraft interior scalar value at the first model level.
+- `ρ_flux_χ`: Grid-mean surface flux of the scalar (mass flux form).
+- `ustar`: Friction velocity.
+- `obukhov_length`: Obukhov length for surface stability.
+- `sfc_local_geometry`: Local geometric factors for the surface exchange.
+
+# Returns
+- Updraft surface flux for scalar `χ` (same units as `ρ_flux_χ`).
+"""
+function sgs_scalar_flux_bc(
+    z_sfc::FT,
+    ρ_sfc,
+    χ_sfc,
+    ᶜχ_int,
+    ᶜaʲ_int,
+    ᶜχʲ_int,
+    ρ_flux_χ,
+    ustar,
+    obukhov_length,
+    sfc_local_geometry,
+) where {FT}
+
+    kinematic_sfc_flux_χ = ρ_flux_χ / ρ_sfc # Convert to kinematic flux [K m/s]
+
+    χ_sfc_var = get_scalar_variance(
+        kinematic_sfc_flux_χ,
+        ustar,
+        z_sfc,
+        obukhov_length,
+        sfc_local_geometry,
+    )
+
+    # TODO: This assumes that there is only one updraft, or that ᶜaʲ_int
+    #       is the specific area fraction for the updraft being considered when
+    #       sampling from the tail of the combined subgrid + grid-mean distribution.
+    #       The percentile range [1 - ᶜaʲ_int, 1] samples the top ᶜaʲ_int fraction.
+    χʲ_sfc_coeff = percentile_bounds_mean_norm(1 - ᶜaʲ_int, FT(1))
+    χʲ_sfc = χ_sfc + χʲ_sfc_coeff * sqrt(χ_sfc_var)
+
+    return (χʲ_sfc - ᶜχʲ_int) / (χ_sfc - ᶜχ_int) * ρ_flux_χ
+end
+
 """
     sgs_scalar_first_interior_bc(
         ᶜz_int::FT,
@@ -66,7 +134,7 @@ function sgs_scalar_first_interior_bc(
 
     kinematic_sfc_flux_scalar = sfc_ρ_flux_scalar / ᶜρ_int # Convert to kinematic flux [K m/s]
 
-    scalar_var = get_first_interior_variance(
+    scalar_var = get_scalar_variance(
         kinematic_sfc_flux_scalar,
         ustar,
         ᶜz_int,
@@ -86,7 +154,7 @@ function sgs_scalar_first_interior_bc(
 end
 
 """
-    get_first_interior_variance(
+    get_scalar_variance(
         kinematic_scalar_flux,
         ustar::FT,
         z,
@@ -115,7 +183,7 @@ Arguments:
 Returns:
 - The estimated variance of the scalar quantity [K² or (kg/kg)²].
 """
-function get_first_interior_variance(
+function get_scalar_variance(
     kinematic_scalar_flux,
     ustar::FT,
     z,

src/prognostic_equations/remaining_tendency.jl

@@ -279,6 +279,7 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
     end
 
     surface_flux_tendency!(Yₜ, Y, p, t)
+    edmfx_surface_flux_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
 
     radiation_tendency!(Yₜ, Y, p, t, p.atmos.radiation_mode)
     if p.atmos.sgs_entr_detr_mode == Explicit()

src/prognostic_equations/surface_flux.jl

@@ -123,6 +123,7 @@ function surface_flux_tendency!(Yₜ, Y, p, t)
     )
     btt = boundary_tendency_scalar(ᶜh_tot, sfc_conditions.ρ_flux_h_tot)
     @. Yₜ.c.ρe_tot -= btt
+
     ρ_flux_χ = p.scratch.sfc_temp_C3
     foreach_gs_tracer(Yₜ, Y) do ᶜρχₜ, ᶜρχ, ρχ_name
         ᶜχ = @. lazy(specific(ᶜρχ, Y.c.ρ))
@@ -138,3 +139,109 @@ function surface_flux_tendency!(Yₜ, Y, p, t)
         end
     end
 end
+
+"""
+    edmfx_surface_flux_tendency!(Yₜ, Y, p, t, turbconv_model::PrognosticEDMFX)
+
+Applies surface–flux tendencies to EDMF updraft prognostic variables when the
+turbulent convection scheme is the `PrognosticEDMFX` model.
+
+This function computes and adds contributions from surface fluxes of moist
+static energy (`mse`) and total specific humidity (`q_tot`) to the corresponding
+updraft subdomain tendencies in `Yₜ`. For each EDMF subdomain, it evaluates
+subgrid scalar fluxes using `sgs_scalar_flux_bc` and converts these fluxes into
+tendencies localized to the surface-adjacent grid cell via
+`boundary_tendency_scalar`.
+
+# Arguments
+- `Yₜ`: Tendency state vector, modified in place.
+- `Y`: Current state vector.
+- `p`: Cache containing parameters, thermodynamic settings, precomputed fields,
+       and EDMF configuration information.
+- `t`: Current simulation time.
+- `turbconv_model::PrognosticEDMFX`: Dispatch selector specifying the EDMF scheme.
+"""
+edmfx_surface_flux_tendency!(Yₜ, Y, p, t, turbconv_model) = nothing
+function edmfx_surface_flux_tendency!(Yₜ, Y, p, t, turbconv_model::PrognosticEDMFX)
+    p.atmos.disable_surface_flux_tendency && return
+
+    (; params) = p
+    (; ᶜρʲs, ᶜK, ᶜts) = p.precomputed
+    thermo_params = CAP.thermodynamics_params(params)
+    n = n_mass_flux_subdomains(p.atmos.turbconv_model)
+
+    ρ_flux_χʲ = p.scratch.sfc_temp_C3
+    ᶜaʲ = p.scratch.ᶜtemp_scalar
+
+    (; ρ_flux_h_tot, ρ_flux_q_tot, ustar, obukhov_length, ts) =
+        p.precomputed.sfc_conditions
+
+    ρ_flux_h_tot_val = Fields.field_values(ρ_flux_h_tot)
+    ρ_flux_q_tot_val = Fields.field_values(ρ_flux_q_tot)
+    ustar_val = Fields.field_values(ustar)
+    obukhov_length_val = Fields.field_values(obukhov_length)
+    sfc_local_geometry_val = Fields.field_values(
+        Fields.local_geometry_field(Fields.level(Y.f, Fields.half)),
+    )
+
+    ρ_sfc = @. lazy(TD.air_density(thermo_params, ts))
+    h_sfc = @. lazy(TD.specific_enthalpy(thermo_params, ts))
+    q_tot_sfc = @. lazy(TD.total_specific_humidity(thermo_params, ts))
+    z_sfc_val =
+        Fields.field_values(Fields.level(Fields.coordinate_field(Y.f).z, Fields.half))
+    ρ_sfc_val = Fields.field_values(ρ_sfc)
+    h_sfc_val = Fields.field_values(h_sfc)
+    q_tot_sfc_val = Fields.field_values(q_tot_sfc)
+
+    ᶜh_tot = @. lazy(
+        TD.total_specific_enthalpy(
+            thermo_params,
+            ᶜts,
+            specific(Y.c.ρe_tot, Y.c.ρ),
+        ),
+    )
+    ᶜh_tot_int_val = Fields.field_values(Fields.level(ᶜh_tot, 1))
+    ᶜK_int_val = Fields.field_values(Fields.level(ᶜK, 1))
+    ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
+    ᶜq_tot_int_val = Fields.field_values(Fields.level(ᶜq_tot, 1))
+
+    ρ_flux_χʲ_val = Fields.field_values(ρ_flux_χʲ)
+
+    for j in 1:n
+        @. ᶜaʲ = draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j))
+        ᶜaʲ_int_val = Fields.field_values(Fields.level(ᶜaʲ, 1))
+
+        ᶜmseʲ_int_val = Fields.field_values(Fields.level(Y.c.sgsʲs.:($j).mse, 1))
+        @. ρ_flux_χʲ_val = sgs_scalar_flux_bc(
+            z_sfc_val,
+            ρ_sfc_val,
+            h_sfc_val,
+            ᶜh_tot_int_val - ᶜK_int_val,
+            ᶜaʲ_int_val,
+            ᶜmseʲ_int_val,
+            ρ_flux_h_tot_val,
+            ustar_val,
+            obukhov_length_val,
+            sfc_local_geometry_val,
+        )
+        btt = boundary_tendency_scalar(Y.c.sgsʲs.:(1).mse, ρ_flux_χʲ)
+        @. Yₜ.c.sgsʲs.:($$j).mse -= btt / p.precomputed.ᶜρʲs.:($$j)
+
+        ᶜq_totʲ_int_val = Fields.field_values(Fields.level(Y.c.sgsʲs.:($j).q_tot, 1))
+        @. ρ_flux_χʲ_val = sgs_scalar_flux_bc(
+            z_sfc_val,
+            ρ_sfc_val,
+            q_tot_sfc_val,
+            ᶜq_tot_int_val,
+            ᶜaʲ_int_val,
+            ᶜq_totʲ_int_val,
+            ρ_flux_q_tot_val,
+            ustar_val,
+            obukhov_length_val,
+            sfc_local_geometry_val,
+        )
+        btt = boundary_tendency_scalar(Y.c.sgsʲs.:(1).q_tot, ρ_flux_χʲ)
+        @. Yₜ.c.sgsʲs.:($$j).q_tot -= btt / p.precomputed.ᶜρʲs.:($$j)
+    end
+
+end