Skip to content

Commit

Permalink
Merge #1358
Browse files Browse the repository at this point in the history 
1358: Explicitly declare InterpolateC2F kwargs r=costachris a=costachris

Explicitly declare InterpolateC2F kwargs for bottom and top boundaries on area.

Co-authored-by: costachris <christopouloscosta@gmail.com>
  • Loading branch information
@bors @costachris
bors[bot] and costachris authored Oct 19, 2023
2 parents 2f13254 + f582c2b commit e7c0a4f
Show file tree
Hide file tree
Showing 4 changed files with 22 additions and 26 deletions.
6 changes: 2 additions & 4 deletions driver/compute_diagnostics.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -303,16 +303,14 @@ function compute_diagnostics!(
end
end

Ifabulk = CCO.InterpolateC2F()
a_up_bulk_f = TC.face_aux_turbconv(state).bulk.a_up
@. a_up_bulk_f = Ifabulk(a_up_bulk)
@. a_up_bulk_f = TC.ᶠinterp_a(a_up_bulk)

RB_precip = CCO.RightBiasedC2F(; top = CCO.SetValue(FT(0)))

@inbounds for i in 1:N_up
Ifaup = CCO.InterpolateC2F()
a_up_f = aux_up_f[i].area
@. a_up_f = Ifaup(aux_up[i].area)
@. a_up_f = TC.ᶠinterp_a(aux_up[i].area)
@inbounds for k in TC.real_face_indices(grid)
diag_tc_f.nh_pressure[k] = 0
diag_tc_f.nh_pressure_b[k] = 0
Expand Down
24 changes: 12 additions & 12 deletions src/EDMF_functions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,23 +63,22 @@ 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
Ifae = CCO.InterpolateC2F()
If = CCO.InterpolateC2F(; bottom = CCO.SetValue(FT(0)), top = CCO.SetValue(FT(0)))

# Compute the mass flux and associated scalar fluxes
@. massflux = ρ_f * Ifae(a_en) * (w_en - toscalar(w_gm))
@. massflux_h = ρ_f * Ifae(a_en) * (w_en - toscalar(w_gm)) * (If(θ_liq_ice_en) - If(θ_liq_ice_gm))
@. massflux_qt = ρ_f * Ifae(a_en) * (w_en - toscalar(w_gm)) * (If(q_tot_en) - If(q_tot_gm))
@. massflux = ρ_f * ᶠinterp_a(a_en) * (w_en - toscalar(w_gm))
@. massflux_h = ρ_f * ᶠinterp_a(a_en) * (w_en - toscalar(w_gm)) * (If(θ_liq_ice_en) - If(θ_liq_ice_gm))
@. massflux_qt = ρ_f * ᶠinterp_a(a_en) * (w_en - toscalar(w_gm)) * (If(q_tot_en) - If(q_tot_gm))
@inbounds for i in 1:N_up
aux_up_f_i = aux_up_f[i]
aux_up_i = aux_up[i]
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 * 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)))
@. aux_up_f[i].massflux = ρ_f * ᶠinterp_a(a_up) * (w_up_i - toscalar(w_gm))
@. massflux_h += ρ_f * (ᶠinterp_a(a_up) * (w_up_i - toscalar(w_gm)) * (If(θ_liq_ice_up) - If(θ_liq_ice_gm)))
@. massflux_qt += ρ_f * (ᶠinterp_a(a_up) * (w_up_i - toscalar(w_gm)) * (If(q_tot_up) - If(q_tot_gm)))
end

if edmf.moisture_model isa NonEquilibriumMoisture
Expand All @@ -90,7 +89,7 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
q_ice_en = aux_en.q_ice
q_liq_gm = prog_gm.q_liq
q_ice_gm = prog_gm.q_ice
@. massflux_en = ρ_f * Ifae(a_en) * (w_en - toscalar(w_gm))
@. massflux_en = ρ_f * ᶠinterp_a(a_en) * (w_en - toscalar(w_gm))
@. massflux_ql = massflux_en * (If(q_liq_en) - If(q_liq_gm))
@. massflux_qi = massflux_en * (If(q_ice_en) - If(q_ice_gm))
@inbounds for i in 1:N_up
Expand Down Expand Up @@ -316,6 +315,8 @@ function surface_helper(surf::SurfaceBase, grid::Grid, state::State)
return (; ustar, zLL, oblength, ρLL)
end

const ᶠinterp_a = CCO.InterpolateC2F(bottom = CCO.Extrapolate(), top = CCO.Extrapolate())

function area_surface_bc(surf::SurfaceBase{FT}, edmf::EDMFModel, i::Int, bc::FixedSurfaceAreaBC)::FT where {FT}
N_up = n_updrafts(edmf)
return surf.bflux > 0 ? edmf.surface_area / N_up : FT(0)
Expand Down Expand Up @@ -576,7 +577,6 @@ 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...)
Expand All @@ -593,7 +593,8 @@ function compute_up_tendencies!(edmf::EDMFModel, grid::Grid, state::State, param
buoy = aux_up[i].buoy

@. tends_ρaw = -(∇f(wvec(LBC(ρaw * w_up))))
@. tends_ρaw += (ρaw * (I0f(entr_w) * w_en - I0f(detr_w) * w_up)) + (ρ_f * Iaf(a_up) * I0f(buoy)) + nh_pressure
@. tends_ρaw +=
(ρaw * (I0f(entr_w) * w_en - I0f(detr_w) * w_up)) + (ρ_f * ᶠinterp_a(a_up) * I0f(buoy)) + nh_pressure
tends_ρaw[kf_surf] = 0
end

Expand Down Expand Up @@ -637,10 +638,9 @@ function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::Su
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)
@. prog_up_f[i].ρaw = Int(If(prog_up[i].ρarea) >= ρ_f * a_min) * prog_up_f[i].ρaw
@. prog_up_f[i].ρaw = Int(ᶠinterp_a(prog_up[i].ρarea) >= ρ_f * a_min) * prog_up_f[i].ρaw
end

@inbounds for k in real_center_indices(grid)
Expand Down
8 changes: 4 additions & 4 deletions src/closures/perturbation_pressure.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,17 +55,17 @@ function compute_nh_pressure!(state::State, grid::Grid, edmf::EDMFModel, surf)

b_bcs = (; bottom = CCO.SetValue(b_up[kc_surf]), top = CCO.SetValue(b_up[kc_toa]))
Ifb = CCO.InterpolateC2F(; b_bcs...)
Ifa = CCO.InterpolateC2F()

nh_press_buoy = aux_up_f[i].nh_pressure_b
nh_press_adv = aux_up_f[i].nh_pressure_adv
nh_press_drag = aux_up_f[i].nh_pressure_drag
nh_pressure = aux_up_f[i].nh_pressure

@. nh_press_buoy = Int(Ifa(a_up) > 0) * -α_b / (1 + α₂_asp_ratio²) * ρ_f * Ifa(a_up) * Ifb(b_up)
@. nh_press_adv = Int(Ifa(a_up) > 0) * ρ_f * Ifa(a_up) * α_a * w_up * ∇(wvec(Ifc(w_up)))
@. nh_press_buoy = Int(ᶠinterp_a(a_up) > 0) * -α_b / (1 + α₂_asp_ratio²) * ρ_f * ᶠinterp_a(a_up) * Ifb(b_up)
@. nh_press_adv = Int(ᶠinterp_a(a_up) > 0) * ρ_f * ᶠinterp_a(a_up) * α_a * w_up * ∇(wvec(Ifc(w_up)))
# drag as w_dif and account for downdrafts
@. nh_press_drag = Int(Ifa(a_up) > 0) * -1 * ρ_f * Ifa(a_up) * α_d * (w_up - w_en) * abs(w_up - w_en) / H_up
@. nh_press_drag =
Int(ᶠinterp_a(a_up) > 0) * -1 * ρ_f * ᶠinterp_a(a_up) * α_d * (w_up - w_en) * abs(w_up - w_en) / H_up
@. nh_pressure = nh_press_buoy + nh_press_adv + nh_press_drag
end
return nothing
Expand Down
10 changes: 4 additions & 6 deletions src/update_aux.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -242,24 +242,22 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
# TODO: figure out why `ifelse` is allocating
# clip updraft w below minimum area threshold
@inbounds for i in 1:N_up
If = CCO.InterpolateC2F()
a_min = edmf.minimum_area
a_up = aux_up[i].area
@. aux_up_f[i].w = ifelse(If(a_up) >= a_min, max(prog_up_f[i].ρaw / (ρ_f * If(a_up)), 0), FT(0))
@. aux_up_f[i].w = ifelse(ᶠinterp_a(a_up) >= a_min, max(prog_up_f[i].ρaw / (ρ_f * ᶠinterp_a(a_up)), 0), FT(0))
end
@inbounds for i in 1:N_up
aux_up_f[i].w[kf_surf] = w_surface_bc(surf)
end

parent(aux_tc_f.bulk.w) .= 0
Ifb = CCO.InterpolateC2F()
@inbounds for i in 1:N_up
a_up = aux_up[i].area
Ifu = CCO.InterpolateC2F()
@. aux_tc_f.bulk.w += ifelse(Ifb(aux_bulk.area) > 0, Ifu(a_up) * aux_up_f[i].w / Ifb(aux_bulk.area), FT(0))
@. aux_tc_f.bulk.w +=
ifelse(ᶠinterp_a(aux_bulk.area) > 0, ᶠinterp_a(a_up) * aux_up_f[i].w / ᶠinterp_a(aux_bulk.area), FT(0))
end
# Assuming w_gm = 0!
@. aux_en_f.w = -Ifb(aux_bulk.area) / (1 - Ifb(aux_bulk.area)) * aux_tc_f.bulk.w
@. aux_en_f.w = -1 * ᶠinterp_a(aux_bulk.area) / (1 - ᶠinterp_a(aux_bulk.area)) * aux_tc_f.bulk.w

#####
##### diagnose_GMV_moments
Expand Down

0 comments on commit e7c0a4f

Please sign in to comment.