Skip to content

Commit

Permalink
Make env entr-detr src terms explicit
Browse files Browse the repository at this point in the history
  • Loading branch information
@charleskawczynski
charleskawczynski committed Oct 30, 2021
1 parent b8d94f8 commit 69959d4
Show file tree
Hide file tree
Showing 2 changed files with 28 additions and 47 deletions.
32 changes: 3 additions & 29 deletions src/Operators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -538,31 +538,19 @@ some terms are treated implicitly. Here is a list of all
the terms in the matrix (and the overall equation solved):
N_a terms = 1 (diffusion)
N_diagonal terms = 6 (unsteady+upwind_advection+2diffusion+entr_detr+dissipation)
N_diagonal terms = 5 (unsteady+upwind_advection+2diffusion+dissipation)
N_c terms = 2 (diffusion+upwind_advection)
∂_t ρa*tke ✓ unsteady
+ ∂_z(ρaw*tke) = ✓ advection
+ ρaK*(∂²_z(u)+∂²_z(v)+∂²_z(w̄))
+ ρawΣᵢ(εᵢ(wⱼ-w₀)²-δ₀*tke) ✓ entr_detr
+ ρawΣᵢ(εᵢ(wⱼ-w₀)²-δ₀*tke)
+ ∂_z(ρa₀K ∂_z(tke)) ✓ diffusion
+ ρa₀*w̅₀b̅₀
- a₀(u⋅∇p)
+ ρa₀D ✓ dissipation
=#
function construct_tridiag_diffusion_en(
grid::Grid,
param_set::APS,
state,
TS,
n_updrafts::Int,
minimum_area::Float64,
pressure_plume_spacing::Vector,
frac_turb_entr,
entr_sc,
mixing_length,
is_tke,
)
function construct_tridiag_diffusion_en(grid::Grid, param_set::APS, state, TS, n_updrafts::Int, mixing_length, is_tke)

c_d = CPEDMF.c_d(param_set)
kc_surf = kc_surface(grid)
Expand Down Expand Up @@ -599,19 +587,6 @@ function construct_tridiag_diffusion_en(
end

@inbounds for k in real_center_indices(grid)
D_env = 0.0

@inbounds for i in 1:n_updrafts
if aux_up[i].area[k] > minimum_area
turb_entr = frac_turb_entr[i, k]
R_up = pressure_plume_spacing[i]
w_up_c = interpf2c(aux_up_f[i].w, grid, k)
D_env += ρ0_c[k] * aux_up[i].area[k] * w_up_c * (entr_sc[i, k] + turb_entr)
else
D_env = 0.0
end
end

# TODO: this tridiagonal matrix needs to be re-verified, as it's been pragmatically
# modified to not depend on ghost points, and these changes have not been
# carefully verified.
Expand All @@ -625,7 +600,6 @@ function construct_tridiag_diffusion_en(
ρ0_c[k] * ae[k] * dti - ρ0_c[k] * ae[k] * w_en_c[k] * Δzi +
rho_ae_K_m[k + 1] * Δzi * Δzi +
rho_ae_K_m[k] * Δzi * Δzi +
D_env +
ρ0_c[k] * ae[k] * c_d * sqrt(max(prog_en.tke[k], 0)) / max(mixing_length[k], 1)
)
if is_toa_center(grid, k)
Expand Down
43 changes: 25 additions & 18 deletions src/Turbulence_PrognosticTKE.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -279,28 +279,17 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca
implicit_eqs = edmf.implicit_eqs
# Matrix is the same for all variables that use the same eddy diffusivity, we can construct once and reuse

common_args = (
grid,
param_set,
state,
TS,
up.n_updrafts,
edmf.minimum_area,
edmf.pressure_plume_spacing,
edmf.frac_turb_entr,
edmf.entr_sc,
edmf.mixing_length,
)
common_args = (grid, param_set, state, TS, up.n_updrafts, edmf.mixing_length)

implicit_eqs.A_TKE .= construct_tridiag_diffusion_en(common_args..., true)
implicit_eqs.A_Hvar .= construct_tridiag_diffusion_en(common_args..., false)
implicit_eqs.A_QTvar .= construct_tridiag_diffusion_en(common_args..., false)
implicit_eqs.A_HQTcov .= construct_tridiag_diffusion_en(common_args..., false)

implicit_eqs.b_TKE .= en_diffusion_tendencies(grid, state, TS, :tke, n_updrafts)
implicit_eqs.b_Hvar .= en_diffusion_tendencies(grid, state, TS, :Hvar, n_updrafts)
implicit_eqs.b_QTvar .= en_diffusion_tendencies(grid, state, TS, :QTvar, n_updrafts)
implicit_eqs.b_HQTcov .= en_diffusion_tendencies(grid, state, TS, :HQTcov, n_updrafts)
implicit_eqs.b_TKE .= en_diffusion_tendencies(edmf, grid, state, TS, :tke, n_updrafts)
implicit_eqs.b_Hvar .= en_diffusion_tendencies(edmf, grid, state, TS, :Hvar, n_updrafts)
implicit_eqs.b_QTvar .= en_diffusion_tendencies(edmf, grid, state, TS, :QTvar, n_updrafts)
implicit_eqs.b_HQTcov .= en_diffusion_tendencies(edmf, grid, state, TS, :HQTcov, n_updrafts)
# -----------

###
Expand Down Expand Up @@ -868,17 +857,24 @@ function compute_covariance_dissipation(edmf::EDMF_PrognosticTKE, grid, state, c
return
end

function en_diffusion_tendencies(grid::Grid, state, TS, covar_sym::Symbol, n_updrafts)
function en_diffusion_tendencies(edmf, grid::Grid, state, TS, covar_sym::Symbol, n_updrafts)
dti = TS.dti
b = center_field(grid)
ρ0_c = center_ref_state(state).ρ0
prog_en = center_prog_environment(state)
aux_en_2m = center_aux_environment_2m(state)
aux_up_f = face_aux_updrafts(state)
prog_covar = getproperty(prog_en, covar_sym)
aux_covar = getproperty(aux_en_2m, covar_sym)
aux_up = center_aux_updrafts(state)

minimum_area = edmf.minimum_area
frac_turb_entr = edmf.frac_turb_entr
pressure_plume_spacing = edmf.pressure_plume_spacing
entr_sc = edmf.entr_sc

ae = center_field(grid)
FT = eltype(grid)

@inbounds for k in real_center_indices(grid)
ae[k] = 1 .- sum(ntuple(i -> aux_up[i].area[k], n_updrafts))
Expand All @@ -888,6 +884,17 @@ function en_diffusion_tendencies(grid::Grid, state, TS, covar_sym::Symbol, n_upd
covar_surf = prog_covar[kc_surf]

@inbounds for k in real_center_indices(grid)
D_env = sum(1:n_updrafts) do i
if aux_up[i].area[k] > minimum_area
turb_entr = frac_turb_entr[i, k]
R_up = pressure_plume_spacing[i]
w_up_c = interpf2c(aux_up_f[i].w, grid, k)
D_env_i = ρ0_c[k] * aux_up[i].area[k] * w_up_c * (entr_sc[i, k] + turb_entr)
else
D_env_i = FT(0)
end
D_env_i
end
if is_surface_center(grid, k)
b[k] = covar_surf
else
Expand All @@ -897,7 +904,7 @@ function en_diffusion_tendencies(grid::Grid, state, TS, covar_sym::Symbol, n_upd
aux_covar.buoy[k] +
aux_covar.shear[k] +
aux_covar.entr_gain[k] +
aux_covar.rain_src[k]
aux_covar.rain_src[k] - D_env * prog_covar[k]
)
end
end
Expand Down

0 comments on commit 69959d4

Please sign in to comment.