Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Include prognostic and closure options for surface area BC #1354

Merged
merged 1 commit into from
Oct 12, 2023
Merged

Include prognostic and closure options for surface area BC #1354

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 2 additions & 4 deletions driver/compute_diagnostics.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 2 additions & 0 deletions driver/generate_namelist.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
33 changes: 28 additions & 5 deletions driver/initial_conditions.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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
2 changes: 1 addition & 1 deletion integration_tests/Artifacts.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,4 +5,4 @@ git-tree-sha1 = "3d6bb62d803fb01908ffe85469a89933f28ced45"
git-tree-sha1 = "c2d323fb412d1250182be3aeadf9d94e816550a0"

[SCAMPy_output]
git-tree-sha1 = "6b82611f969d3cf99102baa7d898e3c759cf014c"
git-tree-sha1 = "f725942089ff0d0fa8138dcbd1349996c746a4f9"
86 changes: 39 additions & 47 deletions src/EDMF_functions.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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}
Expand All @@ -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_)
Expand All @@ -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)
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
3 changes: 1 addition & 2 deletions src/closures/perturbation_pressure.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading
Loading