diff --git a/integration_tests/utils/main.jl b/integration_tests/utils/main.jl index d418090672..dca6d77fe3 100644 --- a/integration_tests/utils/main.jl +++ b/integration_tests/utils/main.jl @@ -175,7 +175,11 @@ cent_aux_vars_edmf(FT, n_up) = (; KM = FT(0), KH = FT(0), mixing_length = FT(0), + m_entr_detr = FT(0), + ∇m_entr_detr = FT(0), θ_virt = FT(0), + ∇θ_virt = FT(0), + ∇Ri_bulk = FT(0), massflux_tendency_h = FT(0), massflux_tendency_qt = FT(0), diffusive_tendency_h = FT(0), diff --git a/src/Operators.jl b/src/Operators.jl index 9ec766a3c1..ecf53e1af5 100644 --- a/src/Operators.jl +++ b/src/Operators.jl @@ -36,24 +36,6 @@ end ∇c2f(f::SA.SVector, grid::Grid, ::BottomBCTag, bc::SetValue) = (f[1] - bc.value) * grid.Δzi * 2.0 ∇c2f(f::SA.SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value -c∇_upwind(f, grid::Grid, k; bottom = NoBCGivenError(), top = NoBCGivenError()) = - c∇_upwind(ccut_upwind(f, grid, k), grid, k; bottom, top) - -function c∇_upwind(f_dual::SA.SVector, grid::Grid, k; bottom = NoBCGivenError(), top = NoBCGivenError()) - if is_surface_center(grid, k) - return c∇_upwind(f_dual, grid, BottomBCTag(), bottom) - elseif is_toa_center(grid, k) - return c∇_upwind(f_dual, grid, TopBCTag(), top) - else - return c∇_upwind(f_dual, grid, InteriorTag()) - end -end -c∇_upwind(f::SA.SVector, grid::Grid, ::InteriorTag) = (f[2] - f[1]) * grid.Δzi -c∇_upwind(f::SA.SVector, grid::Grid, ::TopBCTag, bc::FreeBoundary) = (f[2] - f[1]) * grid.Δzi -c∇_upwind(f::SA.SVector, grid::Grid, ::TopBCTag, bc::SetGradient) = bc.value -c∇_upwind(f::SA.SVector, grid::Grid, ::BottomBCTag, bc::SetValue) = (f[1] - bc.value) * (grid.Δzi * 2) -c∇_upwind(f::SA.SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value - # Used when traversing cell faces interpc2f(f, grid::Grid, k::CCO.PlusHalf; bottom = NoBCGivenError(), top = NoBCGivenError()) = diff --git a/src/turbulence_functions.jl b/src/turbulence_functions.jl index 7e70b9adb4..67f67eb9eb 100644 --- a/src/turbulence_functions.jl +++ b/src/turbulence_functions.jl @@ -16,15 +16,20 @@ function get_inversion(grid::Grid, state, param_set, Ri_bulk_crit) u = center_prog_grid_mean(state).u v = center_prog_grid_mean(state).v Ri_bulk = center_aux_bulk(state).Ri + ∇θ_virt = center_aux_turbconv(state).∇θ_virt + ∇Ri_bulk = center_aux_turbconv(state).∇Ri_bulk θ_virt_b = θ_virt[kc_surf] z_c = grid.zc + ∇c = CCO.DivergenceF2C() + wvec = CC.Geometry.WVector # test if we need to look at the free convective limit if (u[kc_surf]^2 + v[kc_surf]^2) <= 0.01 kmask = map(k -> (k, θ_virt[k] > θ_virt_b), real_center_indices(grid)) k_star = first(kmask[findlast(km -> km[2], kmask)]) - ∇θ_virt = c∇_upwind(θ_virt, grid, k_star; bottom = SetGradient(0), top = FreeBoundary()) - h = (θ_virt_b - θ_virt[k_star - 1]) / ∇θ_virt + z_c[k_star - 1] + LB = CCO.LeftBiasedC2F(; bottom = CCO.SetValue(θ_virt[kc_surf])) + @. ∇θ_virt = ∇c(wvec(LB(θ_virt))) + h = (θ_virt_b - θ_virt[k_star - 1]) / ∇θ_virt[k_star] + z_c[k_star - 1] else Ri_bulk_fn(k) = g * (θ_virt[k] - θ_virt_b) * z_c[k] / θ_virt_b / (u[k] * u[k] + v[k] * v[k]) @@ -33,8 +38,9 @@ function get_inversion(grid::Grid, state, param_set, Ri_bulk_crit) end kmask = map(k -> (k, Ri_bulk_fn(k) > Ri_bulk_crit), real_center_indices(grid)) k_star = first(kmask[findlast(km -> km[2], kmask)]) - ∇Ri_bulk = c∇_upwind(Ri_bulk, grid, k_star; bottom = SetGradient(0), top = FreeBoundary()) - h = (Ri_bulk_crit - Ri_bulk[k_star - 1]) / ∇Ri_bulk + z_c[k_star - 1] + LB = CCO.LeftBiasedC2F(; bottom = CCO.SetValue(Ri_bulk[kc_surf])) + @. ∇Ri_bulk = ∇c(wvec(LB(Ri_bulk))) + h = (Ri_bulk_crit - Ri_bulk[k_star - 1]) / ∇Ri_bulk[k_star] + z_c[k_star - 1] end return h diff --git a/src/update_aux.jl b/src/update_aux.jl index 093d27573b..749cf3dd40 100644 --- a/src/update_aux.jl +++ b/src/update_aux.jl @@ -35,6 +35,9 @@ function update_aux!(edmf::EDMF_PrognosticTKE{N_up}, gm, grid, state, Case, para prog_up_f = face_prog_updrafts(state) aux_en_unsat = aux_en.unsat aux_en_sat = aux_en.sat + m_entr_detr = aux_tc.m_entr_detr + ∇m_entr_detr = aux_tc.∇m_entr_detr + wvec = CC.Geometry.WVector ##### ##### center variables @@ -238,20 +241,20 @@ function update_aux!(edmf::EDMF_PrognosticTKE{N_up}, gm, grid, state, Case, para ##### compute_updraft_closures ##### - @inbounds for k in real_center_indices(grid) - @inbounds for i in 1:N_up + Ic = CCO.InterpolateF2C() + ∇c = CCO.DivergenceF2C() + LB = CCO.LeftBiasedC2F(; bottom = CCO.SetValue(FT(0))) + + @inbounds for i in 1:N_up + # compute ∇m at cell centers + a_up = aux_up[i].area + w_up = aux_up_f[i].w + w_gm = prog_gm_f.w + @. m_entr_detr = a_up * (Ic(w_up) - Ic(w_gm)) + @. ∇m_entr_detr = ∇c(wvec(LB(m_entr_detr))) + @inbounds for k in real_center_indices(grid) # entrainment if aux_up[i].area[k] > 0.0 - # compute ∇m at cell centers - a_up_c = aux_up[i].area[k] - w_up_c = interpf2c(aux_up_f[i].w, grid, k) - w_gm_c = interpf2c(prog_gm_f.w, grid, k) - m = a_up_c * (w_up_c - w_gm_c) - a_up_cut = ccut_upwind(aux_up[i].area, grid, k) - w_up_cut = daul_f2c_upwind(aux_up_f[i].w, grid, k) - w_gm_cut = daul_f2c_upwind(prog_gm_f.w, grid, k) - m_cut = a_up_cut .* (w_up_cut .- w_gm_cut) - ∇m = FT(c∇_upwind(m_cut, grid, k; bottom = SetValue(0), top = FreeBoundary())) w_min = 0.001 @@ -267,8 +270,8 @@ function update_aux!(edmf::EDMF_PrognosticTKE{N_up}, gm, grid, state, Case, para b_up = aux_up[i].buoy[k], b_en = aux_en.buoy[k], tke = aux_en.tke[k], - dMdz = ∇m, - M = m, + dMdz = ∇m_entr_detr[k], + M = m_entr_detr[k], a_up = aux_up[i].area[k], a_en = aux_en.area[k], R_up = edmf.pressure_plume_spacing[i], @@ -298,7 +301,6 @@ function update_aux!(edmf::EDMF_PrognosticTKE{N_up}, gm, grid, state, Case, para end end - wvec = CC.Geometry.WVector @inbounds for i in 1:N_up # pressure