pedmf bcs; updraft fluxes at the bottom boundary are computed based on surface-interior tracer contrasts and grid mean fluxes

Author: sajjadazimi

reproducibility_tests/ref_counter.jl

@@ -1,4 +1,4 @@
-274
+275
 
 # **README**
 #
@@ -20,6 +20,10 @@
 
 
 #=
+275
+- Change boundary conditions for prognostic EDMF; updraft surface fluxes are
+  computed based on surface-interior tracer contrasts and grid-mean fluxes
+
 274
 - Remove unused calculation of TKE exchange in mixing length
 

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,54 @@
 ##### EDMFX SGS boundary condition
 #####
 
+"""
+    sgs_scalar_flux_bc(
+        χ_sfc, ᶜχ_int, ᶜχʲ_int, ᶜaʲ_int, ρ_flux_χ
+    )
+
+Computes the surface scalar flux for an EDMF updraft subdomain by scaling the
+grid-mean surface flux according to the relative surface–interior scalar
+contrasts in the updraft and the grid mean.
+
+In this simplified formulation, the updraft scalar value at the surface is
+assumed to equal the grid-mean surface value (`χ_sfc`). The updraft flux is
+obtained by multiplying the grid-mean scalar flux (`ρ_flux_χ`) by the ratio
+(χ_sfc - ᶜχʲ_int)/(χ_sfc - ᶜχ_int), with appropriate limiting for numerical 
+stability. If the surface–interior contrast in the grid mean is negligible, 
+the grid-mean flux is returned to avoid division by zero.
+
+# Arguments
+- `χ_sfc`: Scalar value at the surface (same for grid-mean and updraft).
+- `ᶜχ_int`: Grid-mean interior scalar at the first model level.
+- `ᶜχʲ_int`: Updraft interior scalar at the first model level.
+- `ᶜaʲ_int`: Updraft fractional area at the first model level.
+- `ρ_flux_χ`: Grid-mean surface scalar flux (mass-flux form).
+
+# Returns
+- Updraft surface scalar flux for `χ` (same units as `ρ_flux_χ`), scaled by the
+  limited ratio of updraft to grid-mean scalar contrasts.
+"""
+
+function sgs_scalar_flux_bc(
+    χ_sfc::FT,
+    ᶜχ_int,
+    ᶜχʲ_int,
+    ᶜaʲ_int,
+    ρ_flux_χ,
+) where {FT}
+
+    # when surface-interior difference on the grid mean is zero (negligible), 
+    # return grid mean flux (which is zero or negligible) to avoid division by zero
+    if abs(χ_sfc - ᶜχ_int) < sqrt(floatmin(FT))
+        return ρ_flux_χ
+    end
+
+    # we limit the ratio of sgs to gs scalar flux for numerical stability; physically we don't need a limit
+    limit = max(0, 1 / ᶜaʲ_int)
+    sgs_to_gs_flux_ratio = max(-limit, min(limit, (χ_sfc - ᶜχʲ_int) / (χ_sfc - ᶜχ_int)))
+    return sgs_to_gs_flux_ratio * ρ_flux_χ
+end
+
 """
     sgs_scalar_first_interior_bc(
         ᶜz_int::FT,

src/prognostic_equations/mass_flux_closures.jl

@@ -206,6 +206,7 @@ end
 Apply EDMF filters:
  - Relax u_3 to zero when it is negative
  - Relax ρa to zero when it is negative
+ - Relax ρa to ρ when a is larger than one
 
 This function modifies the tendency `Yₜ` in place based on the current state `Y`,
 parameters `p`, time `t`, and the turbulence convection model `turbconv_model`.
@@ -225,6 +226,8 @@ function edmfx_filter_tendency!(Yₜ, Y, p, t, turbconv_model::PrognosticEDMFX)
             @. Yₜ.f.sgsʲs.:($$j).u₃ -=
                 C3(min(Y.f.sgsʲs.:($$j).u₃.components.data.:1, 0)) / FT(dt)
             @. Yₜ.c.sgsʲs.:($$j).ρa -= min(Y.c.sgsʲs.:($$j).ρa, 0) / FT(dt)
+            @. Yₜ.c.sgsʲs.:($$j).ρa -=
+                max(Y.c.sgsʲs.:($$j).ρa - p.precomputed.ᶜρʲs.:($$j), 0) / FT(dt)
         end
     end
 end

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,101 @@ 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_h_tot, ts) =
+        p.precomputed.sfc_conditions
+
+    ᶜaʲ = p.scratch.ᶜtemp_scalar
+    ρ_flux_χʲ = p.scratch.sfc_temp_C3
+    # We need field_values everywhere because we are mixing
+    # information from surface and first interior inside the
+    # sgs_scalar_flux_bc call.
+    ρ_flux_χʲ_val = Fields.field_values(ρ_flux_χʲ)
+
+    # Based on boundary conditions for updrafts we compute
+    # the tendency due to the surface flux for EDMFX ᶜmseʲ...
+    ρ_flux_h_tot_val = Fields.field_values(ρ_flux_h_tot)
+    h_tot_sfc = @. lazy(TD.specific_enthalpy(thermo_params, ts))
+    h_tot_sfc_val = Fields.field_values(h_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))
+
+    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(
+            h_tot_sfc_val,
+            ᶜh_tot_int_val - ᶜK_int_val,
+            ᶜmseʲ_int_val,
+            ᶜaʲ_int_val,
+            ρ_flux_h_tot_val,
+        )
+        btt = boundary_tendency_scalar(Y.c.sgsʲs.:(1).mse, ρ_flux_χʲ)
+        @. Yₜ.c.sgsʲs.:($$j).mse -= btt / p.precomputed.ᶜρʲs.:($$j)
+    end
+
+    # ... and the tendency for EDMFX ᶜq_totʲ.
+    if !(p.atmos.moisture_model isa DryModel)
+        ρ_flux_q_tot_val = Fields.field_values(p.precomputed.sfc_conditions.ρ_flux_q_tot)
+        q_tot_sfc = @. lazy(TD.total_specific_humidity(thermo_params, ts))
+        q_tot_sfc_val = Fields.field_values(q_tot_sfc)
+        ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
+        ᶜq_tot_int_val = Fields.field_values(Fields.level(ᶜq_tot, 1))
+
+        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))
+            ᶜq_totʲ_int_val = Fields.field_values(Fields.level(Y.c.sgsʲs.:($j).q_tot, 1))
+
+            @. ρ_flux_χʲ_val = sgs_scalar_flux_bc(
+                q_tot_sfc_val,
+                ᶜq_tot_int_val,
+                ᶜq_totʲ_int_val,
+                ᶜaʲ_int_val,
+                ρ_flux_q_tot_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
+
+end