Merge #1354

1354: Include prognostic and closure options for surface area BC r=costachris a=costachris

- Adds options for prognostic and closure-based determination of bottom boundary condition on area. For prognostic option, surface area is prognostically determined in the bottom cell center, instead of prescribed as surface_area parameter. For closure option, bottom area determined as a linear function of surface heat fluxes.
-  Initial conditions for prognostic surface area: area fraction and updraft w set to nonzero values in bottom cell center and top face, along with the associated prognostic variables. This initialization is to break symmetry and prevent trivial a = 0, w = 0 solution.
- Swaps order of set_edmf_surface_bc and filter_updraft_vars, since updraft variables are used to compute surface theta and qt fluxes.

Co-authored-by: costachris <christopouloscosta@gmail.com>

driver/compute_diagnostics.jl

@@ -303,16 +303,14 @@ function compute_diagnostics!(
         end
     end
 
-    a_bulk_bcs = TC.a_bulk_boundary_conditions(surf, edmf)
-    Ifabulk = CCO.InterpolateC2F(; a_bulk_bcs...)
+    Ifabulk = CCO.InterpolateC2F()
     a_up_bulk_f = TC.face_aux_turbconv(state).bulk.a_up
     @. a_up_bulk_f = Ifabulk(a_up_bulk)
 
     RB_precip = CCO.RightBiasedC2F(; top = CCO.SetValue(FT(0)))
 
     @inbounds for i in 1:N_up
-        a_up_bcs = TC.a_up_boundary_conditions(surf, edmf, i)
-        Ifaup = CCO.InterpolateC2F(; a_up_bcs...)
+        Ifaup = CCO.InterpolateC2F()
         a_up_f = aux_up_f[i].area
         @. a_up_f = Ifaup(aux_up[i].area)
         @inbounds for k in TC.real_face_indices(grid)

driver/generate_namelist.jl

@@ -88,6 +88,8 @@ function default_namelist(
 
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"] = Dict()
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["surface_area"] = 0.1
+    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["surface_area_bc"] = "Fixed" #{"Fixed", "Prognostic", "Closure"}
+    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["surface_area_bc_params"] = [1e-3, 1e-3, 1e-3]
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["max_area"] = 0.9
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["min_area"] = 1e-5
 

driver/initial_conditions.jl

@@ -12,7 +12,7 @@ function initialize_edmf(edmf::TC.EDMFModel, grid::TC.Grid, state::TC.State, sur
     @. aux_gm.θ_virt = TD.virtual_pottemp(thermo_params, ts_gm)
     surf = get_surface(surf_params, grid, state, t, param_set)
     if case isa Cases.DryBubble
-        initialize_updrafts_DryBubble(edmf, grid, state)
+        initialize_updrafts_DryBubble(edmf, grid, state, surf)
     else
         initialize_updrafts(edmf, grid, state, surf)
     end
@@ -48,14 +48,18 @@ end
 function initialize_updrafts(edmf, grid, state, surf)
     N_up = TC.n_updrafts(edmf)
     kc_surf = TC.kc_surface(grid)
+    kf_surf = TC.kf_surface(grid)
     aux_up = TC.center_aux_updrafts(state)
     prog_gm = TC.center_prog_grid_mean(state)
     aux_up = TC.center_aux_updrafts(state)
     aux_up_f = TC.face_aux_updrafts(state)
     aux_gm = TC.center_aux_grid_mean(state)
+    aux_gm_f = TC.face_aux_grid_mean(state)
     prog_up = TC.center_prog_updrafts(state)
     prog_up_f = TC.face_prog_updrafts(state)
     ρ_c = prog_gm.ρ
+    ρ_f = aux_gm_f.ρ
+    FT = TC.float_type(state)
     @inbounds for i in 1:N_up
         @inbounds for k in TC.real_face_indices(grid)
             aux_up_f[i].w[k] = 0
@@ -81,16 +85,31 @@ function initialize_updrafts(edmf, grid, state, surf)
             @. prog_up[i].δ_nondim = 0
         end
 
-        a_surf = TC.area_surface_bc(surf, edmf, i)
-        aux_up[i].area[kc_surf] = a_surf
-        prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_surf
+        if edmf.surface_area_bc isa TC.FixedSurfaceAreaBC || edmf.surface_area_bc isa TC.ClosureSurfaceAreaBC
+            a_surf = TC.area_surface_bc(surf, edmf, i, edmf.surface_area_bc)
+            aux_up[i].area[kc_surf] = a_surf
+            prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_surf
+
+            # To avoid degeneracy from 0 BC on surface updraft area fraction, add small perturbation
+            # to relevant prognostic variables in bottom cell -> This avoids trivial solution with no updrafts.
+        elseif edmf.surface_area_bc isa TC.PrognosticSurfaceAreaBC
+            a_up_initial = FT(0.1)
+            aux_up[i].area[kc_surf] = a_up_initial
+            prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_up_initial
+
+            aux_up_f[i].w[kf_surf + 1] = FT(0.01) # small velocity perturbation
+            prog_up_f[i].ρaw[kf_surf + 1] = ρ_f[kf_surf + 1] * a_up_initial * FT(0.01) # small ρaw perturbation
+
+            prog_up[i].ρaq_tot[kc_surf] = ρ_c[kc_surf] * a_up_initial * aux_gm.q_tot[kc_surf]
+            prog_up[i].ρaθ_liq_ice[kc_surf] = ρ_c[kc_surf] * a_up_initial * aux_gm.θ_liq_ice[kc_surf]
+        end
     end
     return
 end
 
 import AtmosphericProfilesLibrary
 const APL = AtmosphericProfilesLibrary
-function initialize_updrafts_DryBubble(edmf, grid, state)
+function initialize_updrafts_DryBubble(edmf, grid, state, surf)
 
     # criterion 2: b>1e-4
     aux_up = TC.center_aux_updrafts(state)
@@ -100,6 +119,7 @@ function initialize_updrafts_DryBubble(edmf, grid, state)
     prog_gm = TC.center_prog_grid_mean(state)
     prog_up = TC.center_prog_updrafts(state)
     prog_up_f = TC.face_prog_updrafts(state)
+    kc_surf = TC.kc_surface(grid)
     ρ_0_c = prog_gm.ρ
     ρ_0_f = aux_gm_f.ρ
     N_up = TC.n_updrafts(edmf)
@@ -143,6 +163,9 @@ function initialize_updrafts_DryBubble(edmf, grid, state)
             end
         end
         @. prog_up_f[i].ρaw = If(prog_up[i].ρarea) * aux_up_f[i].w
+        a_surf = TC.area_surface_bc(surf, edmf, i, edmf.surface_area_bc)
+        aux_up[i].area[kc_surf] = a_surf
+        prog_up[i].ρarea[kc_surf] = ρ_0_c[kc_surf] * a_surf
     end
     return nothing
 end

integration_tests/Artifacts.toml

@@ -5,4 +5,4 @@ git-tree-sha1 = "3d6bb62d803fb01908ffe85469a89933f28ced45"
 git-tree-sha1 = "c2d323fb412d1250182be3aeadf9d94e816550a0"
 
 [SCAMPy_output]
-git-tree-sha1 = "6b82611f969d3cf99102baa7d898e3c759cf014c"
+git-tree-sha1 = "f725942089ff0d0fa8138dcbd1349996c746a4f9"

src/EDMF_functions.jl

@@ -63,8 +63,7 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
     θ_liq_ice_gm = aux_gm.θ_liq_ice
     q_tot_gm = aux_gm.q_tot
     q_tot_en = aux_en.q_tot
-    a_en_bcs = a_en_boundary_conditions(surf, edmf)
-    Ifae = CCO.InterpolateC2F(; a_en_bcs...)
+    Ifae = CCO.InterpolateC2F()
     If = CCO.InterpolateC2F(; bottom = CCO.SetValue(FT(0)), top = CCO.SetValue(FT(0)))
 
     # Compute the mass flux and associated scalar fluxes
@@ -74,15 +73,13 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
     @inbounds for i in 1:N_up
         aux_up_f_i = aux_up_f[i]
         aux_up_i = aux_up[i]
-        a_up_bcs = a_up_boundary_conditions(surf, edmf, i)
-        Ifau = CCO.InterpolateC2F(; a_up_bcs...)
         a_up = aux_up[i].area
         w_up_i = aux_up_f[i].w
         q_tot_up = aux_up_i.q_tot
         θ_liq_ice_up = aux_up_i.θ_liq_ice
-        @. aux_up_f[i].massflux = ρ_f * Ifau(a_up) * (w_up_i - toscalar(w_gm))
-        @. massflux_h += ρ_f * (Ifau(a_up) * (w_up_i - toscalar(w_gm)) * (If(θ_liq_ice_up) - If(θ_liq_ice_gm)))
-        @. massflux_qt += ρ_f * (Ifau(a_up) * (w_up_i - toscalar(w_gm)) * (If(q_tot_up) - If(q_tot_gm)))
+        @. aux_up_f[i].massflux = ρ_f * Ifae(a_up) * (w_up_i - toscalar(w_gm))
+        @. massflux_h += ρ_f * (Ifae(a_up) * (w_up_i - toscalar(w_gm)) * (If(θ_liq_ice_up) - If(θ_liq_ice_gm)))
+        @. massflux_qt += ρ_f * (Ifae(a_up) * (w_up_i - toscalar(w_gm)) * (If(q_tot_up) - If(q_tot_gm)))
     end
 
     if edmf.moisture_model isa NonEquilibriumMoisture
@@ -242,8 +239,8 @@ function affect_filter!(edmf::EDMFModel, grid::Grid, state::State, param_set::AP
     ###
     ### Filters
     ###
-    set_edmf_surface_bc(edmf, grid, state, surf, param_set)
     filter_updraft_vars(edmf, grid, state, surf)
+    set_edmf_surface_bc(edmf, grid, state, surf, param_set)
 
     @inbounds for k in real_center_indices(grid)
         prog_en.ρatke[k] = max(prog_en.ρatke[k], 0.0)
@@ -274,12 +271,12 @@ function set_edmf_surface_bc(edmf::EDMFModel, grid::Grid, state::State, surf::Su
     cp = TD.cp_m(thermo_params, ts_gm[kc_surf])
     ρ_c = prog_gm.ρ
     ρ_f = aux_gm_f.ρ
-    ae_surf::FT = 1
+    ρa_env_surf::FT = ρ_c[kc_surf]
     @inbounds for i in 1:N_up
         θ_surf = θ_surface_bc(surf, grid, state, edmf, i)
         q_surf = q_surface_bc(surf, grid, state, edmf, i)
-        a_surf = area_surface_bc(surf, edmf, i)
-        prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_surf
+
+        ρa_surf = prog_up[i].ρarea[kc_surf]
         prog_up[i].ρaθ_liq_ice[kc_surf] = prog_up[i].ρarea[kc_surf] * θ_surf
         prog_up[i].ρaq_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * q_surf
         if edmf.moisture_model isa NonEquilibriumMoisture
@@ -289,7 +286,7 @@ function set_edmf_surface_bc(edmf::EDMFModel, grid::Grid, state::State, surf::Su
             prog_up[i].ρaq_ice[kc_surf] = prog_up[i].ρarea[kc_surf] * q_ice_surf
         end
         prog_up_f[i].ρaw[kf_surf] = ρ_f[kf_surf] * w_surface_bc(surf)
-        ae_surf -= a_surf
+        ρa_env_surf -= ρa_surf
     end
 
     flux1 = surf.ρθ_liq_ice_flux
@@ -298,13 +295,12 @@ function set_edmf_surface_bc(edmf::EDMFModel, grid::Grid, state::State, surf::Su
     ustar = surf.ustar
     oblength = surf.obukhov_length
     ρLL = prog_gm.ρ[kc_surf]
-    ρ_ae = ρ_c[kc_surf] * ae_surf
     mix_len_params = mixing_length_params(edmf)
-    prog_en.ρatke[kc_surf] = ρ_ae * get_surface_tke(mix_len_params, surf.ustar, zLL, surf.obukhov_length)
+    prog_en.ρatke[kc_surf] = ρa_env_surf * get_surface_tke(mix_len_params, surf.ustar, zLL, surf.obukhov_length)
     if edmf.thermo_covariance_model isa PrognosticThermoCovariances
-        prog_en.ρaHvar[kc_surf] = ρ_ae * get_surface_variance(flux1 / ρLL, flux1 / ρLL, ustar, zLL, oblength)
-        prog_en.ρaQTvar[kc_surf] = ρ_ae * get_surface_variance(flux2 / ρLL, flux2 / ρLL, ustar, zLL, oblength)
-        prog_en.ρaHQTcov[kc_surf] = ρ_ae * get_surface_variance(flux1 / ρLL, flux2 / ρLL, ustar, zLL, oblength)
+        prog_en.ρaHvar[kc_surf] = ρa_env_surf * get_surface_variance(flux1 / ρLL, flux1 / ρLL, ustar, zLL, oblength)
+        prog_en.ρaQTvar[kc_surf] = ρa_env_surf * get_surface_variance(flux2 / ρLL, flux2 / ρLL, ustar, zLL, oblength)
+        prog_en.ρaHQTcov[kc_surf] = ρa_env_surf * get_surface_variance(flux1 / ρLL, flux2 / ρLL, ustar, zLL, oblength)
     end
     return nothing
 end
@@ -320,26 +316,15 @@ function surface_helper(surf::SurfaceBase, grid::Grid, state::State)
     return (; ustar, zLL, oblength, ρLL)
 end
 
-function a_up_boundary_conditions(surf::SurfaceBase, edmf::EDMFModel, i::Int)
-    a_surf = area_surface_bc(surf, edmf, i)
-    return (; bottom = CCO.SetValue(a_surf), top = CCO.Extrapolate())
-end
-
-function a_bulk_boundary_conditions(surf::SurfaceBase, edmf::EDMFModel)
-    N_up = n_updrafts(edmf)
-    a_surf = sum(i -> area_surface_bc(surf, edmf, i), 1:N_up)
-    return (; bottom = CCO.SetValue(a_surf), top = CCO.Extrapolate())
-end
-
-function a_en_boundary_conditions(surf::SurfaceBase, edmf::EDMFModel)
+function area_surface_bc(surf::SurfaceBase{FT}, edmf::EDMFModel, i::Int, bc::FixedSurfaceAreaBC)::FT where {FT}
     N_up = n_updrafts(edmf)
-    a_surf = 1 - sum(i -> area_surface_bc(surf, edmf, i), 1:N_up)
-    return (; bottom = CCO.SetValue(a_surf), top = CCO.Extrapolate())
+    return surf.bflux > 0 ? edmf.surface_area / N_up : FT(0)
 end
 
-function area_surface_bc(surf::SurfaceBase{FT}, edmf::EDMFModel, i::Int)::FT where {FT}
-    N_up = n_updrafts(edmf)
-    return surf.bflux > 0 ? edmf.surface_area / N_up : FT(0)
+function area_surface_bc(surf::SurfaceBase{FT}, edmf::EDMFModel, i::Int, bc::ClosureSurfaceAreaBC)::FT where {FT}
+    params = bc.params
+    surf_area_bc_pred = params[1] + params[2] * surf.lhf + params[3] * surf.shf
+    return min(max(0, FT(surf_area_bc_pred)), edmf.max_area)
 end
 
 function w_surface_bc(::SurfaceBase{FT})::FT where {FT}
@@ -351,14 +336,18 @@ end
 function θ_surface_bc(surf::SurfaceBase{FT}, grid::Grid, state::State, edmf::EDMFModel, i::Int)::FT where {FT}
     aux_gm = center_aux_grid_mean(state)
     prog_gm = center_prog_grid_mean(state)
+    prog_up = center_prog_updrafts(state)
     ρ_c = prog_gm.ρ
     kc_surf = kc_surface(grid)
     ts_gm = aux_gm.ts
     UnPack.@unpack ustar, zLL, oblength, ρLL = surface_helper(surf, grid, state)
+    N_up = n_updrafts(edmf)
 
     surf.bflux > 0 || return FT(0)
-    a_total = edmf.surface_area
-    a_ = area_surface_bc(surf, edmf, i)
+
+    a_total = sum(i -> prog_up[i].ρarea[kc_surf] / ρ_c[kc_surf], 1:N_up)
+    a_ = prog_up[i].ρarea[kc_surf] / ρ_c[kc_surf]
+
     ρθ_liq_ice_flux = surf.ρθ_liq_ice_flux # assuming no ql,qi flux
     h_var = get_surface_variance(ρθ_liq_ice_flux / ρLL, ρθ_liq_ice_flux / ρLL, ustar, zLL, oblength)
     surface_scalar_coeff = percentile_bounds_mean_norm(1 - a_total + (i - 1) * a_, 1 - a_total + i * a_)
@@ -367,11 +356,16 @@ end
 function q_surface_bc(surf::SurfaceBase{FT}, grid::Grid, state::State, edmf::EDMFModel, i::Int)::FT where {FT}
     aux_gm = center_aux_grid_mean(state)
     prog_gm = center_prog_grid_mean(state)
+    prog_up = center_prog_updrafts(state)
     ρ_c = prog_gm.ρ
     kc_surf = kc_surface(grid)
+    N_up = n_updrafts(edmf)
+
     surf.bflux > 0 || return aux_gm.q_tot[kc_surf]
-    a_total = edmf.surface_area
-    a_ = area_surface_bc(surf, edmf, i)
+
+    a_total = sum(i -> prog_up[i].ρarea[kc_surf] / ρ_c[kc_surf], 1:N_up)
+    a_ = prog_up[i].ρarea[kc_surf] / ρ_c[kc_surf]
+
     UnPack.@unpack ustar, zLL, oblength, ρLL = surface_helper(surf, grid, state)
     ρq_tot_flux = surf.ρq_tot_flux
     qt_var = get_surface_variance(ρq_tot_flux / ρLL, ρq_tot_flux / ρLL, ustar, zLL, oblength)
@@ -516,7 +510,6 @@ function compute_up_tendencies!(edmf::EDMFModel, grid::Grid, state::State, param
         tends_ρarea = tendencies_up[i].ρarea
         tends_ρaθ_liq_ice = tendencies_up[i].ρaθ_liq_ice
         tends_ρaq_tot = tendencies_up[i].ρaq_tot
-
         @. tends_ρarea =
             -∇c(wvec(LBF(Ic(w_up) * ρarea))) + (ρarea * Ic(w_up) * entr_turb_dyn) - (ρarea * Ic(w_up) * detr_turb_dyn)
 
@@ -573,7 +566,6 @@ function compute_up_tendencies!(edmf::EDMFModel, grid::Grid, state::State, param
             @. tends_δ_nondim = δ_λ * (mean_detr - δ_nondim)
         end
 
-        tends_ρarea[kc_surf] = 0
         tends_ρaθ_liq_ice[kc_surf] = 0
         tends_ρaq_tot[kc_surf] = 0
     end
@@ -584,13 +576,12 @@ function compute_up_tendencies!(edmf::EDMFModel, grid::Grid, state::State, param
     # and buoyancy should not matter in the end
     zero_bcs = (; bottom = CCO.SetValue(FT(0)), top = CCO.SetValue(FT(0)))
     I0f = CCO.InterpolateC2F(; zero_bcs...)
+    Iaf = CCO.InterpolateC2F()
     adv_bcs = (; bottom = CCO.SetValue(wvec(FT(0))), top = CCO.SetValue(wvec(FT(0))))
     LBC = CCO.LeftBiasedF2C(; bottom = CCO.SetValue(FT(0)))
     ∇f = CCO.DivergenceC2F(; adv_bcs...)
 
     @inbounds for i in 1:N_up
-        a_up_bcs = a_up_boundary_conditions(surf, edmf, i)
-        Iaf = CCO.InterpolateC2F(; a_up_bcs...)
         ρaw = prog_up_f[i].ρaw
         tends_ρaw = tendencies_up_f[i].ρaw
         nh_pressure = aux_up_f[i].nh_pressure
@@ -645,11 +636,10 @@ function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::Su
             prog_up[i].ρaq_ice .= max.(prog_up[i].ρaq_ice, 0)
         end
     end
-
+    # apply clipping at 0 and minimum area to ρaw
+    If = CCO.InterpolateC2F()
     @inbounds for i in 1:N_up
         @. prog_up_f[i].ρaw = max.(prog_up_f[i].ρaw, 0)
-        a_up_bcs = a_up_boundary_conditions(surf, edmf, i)
-        If = CCO.InterpolateC2F(; a_up_bcs...)
         @. prog_up_f[i].ρaw = Int(If(prog_up[i].ρarea) >= ρ_f * a_min) * prog_up_f[i].ρaw
     end
 
@@ -687,10 +677,12 @@ function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::Su
         @. prog_up[i].ρarea = ifelse(Ic(prog_up_f[i].ρaw) <= 0, FT(0), prog_up[i].ρarea)
         @. prog_up[i].ρaθ_liq_ice = ifelse(Ic(prog_up_f[i].ρaw) <= 0, FT(0), prog_up[i].ρaθ_liq_ice)
         @. prog_up[i].ρaq_tot = ifelse(Ic(prog_up_f[i].ρaw) <= 0, FT(0), prog_up[i].ρaq_tot)
+        if edmf.surface_area_bc isa FixedSurfaceAreaBC || edmf.surface_area_bc isa ClosureSurfaceAreaBC
+            a_surf = area_surface_bc(surf, edmf, i, edmf.surface_area_bc)
+            prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_surf
+        end
         θ_surf = θ_surface_bc(surf, grid, state, edmf, i)
         q_surf = q_surface_bc(surf, grid, state, edmf, i)
-        a_surf = area_surface_bc(surf, edmf, i)
-        prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_surf
         prog_up[i].ρaθ_liq_ice[kc_surf] = prog_up[i].ρarea[kc_surf] * θ_surf
         prog_up[i].ρaq_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * q_surf
     end

src/closures/perturbation_pressure.jl

@@ -54,9 +54,8 @@ function compute_nh_pressure!(state::State, grid::Grid, edmf::EDMFModel, surf)
         w_en = aux_en_f.w
 
         b_bcs = (; bottom = CCO.SetValue(b_up[kc_surf]), top = CCO.SetValue(b_up[kc_toa]))
-        a_bcs = a_up_boundary_conditions(surf, edmf, i)
         Ifb = CCO.InterpolateC2F(; b_bcs...)
-        Ifa = CCO.InterpolateC2F(; a_bcs...)
+        Ifa = CCO.InterpolateC2F()
 
         nh_press_buoy = aux_up_f[i].nh_pressure_b
         nh_press_adv = aux_up_f[i].nh_pressure_adv